Musaline library

Introduction

Musaline is a C++ library for alignment of sequences, with a focus on musical sequences.

In order to compare two sequences, we can align the elements of the sequences to each other. Given a cost function for aligning elements, this gives a similarity measure between sequences. The similarity is computed with the Needleman-Wunsch algorithm and variants.

This library can align arbitrary sequences of arbitrary elements. The user must supply a class for computing the cost of matching elements. There are several cost functions supplied for the music domain, e.g. for aligning melodies (sequences of notes).

Building and using the library

The library itself is contained in header files, so no building is needed. There are some test files that can be built. The library comes with CMake files that allow for easy building and installing, if you know how to use CMake.

If you build a program that uses the library, and you want maximum efficiency, you should supply GJMATRIX_NO_CHECKS to the preprocessor, usually by giving -DGJMATRIX_NO_CHECKS on the command line, in a makefile or other build system.

When you use g++, you may need to give the option -std=c++0x.

The library should be usable in multithreaded code.

Framework overview

The classes musaline::LinearAligner, musaline::AffineAligner and musaline::ConsolidatingAligner are the main classes.

This file explains the interface. It is a good idea to first study the example program test/string_example.cpp. This program is well commented and explains most of the library. The output of the program (or actually, of a similar program that outputs html), can be found in test/text_result.html.

Framework requirements

The library is a very general framework. That means that the user of the framework must supply some input to get a concrete aligner.

We talked about sequences of elements that are to be aligned. Question that arise are:

• what are the elements?
• how should we represent a sequence?
• how do you decide what alignment is good?

The answer to all those questions is: it is up to you, the user of the framework. But you should specify it in a certain way. So, if you want your elements to be integers and your sequences to be arrays of integers? OK. If you want the sequences to be standard vectors or lists of a struct defined by yourself? Fine.

What and how you should provide depends a bit on the particular aligner. We start with the most simple case.

Linear Aligner requirements

The Linear Aligner aligns an element of a sequence either with an element of the other sequence or with a gap.

The alignment algorithms are templatized by classes CostCalculator, Sequence0 and Sequence1. Usually, the two sequence classes will be the same, or at least they will have the same value type.

Here is a complete example of a program that uses the Linear Aligner.

#include "musaline.hpp"
#include <iostream>

class MyCostCalculator1 {
  public:
    double match(int a, int b) const
      { return a==b? 1.2 : -1.0; }
    double gap_bft() const
      { return -2.0; }
    void preprocess(int const *, int const *) const
      {}
};

int main()
{
    musaline::LinearAligner<MyCostCalculator1> aligner;
    int seq0[] = { 5,3,78,2,1,1,2};
    int seq1[] = { 3,7,2,1,2};
    double score=aligner.global_align(seq0, seq1).score;
    std::cout<<"Score: "<<score<<"\n";
}

In this example we align two arrays (the sequence type) of integers (the element type). The aligner object is an instance of the musaline::LinearAligner class, with template parameter MyCostCalculator1. The latter is the class through which the user supplies the required functionality. Three functions are required in this class. The function 'match' returns the benefit when aligning two elements -ints, in this case. Good matches should return a high value, bad matches a low, possibly negative, value. The function 'gap_bft' returns the benefit for aligning an element with no element. The function 'preprocess' is called once with both sequences as parameter, before the actual alignment starts. In this example, the function does nothing.

#include "musaline.hpp"
#include <iostream>
#include <vector>

struct MyElement {
    int value;
};

class MyCostCalculator1 {
  public:
    double match(MyElement const &a, MyElement const &b) const
      { return a.value==b.value? 1.2 : -1.0; }
    double gap_bft() const
      { return -2.0; }
    void preprocess(std::vector<MyElement> const &, std::vector<MyElement> const &) const
      {}
};

int main()
{
    musaline::LinearAligner<MyCostCalculator1> aligner;
    std::vector<MyElement> seq0;
    std::vector<MyElement> seq1;
    // initialise seq0 and seq1
    double score=aligner.global_align(seq0, seq1).score;
    std::cout<<"Score: "<<score<<"\n";
}

Finally we show that it is also possible to use different sequences with the same aligner. Here, we use arrays, standard vectors and lists. The only adaptation we have to make to the CostCalculator class is to make preprocess a templated function.

#include <iostream>
#include <vector>
#include <list>

struct MyElement {
    int value;
};

class MyCostCalculator1 {
    public:
        double match(MyElement const &a, MyElement const &b) const
        { return a.value==b.value? 1.2 : -1.0; }
        double gap_bft() const
        { return -2.0; }
        template <class Seq0, class Seq1>
        void preprocess(Seq0 const &, Seq1 const &) const
        {}
};

