Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm looking for robust code to solve the "Longest Common Substring" problem. I can just code it up from that description, but I'd thought I'd ask here, first, in case someone knows of an implementation either distributed with Mathematica or available from an open source. I found a hint here that a solution might be part of the (huge) Combinatorica package, but a quick search of the documentation did not disclose it.

share|improve this question
4  
You've seen LongestCommonSubsequence[]? –  J. M. May 28 '12 at 16:51
    
@J.M. Beat me by 10 secs :D –  belisarius May 28 '12 at 16:52
    
As far as I know, LongestCommonSubsequence only returns the first hit. That might not be robust enough as there could be multiple distinct and different hits. –  István Zachar May 28 '12 at 17:14
3  
@J.M. -- I searched the MMA documentation for "Longest Common Substring" and got tutorial/StringPatterns, tutorial/WorkingWithStringPatterns, and guide/SummaryOfNewFeaturesIn60. I didn't search for "Longest Common Subsequence," which would have found LongestCommonSubsequence[], so thanks for the lesson :) –  Reb.Cabin May 28 '12 at 20:20
    
amusingly enough, the wikipedia article on the longest common substring problem says that "it must not be confused with the longest common subsequence problem." IOW, wikipedia will send you to the dark corner of the MMA docs where you will have a hard time finding the function LongestCommonSubsequence, which is, in fact, a solution to Wikipedia's dijoint "longest common substring" problem. –  Reb.Cabin May 29 '12 at 23:01
add comment

4 Answers 4

up vote 19 down vote accepted

Mathematica supports two related functions, LongestCommonSequence[] and LongestCommonSubsequence[]. The first one finds the longest (contiguous or non-contiguous) sequence common to the two strings given as arguments to it:

LongestCommonSequence["AAABBBBCCCCC", "CCCBBBAAABABA"]
"AAABB"

while the second function is constrained to give the longest contiguous sequence:

LongestCommonSubsequence["AAABBBBCCCCC", "CCCBBBAAABABA"]
"AAAB"

These functions became available only in version seven; if you need to do this in an earlier version, István's routine is useful.

share|improve this answer
6  
It's worth mentioning the (undocumented, but System-context) LongestCommonSubsequencePositions and LongestCommonSequencePositions too. –  Szabolcs May 28 '12 at 17:01
6  
@Szabolcs LongestAscendingSequence[] is another interesting undocumented thingie –  belisarius May 28 '12 at 17:05
add comment

Preamble and motivation

While I am much late to the party here, I hope this answer will not be totally useless. This is a first in a series of posts where I will advocate a wider use of Java in our workflow, and present/describe certain toolset to reduce the mental overhead of this. So, my motivation here is not to provide a faster or more elegant solution, but to show that often, we can mindlessly reuse existing (found on the web or elsewhere) Java code, and the process can be made easy and painless.

Simplistic Java reloader

Here I will present a simplistic Java class reloader, which takes a string of Java code, attempts to compile it, and, upon success, load the resulting class into Mathematica via JLink. Note that I only so far tested it on Windows, but hopefully soon will test on other platforms as well.

Note that it is not my intention to present and describe a full workflow involving the reloader, in this post - I will save this for a future one. Here, I just present the code and an example of how it is useful for the case at hand.

Code

BeginPackage["SimpleJavaReloader`", {"JLink`"}];

JCompileLoad::usage = 
"JCompileLoad[javacode_,addToClassPath_] attempts to compile a Java \
class defined by  a string javacode, optionally adding to Java compiler classpath \
files and folders from addToClassPath, and load the resulting class into 
Mathematica";

Begin["`Private`"]

JCompileLoad::dirmakeerr = "Can not create directory `1`";

$stateValid = True;

