On Sunday, June 30, 2013, Erik Erlandson wrote:

...

I have had a quick look. It seems that you've implemented the Needleman-Wunsch algorithm, which takes quadratic time and space. There are less complex algorithms, at least for similar input sequences, e.g. Myers - "An O(ND) Difference Algorithm and Its Variations"; I have an implementation of that here: http://svn.chezphil.org/anyterm/trunk/src/diff.cc & .hh.

What is your motivation for choosing Needleman-Wunsch?

It was the simplest algorithm to get going. I was interested in developing a nice boost-ish programmatic interface, and developing unit tests, so for the implementation I stuck with the easy algorithm that I was familiar with. Now I (or anybody else) can drop in better algorithms with some unit test assurance.

Plus, Needleman-Wunsch is pretty effective for applications involving relatively small strings or sequences. Not so much for applications like large file diffs, or DNA sequences, etc.

At any rate, I view the particular algorithm as somewhat fungible, provided that the interface design is appealing to people. If there is disagreement on that point, I'm interested to know that too.

In the current interface I don't see how one could implement a replacement algorithm. -- -- -- Grafik - Don't Assume Anything -- Redshift Software, Inc. - http://redshift-software.com -- rrivera/acm.org - grafik/redshift-software.com -- 102708583/icq - grafikrobot/aim - grafikrobot/yahoo

Erik Erlandson wrote:

I am grappling with how best to represent the returned "edit script".