int main()
{
    musaline::LinearAligner<MyCostCalculator1> aligner;
    MyElement arr0[] = { 5,3,78,2,1,1,2};
    MyElement arr1[] = { 3,7,2,1,2};
    std::vector<MyElement> vec1(arr1,arr1+5);
    std::list<MyElement> lis1(vec1.begin(),vec1.end());
    double score=aligner.global_align(arr0, arr1).score;
    std::cout<<"Score a-a: "<<score<<"\n";
    score=aligner.global_align(arr0, arr1).score;
    std::cout<<"Score a-v: "<<score<<"\n";
    score=aligner.global_align(arr0, lis1).score;
    std::cout<<"Score a-l: "<<score<<"\n";
}

Consolidating Aligner requirements

The consolidating aligner can align one element of a sequence with several elements of another sequence. We will illustrate this with an example where we will align two melodies. A melody is a sequence of notes and a note has a pitch and a duration. In a melody, e.g., two quarter notes have the same duration as a half note. We take these durations into account during the alignment.

The basic idea of the aligner example here is that while two melodies have the same pitch, they gain 1.0 per duration unit. When the pitches differ, or a note is not matched, they lose 1.0 per duration unit.

Below we see the implementation of the required functions match_gap_with1, match_gap_with2, match_range_with2 and match_range. The implementation of the most complicated function, match_range_with1, is shown later.

#include <vector>
using std::vector;

struct Note
{
    char pitch;
    int duration;
};

struct NoteCostCalculator {
    void preprocess(vector<Note> const &seq1, vector<Note> const &seq2) {}
    template <class Seq1, class Seq2>
    void preprocess(Seq1 const &seq1, Seq2 const &seq2) {}

    double match_gap_with1(Note const &n)
    {
        return -n.duration;
    }

    double match_gap_with2(Note const &n)
    {
        return -n.duration;
    }

    template <class Iter>
    std::vector<double> match_range_with1(Note const &n,
            Iter begin,
            Iter cur);

    template <class Iter>
    std::vector<double> match_range_with2(
            Iter begin,
            Iter cur, Note const &n)
    {
        return match_range_with1(n,begin,cur);
    }

    double match(Note const &n1, Note const &n2)
    {
        auto minmaxpair = std::minmax (n1.duration, n2.duration);
        if (n1.pitch == n2.pitch) {

            return 2.0*minmaxpair.first - minmaxpair.second;
        } else {
            return -1.0*minmaxpair.second;
        }
    }
};

The required member function match_range_with1 is a lot more subtle. The framework will call this function with a note of the first sequence as a first parameter and an initial subsequence of the second sequence, represented as a begin and end iterator as second and third parameter. The function should return the benefit of aligning the note with the last two notes, the last three notes, the last four notes, etcetera. The case for one note or no notes (a gap) is taken care of in the functions match and match_gap_with1. The function should return those alignment benefits in a vector of doubles. It is permitted to cut off the vector. In this case, we do that if the total duration of the sequence of notes exceeds the duration of the single note. By cutting off the number of returned elements, the algorithm will also take a lot less time.

template <class Iter>
std::vector<double>
NoteCostCalculator::match_range_with1(Note const &n,
            Iter begin, Iter cur)
{
    std::vector<double> result;
    int total_duration = 0;
    int matching_duration = 0;
    int i=0;
    while (cur != begin) {
        --cur;
        ++i;
        if (cur->pitch ==n.pitch) {
            matching_duration += cur->duration;
        }
        total_duration += cur->duration;
        if (i>1) {
            int max_duration=std::max (total_duration,n.duration);
            result.push_back(2.0*matching_duration - max_duration);
        }
        if (total_duration>=n.duration)
            break;
    }
    return result;
}

The required function match_range_with2 is the same as match_range_with1, with the two sequences exchanged.In this case, we just forward the call, because the situation is symmetrical. This need not be the case. For instance, we may want to transpose one of the melodies in the preprocess phase.

As an example, consider the two sequences a4-b1-c1-b1-a1-b4 and d4-b4-b4. Here, we represent the pitch as a letter and the duration as an integer. At a certain moment, the framework will call the match_range_with2 function with second note of the second sequence (b4) and the first 5 notes of the first sequence (a4-b1-c1-b1-a1).

During the while loop, the following values will be assumed for *cur, matching_duration, total_duration and the benefit value.

Note	Match	Tot	Value
a1	0	1	-1
b1	1	2	0
c1	1	3	-1
b1	2	4	0

The function will return a vector with contents 0.0, -1.0, 0.0. The first value -1 is skipped, and the loop ends when total_duration is 4.

The algorithmic complexity is dependent on match_range. It is called for every pair of elements from the two sequences. The amount of work the framework does for one call is linear in the length of the resulting vector. The amount of work that is done in the function will also be at least linear in the output size. That means that if the amount of work that is done in match_range is linear in the length of the input sequence, the total running time is cubic in the sequence length. But if we can guarantee that the amount of work done in the function is constant, the running time is quadratic.

In the example of matching melodies, it might be possible to guarantee that the shortest duration is at most -say- 32 times smaller than the largest duration. That would mean that we only spend a constant time (at most 32 loop iterations) in match_range .