$tempJavaDirectory =  FileNameJoin[{$UserBaseDirectory, "Temp", "Java"}];
    $tempClassDirectory = FileNameJoin[{$tempJavaDirectory, "Classes"}];
    $tempJavaLogDirectory = FileNameJoin[{$tempJavaDirectory, "Logs"}];
    $tempCompileLogFile =   FileNameJoin[{$tempJavaLogDirectory, "javac.log"}];
    $jrePath =   
     FileNameJoin[{$InstallationDirectory, "SystemFiles", "Java", $SystemID}];
$javaPath = FileNameJoin[{$jrePath, "bin"}];
    $jlibPath = FileNameJoin[{$jrePath, "lib"}];
    $classPath = {$tempClassDirectory, $jlibPath};


Scan[
   If[! FileExistsQ[#] && CreateDirectory[#] === $Failed,
          Message[JCompileLoad::dirmakeerr, #];
          $stateValid = False
   ] &,
   {
     $tempJavaDirectory,
         $tempClassDirectory,
     $tempJavaLogDirectory
   }];



(* determine a short name of the class from code *)
Clear[getClass];
getClass[classCode_String] :=
  With[{cl =
     StringCases[classCode, 
       "public" ~~ Whitespace ~~ "class"|"interface" ~~ Whitespace ~~ 
         name : (WordCharacter ..) :> name
     ]},
    First@cl /; cl =!= {}];

getClass[__] := Throw[$Failed, error[getClass]];

(* Determine the name of the package for the class *) 
Clear[getPackage];
getPackage[classCode_String] :=
  With[{pk = 
      StringCases[classCode, 
          ShortestMatch["package" ~~ Whitespace ~~ p__ ~~ ";"] :> p
      ]},
    First@pk /; pk =!= {}];

getPackage[classCode_String] := None;

getPackage[__] := Throw[$Failed, error[getPackage]];


ClearAll[getFullClass];
getFullClass[classCode_String] :=
   StringJoin[If[# === None, "", # <> "."] &@
      getPackage[classCode], getClass[classCode]];

(* Note: So far, tested on Windows only. Some specifics of quoting are 
   tuned to work around some bugs in Windows command line processor *)
Clear[makeCompileScript];
makeCompileScript[sourceFile_String] :=
  StringJoin[
    "\"",
    "\"", FileNameJoin[{$javaPath, "javac"}] , "\"",
        " -g ", sourceFile,
        " -d ", $tempClassDirectory,
    " -classpath ", "\"", Sequence @@ Riffle[$classPath, ";"], "\"",
        " 2> ", $tempCompileLogFile,
    "\""
  ];

Clear[getSourceFile];
getSourceFile[javacode_String] :=
   FileNameJoin[{$tempClassDirectory, getClass[javacode] <> ".java"}];

Clear[JCompileLoad];

JCompileLoad::invst =  "The loader is not on valid state. Perhaps some temporary \
     directories do not exist";

JCompileLoad::cmperr = "The following compilation errors were encountered: `1`";

JCompileLoad[javacode_String, addToClassPath_: {}]/; $stateValid :=
      Module[{sourceFile, fullClassName = getFullClass[javacode]},
         sourceFile = getSourceFile[javacode];
         With[{script =
            Block[{$classPath = Join[$classPath, addToClassPath]},
               makeCompileScript[sourceFile]
            ]},
           Export[sourceFile, javacode, "String"];
           If[Run[script] =!= 0,
             Message[
                JCompileLoad::cmperr, 
                Style[#, Red] &@Import[$tempCompileLogFile, "String"]
         ];
         $Failed,
             (*else*)
             ReinstallJava[];
             AddToClassPath @@ Join[$classPath, addToClassPath];
         LoadJavaClass[fullClassName]
       ]
     ]
  ];

JCompileLoad[_String, addToClassPath_: {}] :=
  (
     Message[JCompileLoad::invst];
     $Failed
  )

End[]

EndPackage[]

Notes

Note that you can either put this into a separate package file, or simply copy and paste into the FrontEnd, and run from there, for a quick test.

The package works by saving the string with your Java code into a temporary file, and then invoking Java compiler which comes with the JRE bundled with Mathematica, to compile this class. The compiled class is stored in another temporary location, from where it is then loaded by JLink. In case if compilation errors were encountered, the message generated by Java compiler is printed, and $Failed is returned.

One important limitation is that ReinstallJava is called to recompile / reload any class.

The case at hand

We will now apply the above to our case. First, we need the solution for the longest common substring problem in Java.

Stealing code from the web

I won't stay noble and code that myself. The whole point is - why doing so if we can steal it from someone :)? I get the one from here (the second one). In the simplified workflow I am presenting now, we need a string of Java code, so we define:

jlcsCode = 

"import java.lang.*;

 public class LCS{  
    public static String longestCommonSubstring(String S1, String S2) {
         int Start = 0;
         int Max = 0;
         for (int i = 0; i < S1.length(); i++){
            for (int j = 0; j < S2.length(); j++){
               int x = 0;
               while (S1.charAt(i + x) == S2.charAt(j + x)){
                  x++;
                  if (((i + x) >= S1.length()) || ((j + x) >= S2.length())) 
                     break;
               }
               if (x > Max) {
                  Max = x;
                  Start = i;
               }
            }
         }
         return S1.substring(Start, (Start + Max));
    }
 }";

Compiling and running

First, we have to compile and load this code:

JCompileLoad[jlcsCode]
 JavaClass[LCS,<>]

We are now ready to use the function, no other preparation needed! For example:

LCS`longestCommonSubstring["AAABBBBCCCCC","CCCBBBAAABABA"]
  "AAAB"

Note that, since the function longestCommonSubstring is a static method of the LCS class (meaning that the method belongs to the class rather than specific instance of it), we have to use the syntax className`method[args], and we don't have to create a class instance with JavaNew to use this. The class itself has to be loaded prior to being used, by JCompileLoad does that for us.

Now, some benchmarks:

s = StringJoin@RandomChoice[{"A", "C", "T", "G"}, 10000];
t = StringJoin@RandomChoice[{"A", "C", "T", "G"}, 10000];

Here we use the Mathematica's built-in function:

LongestCommonSubsequence[s,t]//AbsoluteTiming
{0.3232421,TCCACACGGGTAG}

Now our function:

LCS`longestCommonSubstring[s,t]//AbsoluteTiming
{1.1269531,TCCACACGGGTAG}

We see that our function is about 4 times slower, but, given that the Mathematica's built-in function was written in C and heavily optimized, while I just picked the first code snippet on the web I found, I think that the pain/gain ratio is pretty good.

Conclusions

I tried here to make a case for using Java in our workflow more frequently. The good thing about Java is that, in contrast to MathLink/LibraryLink, the JLink interface brings us pretty much all the way there, so there is no preparation at all necessary. The Java class reloader I presented here is very simplistic, but it nevertheless "closes the circle", and now we can protptype everything exclusively from Mathematica. I will expand on this in some future posts, and illustrate the workflow more fully. Note that I don't consider the reloader as anything complete - this is rather a proof of concept at this point.

For the case at hand, it took me literally 5 minutes from start to finish to get this working (the Java reloader I already had), and that includes finding the code on the web, pasting to Mathematica, compiling and using it. Given that there are many cases when Mathematica built-in functions are not available, while Java implementations are ubiquitous, I think this option can significantly expand our possibilities. Of course, to use this one needs some knowledge of Java, but let this not put you off: the things about Java you really need to know for such cases can be picked up in a day or two, especially if you have any experience in C/C++ (but even if not).

share|improve this answer
1  
I would love to bookmark this particular answer in SE. But I can't. +1+1+1... –  István Zachar Jun 3 '12 at 14:32
1  
Wouldn't that make perfect blog material? –  celtschk Jun 3 '12 at 17:12
1  
@celtschk Yes, I also thought about it, but I did not want to wait until we get a blog. I will in any case be able to refactor some of my posts later into a blog post. –  Leonid Shifrin Jun 3 '12 at 17:20
1  
@Leonid brilliant! Similar pattern could enable C# interop as easily, and I could imagine JavaScript interop via a service in Node.JS wouldn't be too far off, either. There is an incredible amount of JavaScript boilerplate out there, and much of it is easier on the mind of a Mathematica dev because JavaScript is fundamentally functional as opposed to fundamentally object-oriented like Java and C#. –  Reb.Cabin Jun 4 '12 at 17:03
1  
Clear[makeCompileScript]; makeCompileScript[sourceFile_String] := StringJoin["\"", "\"", FileNameJoin[{$javaPath, "javac"}], "\"", " -g ", "\"", sourceFile, "\"", " -d ", "\"", $tempClassDirectory, "\"", " -classpath ", "\"", Sequence @@ Riffle[$classPath, ";"], "\"", " 2> ", "\"", $tempCompileLogFile, "\"", "\"", "\""]; –  belisarius Aug 21 '12 at 4:34
show 42 more comments

Also, using pattern matching,just in case:

{{a, b, c, d, e, f, g}, {x, a, r, b, c, j}} /. {{___, Longest[y__], ___}, {___, y__, ___}} -> {y}
(*
-> {b, c}
*)

Edit

With this approach you can do one thing that seems not trivial by using the faster LongestCommonSequence[] function: finding the maximal common subsequence among several lists:

{{1, 2, 3, 4, 7, 8}, {1, 2, 3, 5, 7, 8}, {3, 4, 7, 5, 7, 8}} /. 
                {{___, Longest[y__], ___}, {___, y__, ___} ...} -> {y}
(*
->{7, 8}
*)
share|improve this answer
2  
Not intended for being used with really long sequences –  belisarius May 28 '12 at 17:25
add comment

These things I have coded up before Mathematica 7 and the introduction of the built in function LongestCommonSubsequence. The built-in version is of course faster, though still this implementation might be of interest as it has a bit wider functionality. Also, with some fine-tuning and compilation the performance can be surely increased.

longestCommonSubsequence[s, t] returns a set of the longest common continuous subsequences that can be found in lists s and t. longestCommonSubsequence lists all distinct continuous subsequences (one thing the built-in LongestCommonSubsequence is not capable of).

longestCommonSubsequence[s_List, t_List] := Module[
   {m = Length@s, n = Length@t, longest = 0, l, set = {}},
   l = Table[0, {m + 1}, {n + 1}];
   Do[
    If[s[[i]] === t[[j]],
      l[[i + 1, j + 1]] = l[[i, j]] + 1;
      If[l[[i + 1, j + 1]] > longest, longest = l[[i + 1, j + 1]]; 
       set = {}];
      If[l[[i + 1, j + 1]] === longest, 
       set = Union[set, {Take[s, {i - longest + 1, i}]}]];
      ];
    , {i, m}, {j, n}];
   set
   ];

Or you can just calculate the length, which is faster:

longestCommonSubsequenceLength[s, t] returns the length of the longest common continuous subsequence that can be found in lists s and t.

longestCommonSubsequenceLength[s_List, t_List] := 
  Module[{m = Length@s, n = Length@t, l, longest = 0},
   l = Table[0, {m + 1}, {n + 1}];
   Do[
    If[s[[i]] === t[[j]], l[[i + 1, j + 1]] = l[[i, j]] + 1];
    longest = Max[longest, l[[i + 1, j + 1]]];
    , {i, m}, {j, n}];
   longest
   ];

Example usage:

s = RandomChoice[{"A", "C", "T", "G"}, 200];
t = RandomChoice[{"A", "C", "T", "G"}, 200];

LongestCommonSubsequence[StringJoin@s, StringJoin@t]

"CATATTG"

longestCommonSubsequence[s, t]

{{"C", "A", "T", "A", "T", "T", "G"}, {"G", "T", "C", "A", "A", "T", "G"}}

Note that longestCommonSubsequence has found all instances of common subsequences.

longestCommonSubsequenceLength[s, t]

7

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.