"Best" obviously depends on what the caller is going to do with it. A good design for the interface should emerge once some experience has been gained with using the algorithm in actual applications. Attempting to design a general-purpose interface too soon can be a mistake. Having said that, the most generic way to return the edit script would be to template the algorithm on an output-handling class: template <typename ITER1, typename ITER2, typename OUTPUT> void diff(ITER1 begin1, ITER1 end1, ITER2 begin2, ITER2 end2, OUTPUT& output) { ....... // eventually calls object's methods, something like this: output.from_1(i,j); // "deletion", present in 1 but not in 2 output.from_2(p,q); // "insertion", present in 2 but not in 1 output.from_both(w,x, y,z); // common subsequence present in both } Used as follows: class diff_output // Write edit script to cout in the style of the diff program { public: template <typename ITER> void from_1(ITER a, ITER b) { std::cout << "< " << string(a,b) << "\n"; } template <typename ITER> void from_2(ITER a, ITER b) { std::cout << "> " << string(a,b) << "\n"; } template <typename ITER1, typename ITER2> void from_both(ITER1, ITER1, ITER2, ITER2) { } }; An important feature of this is that it doesn't store the output. If the caller wants to store the output, they can supply an object that does that. Regards, Phil.

Erik Erlandson wrote:

...

...

I am grappling with how best to represent the returned "edit script".

the most generic way to return the edit script would be to template the algorithm on an output-handling class:

template <typename ITER1, typename ITER2, typename OUTPUT> void diff(ITER1 begin1, ITER1 end1, ITER2 begin2, ITER2 end2, OUTPUT& output)

I had a design using this approach, which has obvious advantages. The one thing I wasn't crazy about is that to do this, the implementation layer is committed to always computing and/or storing start+length information and edit-op cost information, so that the output handler has its space of choices about which information it cares about.

Not necessarily; you can allow a variety of methods in the output handler, with different names or signatures, and use "enable-if" SFINAE techniques to enable appropriate parts of the implementation algorithm. (I don't think I've ever seen anything like this done in a std:: or Boost library. Maybe someone else can think of concrete examples?) Regards, Phil.

I reworked both the algorithms and the interface. https://github.com/erikerlandson/edit_distance I devised a variation on the Dijkstra Single Source Shortest Path algorithm, optimized for the particular structure of an edit graph. Like the Meyers algorithm, it is very efficient on long sequences with localized differences (e.g. computing the diff of two similar files, or comparing two similar DNA sequences). Under such conditions, I've been comparing two similar sequences with 100 million elements in under a second. I went with an SSSP-based approach because the Meyers algorithm makes assumptions about the edit operation cost function, and a general-purpose library should allow any cost function (think std::sort() or std::map<> allowing any definition of '<'). I don't feel like I sacrificed a lot of performance. It's quite fast. I doubt it would be hard to allow the user to provide a specialized 'meyers_cost' object, and have it trigger a specific Meyers algorithm under the hood. I also upgraded the interface to a named-argument style function, using the boost::parameter library, so calls now support some named optional parameters: // apply an optional custom cost, and an optional beam constraint int d = edit_distance(seq1, seq2, _cost = my_cost(), _beam=100); In the case of edit_alignment(), the actual sequence of edit operations is received by a user-supplied output handling object: // output_handler provides methods for handling insertion, deletion, substitution and equal elements int d = edit_alignment(seq1, seq2, output_handler, _cost=my_cost(), _beam=1000);

project repo: https://github.com/erikerlandson/edit_distance Since last time: I implemented specializations for the Myers algorithm, which matches whenever the necessary conditions are met: *) sequences are random access *) no substitution (i.e. _allow_sub = boost::false_type(), which is now the default) *) unit cost for insertion and deletion (also the default) *) no optional pruning behaviors (excepting the new _max_cost, which is supported) Default behavior is now such that the Myers algorithm applies, unless options are specifically chosen such that it cannot be used. I altered the semantic for cost function: a separate equality relation is now applied explicitly to determine equality, which takes precedence over substitution cost function. optional behaviors are now: _cost = <cost-function-object> // customize the cost functions _equal = <equal-object> // customize the element equality relation _allow_sub = <true/false/true_type/false_type> // enable or disable substitution, at run time or compile time _max_cost <cost value> // if edit cost grows too large, halt best-path computation and fall back to fast/dumb linear completion _max_cost_exception // if true, throw an exception when max cost is exceeded instead of completing _cost_beam = <cost_type value> // prune "unpromising" paths based on high cost per position _edit_beam = <unsigned int> // prune paths that are too far off the diagonal If substitution is disabled at compile time (false_type()), then cost functions and the output functions are not required to define substitution methods. Now that there are two very different algorithms (Myers and SSSP), there is improved unit-testing. Obtaining the same results on both algorithms is an excellent cross-check. I've had interest in supporting an application related to: https://github.com/erikerlandson/edit_distance/issues/5 I have a plan for supporting that, which will probably take a bit of time to work out. I've been considering possible alternate names for the functions: edit_distance / edit_alignment or: edit_cost / edit_path or: edit_cost / edit_script I've also been considering whether or not to support for _cost_beam and/or _edit_beam. Unsure how much call there is to prune at the expense of correctness.

----- Original Message -----

On 12/14/2013 03:54 AM, Erik Erlandson wrote:

...

I've been considering possible alternate names for the functions: edit_distance / edit_alignment or: edit_cost / edit_path or: edit_cost / edit_script

I suggest that you rename the boost::algorithm::sequence_alignment namespace to boost::algorithm::edit.

My original idea for boost::algorithm::sequence_alignment was that it could also serve as a home for other sequence alignment algorithms, such as BLAST, or HMMs. However, that is pretty speculative.

Then we can drop the "edit_" prefix on the function names, giving boost:algorithm::edit::distance and likewise for alignment.

Just 'distance' or 'alignment' (or 'path') seem like they might be too generic. The likelihood that they'd need to be qualified seems higher. I'll give it some thought.

The edit::distance is reminiscient of std::distance, so I think that name should be retained. I have no strong opinion about edit::alignment. Research papers use "edit alignment" but I find "edit path" slightly more descriptive.

I've been leaning toward 'edit_distance' and 'edit_path'. My survey of the literature hasn't given me any impression of an unambiguous standard for terminology.

Talking about naming, it may be a good idea avoid abbreviations for some of the names. So "allow_sub" could become "allow_substitution" (or just "substitution".)

Agreed, my abbreviations aren't really warranted.

Likewise, the cost function object has names like cost_ins and cost_del, but as these are already located inside the cost function object, they could be renamed to insertion and deletion (or erasure to match the STL namining a bit better.)

That makes sense. I could do the same for the output function object.

On 12/16/2013 10:45 PM, Erik Erlandson wrote:

Come to think of it, now that the interface sits on top of boost::parameter, it ought to be possible to simply present the single 'edit_distance' function, and enable the path-output implementations iff an '_output = outputter_class' is provided.

Very good idea.

Tangentially, I considered 'edit_cost' instead of 'edit_distance' mostly because in the general case, the function may violate the definition of a distance metric, for example if the cost of insertion is different than the cost of deletion. In that case, edit_cost(seq1, seq2) != edit_cost(seq2, seq1).

Which reminds me... If you are running without substitution (_allow_sub set to false), then we may need to have a lower cost for deletions than for insertions, because for substitutions we want a script with delete(A) followed by insert(A), rather than insert(A) followed by delete(A). The latter may result in 'A' being removed rather than replaced.

On 12/17/2013 10:22 PM, Erik Erlandson wrote:

I did a few experiments trying to cajole the algorithms into outputting edit operations in some kind of canonical ordering. Getting the algorithms to do this natively didn't look very promising, and it left me feeling like the best solution would be to collect operations in the output-object, and then sort them after-the-fact into some desired ordering, for cases where that is important.

You are probably right. I can think of two solutions, but both are rather complex. One solution would be to have a path-dependent cost function (as you suggested earlier in issue #5), and the other to prune the insert-delete combinations from the internal search tree. Regarding your suggestion, AFAIK, the contraint is always between adjacent path elements, and it therefore can be handled in O(1) time in the output object via element swapping. So it seems like a feasible solution.