tag:blogger.com,1999:blog-52790024313666850312018-06-16T00:25:45.996-07:00Finding PalindromesWhere do palindromes come from, where do they appear, en how do you find them?Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.comBlogger18125tag:blogger.com,1999:blog-5279002431366685031.post-14330814589126815562012-10-28T03:04:00.000-07:002012-10-28T03:04:18.972-07:00The history of finding palindromes<p>I have given several programs to find palindromes in the blog posts on this blog. Who was the first to construct those? The history of finding palindromes is not as long as the history of palindromes themselves, since, obviously, computers did not exist in the time of the ancient Greeks. Using a computer to find palindromes has probably been done since the 1970s. Before that, palindromes have been used to illustrate the power of computing models. In this blogpost I will go through the history of finding palindromes using a computer. </p><p>The oldest references to papers in computer science I could find mentioning palindromes (or, actually, `strings satisfying the symmetry predicate') were published in 1965 in Russian by Freivald and Bārzdiņš. Since I cannot read Russian, I depend on the description of these papers in a review from Albert Mullin in the Journal of Symbolic Logic in 1970. Both Freivald and Bārzdiņš show how to determine whether or not a string is a palindrome using different variants of a Turing machine. Turing machines are named after Alan Turing, who was born exactly 100 years ago. <div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-x2KVUnQr9K0/UENPPZbIuBI/AAAAAAAAAMU/cOSn4zcItg0/s1600/Alan_Turing-728x409.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="180" width="320" src="http://3.bp.blogspot.com/-x2KVUnQr9K0/UENPPZbIuBI/AAAAAAAAAMU/cOSn4zcItg0/s320/Alan_Turing-728x409.jpg" /></a></div> In the 1930s, Turing devised a model of a machine to prove a fundamental mathematical theorem about the <i>Entscheidungsproblem</i>. The gist of this theorem is that you cannot solve all problems by means of a computer. In particular, you cannot use a computer to determine if an arbitrary program returns it result in a reasonable amount of time. The Turing machine model has had a profound influence on real computing machines built in the 1940s and 1950s. It has also been used a lot in Computer Science, the scientific field that arose since the construction of the first computers. For a long time, the fields of algorithms and complexity theory, which study algorithms and their efficiency, have used Turing machines to show the efficiency of particular algorithms. Modern computers and Turing machines are very dissimilar, and the complexity of software is often determined by factors not present in Turing machines. Probably for that reason, Turing machines are not used much anymore in Computer Science, although they still appear as examples of computing models. </p><p>The most basic form of a Turing machine has a <i>tape</i> on which it can write symbols, a <i>table</i> that tells which actions to take, using a <i>head</i> that can move along the tape, and a <i>state register</i>, which contains the different states the machine can be in, such as a start state, an end state, or a state in between. The machine can read a symbol from the type, write a symbol on the tape, and move the head along the tape. Freivald and Bārzdiņš show that to test whether an input word is a palindrome requires a number of steps that is quadratic in the length of the input string on such a Turing machine. Eight years later in 1973, Slisenko showed that if you are allowed to use more than one head on a tape, or multiple tapes, than a palindrome can be determined in a number of steps linear in the length of the input string. The English translation of the paper `Recognizing a symmetry predicate by multihead Turing machines with input' is a 183 page, dense mathematical text. <div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-hoOIgNc1nlI/UENqrKYvlyI/AAAAAAAAAMo/UbBFylPyvms/s1600/slisenko.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="200" width="131" src="http://1.bp.blogspot.com/-hoOIgNc1nlI/UENqrKYvlyI/AAAAAAAAAMo/UbBFylPyvms/s200/slisenko.jpg" /></a></div>Slisenko announced the result already in 1969, and given the form and the length of his paper (with 183 pages this is more a book than a paper), I find it highly likely that it took him a couple of years to write it. Slisenko's result was a surprisingly strong result, and people started to study and trying to improve upon it. One of these people was Zvi Galil, then a postdoc IBM Yorktown Heights. <div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-bY0zqicpfj8/UENyc6CBkoI/AAAAAAAAAM8/Dwy-h_XEi-Y/s1600/galil.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="200" width="133" src="http://2.bp.blogspot.com/-bY0zqicpfj8/UENyc6CBkoI/AAAAAAAAAM8/Dwy-h_XEi-Y/s200/galil.jpg" /></a></div>He had a hard time understanding Slisenko's paper, which I fully understand, but in 1978 he obtained the same result, only he needed just 18 pages to describe it. The problem of finding palindromes efficiently started off an area within Computer Science now known as stringology, which studies strings and their properties, and algorithms on strings. </p><p>Freivald, Bārzdiņš, Slisenko, and Galil describe algorithms for recognizing palindromes on Turing machines with a varying number of tapes. On other machine models, palindromes have been a popular example too. Stephen Cole showed how to recognize palindromes on iterative arrays of finite-state machines in 1964. Alvy Ray Smith III used cellular automata to recognize palindromes in 1972. Also in the 1970s, Freivald and Yao, independently, looked at recognizing palindromes by means of probabilistic Turing machines. Looking at these results from a distance of around 40 years, it seems like it was a sport to study programming with various restrictions, like not using your left hand, or tying your legs together. </p><p>Some of the machine models on which algorithms for palindromes were developed have rather artificial restrictions, which gives rather artificial algorithms for finding palindromes. But the algorithms on some of the more realistic machine models, such as the Random Access Machine (RAM) model, contain the essential components of the algorithm I give in the blog post on <a href="http://finding-palindromes.blogspot.nl/2012/05/finding-palindromes-efficiently.html">finding palindromes efficiently</a>. Manacher gives a linear-time algorithm on the RAM computing model finding the smallest initial palindrome of even length. He also describes how to adjust his algorithm in order to find the smallest initial palindrome of odd length longer than 3. He did not realize that his approach could be used to find all maximal palindromes in a string in linear time. Zvi Galil and Joel Seiferas did, in 1976. They wrote a paper, titled `Recognizing certain repetitions and reversals within strings', in which they develop an algorithm that finds all maximal palindromes in a string in linear time. As far as I know, this is the first description of the algorithm I give in the blog post on <a href="http://finding-palindromes.blogspot.nl/2012/05/finding-palindromes-efficiently.html">finding palindromes efficiently</a>. </p><p>Twelve years later I rediscovered this algorithm. I published a paper on `The derivation of on-line algorithms, with an application to finding palindromes', in which I show how to obtain the efficient algorithm for finding palindromes from <a href="http://finding-palindromes.blogspot.nl/2012/04/naive-algorithm-for-finding-palindromes.html">the naive algorithm for finding palindromes</a> using algebraic reasoning. The method at which I arrive at the algorithm is completely different from the way Galil and Seiferas present their version. I presented the algorithm and the method I use to construct it to several audiences. I was only 22 years old at the time, and I am afraid my presentation skills were limited. Several people that saw me presenting the efficient algorithm for finding palindromes thought they could do better. An example can be found in the book `Beauty is our business' (which isn't about models, or escort services, but about Computer Science, and dedicated to the Dutch Turing award winning computer scientist Edsger Dijkstra on his sixtieth birthday), in which Jan Tijmen Udding derives the same algorithm using Dijkstra's calculus for calculating programs. </p><p>Zvi Galil developed algorithms for almost all problems related to palindromes in his series of papers on palindromes. In 1985, he even developed an algorithm for finding palindromes using a parallel machine. Turing machines, and the other machine models mentioned above, are sequential machines: computations are performed in sequence, and not at the same time. Parallel machines allow computations to be performed at the same time, in parallel. In the previous century few parallel machines were available, but nowadays, even the laptop on which I am typing this text has four cores, and can do many things in parallel. The importance of parallel machines has increased considerably, also because it is expected that increasing the speed of computers by increasing the clock speed is not going to work anymore in the next couple of years. Extra power of computers has to come from the possibilities offered by using multiple cores for computations. To put these cores to good use, we need efficient parallel algorithms for the problems we want to solve. Apostolico, Breslauer, and Galil improved upon Galil's first parallel algorithm for finding palindromes in their paper on `Parallel detection of all palindromes in a string' in 1995. When finding palindromes in DNA, such parallel algorithms might be very useful. I am not aware of any experience report on using parallel algorithms to find palindromes in DNA, but I hope to try this out myself soon, and report on it on this blog. </p><p>Since the discovery that palindromes appear in DNA, people working on bioinformatics have developed algorithms for finding palindromes in DNA. The first description of these algorithms I could find are in the book Algorithms on Strings, Trees and Sequences: Computer Science and Computational Biology, by Dan Gusfield. <div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-qalX8doH05c/UIw73BotVbI/AAAAAAAAANo/7qJBzOm81hA/s1600/gusfield.jpg" imageanchor="1" style="clear:left; float:left; margin-right:1em; margin-bottom:1em"><img border="0" height="151" width="110" src="http://4.bp.blogspot.com/-qalX8doH05c/UIw73BotVbI/AAAAAAAAANo/7qJBzOm81hA/s200/gusfield.jpg" /></a></div>This is one of the standard works in computational biology. Gusfield shows how to compute all maximal palindromes in a string using a suffix tree. He then moves on to discuss gapped palindromes (called separate palindromes in Gusfield's book), and approximate palindromes, for which he leaves it to the reader to construct an algorithm that returns approximate palindrome in time linear in the product of the maximum number of errors and the length of the input. After Gusfield, several other scientist have worked on finding gapped and approximate palindromes, sometimes improving on or generalizing Gusfield's results. For example, as I already mentioned in the blog post on <a href="http://finding-palindromes.blogspot.nl/2012/09/palindromes-with-errors-or-gaps.html">palindromes with errors or gaps</a>, Porto and Barbosa have developed an efficient algorithm that finds approximate palindromes in which also symbols may be inserted or deleted. </p><p>Are all problems related to finding palindromes solved? I think many are, and for most practical purposes, efficient algorithms for finding palindromes are available. From an algorithm-design perspective I think there are still some open questions. The central concept I use in the design of the efficient algorithm for finding palindromes is the <a href="http://finding-palindromes.blogspot.nl/2012/05/palindromes-in-palindromes.html">palindromes in palindromes</a> property. I think this property should also play a central role in efficient algorithms for finding gapped and approximate palindromes. I have tried for quite a while to design algorithms for these problems using the palindromes in palindromes concept, but failed. Anyone? </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com1tag:blogger.com,1999:blog-5279002431366685031.post-55686860497568941962012-09-25T04:51:00.000-07:002012-10-23T11:54:28.132-07:00Finding palindromes in the Y chromosome<p>To test the program I obtain by combining palindrome finding with errors and gaps, I downloaded a copy of the Y chromosome. The first paper that described the occurrences of huge palindromes in the male chromosome referred to in my <a href="http://finding-palindromes.blogspot.nl/2012/07/palindromes-in-dna.html">blog post on palindromes in DNA</a> is `The male-specific region of the human Y chromosome is a mosaic of discrete sequence classes'. This paper appeared in Nature in 2003, and the Nature web pages provide a link to the Y chromosome studied in the paper. The paper caused quite a stir. It has been cited more than a thousand times in other scientific papers. A thousand citations is a lot: none of my direct colleagues in computer science ever wrote an article with so many citations. </p><p>Since it is computationally infeasible to test the program for finding palindromes with errors and gaps with large numbers of allowed errors, I looked up the information about palindromes in the paper. The paper contains a table that gives information about the eight palindromes the authors found in chromosome Y. </p><table><tr><td>Palindrome </td><td>Arm length (kb) </td><td>Arm-to-arm identity (%) </td><td>Spacer length (kb) </td><td>Palindrome span (kb) </td></tr><tr><td>P1 </td><td>1,450 </td><td>99.97 </td><td>2.1 </td><td>2,902 </td></tr><tr><td>P1.1 </td><td>9,9 </td><td>99.95 </td><td>3.9 </td><td>24 </td></tr><tr><td>P1.2 </td><td>9,9 </td><td>99.95 </td><td>3.9 </td><td>24 </td></tr><tr><td>P2 </td><td>122 </td><td>99.97 </td><td>2.1 </td><td>246 </td></tr><tr><td>P3 </td><td>283 </td><td>99.94 </td><td>170 </td><td>736 </td></tr><tr><td>P4 </td><td>190 </td><td>99.98 </td><td>40 </td><td>419 </td></tr><tr><td>P5 </td><td>496 </td><td>99.98 </td><td>3.5 </td><td>996 </td></tr><tr><td>P6 </td><td>110 </td><td>99.97 </td><td>46 </td><td>266 </td></tr><tr><td>P7 </td><td>8.7 </td><td>99.97 </td><td>12.6 </td><td>30 </td></tr><tr><td>P8 </td><td>36 </td><td>99.997 </td><td>3.4 </td><td>75 </td></tr></table><p>Palindrome P1 consists of almost 3 million symbols, and the arm-to-arm identity is 99.97 percent. An arm has length almost one and a half million, so the number of errors is around 450. I assume that the errors are more or less evenly distributed in the palindrome, which implies that in the central 10,000 symbols of an arm of a palindrome, after the gap, about 3 errors would occur. So instead of trying to find the long palindrome with possibly 450 errors and a gap of size around 2100 symbols, I instead try to find palindromes with a gap of size 2200 (to be on the safe side) of length at least 5000 (to be on the safe side) with at most 5 errors (to be on the safe side). To my surprise, I do not find a single palindromes that satisfies these requirements. This long palindrome with a gap of size 2100 does not seem to appear in the Y chromosome of the authors. Since P2 has the same gap size, it does not appear either. Hmm, what's wrong here? Have I made a mistake? </p><p>Lets have a look at the other palindromes reported in the table. The next smallest gap size reported in the above table is 3.4 kb. If I try to find palindromes with a gap of size 3400, I get two hits. Around position 18905096 I find a palindrome of length 32056 (including the gap) if I allow for 8 errors. If I try to reduce the gap length, I find that I can take a gap length of 2320, and still find the same palindrome. So maybe this is one of the palindromes with a gap of around 2.1 kb? The arm length is (32056-2320)/2 = 14868, which doesn't correspond at all with the reported arm lengths. I also find a palindrome around position 9806955. This palindrome has length 17304, and 5 errors. This palindrome indeed has a gap of size 3400: reducing the gap size leads to much shorter palindromes. But the palindrome around this position is much shorter than the reported length in Nature. I experimented with different gap sizes and numbers of allowed errors, and found the following palindromes in the Y chromosome: <table><tr><td>Center </td><td>Length </td><td>Errors </td><td>Gap length </td></tr><tr><td>9806955 </td><td>13904 </td><td>5 </td><td>3400 </td></tr><tr><td>12078263 </td><td>43356 </td><td>1 </td><td>151498 </td></tr><tr><td>13738581 </td><td>34452 </td><td>0 </td><td>39262 </td></tr><tr><td>14445793 </td><td>64738 </td><td>5 </td><td>132950 </td></tr><tr><td>17989037 </td><td>14934 </td><td>5 </td><td>297114 </td></tr><tr><td>18899672 </td><td>3364 </td><td>7 </td><td>39540 </td></tr><tr><td>18905096 </td><td>29736 </td><td>8 </td><td>2320 </td></tr><tr><td>20530881 </td><td>24374 </td><td>2 </td><td>157760 </td></tr></table></p><p>Since my findings are very different from the results reported in Nature, I contacted the group responsible for the paper, the Whitehead Institute at the Department of Biology at the MIT. The research scientist that answered my questions was kind enough to provide me with the positions at which the palindromes occur in the Y chromosome. Given the positions, I tried to compare the arms by hand (or eye, actually). I got nowhere: the start and end of the arms did not resemble each other at all. </p><p>The final clue I needed was that it is not enough to discard errors constituted by two symbols that do not match in a palindrome. Sometimes I also have to delete symbols in one arm of a palindrome, or, equivalently, add symbols to the other arm. With this last clue, I selected the reported palindromic arms in the Y chromosome, and compared them. It turns out that in each case I only have to delete or insert a couple of DNA symbols to obtain a palindrome. The table of palindromes I find are the palindromes P8 to P1 reported in Nature, with P7 missing. My seventh palindrome does not appear in the list of positions I received, but this is probably one of P1.1 or P1.2, of which I didn't receive the positions. </p><p>After months of playing with the Y chromosome in some of my spare time, I now finally find the palindromes I was after. If I would have to do the same thing again for another DNA string, I would follow the same approach: <ul><li> use the program that finds palindromes with gaps and errors with various gap sizes and numbers of allowed errors to find centers of palindrome candidates</li><li> hand align the sequences around the centers of palindrome candidates to determine the length of the palindromes around these centers, given some number of allowed deletions or insertions</li></ul>I don't think another approach is feasible given the current speed of computers. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com1tag:blogger.com,1999:blog-5279002431366685031.post-66974320584680982402012-09-12T13:19:00.000-07:002012-09-12T13:19:30.515-07:00Palindromes with errors or gaps<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script><!-- > import Debug.Trace > import qualified Data.ByteString as B > import Data.ByteString.Internal as BI > import Data.Word (Word8) > import Data.Char (toUpper) > import Data.List (intercalate) > (=:=) :: Word8 -> Word8 -> Bool > l =:= r = let cl = toUpper (BI.w2c l) > cr = toUpper (BI.w2c r) > in if cl `elem` "ATCG" && cr `elem` "ATCG" > then cl == negateDNA cr > else False > negateDNA :: Char -> Char > negateDNA 'A' = 'T' > negateDNA 'T' = 'A' > negateDNA 'C' = 'G' > negateDNA 'G' = 'C' > negateDNA _ = error "negateDNA: not a DNA character" > lengthPalindromeAroundDNA input lastPos start end > | start < 0 || end > lastPos = end-start-1 > | B.index input start =:= B.index input end = > lengthPalindromeAroundDNA > input > lastPos > (start-1) > (end+1) > | otherwise = end-start-1 > approximatePalindromesDNA :: Int -> B.ByteString -> String > approximatePalindromesDNA k input = > intercalate "\n" > $ map (showPalindromeDNA input) > $ zip (dnaApproximatePalindromesAroundCentres input k) [0..] > dnaApproximatePalindromesAroundCentres :: ByteString -> Int -> [Int] > dnaApproximatePalindromesAroundCentres input k = > let centers = [0 .. B.length input] > lastPos = B.length input - 1 > in map (dnaLengthApproximatePalindromeAround input lastPos k) centers > dnaLengthApproximatePalindromeAround > :: B.ByteString -> Int -> Int -> Int -> Int > dnaLengthApproximatePalindromeAround input lastPos k center > = lengthApproximatePalindromeAroundDNA input (B.length input -1) k (center-1) center > showPalindromeDNA :: B.ByteString -> (Int,Int) -> String > showPalindromeDNA input (len,pos) = > let startpos = pos - len `div` 2 > in (show startpos ++) > . (" to " ++) > . (show (startpos+len) ++) > . ("\t" ++) > . (show (B.take len $ B.drop startpos input) ++) > . ("\t" ++) > $ show len > palindrome s = s == reverse s > input = B.pack (map c2w "AATAATTCGAATTCGCGAAAA") --><p>Sometimes I not only want to find perfect palindromes, but also palindromes that contain a limited number of errors. A palindrome with a limited number of errors is often called an <i>approximate</i> palindrome. For example, in the book Judges in the King James Bible, verse 19:9 reads: <pre><br /> And when the man rose up to depart, he, and his <br /> concubine, and his servant, his father in law, the <br /> damsel's father, said unto him, Behold, now the day <br /> draweth toward evening, I pray you tarry all night: <br /> behold, the day groweth to an end, lodge here, that <br /> thine heart may be merry; and to morrow get you early <br /> on your way, that thou mayest go home.<br /></pre>The substring "draweth toward" is a text palindrome with one error: the `e' and the `o' don't match. To be precise, a string <code>s</code> is an approximate palindrome with <code>k</code> errors, if it satisfies the predicate <code>approximatePalindrome k s</code>. I define <code>approximatePalindrome</code> as a program by taking the first half of both the input string and its reverse, and comparing each of the elements of these strings pairwise. I don't have to compare the second halves of these strings, since that would give the same result as the comparison of the first halves, and mismatches would be counted twice. <pre><br /><br />> approximatePalindrome k s = <br />> let half = div (length s) 2<br />> approximatePalindrome' k [] [] = k>=0<br />> approximatePalindrome' k (x:xs) (y:ys) <br />> | x == y = approximatePalindrome' k xs ys<br />> | k >= 0 = approximatePalindrome' (k-1) xs ys<br />> | otherwise = False<br />> in approximatePalindrome' k <br />> (take half s) <br />> (take half (reverse s)) <br /><br /></pre></p><p>A palindrome with a gap is a palindrome in which a gap of a particular size in the middle is ignored. An example of a palindrome with a gap is found in Revelations, where verses 20:7-8 read: <pre><br /> And when the thousand years are expired, Satan shall <br /> be loosed out of his prison, And shall go out to <br /> deceive the nations which are in the four quarters of <br /> the earth, Gog, and Magog, to gather them together to <br /> battle: the number of whom is as the sand of the sea.<br /></pre>Here "Gog, and Magog" is a text palindrome with a gap of length three in the middle: the `n' and the `M' around the central `d' don't match. Since the gap appears in the middle of the string, the length of the gap is odd if the length of the palindrome is odd, and even if the length of the palindrome is even. To be precise, a string <code>s</code> is a palindrome with a gap of length <code>g</code> in the middle, if it satisfies the predicate <code>gappedPalindrome g s</code>: <pre><br /><br />> gappedPalindrome g s = <br />> let ls = length s<br />> removeGap g s = <br />> let armLength = div (ls - g) 2<br />> (before,rest) = splitAt armLength s<br />> (gap,after) = splitAt g rest<br />> in before ++ after <br />> in if g < ls && odd g == odd ls <br />> then palindrome (removeGap g s)<br />> else error "gappedPalindrome: incorrect gap length"<br /> <br /></pre>If the gap is shorter than the length of the input string, and both the gap and input string have even or odd length, I remove the gap from the middle of the string, and check that the remaining string is a palindrome. I remove the gap by taking half of the elements of the input string, minus half of the gap elements, then dropping the gap, and then taking the rest of the string. </p><p>To find palindromes with errors or gaps, I adapt my software for finding palindromes. </p><p>It is quite easy to adapt the quadratic-time algorithm for finding palindromes to also find approximate palindromes, and palindromes with gaps. In this blog post I show how to adapt the quadratic-time algorithm that finds palindromes in DNA. Adapting the quadratic-time algorithms for finding exact or text palindromes is done in exactly the same way. To find approximate palindromes in DNA, I pass an extra argument <code>k</code>, representing the maximum number of errors allowed, to the function <code>qMaximalPalindromesDNA</code>, which I now call <code>approximatePalindromesDNA</code>. The only function I have to adapt is the function <code>lengthPalindromeAroundDNA</code>: <pre><br /><br />> lengthApproximatePalindromeAroundDNA <br />> input lastPos k start end <br />> | start < 0 || end > lastPos = end-start-1<br />> | B.index input start =:= B.index input end = <br />> lengthApproximatePalindromeAroundDNA <br />> input <br />> lastPos <br />> k <br />> (start-1) <br />> (end+1) <br />> | k > 0 = <br />> lengthApproximatePalindromeAroundDNA <br />> input <br />> lastPos <br />> (k-1) <br />> (start-1) <br />> (end+1) <br />> | otherwise = end-start-1<br /><br /></pre>The time complexity of <code>approximatePalindromesDNA</code> is related to the sum of the lengths of the approximate palindromes found in the input. As I described in the blog post on <a href="http://finding-palindromes.blogspot.nl/2012/08/finding-palindromes-in-dna.html">finding palindromes in DNA</a>, this is not a problem for perfect palindromes in DNA. In the case of approximate palindromes, the number of allowed errors has a big influence on the time complexity. If I increase the number of allowed errors by one, the length of all palindromes increases with at least two, since a single error represents two symbols, and since after the error there might be matching symbols. For example, if I allow for 500 errors to find a palindrome consisting of almost $3$ million symbols, each approximate palindrome in the input has at least length $1000$, except for the palindromes at the start and the end of the string. The average approximate palindrome is quite a bit longer, since the chance for an accidental match is one in four. Since each palindrome has length on average more than $1250$, applying the function <code>approximatePalindromesDNA</code> to a DNA string of length $25$ million requires more than $30$ billion steps only to compute the lengths of the approximate palindromes. This number is very challenging for my laptop, and until now I haven't been able to convince my machine to calculate these values for me, despite letting it run for several nights. </p><p>To find palindromes with gaps in DNA, I use an extra argument <code>g</code>, which represents the length of the gap in the middle of the palindrome which I can ignore. A gap may be odd or even, and in the former case I get a palindrome around a center <i>on</i> a symbol in a DNA string, instead of in between two symbols. So where an exact palindromes in a DNA string has its center in between two symbols, a gapped palindrome may have a center in between two symbols or on a symbol, depending on whether or not the gap is even. The function <code>qMaximalPalindromesDNA</code>, which is now called <code>gappedMaximalPalindromes</code>, is the only function I have to adapt. For each center, I have to subtract half of the gap of the center for the start position, and add half of the gap to the end position. When doing so, I have to check that these positions do not point before the beginning of the string or after the end of the string. If they do, I have to replace half of the gap by the appropriate smaller number. <pre><br /><br />> gappedMaximalPalindromesDNA input g =<br />> let halfg = div g 2<br />> lengthInput = B.length input<br />> centers = if even g <br />> then [0 .. lengthInput]<br />> else [0 .. lengthInput-1]<br />> halfg' c | c < halfg = c<br />> | c+halfg > lengthInput = lengthInput-c<br />> | otherwise = halfg<br />> left c = c-1-halfg' c<br />> right c = if even g<br />> then c+halfg' c<br />> else c+1+halfg' c<br />> in map <br />> (\c -> lengthPalindromeAroundDNA <br />> input <br />> (lengthInput-1)<br />> (left c)<br />> (right c)<br />> ) <br />> centers<br /><br /></pre></p><p>By combining the functions <code>lengthApproximatePalindromeAroundDNA</code> and <code>gappedMaximalPalindromesDNA</code> I obtain a function for determining palindromes with both gaps in the middle, and errors in the palindromic arms. </p><p>Adapting the quadratic-time algorithm for finding palindromes to also allow for errors and gaps in palindromes is not very hard. How about adapting the linear-time algorithm for finding palindromes? </p><p>Adapting the linear-time algorithm to account for errors and gaps is <i>very</i> hard. I spent a long time studying the problem of finding palindromes with a given maximum number of errors in time linear in some combination of the number of errors and the length of the input, but failed to find a satisfying solution. During a couple of long trips I drew many ideas in my notebook to adapt the main components of the linear-time algorithm, the longest tail palindrome and the list of maximal palindromes with centers before the longest tail palindrome, to allow for errors and gaps. In vain. </p><p>I tried to adapt the linear-time algorithm for finding palindromes to also find approximate palindromes before I had realized that the quadratic-time algorithm returns its results much faster on DNA strings than the linear-time algorithm. My frustrations about not being able to find a solution have almost disappeared now. </p><p>In the scientific literature I found two solutions to the problem of finding approximate palindromes in time linear in some combination of the input and the number of allowed errors. The first solution uses so-called suffix trees. A suffix tree is an advanced data structure that can be used to find all kinds of structures in strings, among which palindromes. I compared the efficiency of my linear-time algorithm for finding exact palindromes with an algorithm using suffix trees, and my linear-time algorithm was far superior. The solutions for finding approximate palindromes using suffix trees use even more arguments, and they will run even slower than the algorithm for finding exact palindromes in strings using suffix trees. The quadratic-time solution for finding approximate palindromes with gaps is no more complicated than the above programs, so although I did not perform any tests, I am pretty confident that the above quadratic-time algorithm for finding approximate palindromes runs much faster than the solutions using suffix trees when applied to DNA strings. The second solution uses advanced algorithms on graphs to determine edit actions that turn one string into another string. I'm afraid I don't understand the details sufficiently well to discuss the algorithm here. The time-complexity of this algorithm is related to $k^2 \times n$, where $k$ is the number of allowed errors, and $n$ the length of the input string. For the DNA string of length $25$ million and the $500$ errors used as examples above, this would amount to more than $6$ trillion steps. My machine would just laugh at it. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com2tag:blogger.com,1999:blog-5279002431366685031.post-29381784773551704432012-08-16T02:40:00.001-07:002012-08-16T02:40:49.306-07:00Finding palindromes in DNA<p>Palindromes play an important role in DNA, and many geneticists want to find palindromes in the DNA they have. To find palindromes in DNA, I cannot directly apply the algorithm for finding palindromes described in a previous blog post, since the equality used for palindromes in DNA differs from the normal equality. A string is a palindrome in DNA if each symbol in the reversed string is the negation of the symbol in the original string, where `A' is the negation of `T' and vice versa, and `C' is the negation of `G', and vice versa. The operator <code>=:=</code> introduced in my blog post on <a href="http://finding-palindromes.blogspot.no/2012/04/what-is-palindrome-ii-specification-in.html">What is a palindrome II? A specification in Haskell</a> implements this equality. By negating the result of <code>=:=</code> I get the inequality operator <code>=/=</code> on DNA symbols. </p><p>A first attempt to find palindromes in DNA is to adapt the various components to deal with DNA in the functions <code>extendPalindrome</code>, <code>moveCenter</code>, and <code>finalPalindromes</code> introduced in the blog post on <a href="http://finding-palindromes.blogspot.no/2012/05/finding-palindromes-efficiently.html">finding palindromes efficiently</a>. In the function <code>extendPalindrome</code> the expression <code>a!... /= a!...</code>, where the dots should be replaced by the appropriate indices in the input, is replaced by the expression <code>a!... =/= a!...</code> using the DNA inequality operator. This is sufficient to turn the function <code>extendPalindrome</code> into a function that finds palindromes in DNA. The functions <code>moveCenter</code> and <code>finalPalindromes</code> do not use the (in)equality operator, so you might think that this is all that is needed to obtain a program for finding palindromes in DNA. </p><p>There is one more location in the algorithm for finding palindromes where I use a property of palindromes that does not hold for DNA. In the function <code>moveCenter</code> I call function <code>extendPalindrome</code> with a new "current palindrome" of length one. But a DNA string of length one is not a palindrome, since if I reverse it, I get the same string, and to be a palindrome, the string has to contain the complement symbol. Every DNA string of odd length is not a palindrome, since when such a string is reversed, the symbol in the middle is compared against itself, and DNA equality fails. To find palindromes in a DNA string, I only have to look at the even positions in the string, which correspond to the positions in between symbols. The palindromes around odd positions on the symbols always have length zero. </p><p>I adapt the algorithm for efficiently finding palindromes to only find palindromes around even positions. Besides adapting the function <code>extendPalindrome</code> with a different inequality operator, I adapt the function <code>moveCenter</code> by calling <code>extendPalindrome</code> with a new longest current palindrome of length zero, and by counting the number of centers down with steps of size two instead of one, since I am not interested in the odd center positions anymore. These are all changes needed, and I will not give the resulting algorithm. The interested reader can find them in <a href="http://www.staff.science.uu.nl/~jeuri101/homepage/palindromes/index.html">my software package for finding palindromes</a>. </p><p>In spring 2012 I was contacted by an Indian scientist who was developing an algorithm for finding palindromes in DNA. I soon found out that her algorithm was similar to my <a href="http://finding-palindromes.blogspot.nl/2012/05/naive-algorithm-for-finding-maximal.html">naive, quadratic-time, algorithm for finding maximal palindromes</a>. What worried me was that her algorithm seemed to run faster than my efficient algorithm for finding palindromes, adapted to finding palindromes in DNA. By thinking that my efficient algorithm would run faster, I had made the classical mistake of confusing worst-case time complexity with practice. </p><p>The worst-case time complexity of my efficient algorithm for finding palindromes is linear. The number of steps the algorithm takes is linear in the length of the input string. The worst case time-complexity of the naive algorithm is quadratic. A string in which all substrings are palindromes leads to "worst case" time behaviour. An example of a DNA string satisfying this property is a string looking like "ATATATAT...". For a string consisting of three million copies of "AT", the efficient algorithm returns the length of all maximal palindromes about equally fast as for the the Amycolatopsis marina, which is about equally long, but the quadratic time algorithm never finishes on the repetitions of ATs. Strings like "ATATATAT..." appear in DNA, but they are short. For example, the longest maximal palindrome in chromosome 18 consisting of just "ATAT..." has length 66. The longest overall maximal palindrome has length 76. The majority of the maximal palindromes is much shorter. The time-complexity of the naive algorithm for finding palindromes is the sum of the length of the palindromes found in the argument string. If the maximum length of a palindrome is 76, then the number of steps performed by the naive algorithm is at most 76 times the length of the input string. Since chromosome 18 consists of almost 75 million symbols, this is more than five and a half billion, which is quite a lot, but doable for a modern laptop. Indeed, if I measure the running time of the quadratic and linear-time solutions on chromosome 18, I find that the quadratic solution is about four to five times faster than the linear solution. To be precise, on my MacBook Pro with 2.4 GHz Intel Core i5 with 4GB memory, with very few other processes running, the linear-time solution requires slightly more than 30 minutes to analyse the almost 75 megabytes of chromosome 18, whereas the quadratic solution uses around 6 minutes. Why is it faster? </p><p>It is hard to exactly determine what one program makes faster than another. The linear time solution uses quite a few arguments, in particular the list of maximal palindromes found until now. This is the central idea to obtain linear-time behaviour, but it also requires constructing, passing on, inspecting, and deconstructing this list to determine maximal palindromes with centers to the right of the current longest tail palindrome. I have experimented with several kinds of representations for the list of maximal palindromes used as argument to the functions <code>extendPalindrome</code> and <code>moveCenter</code>, but none of these solutions can compare to the quadratic solution when looking at running times on large DNA strings. An important motivation for developing a linear-time algorithm for finding palindromes turns out not to be valid at all... </p><p>The quadratic-time solution for finding palindromes in DNA is very similar to the naive algorithm for finding maximal palindromes. There are three differences. First, the input to the function <code>maximalPalindromes</code>, which I now call <code>qMaximalPalindromesDNA</code>, is of type <code>ByteString</code> instead of <code>Array Int Char</code>. The type <code>ByteString</code> is a data structure specifically designed to efficiently deal with strings of which the length is known. It is very similar to the type <code>String</code>, but functions like <code>read</code> for reading strings, <code>map</code> for mapping a function over a string, <code>length</code> for determining the length of a stringdefined in the module <code>ByteString</code>, all are usually much faster than their equivalents on <code>String</code>, or even on <code>Array Int Char</code>. I write <code>B.ByteString</code>, <code>B.length</code>, and so on, to make explicit that these components come from the module <code>ByteString</code>: <pre><br /><br />> import qualified Data.ByteString as B<br /><br /></pre><!-- > import Data.ByteString.Internal as BI > import Data.Word (Word8) > import Data.Char (toUpper) > (=:=) :: Word8 -> Word8 -> Bool > l =:= r = let cl = toUpper (BI.w2c l) > cr = toUpper (BI.w2c r) > in if cl `elem` "ATCG" && cr `elem` "ATCG" > then cl == negateDNA cr > else False > negateDNA :: Char -> Char > negateDNA 'A' = 'T' > negateDNA 'T' = 'A' > negateDNA 'C' = 'G' > negateDNA 'G' = 'C' > negateDNA _ = error "negateDNA: not a DNA character" -->Second, in the function <code>qMaximalPalindromesDNA</code> I only consider positions in between symbols around which to find maximal palindromes, since palindromes in DNA have even length. The length of the list of <code>centers</code> has length one plus the length of the input, to account for the center at the end of the input. <pre><br /><br />> qMaximalPalindromesDNA :: B.ByteString -> [Int]<br />> qMaximalPalindromesDNA input = <br />> let lengthInput = B.length input<br />> centers = [0 .. lengthInput]<br />> lastPos = lengthInput - 1<br />> in map <br />> (\c -> lengthPalindromeAroundDNA <br />> input <br />> lastPos <br />> (c-1) <br />> c<br />> ) <br />> centers<br /><br /></pre>The third difference is that function <code>lengthPalindromeAroundDNA</code> doesn't distinguish between even-length and odd-length palindromes anymore, since all palindromes have even length, and that it needs to use the equality operator for DNA. <pre><br /><br />> lengthPalindromeAroundDNA input lastPos start end <br />> | start < 0 || end > lastPos = end-start-1<br />> | B.index input start =:= B.index input end = <br />> lengthPalindromeAroundDNA <br />> input <br />> lastPos <br />> (start-1) <br />> (end+1) <br />> | otherwise = end-start-1<br /><br /></pre></p><p>This solves the problem of finding palindromes in DNA. But I cannot yet use this algorithm to find the huge palindromes in the Y chromosome, since those palindromes have some noise in the middle, and the palindromic arms contain a few errors. I will look into this in a following blog post. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com2tag:blogger.com,1999:blog-5279002431366685031.post-48762835902694525762012-07-22T14:38:00.000-07:002012-07-22T14:38:24.642-07:00Palindromes in DNA<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>DNA strings are often millions of letters long. For example, the copy I have of the DNA of the Amycolatopsis marina, a bacteria discovered in 2009, consists of $6,503,724$ characters. <div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-dWMSDZBWAMk/T_fgQYG4SRI/AAAAAAAAALo/WWCNJakV6mI/s1600/Actinomyces_israelii.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="130" width="200" src="http://3.bp.blogspot.com/-dWMSDZBWAMk/T_fgQYG4SRI/AAAAAAAAALo/WWCNJakV6mI/s200/Actinomyces_israelii.jpg" /></a></div>Bacteria are of course some of the smallest living beings. The DNA in the copy I have of the human chromosome 18 is an astounding $74,657,229$ characters long. Chromosome 18 is one of the 23 chromosomes of human beings, and is estimated to contain in between $300$ and $500$ genes. A gene is responsible for part of our functioning as human beings. An example of a gene contained within chromosome 18 is the NPC1 gene, named after the Niemann-Pick disease, type C1. This gene is named after the disease you get when you have an error on the gene. Experiments with mice show that this gene probably controls the appetite, and people with a mutation on the NPC1 gene often suffer from obesitas. Interestingly, the same mutation on NPC1 might make you immune for the ebola virus, one of the deadliest known viruses. So the mutation is partly a blessing in disguise. <br/><br/><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-TbiQ5h9vOOA/T_f3pNz6PII/AAAAAAAAAL4/2EGSMpaFIe8/s1600/Protein_NPC1_PDB_3GKH.png" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="166" width="200" src="http://4.bp.blogspot.com/-TbiQ5h9vOOA/T_f3pNz6PII/AAAAAAAAAL4/2EGSMpaFIe8/s200/Protein_NPC1_PDB_3GKH.png" /></a></div> If I search for the keyword palindrome in the electronic publications available at the library of my university, I get more than $500$ hits. The first ten of these hits are all about palindromes in DNA. My guess is that at least $90$% of the list of publications are about palindromes in DNA. So what is the interest in palindromes in DNA? Surely, if your strings only contain `A', `T', `C', and `G' you are likely to get many palindromes? </p><p>Let us look at how many palindromes we expect to find in the DNA of the Amycolatopsis marina. Suppose I want to find palindromes of length twelve. I calculate the chance that an arbitrary DNA string of length twelve is a palindrome as follows. The first six characters of the string don't matter. The seventh character needs to match with the sixth character, for which we have a chance of one in four. Remember that in DNA, `A' matches with `T' and vice versa, but both `A' and `T' do not match with themselves, or with `C' and `G'. The eighth character needs to match with the fifth character, for which we also have a chance of one in four. This goes on until the twelfth character, so I get a chance of \[ \frac{1}{4} \times \frac{1}{4} \times \frac{1}{4} \times \frac{1}{4} \times \frac{1}{4} \times \frac{1}{4} = (\frac{1}{4})^6 = \frac{1}{4^6} = \frac{1}{4096} \] Since the Amycolatopsis marina is $6,503,724$ characters long, there are $6,503,713$ substrings of length twelve. Multiplying this number with the chance that it is a palindrome, I expect to get $1,589$ palindromes. Using my algorithm for finding palindromes I get $1,784$ palindromes. This is slightly above $10$% more than expected, but that might be an accident. If I look at the palindromes of length fourteen in the $18$th human chromosome, using a similar calculation I expect to find $4,556$ palindromes. My algorithm finds $25,323$. More than five times as many palindromes as expected! This is not an accident anymore. Palindromes play a role in DNA. But what role? </p><p>Palindromes perform several tasks in human DNA. I will discuss one particularly intriguing task in this blog post. </p><p>Why do we humans have sex? I can answer this question from several perspectives. The most common answer will probably mention love, pleasure, or children. Looking at the question from a biological perspective, the answer discusses the biological advantages of having sex. Our children get their genes from two parents. They get two complete copies of human DNA, and merge these copies in some way to make them the way they are. In this merging process, damaged DNA from one parent can be replaced by functioning DNA from the other parent. There is a lot of DNA to repair: every day, each cell may be damaged at one million different places. Most of this damage is harmless, since a lot of DNA does not seem to have any function at all, but some kinds of damage may cause a lot of problems. Being able to repair non-functioning DNA when combining DNA from parents is essential for keeping the human race in a good state. The American geneticist Hermann Joseph Muller used this argument to explain why sexual reproduction is favored over asexual reproduction in organisms. When an organism reproduces asexually it passes all its DNA errors on to its offspring, without a possibility to repair them, eventually leading to the extinction of the population. The process has been dubbed Muller's ratchet, after the ratchet device, which can only turn in one direction. This is the theory, practice is slightly more complicated. </p><p>Muller's ratchet should already be at work in humans, since there is one human chromosome in which no combination of parental genes takes place: chromosome 23. Chromosome 23 determines whether we become a man or a woman. A man gets a copy of the chromosome called X from his mother and a copy called Y from his father. A woman gets an X from her mother and an X from her father. A woman can merge both X's and repair possible errors, but a man has two different copies of the chromosome, and has no possibility to combine, let alone repair, them. The Y is passed on from father to son, with no involvement of women. Muller's ratchet theory says the genes on the chromosome that make a man should deteriorate quickly, and men should soon become extinct. </p><p>There is no sign of men becoming extinct. Apparently there are some other mechanisms at work to repair DNA on the Y chromosome. If I cannot obtain a copy of some piece of DNA from my mother, maybe I can store copies of the DNA in the Y chromosome itself? If I maintain two copies, I can always check if a piece of DNA is correct by comparing it against the copy of this piece of DNA. This is the mechanism used by the Y chromosome, where the copies are stored as palindromes, with some noise in the middle of these palindromes. Such palindromes with gaps in the middle are often called inverted repeats in biology. The Y chromosome contains eight huge palindromes, the longest of which consists of almost three million characters. Around 25% of the Y chromosome consists of palindromes. The DNA in the palindromes carries genes for describing the male testes. So the mechanism by means of which men survive is called palindrome... </p><p> </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-2404226628012499362012-07-01T12:18:00.001-07:002014-02-05T13:01:30.472-08:00Other implementations and solutions<p>Next year it is 25 years ago since I constructed an algorithm for finding palindromes efficiently. I was 22 years old then, and had just started as a PhD student. I was quite excited about having solved the problem of finding palindromes efficiently, and wrote a <a href="http://www.staff.science.uu.nl/~jeuri101/homepage/Publications/algorithmicapalindromes.pdf">scientific paper</a> about it. My excitement wasn't shared by the scientific community though. In the last twenty-five years this paper has been cited less than ten times, and appears at the bottom end of my most cited papers list. </p><p>In 2007 we developed the <a href="http://en.wikipedia.org/wiki/ICFP_Programming_Contest">ICFP programming contest</a>. The ICFP programing contest is a very challenging contest, in which thousands of programmers try to show off their programming talents in their programming language of choice. Our contest asked participants to transform the left picture below to the right picture using as few commands as possible. We included a problem related to palindromes in our contest, since this was my pet-problem ever since 1988. After the contest I explained how to solve the palindrome problem efficiently in a <a href="http://johanjeuring.blogspot.se/2007/08/finding-palindromes.html">blog post</a>. </p><a href="http://2.bp.blogspot.com/-ljuPcSCihpA/T-ohj_MK7ZI/AAAAAAAAAK8/Ebr4ufBk94o/s1600/source.png" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="200" width="200" src="http://2.bp.blogspot.com/-ljuPcSCihpA/T-ohj_MK7ZI/AAAAAAAAAK8/Ebr4ufBk94o/s200/source.png" /></a><a href="http://2.bp.blogspot.com/-jRmjZpKkaTk/T-ohsXAMEaI/AAAAAAAAALI/8hKqooTf5ts/s1600/target.png" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="200" width="200" src="http://2.bp.blogspot.com/-jRmjZpKkaTk/T-ohsXAMEaI/AAAAAAAAALI/8hKqooTf5ts/s200/target.png" /></a><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><br/><p>When some years later I looked at the number of hits the contest pages received, I found that each month, thousands of people are reading the blog message about finding palindromes. Why does this page attract so many visitors? </p><p>The palindrome problem is a common question in interviews for jobs for software developers, so the blog post attracts software developers looking for a new job, and preparing themselves for an interview. Another reason the blog post attracts visitors is that I think quite a few people are genuinely interested in the problem and its solution. The concept of palindromes appears in almost all languages, which means that the question of finding palindromes is asked all over the world. The blog post indeed attracts visitors from all over the world. The last 100 (July 1, 2012) visitors come from all continents except Australia. </p><p>Many people ask the question of how to find palindromes, but also many people try to answer the question. You can find hundreds of solutions for finding palindromes on the internet. Some of these are variants of my linear-time solution, others are more naive quadratic or sometimes even cubic-time solutions. Below I give the list I found, ordered on the programming language used for the implementation. If there exists a linear-time implementation, I don't list less efficient solutions. <center><table border="0" cellspacing="15"> <tr><td><a href="http://www.careercup.com/question?id=245679">C</a></td><td>The same linear-time solution as my blog post. </td></tr> <tr><td><a href="http://www.leetcode.com/2011/11/longest-palindromic-substring-part-ii.html">C++ 1</a></td><td>The same linear-time solution as my blog post. This post has an extensive description of the palindromes in palindromes property, including some pictures which try to explain it. </td></tr> <tr><td><a href="http://wcipeg.com/wiki/index.php?title=Longest_palindromic_substring">C++ 2</a></td><td>A C++ implementation of Manacher's algorithm. </td></tr> <tr><td><a href="https://gist.github.com/fpelliccioni/8828315">C++ 3</a></td><td>The same linear-time solution as my blog post in C++ by Fernando Pelliccioni. </td></tr> <tr><td><a href="http://efreedom.com/Question/1-2677000/Find-Longest-Palindrome-Given-String">C#</a></td><td>As far as I can see, this is a quadratic-time solution. </td></tr> <tr><td><a href="http://re-factor.blogspot.com/2010/10/longest-palindrome.html">Factor</a></td><td>A cubic-time solution. </td></tr> <tr> <td><a href="http://programmingpraxis.com/2010/10/15/find-the-longest-palindrome-in-a-string/">F#</a></td><td>A quadratic-time solution. </td></tr> <tr> <td><a href="http://programmingpraxis.com/2010/10/15/find-the-longest-palindrome-in-a-string/">Go</a></td><td>A quadratic-time solution, by Lars Björnfot. </td></tr> <tr><td><a href="http://rampion.blogspot.se/2012/01/what-i-did-on-my-winter-vacation-or-ro.html">Haskell</a></td><td>Just as the program for finding palindromes described in my blog post, this blog post describes a linear-time program for finding palindromes in Haskell. The interesting aspect of this solution is that it returns results lazily: as soon as it finds a maximal palindrome, it writes it to the output. </td></tr> <tr><td><a href="http://alg-code.blogspot.com/2011/04/find-longest-palindromes.html">Java 1</a></td><td>A quadratic-time solution. </td></tr> <tr><td><a href="http://stackoverflow.com/questions/1115001/write-a-function-that-returns-the-longest-palindrome-in-a-given-string">Java 2</a></td><td>Another quadratic-time solution in Java, by Marcello de Sales. </td></tr> <tr><td><a href="http://stackoverflow.com/questions/1115001/write-a-function-that-returns-the-longest-palindrome-in-a-given-string">Java 3</a></td><td>Yet another quadratic-time solution in Java, by Mohit Bhandari. </td></tr> <tr> <td><a href="http://programmingpraxis.com/2010/10/15/find-the-longest-palindrome-in-a-string/">Matlab</a></td><td>A quadratic-time solution, by Lalit Mohan. </td></tr> <tr> <td><a href="http://stackoverflow.com/questions/1115001/write-a-function-that-returns-the-longest-palindrome-in-a-given-string">PHP 1</a></td><td>As far as I can see this is a quadratic-time solution, by Marc Donaldson. </td></tr> <tr> <td><a href="http://www.biophp.org/minitools/find_palindromes/">PHP 2</a></td><td>This is a quadratic- or cubic-time solution, by Joseba Bikandi. You can try it <a href="http://www.biophp.org/minitools/find_palindromes/demo.php">on-line</a>. </td></tr> <tr> <td><a href="http://www.akalin.cx/longest-palindrome-linear-time">Python</a></td><td>The same linear-time solution as my blog post, developed by Fred Akalin. This post contains an extensive description of the palindromes in palindromes property too. </td></tr> <tr> <td><a href="http://svitsrv25.epfl.ch/R-doc/library/Biostrings/html/findPalindromes.html">R</a></td><td>This page describes the interface of functionality for finding palindromes in DNA strings. It doesn't say how efficient the software is. </td></tr> <tr> <td><a href="http://matthewkerry.blogspot.co.uk/2012/02/ruby-longest-palindrome-in-string.html">Ruby 1</a></td><td>A quadratic-time solution, by Matthew Kerry. </td></tr> <tr> <td><a href="http://talklikeaduck.denhaven2.com/2011/03/13/algorithm-choice-trumps-programming-language-choice-for">Ruby 2</a></td><td>Another quadratic-time solution in Ruby, by Rick DeNatale. The post and the comments mention some more Ruby implementations. </td></tr> <tr> <td><a href="http://mfcoding.wordpress.com/2011/11/14/longest-palindrome-in-a-string/">Ruby 3</a></td><td>Yet another quadratic-time solution in Ruby, by Mitchell Fang. </td></tr> <tr> <td><a href="http://programmingpraxis.com/2010/10/15/find-the-longest-palindrome-in-a-string/">Scheme</a></td><td>A quadratic-time solution. </td></tr> </table></center> I could only find linear-time implementations for finding palindromes in C, C++, Haskell, and Python. Please let me know if you find better or alternative implementations for finding palindromes. This list easily gets outdated, so please keep me informed. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com9tag:blogger.com,1999:blog-5279002431366685031.post-29554620492123481372012-06-29T04:57:00.000-07:002013-07-16T23:06:04.100-07:00The Sator square<p>The Sator square is a square that reads the same forwards, backwards, upwards and downwards: <center><table><tr><td>S</td><td>a</td><td>t</td><td>o</td><td>r</td></tr><tr><td>A</td><td>r</td><td>e</td><td>p</td><td>o</td></tr><tr><td>T</td><td>e</td><td>n</td><td>e</td><td>t</td></tr><tr><td>O</td><td>p</td><td>e</td><td>r</td><td>a</td></tr><tr><td>R</td><td>o</td><td>t</td><td>a</td><td>s</td></tr></table></center><br/>The square has been found at many places around the world: on amulets, carved into walls or stones, written on disks used for extinguishing fire, and so on. Its first known appearance is in the ruins of Pompeii. The square was found in a building that was decorated in a style that became popular since the year 50. Pompeii was covered in ash from the Vesuvius in the year 79, and it follows that this example is almost 2000 years old. This makes the square one of the oldest known palindromes.<br/> <a href="http://4.bp.blogspot.com/-8NDeUfS8w-w/T-moZ6YxcvI/AAAAAAAAAKs/ZGgp7vFwegc/s1600/Sator_Square_at_Opp%25C3%25A8de.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="198" width="200" src="http://4.bp.blogspot.com/-8NDeUfS8w-w/T-moZ6YxcvI/AAAAAAAAAKs/ZGgp7vFwegc/s200/Sator_Square_at_Opp%25C3%25A8de.jpg" /></a></p><p>What does the Sator square mean? Except for `Arepo', scholars mostly agree on the meaning of the individual words. `Sator' is the seeder, sower, or begetter. It is also used as a metaphor for God, not necessarily Christian, but also Roman. `Tenet' is a verb, and means it, he or she holds. `Opera' is often considered to mean with effort. It is related to opus or opera, which means work. Finally, `Rotas' is generally considered a noun, meaning the wheels. But it might also be a verb meaning turn. Nobody has ever found a word in Roman, Greek, Etruscan, or any Indo-European language that explains `Arepo'. Some people think that it is a name, others that it stands for plough, misspelled to fit the palindrome, and yet others that it is just nonsense. For a long period, people thought the square had a Christian meaning by rearranging the letters into two Pater Nosters in the form of a cross. The Pre-Christian Pompeii find crushed this theory. If Arepo is a name, the Sator square might read: `The sower Arepo holds the wheels with effort', and if it stands for plough, you might get: `God holds the plough, but you turn the furrows'. According to <a href="http://www.johntcullen.com/">John Cullen</a>, this could have been a motto for farmers in Roman times. Most explanations I have read call the Sator square an ancient meme, with as much meaning as the sentence `all your base are belong to us'. In 1937, the Italian Antonio Ferrua probably gives the best explanation for what the square means: `Esattamente quello che si vuole'(!) E basta di questo argumento (it means exactly what you want it to mean. And so much for that argument!). The Sator square is an example of a meme that went viral long before the internet. <div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-McXoBJVaJAg/T-yVt3puorI/AAAAAAAAALY/Gb4SHzK9xqw/s1600/SatorSkara.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="168" width="200" src="http://4.bp.blogspot.com/-McXoBJVaJAg/T-yVt3puorI/AAAAAAAAALY/Gb4SHzK9xqw/s200/SatorSkara.jpg" /></a></div>Here is an instance of the square I found in the library in Skara, Sweden, written by an unknown author probably in Stockholm, Sweden, in the year 1722. It starts with Rotas instead of Sator, just as the first appearances of the square. For some reason the version starting with Sator has become more popular. </p><p>Because it is a palindrome, the Sator square was thought to have many healing effects, curing snake-bites, headaches, jaundice, and many other illnesses. Medieval books mention the square as a cure for fever and insanity. Interestingly, although we now know that using words to cure an illness is of little help, we humans do use palindromes to repair our body. I will say more about this in a later blog post on palindromes in DNA. </p><p>The amount of human effort gone into explaining the meaning of the Sator square is unbelievable. Since 1889, when Francis John Haverfield described a find of the square in a Romano-Bristish building in Cirencester, there has been a steady stream of articles on the Sator square, and the total number of articles easily surpasses a hundred. Rose Mary Sheldon recently published a 54-page annotated biography of the literature on the Sator square: The Sator rebus: An unsolved cryptogram? Charles Douglas Gunn wrote his PhD thesis on the square at Yale in 1969. He suggests that the square was written by a Roman who wanted to take palindromic squares one step further from the misformed four-letter word square Roma tibi subito montibus ibit amor, meaning `For by my efforts you are about to reach Rome, the object of your travel'. He wrote software to generate all possible five-letter Latin word squares. These squares take up more than a hundred pages in his thesis. He concludes that the Sator square is the best. </p><p>There are not many other examples of squares like the Sator square, consisting of five words that together constitute a meaningful palindrome. Piet Burger discovered the following Dutch square in 1990. <center><table><tr><td>S</td><td>n</td><td>e</td><td>e</td><td>d</td></tr><tr><td>n</td><td>e</td><td>t</td><td>t</td><td>e</td></tr><tr><td>e</td><td>t</td><td>s</td><td>t</td><td>e</td></tr><tr><td>e</td><td>t</td><td>t</td><td>e</td><td>n</td></tr><tr><td>d</td><td>e</td><td>e</td><td>n</td><td>s</td></tr></table></center><br/>Its meaning is so far-fetched that I won't try to give an English translation. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com9tag:blogger.com,1999:blog-5279002431366685031.post-89494288115584823012012-05-31T23:11:00.000-07:002013-02-11T03:10:21.595-08:00Finding palindromes efficiently<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>Using the <a href="http://finding-palindromes.blogspot.se/2012/05/palindromes-in-palindromes.html">palindromes in palindromes</a> property, I now develop an algorithm for finding palindromes that requires a number of steps approximately equal to the length of its input. This linear-time algorithm can be used to find palindromes in documents of any size, even in DNA strings. Finding palindromes in a string of length $20,000,000$ using this algorithm requires a number of seconds on a modern laptop. It is impossible to find palindromes significantly faster, unless you have a machine with multiple cores, and use a parallel algorithm. </p><p>The program for efficiently finding palindromes is only about 25 lines long, but for formatting purposes it extends to over 80 lines in this blog post. Although the program is short, it is rather intricate. I guess that you need to experiment a bit with to find out how and why it works. </p><p>The <a href="http://finding-palindromes.blogspot.se/2012/05/naive-algorithm-for-finding-maximal.html">naive algorithm for finding palindromes</a> discussed in a previous blog post uses the following top-level definition for finding maximal palindromes: <pre><br /><br />< maximalPalindromes :: Array Int Char -> [Int]<br />< maximalPalindromes a = <br />< let (first,last) = bounds a<br />< centers = [0 .. 2*(last-first+1)]<br />< in map (lengthPalindromeAround a) centers<br /><br /></pre>The reason why this algorithm is naive is that <code>lengthPalindromeAround</code> calculates the maximal palindrome around a center independently of the palindromes calculated previously. I now change this by calculating the maximal palindromes from left to right around the centers of a string. In this calculation I either extend a palindrome around a center, or I move the center around which I determine the maximal palindrome rightwards because I have found a maximal palindrome around a center. So I replace the above definition by <pre><br /><br />> maximalPalindromes :: Array Int Char -> [Int]<br />> maximalPalindromes a = <br />> let (first,last) = bounds a<br />> in reverse (extendPalindrome a first 0 [])<br /><br /></pre>Before I introduce and explain function <code>extendPalindrome</code>, I give an example of how the algorithm works. Suppose I want to find the maximal palindromes in the string "yabadabadoo". The algorithm starts by finding the maximal palindrome around the position in front of the string, which cannot be anything else than the empty string. It moves the position around which to find the maximal palindrome one step to point to the `y'. The maximal palindrome around this position is "y", since there is no character in front of it. It again moves the position around which to find palindromes one step to point to the position in between `y' and `a'. Since `y' and `a' are different, the maximal palindrome around this position is the empty string. Moving the center to `a', it finds that "a" is the maximal palindrome around this center, since `y' and `b' are different. The maximal palindrome around the next center in between `a' and `b' is again the empty string. Moving the center to `b', it can extend the current longest palindrome "b" around this center, since both before and after `b' it finds an `a'. It cannot further extend the palindrome "aba", since `y' and `d' are different. To determine the maximal palindrome around the center in between `b' and `a', the next center position, it uses the fact that "aba" is a palindrome, and that it already knows that the maximal palindrome around the center in between `a' and `b' is the empty string. Using the palindromes in palindromes property, it finds that the maximal palindrome around the position in between `b' and `a' is also the empty string, without having to look at the `b' and the `a'. To determine the maximal palindrome around the next center position on the last `a' of "aba", it has to determine if `d' equals `b', which it doesn't of course. Also here it uses the palindromes in palindromes property. Since "a" is the maximal palindrome around the center of the first `a' in "aba", and it reaches until the start of the palindrome "aba", I have to determine if the palindrome "a" around the second `a' can be extended. I won't describe all steps <code>extendPalindrome</code> takes in detail, but only give one more detail I already described above: the second occurrence of the palindrome "aba" in "yabadabadoo" is not found by extending the palindrome around its center, but by using the palindromes in palindromes property to find "aba" a second time in "abadaba". </p><p>Function <code>extendPalindrome</code> takes four arguments. The first argument is the array <code>a</code> in which we are determining maximal palindromes. The second argument is the position in the array directly after the longest palindrome around the current center (the longest palindrome around the center before the first symbol has length $0$, so the position directly after the empty palindrome around the first center is the <code>first</code> position in the array). I will call this the current rightmost position. The third argument is the length of the current longest palindrome around that center (starting with <code>0</code>), and the fourth and final argument is a list of lengths of longest palindromes around positions before the center of the current longest tail palindrome, in reverse order (starting with the empty list <code>[]</code>). It returns the list of lengths of maximal palindromes around the centers of the array, in reverse order. Applying function <code>reverse</code> to the result gives the maximal palindromes in the right order. The function <code>extendPalindrome</code> maintains the invariant that the current palindrome is the longest palindrome that reaches until the current rightmost position. </p> <p>There are three cases to be considered in function <code>extendPalindrome</code>. <ul><li>If the current position is after the end of the array, so <code>rightmost</code> is greater than <code>last</code>, I cannot extend the current palindrome anymore, and it follows that it is maximal. This corresponds to the <code>null after</code> case in the definition of <a href="http://finding-palindromes.blogspot.se/2012/04/maximal-palindromes.html">maximal palindromes</a>. It only remains to find the maximal palindromes around the centers between the current center and the end of the array, for which I use the function <code>finalPalindromes</code>. </li><li>If the current palindrome extends to the start of the array, or it cannot be extended, it is also maximal, and I add it to the list of maximal palindromes found. This case corresponds to the <code>null after || last before /= head after</code> case in the definition of <a href="http://finding-palindromes.blogspot.se/2012/04/maximal-palindromes.html">maximal palindromes</a>. I then determine the maximal palindrome around the following center by means of the function <code>moveCenter</code>. </li><li>If the element at the current rightmost position in the array equals the element before the current palindrome I extend the current palindrome. </li> </ul>In the following Haskell code, I list the arguments to a function below the function, since they don't fit on a single line anymore. I have also included some comments in the code, which are the pieces of text appearing after <code>--</code>. <pre><br /><br />> extendPalindrome <br />> a <br />> rightmost <br />> currentPalindrome <br />> currentMaximalPalindromes<br />> | rightmost > last = <br />> -- reached the end of the array<br />> finalPalindromes <br />> currentPalindrome <br />> currentMaximalPalindromes <br />> (currentPalindrome:currentMaximalPalindromes)<br />> | rightmost-currentPalindrome == first ||<br />> a!rightmost /= a!(rightmost-currentPalindrome-1) = <br />> -- the current palindrome extends to the start <br />> -- of the array, or it cannot be extended<br />> moveCenter <br />> a <br />> rightmost <br />> (currentPalindrome:currentMaximalPalindromes) <br />> currentMaximalPalindromes <br />> currentPalindrome <br />> | otherwise = <br />> -- the current palindrome can be extended<br />> extendPalindrome <br />> a <br />> (rightmost+1) <br />> (currentPalindrome+2) <br />> currentMaximalPalindromes <br />> where (first,last) = bounds a<br /><br /></pre>In two of the three cases, function <code>extendPalindrome</code> finds a maximal palindrome, and goes on to the next center by means of function <code>finalPalindromes</code> or <code>moveCenter</code>. In the other case it extends the current palindrome, and moves the rightmost position one further to the right. </p><p>Function <code>moveCenter</code> moves the center around which the algorithm determines the maximal palindrome. In this function I make essential use of the <a href="http://finding-palindromes.blogspot.se/2012/05/palindromes-in-palindromes.html">palindromes in palindromes</a> property. It takes the array as argument, the current rightmost position in the array, the list of maximal palindromes to be extended, the list of palindromes around centers before the center of the current palindrome, and the number of centers in between the center of the current palindrome and the rightmost position. It uses the palindromes in palindromes property to calculate the longest palindrome around the center after the center of the current palindrome. <ul><li>If the last center is on the last element, there is no center in between the rightmost position and the center of the current palindrome. I call <code>extendPalindrome</code> with rightmost position one more than the previous position, and a current palindrome of length 1. </li><li>If the previous element in the list of maximal palindromes reaches exactly to the left end of the current palindrome, I use the palindromes in palindromes property of palindromes to find the next current palindrome using <code>extendPalindrome</code>. </li><li>In the other case, I have found the longest palindrome around a center, add that to the list of maximal palindromes, and proceed by moving the center one position, and calling <code>moveCenter</code> again. I only know that the previous element in the list of maximal palindromes does not reach exactly to the left end of the current palindrome, so it might be either shorter or longer. If it is longer, I need to cut off the new maximal palindrome found, so that it reaches exactly to the current rightmost position. </li></ul><pre><br /><br />> moveCenter <br />> a <br />> rightmost <br />> currentMaximalPalindromes <br />> previousMaximalPalindromes <br />> nrOfCenters<br />> | nrOfCenters == 0 = <br />> -- the last center is on the last element: <br />> -- try to extend the palindrome of length 1<br />> extendPalindrome <br />> a <br />> (rightmost+1) <br />> 1 <br />> currentMaximalPalindromes<br />> | nrOfCenters-1 == head previousMaximalPalindromes = <br />> -- the previous maximal palindrome reaches <br />> -- exactly to the end of the last current <br />> -- palindrome. Use the palindromes in palindromes <br />> -- property to extend the current palindrome<br />> extendPalindrome <br />> a <br />> rightmost <br />> (head previousMaximalPalindromes) <br />> currentMaximalPalindromes<br />> | otherwise = <br />> -- move the center one step. Add the length of <br />> -- the longest palindrome to the maximal<br />> -- palindromes<br />> moveCenter <br />> a <br />> rightmost <br />> (min (head previousMaximalPalindromes) <br />> (nrOfCenters-1)<br />> :currentMaximalPalindromes) <br />> (tail previousMaximalPalindromes) <br />> (nrOfCenters-1)<br /><br /></pre>In the first case, function <code>moveCenter</code> moves the rightmost position one to the right. In the second case it calls <code>extendPalindrome</code> to find the maximal palindrome around the next center, and in the third case it adds a maximal palindrome to the list of maximal palindromes, and moves the center of the current palindromes one position to the right. </p><p>Function <code>finalPalindromes</code> calculates the lengths of the longest palindromes around the centers that come after the center of the current palindrome of the array. These palindromes are again obtained by using the palindromes in palindromes property. Function <code>finalPalindromes</code> is called when we have reached the end of the array, so it is impossible to extend a palindrome. We iterate over the list of maximal palindromes, and use the palindromes in palindromes property to find the maximal palindrome at the final centers. As in the function <code>moveCenter</code>, if the previous element in the list of maximal palindromes reaches before the left end of the current palindrome, I need to cut off the new maximal palindrome found, so that it reaches exactly to the end of the array. <pre><br /><br />> finalPalindromes <br />> nrOfCenters <br />> previousMaximalPalindromes <br />> currentMaximalPalindromes <br />> | nrOfCenters == 0 = currentMaximalPalindromes<br />> | otherwise =<br />> finalPalindromes <br />> (nrOfCenters-1) <br />> (tail previousMaximalPalindromes) <br />> (min (head previousMaximalPalindromes) <br />> (nrOfCenters-1)<br />> :currentMaximalPalindromes)<br /><br /></pre>In each step, function <code>finalPalindromes</code> adds a maximal palindrome to the list of maximal palindromes, and moves on to the next center. </p><p>I have discussed the number of steps this algorithm takes for each function. At a global level, this algorithm either extends the current palindrome, and moves the rightmost position in the array, or it extends the list of lengths of maximal palindromes, and moves the center around which we determine the maximal palindrome. If the length of the input array is $n$, the number of steps the algorithm is $n$ for the number of moves of the rightmost position, plus $2n+1$ for the number of center positions. This is a linear-time algorithm. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com6tag:blogger.com,1999:blog-5279002431366685031.post-66853643432000897422012-05-22T23:12:00.000-07:002012-05-22T23:12:11.918-07:00Sotades<p>His name has been mentioned already a couple of times: <i>Sotades</i> reportedly wrote some of the first palindromes. <div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-JB9Y5cXN9-A/T69h0TlF9KI/AAAAAAAAAJY/R0xHWyIaCSg/s1600/sotades2.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="152" width="162" src="http://4.bp.blogspot.com/-JB9Y5cXN9-A/T69h0TlF9KI/AAAAAAAAAJY/R0xHWyIaCSg/s200/sotades2.jpg" /></a></div></p> <p>Sotades was an Ancient Greek poet, living in the third century before Christ. He lived in Alexandria during the reign of Ptolemy II Philadelphus (285 BC-246 BC). Ptolemy II was born in Kos in about 309 BC. As a youth, he enjoyed the best tutors. The practice of getting the best scholars or poets available to educate the crown prince was something that Ptolemy II's father (also called Ptolemy, which explains the II in his son's name), had the occasion to observe in Macedonia, where the young Alexander was taught by no less a figure than Aristotle himself. Ptolemy II and his father set up the Alexandrian University and Library, and attracted the most learned men of the times to Alexandria. Ptolemy II provided these scholars with a life free from want and taxes, allowing them to study, write, perform research, and lecture in their respective disciplines. </p> <p>Ptomely II was married to Arsinoe (I), and got three children with her. He had also a sister called Arsinoe (II), who was married to a Greek king. Sister Arsinoe was a scheming woman. She had the oldest son of her husband, by a different woman, murdered to position her own children for the throne. When her husband was killed in the battle that ensued after the widow of the murdered son fled to a neighboring country, she married her half-brother (also called Ptolemy), who claimed the same throne as her husband had had. When her new husband became too powerful, she and her sons conspired against him. They were found out. Two sons got killed, and she and her eldest son (also called Ptolemy) fled to her full brother Ptolemy II in Egypt. In Eqypt she managed to get rid of Ptolemy's first wife Arsinoe I by accusing her of conspiring against her husband. Arsinoe I was sent in exile, and Arsinoe II married her brother.<div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-ihRfINAKcRw/T6_NNsg_HxI/AAAAAAAAAJo/olOBhzCZfuI/s1600/Oktadrachmon_Ptolemaios_II_Arsinoe_II.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="189" width="200" src="http://3.bp.blogspot.com/-ihRfINAKcRw/T6_NNsg_HxI/AAAAAAAAAJo/olOBhzCZfuI/s200/Oktadrachmon_Ptolemaios_II_Arsinoe_II.jpg" /></a></div></p><p>Alexandria not only attracted scholars, but also poets and satirists, such as Sotades. One of his poems attacked Ptolemy II's marriage to his sister Arsinoe, from which came the infamous line: "You're sticking your prick in an unholy hole." For this, Sotades was imprisoned, but he escaped to the island of Kaunos, where he was captured by Ptolemy's admiral, shut up in a leaden chest, and thrown into the sea. </p> <p>Sotades wrote obscene verses, in particular about the love between an older man and an adolescent boy. He invented his own metre for such verses, the sotadic metre. This metre spread around the world, and was used by other satirists to hint at an obscene or satirical interpretation. For example, Ennius wrote his `Sota' already in the same century in southern Italy. Sotades' verses could be read both forwards and backwards, where the backwards reading was the same, or contained an obscene interpretation. Unfortunately, none of the palindromic verses of Sotades survived. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-29653533795628798662012-05-15T03:54:00.000-07:002012-05-15T03:54:38.365-07:00Palindromes in palindromes<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>In a <a href="http://finding-palindromes.blogspot.se/2012/05/naive-algorithm-for-finding-maximal.html">previous blog post</a>, I develop a quadratic time algorithm for finding the length of the maximal palindromes around all centers of a string. Although it is a lot faster than <a href="http://finding-palindromes.blogspot.se/2012/04/naive-algorithm-for-finding-palindromes.html">the naive algorithm that finds all palindromes</a> in a string, it is still not very useful when trying to find palindromes in DNA, or in a text as long as the Bible. </p> <p>To design an efficient algorithm for finding palindromes, I reuse as much information as possible when calculating the length of maximal palindromes around centers. In this blog post, I will explain the central idea that allows me to reuse previously calculated maximal palindromes. </p> <p>In the naive algorithm for finding the maximal palindromes around all centers of a string, I start with calculating the maximal palindrome around the first, leftmost, center, then the palindrome around the next center on the first character, then the palindrome around the center in between the first two characters, and so on. I calculate the palindromes from left to right. In the illustration below I depict a string with the letters blackened out. But you may assume that it represents the string "yabadabadoo". In this string, <code>p</code> is the maximal palindrome around center <code>b</code>. For example, <code>p</code> might be "abadaba", with <code>b</code> equal to $9$. I assume <code>p</code> is the longest palindrome that reaches to its rightmost position, so that no palindrome around a center before <code>b</code> reaches there. <div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-36UHjp10tJY/T6gmWxU-Y0I/AAAAAAAAAIM/3cFsgUTA6Ak/s1600/pinp.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="185" width="400" src="http://4.bp.blogspot.com/-36UHjp10tJY/T6gmWxU-Y0I/AAAAAAAAAIM/3cFsgUTA6Ak/s400/pinp.jpg" /></a></div>I calculate the maximal palindromes from left to right, so I have already calculated the maximal palindrome around center <code>a</code>, which appears within palindrome <code>p</code>. I call the maximal palindrome around center <code>a</code> <code>q</code>. In the "yabadabadoo" example <code>a</code> might be center $5$, with "aba" as maximal palindrome <code>q</code> around it. In the illustration <code>q</code> is completely within <code>p</code>, but it might also extend to before <code>p</code>. It cannot extend after <code>p</code>, since then there would be a longer palindrome that extends to the rightmost position of <code>p</code>. Center <code>a'</code> is the center I get when mirroring center <code>a</code> with respect to center <code>b</code>: <code>b</code> plus the difference between <code>b</code> and <code>a</code>. This would be $13$ ($=9+(9-5)$) in our example. What is the maximal palindrome around <code>a'</code>? </p> <p>Since <code>p</code> is a palindrome, the string with the same length as <code>q</code> around <code>a'</code> is <code>reverse q</code>, provided <code>q</code> is completely within <code>p</code>. If <code>q</code> extends to before <code>p</code>, we cut it off at the start of <code>p</code>, and remove an equally long part at the end of <code>q</code>. Since <code>q</code> is a palindrome, <code>reverse q</code> equals <code>q</code>, and is hence also a palindrome. I will call this palindrome <code>q'</code> to distinguish it from <code>q</code>. Is it the maximal palindrome around <code>a'</code>? </p> <p>To find out whether or not <code>q'</code> is the maximal palindrome around <code>a'</code>, I determine whether or not the characters before and after <code>q'</code> are the same. If the maximal palindrome around center <code>a</code> does not extend to the start of <code>p</code>, I know that the characters around <code>q</code> are different. Since the same characters appear around the palindrome <code>q'</code>, it is the maximal palindrome around center <code>a'</code>. In this case I don't even have to look at the characters before and after <code>q'</code>: the fact that <code>q</code> does not extend to the start of <code>p</code> gives me enough information to determine that <code>q'</code> is maximal. If the maximal palindrome around center <code>q</code> <i>does</i> extend to the start of <code>p</code>, as in our example, where "abadaba" starts with "aba", or even before, I have to compare the character before <code>q'</code> with the character after <code>p</code>. If these are different, <code>q'</code> is the maximal palindrome around <code>a'</code>, otherwise, I have to investigate further to find out the maximal palindrome around <code>a'</code>. I will do this in the blog post that discusses an efficient algorithm for finding the length of all maximal palindromes in a string. In our example, "aba" is not the maximal palindrome around center <code>a'</code>, since the character 'd' appears both before and after <code>q'</code> </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-26746828191831073292012-05-09T22:11:00.000-07:002012-05-09T22:19:16.055-07:00The word `palindrome'<p>The word `palindrome' comes from the Greek παλινδρόμους. The word παλιν translates to `again', and the word δρόμους (or δρόμος) to `way' or `direction', and the combination is often translated to `running back again'. In the 17th century, the English poet Ben Johnson connected the <i>word</i> palindrome to the <i>concept</i> of a word or a sequence of words that reads, letter for letter, the same backwards as forwards. The palindrome concept existed long before the 17th century. Reportedly, the first recorded palindromes were written by Sotades in (again) Greece in the third century B.C. </p> <p>`Palindrome' is an old word, which already appears in some early Greek texts. For example, it appears in the work Timon of Lucian of Samosata. <div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-fDwBMTkDYN0/T6oL4dYj7xI/AAAAAAAAAIc/14xamNzdZws/s1600/Lucianus.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="200" width="125" src="http://1.bp.blogspot.com/-fDwBMTkDYN0/T6oL4dYj7xI/AAAAAAAAAIc/14xamNzdZws/s200/Lucianus.jpg" /></a></div>Lucian was a rhetorician and satirist living in the 2nd century A.D. He wrote more than eighty works, among which Timon, sometimes also translated as The Misanthrope. Shakespeare's play Timon of Athens has been influenced by the work of Lucian. </p><p>Timon of Athens was a citizen of Athens during the era of the Peloponnesian War between Athens and Sparta (431 BC–404 BC). He was wealthy and lavished his money on flattering friends. When funds ran out, friends deserted and Timon was reduced to working in the fields. One day, he found a pot of gold and soon his unreliable friends were back. This time, he drove them away with dirt clods. Timon became an angry despiser of mankind, and his reputation for misanthropy grew to legendary status. <div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-szFFOaT4a4o/T6oTmYUTwgI/AAAAAAAAAIs/Ivc9rkdoR5k/s1600/Howard-TimonAct4.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="200" width="124" src="http://1.bp.blogspot.com/-szFFOaT4a4o/T6oTmYUTwgI/AAAAAAAAAIs/Ivc9rkdoR5k/s200/Howard-TimonAct4.jpg" /></a></div> He can't have been nice company: a typical sentence of the misanthropic Timon in Lucian's work is `Go your ways, then, Hermes, and take Plutus back to Zeus. I am quite content to let every man of them go hang.' Hermes is the messenger god, en Plutus the god of wealth. In Greek, this sentence reads: ὥστε παλίνδρομος ἄπιθι, ὦ Ἑρμῆ, τὸν Πλοῦτον ἐπανάγων τῷ Διί: ἐμοὶ δὲ τοῦτο ἱκανὸν ἦν, πάντας ἀνθρώπους ἡβηδὸν οἰμώζειν ποιῆσαι. The second word is the occurrence of `palindrome'. It is not translated by `running back again', but clearly something goes back, or returns, here. </p> <p>Another occurrence of `palindrome' is found in Diogenes Laërtius' book `Lives of Eminent Philosophers', which was probably published somewhere in the 3rd century A.D. Laërtius introduces Aristippus in Chapter 8. <div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-adx63E__YtU/T6oZSKxRfbI/AAAAAAAAAI8/hABVYSu1fGA/s1600/Aristippus.jpg" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="200" width="142" src="http://2.bp.blogspot.com/-adx63E__YtU/T6oZSKxRfbI/AAAAAAAAAI8/hABVYSu1fGA/s200/Aristippus.jpg" /></a></div> Aristippus was drawn to Athens by the fame of Socrates. He became a lecturer himself, and was one of the first to both pay and charge fees for lecturing. `And on one occasion the sum of twenty minae which he had sent was returned to him, Socrates declaring that the supernatural sign would not let him take it; the very offer, in fact, annoyed him.' In Greek, this reads: καί ποτε πέμψας αὐτῷ μνᾶς εἴκοσι παλινδρόμους ἀπέλαβεν, εἰπόντος Σωκράτους τὸ δαιμόνιον αὐτῷ μὴ ἐπιτρέπειν: ἐδυσχέραινε γὰρ ἐπὶ τούτῳ. The 7th word of this sentence is `palindrome' in Greek. Again, it is not translated exactly as `running back again', but something is returned. </p><p>Although the roots of the word `palindrome' are in Greek, the Greek themselves describe the palindrome concept by <i>karkinikê epigrafê</i> (καρκινικὴ επιγραφή; "crab inscription"), or simply karkinoi (καρκίνοι; "crabs"), alluding to the movement of crabs. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-48998273087711374522012-05-03T12:37:00.000-07:002012-05-22T23:21:21.919-07:00A naive algorithm for finding maximal palindromes<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>The <a href="http://finding-palindromes.blogspot.se/2012/04/maximal-palindromes.html">previous blog post</a> discusses why it is sufficient to calculate the length of the maximal palindrome around each center position in a string if I want to find palindromes in a string. In this blog post I will describe a naive algorithm for finding the length of all maximal palindromes in a string. This algorithm is slightly less naive than the naive algorithm for finding palindromes given in an <a href="http://finding-palindromes.blogspot.se/2012/04/naive-algorithm-for-finding-palindromes.html">earlier blog post</a>. The latter algorithm requires a supercomputer for a couple of days to find palindromes in a long string, the algorithm I will develop in this blogpost requires several days of computing power of the laptop on which I am typing this post. In numbers, this algorithm is about a milliard times faster on such a long string. </p> <p>Given a string as input, I want to find the list of lengths of maximal palindromes around all centers of the string. I will use the function <code>maximalPalindromes</code> for this purpose. It has the following type <pre><br /><br />< maximalPalindromes :: String -> [Int]<br /><br /></pre>I want to find the length of the maximal palindrome around each center in a string. I will do this by trying to extend the trivial palindromes consisting of either a single letter (for odd centers) or of the empty string (for even centers) around each center. To extend a palindrome, I have to compare the characters before and after the current palindrome. It would be helpful if I had random access into the string, so that looking up the character at a particular position in a string can be done in constant time. Since an array allows for constant time lookup, I change the input type of the <code>maximalPalindromes</code> function to an array. I have to import the module <code>Data.Array</code> to use arrays. <pre><br /><br />> import Data.Array<br />><br />> maximalPalindromes :: Array Int Char -> [Int]<br /><br /></pre>If I change my input type from strings to arrays, I have to convert an input string into an array. The <code>Data.Array</code> module contains a function just for this purpose: the function <code>listArray</code> creates an array from a string together with a pair of indices for the first and last position in the string. Function <code>maximalPalindromes</code> calculates the following lengths of maximal palindromes in the string "abb": <pre><br /><br />?maximalPalindromes (listArray (0,2) "abb")<br />[0,1,0,1,2,1,0]<br /><br /></pre></p> <p>Function <code>maximalPalindromes</code> calculates the length of maximal palindromes by first calculating all center positions of an input array, and then the length of the maximal palindrome around each of these centers. The center positions of an array with <code>first</code> and <code>last</code> as <code>bounds</code> are the elements in the list <code>[0 .. 2*(last-first+1)]</code>. I will always use arrays that start at index <code>0</code>, which implies that the length of an array is <code>last+1</code>. <pre><br /><br />> maximalPalindromes a = <br />> let (first,last) = bounds a<br />> centers = [0 .. 2*(last-first+1)]<br />> in map (lengthPalindromeAround a) centers<br /><br /></pre>Function <code>lengthPalindromeAround</code> takes an array and a center position, and calculates the length of the longest palindrome around that position. Center positions do not exactly correspond to indices in the array, since they not only describe locations of characters, but also locations in between characters. To obtain an array index from a center position I divide it by two. If the center position is even, it is in between two characters, and if I divide it by two, I get an array index that points at the character after the center position. I compare the character pointed to by the index with the character before it to find out if the (empty) palindrome can be extended. If it can, I subtract one from the starting position, and add one to the end position, and try to extend the palindrome with the new indices. If the palindrome cannot be extended, because the elements around it are different, or the current palindrome starts at the start or ends at the end of the array, I return the length of the maximal palindrome found, which is the difference between the end and the start position minus one. For odd center positions I can start with the indices before and after the character pointed to by the index obtained by dividing the center position by two, rounding down. <pre><br /><br />> lengthPalindromeAround :: Array Int Char -> Int -> Int<br />> lengthPalindromeAround a center <br />> | even center = <br />> lengthPalindrome (first+c-1) (first+c) <br />> | odd center = <br />> lengthPalindrome (first+c-1) (first+c+1) <br />> where c = div center 2<br />> (first,last) = bounds a<br />> lengthPalindrome start end = <br />> if start < 0 <br />> || end > last-first <br />> || a!start /= a!end<br />> then end-start-1<br />> else lengthPalindrome (start-1) (end+1) <br /><br /></pre>For each position, this function may take an amount of steps linear in the length of the array, so this is a quadratic-time algorithm. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-55787633971100498122012-04-27T11:48:00.000-07:002012-05-20T05:34:33.784-07:00Maximal palindromes<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>The algorithm for finding palindromic substrings described in my <a href="http://finding-palindromes.blogspot.se/2012/04/naive-algorithm-for-finding-palindromes.html">previous blog post</a> has two major problems. First, if I want to find the palindromes in a string of substantial length, it will take a normal computer many years to calculate the result. I didn't mention the second problem explicitly in my previous blog post, but it is pretty obvious. The number of palindromic substrings returned by the algorithm is at least as large as the length of the string, since all single letters are palindromes, and there might be many more. If the string contains a palindrome of length $1,000,000$, then the algorithm will also return the palindrome of length $999,998$ obtained by removing the first and last element of the million letter palindrome, the palindromes of length $999,996$, $999,994$, and so on. This implies that the total number of palindromes occurring in a long string is huge, and I can easily drown in the palindromes returned. </p><p>Of the series of palindromes of length $1,000,000$, $999,998$, $999,996$ etc, all palindromes have the same <i>center</i>, and all shorter palindromes can be derived from the longest palindrome by removing equally many characters from the front and from the end. In the string "yabadabadoo", the palindromes around the center at the second occurrence of a 'b', are "b", "aba", and "dabad". "dabad" is the <i>maximal</i> palindrome around this center, since its extension "adabado" is not a palindrome. Maximality does not imply it is the longest palindromic substring, it only is the maximal palindrome around a particular center in the string. For example, in "yabadabadoo", "abadaba" is a longer palindrome, but with a different center. A center position is either on a letter, as in "dabad" on 'b', or in between two letters, as in "oo", where the center is in between the two o's. The string "yabadabadoo" has $23$ centers: one on each letter (of which there are eleven), one before each letter (another eleven), and one after the last letter of the string. If I assign $0$ to the center before the first letter, $1$ to the center on the first letter, $2$ to the center after the first letter, and so on, the maximal palindrome "dabad" has center $13$, and the maximal palindrome "abadaba" $9$. For a string of length $n$, there are $2n+1$ center positions in the string. Since we can derive all shorter palindromes around a center from the maximal palindrome around the center, it is sufficient to calculate the maximal palindromes around the centers in a string. So instead of calculating all palindromic substrings of a string, I calculate all maximal palindromic substrings of a string. </p><p>There are just as many maximal palindromes in a string as there are center positions, so there still are quite a few maximal palindromes. But the number of maximal palindromes might be substantially lower than the number of (not necessarily maximal) palindromic substrings in a string. Take for example the string of length $n$ just containing the character 'a'. The number of maximal palindromes is equal to the number of center positions in the string, $2n+1$. The total number of palindromes in the string is the number of substrings of the string, since all substrings consist of only 'a's, and hence all are palindromes. Since <code>substrings</code> is defined as <code>concat . map tails . inits</code>, I can calculate the total number of substrings if I know how many initial and final substrings appear in a string of length $n$. Function <code>inits</code> returns $n+1$ substrings, of length $0...n$, respectively. Similarly, for a string of length $n$, function <code>tails</code> returns $n+1$ substrings, also of length $0...n$. So the total number of substrings is: \[ \sum\limits_{i=0}^{n} (i+1) = \frac{n(n+1)}{2} + n+1 = \frac{(n+2)(n+1)}{2} = \frac{n^2}{2}+\frac{3}{2}n+1\] So in the string "aaaaaaaaaa" of length $10$, there are $21$ maximal palindromes, and $66$ palindromic substrings. </p><p>If I know the center and the length of a palindrome, I can recover the palindrome if I have the original string. For example, the maximal palindrome of length $5$ around center $13$ in the string "yabadabadoo" is the string "dabad", and the maximal palindrome of length $7$ around center $9$ is the string "abadaba". So to find maximal palindromes, it suffices to find their center and length. </p><p>Given a <code>string</code>, a <code>center</code> within the string, and a <code>lengthp</code> denoting the length of the maximal palindrome in the <code>string</code> around the <code>center</code>, I can use a Haskell program to check if the length is indeed the length of the maximal palindrome around the center in the string. To determine the substring denoted by a center and a length, I first split the string into the initial part of length <code>div (center-lengthp) 2</code>, the part before the maximal palindrome, and the following <code>lengthp</code> letters. Function <code>div</code> divides its first argument by its second argument, throwing away the remainder. Function <code>splitAt</code> takes a positive integer $n$ and a list, and splits the list into two lists. The first list contains the first $n$ letters, and the second list contains the rest of the letters. The <code>lengthp</code> letters after <code>before</code> should form a palindrome. The letters around it, the <code>last</code> letter of the part of the string before it and the <code>head</code> of the string after it, should be different to make it maximal. Functions <code>last</code> and <code>head</code> are predefined functions which return the last or head element of a list, respectively. The palindrome is maximal too if at least one of the strings before or after it is empty. <pre><br /><br />> maximalPalindrome string center lengthp =<br />> let (before,rest) = splitAt <br />> (div (center-lengthp) 2) <br />> string<br />> (p,after) = splitAt lengthp rest<br />> in length p == lengthp <br />> && odd center == odd lengthp<br />> && palindrome p <br />> && ( null before <br />> || null after <br />> || last before /= head after<br />> )<br /><br /></pre>The function <code>maximalPalindrome</code> uses two sanity checks. It first checks that there indeed exists a substring of length <code>lengthp</code> around <code>center</code> by means of <code>length p == lengthp</code>, where function <code>length</code> is a predefined function that returns the length of a list. The other sanity check checks that the maximal palindrome around an even center has even length, and a maximal palindrome around an odd center has odd length, by means of <code>odd center == odd lengthp</code>, which equals <code>True</code> only if both <code>center</code> and <code>lengthp</code> are odd or both are even. Function <code>odd</code> is a predefined function that determines whether or not a number is odd. The logic or operator <code>||</code> returns <code>True</code> if one of its arguments is <code>True</code>. If its left-hand side argument is <code>True</code> it doesn't evaluate its right-hand argument, which is just as well for the above definition, since evaluating <code>last before</code> if <code>before</code> is empty (<code>null before</code>) would lead to a undefined value. </p><p> In the next blog post I will show how to calculate the length of all maximal palindromes in a string. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-4233099619027241502012-04-18T11:52:00.000-07:002012-04-25T01:03:50.373-07:00A naive algorithm for finding palindromes<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); </script> <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"> </script> <p>It is easy to check by hand that a sentence like "Madam, I'm Adam" is a palindrome if you ignore whitespace and punctuation characters. But checking that <a href="http://norvig.com/palindrome.html">the world's longest palindrome</a>, consisting of 17,826 words, is indeed a palindrome is another matter. For such a string, I need software to check that it is a palindrome. </p> <p>I have already defined the necessary software to check that a string is a palindrome. The <code>palindrome</code> (or <code>textPalindrome</code> or <code>dnaPalindrome</code>) function defined in the <a href="http://finding-palindromes.blogspot.se/2012/04/what-is-palindrome-ii-specification-in.html">previous blog post</a> does the job. Given an implementation of the function <code>reverse</code> that takes about as many steps as the length of the string it reverses, the <code>palindrome</code> function determines whether or not a string is a palindrome in about the same amount of steps. It is impossible to do this significantly faster on a computer. </p><p>Problem solved? </p><p>Well, it depends on what you want. If you want to check for a complete string whether or not it is a palindrome, then the problem is solved. But often you want to find occurrences of palindromes inside another, possibly non-palindromic, string. For example, you might want to find the about 5,700,000 characters that together form eight palindromes, in a string of about 20,000,000 characters representing the male-specific region of the Y chromosome. Or you night want to find the longest palindrome in the Bible. These are different problems, for which we need more than just the <code>palindrome</code> function. </p><p>To find palindromes in the male-specific region of the Y chromosome, I have to consider <i>all</i> substrings of the DNA string, and check if each of the substrings is a palindrome. How can I find all substrings that are palindromes in a string? A string may contain many substrings that are palindromes. Just the string "abba" contains the palindromes "a" (twice), "b" (twice), "bb", and "abba", and we also consider the empty substring, which appears a lot, a palindrome. </p><p>Finding all palindromic substrings is specified by: <pre><br /><br />> palindromes :: String -> [String]<br />> palindromes = filter palindrome . substrings<br /><br /></pre>The first line gives the <i>type</i> of the function <code>palindromes</code>. Function <code>palindromes</code> takes a value of type <code>String</code> as argument, and returns a list of palindromic strings, so a value of type <code>[String]</code>. The argument type is to the left of the arrow <code>-></code>, the result type to the right. Function <code>substrings</code> calculates all substrings of a string. The function <code>substrings</code> is composed with the function <code>filter palindrome</code> using the composition operator <code>.</code> represented by the dot symbol. <code>filter palindrome</code> removes all substrings that are not palindromes. It remains to define function <code>substrings</code>. Function <code>substrings</code> is defined in terms of two helper functions <code>inits</code> and <code>tails</code>. Function <code>inits</code> returns the list of all initial substring of a string. For example, <code>inits "abba"</code> gives <code>["", "a", "ab", "abb", "abba"]</code>. Similarly, function <code>tails</code> returns the list of all final substrings of a string, so <code>tails "abba"</code> gives <code>["abba", "bba", "ba", "a", ""]</code>. <pre><br /><br />> inits [] = [[]]<br />> inits (x:xs) = [[]] ++ map (x:) (inits xs)<br /><br />> tails [] = [[]]<br />> tails (x:xs) = [x:xs] ++ tails xs<br /><br /></pre>The empty list has just one initial list, namely the empty list. The initial lists of a list consisting of an element <code>x</code> followed by a list <code>xs</code> are calculated by first calculating the initial lists of <code>xs</code>, prepending <code>x</code> to each of these lists by means of the function <code>(x:)</code>, and adding the empty initial list. The empty list also has just one final list, namely the empty list. The final lists of a list consisting of an element <code>x</code> followed by a list <code>xs</code> are calculated by first calculating the final lists of <code>xs</code>, and adding the complete list <code>x:xs</code>. Function <code>substrings</code> is now defined by either taking all initial substrings of all final substrings, or all final substrings of all initial substrings. It doesn't matter much which I take, so I define <pre><br /><br />> substrings = concat . map tails . inits<br /><br /></pre>The standard function <code>concat</code> <i>flattens</i> takes a list of lists, and turns that into a list. So <code>concat [[3,1],[2,4]]</code> results in <code>[3,1,2,4]</code>. It removes one nesting level. Applying function <code>substrings</code> to <code>"abba"</code> gives <code>["", "a", "", "ab", "b", "", "abb", "bb", "b", "", "abba", "bba", "ba", "a", ""]</code>. </p> <p>Although we now have a program to find all palindromes in a string, it is of no use when the string in which you want to find palindromes is very large. You might just as well try to find the palindromes by hand, since waiting for the program to find them is easily going to take years. Just how impossible it is to find palindromes in a string of length 20,000,000 using this program, consider the amount of computation necessary to find them. Since we want to find palindromic substrings of this string, we have to calculate the substrings first. We have a single substring of length 20,000,000, two substrings of length 19,999,999 each, three of length 19,999,998, four of length 19,999,997, and so on. So the length of all substrings we need to consider when looking for palindromes is \[ 1 \times n + 2 \times (n-1) + 3 \times (n-2) + \ldots + n \times 1 \] where $n$ is the length of the string in which we are looking for palindromes. I will use some simple arithmetic to obtain a value that is easier to calculate. The above expression can be reformulated to: \[ \sum\limits_{i=1}^{n} i \times (n+1-i) \] Now I apply a law which says that constant factors can be moved outside a sum, so \[ \sum\limits_{i=1}^{n} i \times k = k \times \sum\limits_{i=1}^{n} i \] where I instantiate $k$ with $n+1$, and I apply laws for sums of consecutive integers, namely \[ \sum\limits_{k=1}^n k = \frac{n(n+1)}{2} \\ \sum\limits_{k=1}^n k^2 = \frac{2n^3 + 3n^2 + n}{6} \] and some simple arithmetic to find that the above expression equals: \[ \frac{n^3+3n^2+2n}{6} \] So given a string of length $n$, the length of all substrings that appear in the string is $\frac{n^3+3n^2+2n}{6}$, which is slightly less than $n^3$. </p><p>Let me return now to the DNA string mentioned at the start of this blog post. Since its length is around 20,000,000, the total length of the substrings that appear inside it amounts to 1,333,333,533,333,340,000,000. This is more than a trilliard characters to be inspected! Suppose, unrealistically, that I can use the fastest supercomputer available on earth as of 2011, the Fujitsi K computer, named after the japanese word for 10 billiard, kei. <div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-Gjq1lejFr8o/T46W4wsQX4I/AAAAAAAAAGw/VkYmIKoKSKs/s1600/k-computer.jpg" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="139" width="200" src="http://4.bp.blogspot.com/-Gjq1lejFr8o/T46W4wsQX4I/AAAAAAAAAGw/VkYmIKoKSKs/s200/k-computer.jpg" /></a></div>The kei computer is actually not so much a computer, but a collection of over 80,000 computers put together. The speed of a supercomputer is specified as the number of basic operations it can perform per second. The speed of the K computer is 10.51 petaflops, where a petaflop is a billiard of such basic operations per second. So the K computer is appropriately but modestly named after its speed. If a single basic operation is enough to compare two characters (which it isn't), then the K computer would need $1,333,333$, the length of the substrings divided by a billiard, divided by $10.51$, the number of billiard operations the K computer can perform per second, seconds, or about one and a half day, to calculate these comparisons. A normal computer would spend its entire lifetime on this problem, and not finish it. </p><p>That is not good enough. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-66520746610912668132012-04-15T11:27:00.001-07:002012-06-29T22:59:15.681-07:00What is a palindrome II? A specification in Haskell<p>In this blog post I will describe the concept of a palindrome more precisely than in my last blog post, by giving the property that determines whether or not a sequence of symbols is a palindrome in the functional programming language Haskell. </p> <p>Why is a specification of a palindrome as a program more precise than the specification given in text in the <a href="http://finding-palindromes.blogspot.se/2012/04/what-is-palindrome.html">previous blog post</a>, which was taken from the Oxford English Dictionary (OED)? And why should someone interested in the concept of palindromes be interested in a program that describes the concept? These are philosophical questions, to which computer science helps to give an answer. The purpose of a definition is to distinguish certain things from other things. The definition of a palindrome distinguishes palindromic words or sequences of words from non-palindromic ones. How do I determine if a sequence of words is a palindrome? The textual definition in the OED says that it should read the same backwards as forwards, letter for letter. The first example in the OED is "Lewd did I live, and evil I did dwel", which doesn't read the same backwards and forwards, unless capital and underscore letters are considered equal, and comma's are ignored. Hmmm. Are there any more exceptions I should know of? </p> <p>A much better way to define the concept of palindromes is to provide a program, which takes a string as input, and which returns either yes or no, determining whether or not the input string is a palindrome. The OED definition leads to three different programs for determining whether or not an input string is a palindrome: a definition which compares strings letter for letter, and only accepts string which are literally the same forwards and backwards, a definition which ignores punctuation and capitalization, and a definition which compares DNA strings. Such a program is a much more precise way to define what it means to be a palindrome. To determine if a string is a palindrome it suffices to run the program on the string. To study the concept of palindromes, it suffices to study the program. No surprises. </p> <p>Not just palindromes profit from being defined by means of a program, many other concepts, ideas, and regulations would be better off being specified as a program. Besides the obvious anagrams, pangrams, ambigrams, and so on, tax laws, testaments, exam regulations, and anything that follows some rules is best described by means of a program. A program is precise about the rules, the exceptions, and when they apply. If a concept cannot be described precisely, using a program becomes harder. I wouldn't know how to construct a Turing test: a program to determine whether or not a species is intelligent. And even a simple concept like a chair is not easily captured in a program. </p> <p><a href = "http://www.haskell.org/">Haskell</a> is a programming language celebrating its 25th anniversary in 2012. It was conceived at an academic conference on functional programming in 1987, in an attempt to combine efforts in developing a purely functional programming language. In the nineties of the last century it was used for research on advanced programming language concepts, and in teaching at many universities. Since a number of years it is becoming more popular in industry, in particular in the investment banking industry, where it is used to specify and implement involved mathematical financial models. I use Haskell because it allows me to describe ideas precisely and concisely. I will introduce the necessary Haskell concepts as I go. </p> <p>This blog post itself is a <i>literate</i> Haskell program. If you save it as a file and make sure its name ends in .lhs, you can load it in <a href = "http://hackage.haskell.org/platform/">ghci</a>, an <i>interpreter</i> for Haskell programs, and use the definitions given in this post. To make this work I need the following two lines. <pre><br /><br />> import Prelude hiding (reverse)<br />> import Data.Char hiding (isLetter)<br /><br /></pre>The first line says that I can use all the standard Haskell functions except the function <code>reverse</code>, which I am going to define in this blog post. The second line says that I can use all basic functions on characters, such as the function <code>isSpace</code>, which tells whether or nor a character is a space. The function <code>isLetter</code> is excluded, since I am going to define that function in this blog post too. </p> <p>I use the following notation for a list of characters. The empty list of characters is denoted by <code>[]</code> and a symbol <code>x</code> followed by a list of characters <code>xs</code> is denoted by <code>x:xs</code>. An individual character is surrounded by quotes to distinguish it from arbitrary variables like <code>x</code>. For example, in this notation, I can write "refer" as <code>'r':'e':'f':'e':'r':[]</code>. Furthermore, I can concatenate two lists of characters by means of the operator <code>++</code>, so that<code> "Madam, " ++ "I'm Adam"</code> is<code> "Madam, I'm Adam"</code>. </p> <p>How can I determine whether or not a string (a list of characters) is a palindrome? The simplest method is to reverse the string and to compare it with itself. So the string <code>xs</code> is a palindrome (<code>palindrome xs</code>), if <code>xs</code> is equal to its reverse: <code>xs == reverse xs</code>, where <code>xs == ys</code> is <code>True</code> only when the strings <code>xs</code> and <code>ys</code> are exactly equal. <pre><br /><br />> palindrome xs = xs == reverse xs<br /><br /></pre>This equation is actually a <i>program</i> in Haskell. We can load this program in ghci and run it. For example, <code>palindrome "abba"</code> evaluates, unsurprisingly, to <code>True</code>, and <code>palindrome "yabadabadoo"</code> evaluates to <code>False</code>. <p>How do I compute the reverse of a string? For this purpose, I define a <i>recursive function</i>: if I know how to calculate the reverse of the empty string, and if I know how to calculate the reverse of <code>x:xs</code>, given that I know how to calculate the reverse of <code>xs</code>, then I can use a computer to calculate the reverse of any string. The two cases for reverse are not hard: the reverse of the empty string is the empty string, and the reverse of <code>x:xs</code> is the reverse of <code>xs</code> followed by the string consisting of the single character <code>x</code>: <pre><br /><br />> reverse [] = []<br />> reverse (x:xs) = reverse xs ++ [x]<br /><br /></pre>Here <code>[x]</code> is the list containing the single element <code>x</code>, which can also be written<code> "x"</code>. This completely defines what it means to be a palindrome. </p> <p>The above definition only qualifies strings as palindromes when they are exactly equal when reversed. So<code> "Madam, I'm Adam"</code> does not pass this test. For this string to also pass the <code>palindrome</code> test, I slightly adapt the definition. I now say that a string is a palindrome if it is equal to its reverse after throwing away all punctuation symbols such as spaces, comma's, periods, etc, and after turning all characters into lower case characters. I call such a palindrome a <code>textPalindrome</code>. <pre><br /><br />> textPalindrome xs = let ys = filter isLetter xs<br />> zs = map toLower ys<br />> in zs == reverse zs<br />><br />> isLetter l = not (isPunctuation l) && not (isSpace l)<br /><br /></pre>Haskell's <code>let ... in ...</code> construct allows me to introduce new definitions after the <code>let</code>, which I can use in the definitions after the <code>in</code>. Here I use two intermediate results <code>ys</code> and <code>zs</code>. By filtering the original string <code>xs</code> with the predicate <code>isLetter</code>, I obtain a new string <code>ys</code> in which no punctuation or space characters appear anymore. So<code> "Madam, I'm Adam"</code> is turned into<code> "MadamImAdam"</code>. Function <code>filter</code> is a standard function in Haskell which takes a predicate <code>p</code> and a list <code>xs</code>, and keeps all the elements of <code>xs</code> that satisfy the predicate <code>p</code>. The function <code>isLetter</code> uses two basic Haskell functions on characters, <code>isPunctuation</code> and <code>isSpace</code>. <code>isPunctuation</code> returns <code>True</code> for all punctuation characters such as the dot, comma, space, and exclamation mark, and <code>False</code> for all non-punctuation characters, such as letters. Function <code>isSpace</code> works similarly on space characters. Function <code>isLetter</code> takes a character, and checks that it is <code>not</code> a punctuation character, and (denoted by <code>&&</code>) <code>not</code> a space character. After filtering out all non-letter characters from <code>xs</code> giving <code>ys</code>, I map each letter in <code>ys</code> to its lower case by means of <code>map toLower ys</code>. So<code> "MadamImAdam"</code> is turned into<code> "madamimadam"</code>. <code>map</code> is also a standard Haskell function, which takes a function as argument, <code>toLower</code> in this example, and applies this function to all characters in a string. The function <code>toLower</code> turns a capital letter into lowercase, and does nothing to lowercase letters. Function <code>textPalindrome</code> accepts all palindromes, irrespective of their punctuation. </p> <p>Both <code>palindrome</code> and <code>textPalindrome</code> cannot be used to check if a DNA sequence is a palindrome, because the symbol <code>'A'</code> should be considered equal to <code>'T'</code>, and <code>'C'</code> to <code>'G'</code>, and the equality <code>==</code> used in the definition of <code>palindrome</code> considers them different. We need to change the equality function to compare characters in DNA strings. We define the DNA character equality function <code>=:=</code> by <pre><br /><br />> 'A' =:= 'T' = True<br />> 'T' =:= 'A' = True<br />> 'C' =:= 'G' = True<br />> 'G' =:= 'C' = True<br />> _ =:= _ = False<br /><br /></pre>We use this new equality function in a definition of <code>dnaPalindrome</code> for sequences of DNA symbols. We pairwise combine the elements of <code>xs</code> and <code>reverse xs</code> with the equality function <code>=:=</code> using the function <code>zipWith</code>. <code>zipWith</code> is another standard function, which takes an operator and two lists, and `zips' the two lists with the operator. Thus we obtain a list of comparisons, which we fold to a single result by requiring that each element in the list is <code>True</code>. Function <code>and</code> takes a list of boolean values, and returns <code>True</code> only if all elements in the list are <code>True</code>. <pre><br /><br />> dnaPalindrome xs = and (zipWith (=:=) xs (reverse xs))<br /><br /></pre>Thus we have three functions for checking whether or not a string is a palindrome: <code>palindrome</code> for strings that are exactly equal when reversed, <code>textPalindrome</code> for strings that are exactly equal when reversed modulo punctuation and space characters, and <code>dnaPalindrome</code> for DNA strings. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-78146785061156510012012-04-08T05:03:00.000-07:002012-04-22T23:32:49.343-07:00What is a palindrome?<p>The Oxford English Dictionary describes a palindrome as "a word or a sequence of words that reads, letter for letter, the same backwards as forwards." In extended use it appears in: <ul> <li><i>Music</i>. A piece of music in which the second half is a retrograde repetition of the first half; the retrograde itself.</li><li><i>Numbers</i>. A number, or a date expressed numerically, that is unchanged when the order of its digits is reversed.</li><li><i>Biology</i>. A nucleic acid sequence that is identical to its complementary sequence when each is read in the same direction (which is usually the direction called 5' - 3').</li></ul>I will return to these definitions in later blog posts. </p> <p>To determine whether a sequence of symbols is a palindrome I need to know what the symbols are from which the palindrome is composed, which symbols I can ignore, and when two symbols are considered equal. </p><p>The above definitions tell me that the symbols may be single characters, notes, digits, or elements of DNA. But even complete words, or sentences, are used in palindromes. J.A. Lindon, who wrote many palindromes, published the following poem in 1955: <blockquote>As I was passing near the jail<br/>I met a man, but hurried by.<br/> His face was ghastly, grimly pale.<br/>He had a gun: I wondered why,<br/>His face was ghastly, grimly pale,<br/>I met a man, but hurried by,<br/> As I was passing near the jail.<br/> </blockquote> Reportedly, the first recorded palindromes, written by Sotades in Greece in the third century BC are also on the level of sentences. </p><p>A palindrome is a simple kind of <a href="http://en.wikipedia.org/wiki/Symmetry_in_mathematics">mathematical symmetry</a>. Where symmetries are exact and the smallest imperfection disqualifies a symmetry, palindromes in which the symbols are characters are often not exact symmetries. The word "refer" reads "refer" when you read it backwards, but the sequence of words in Adam's introduction "Madam, I'm Adam" reads "madA m'I ,madaM" when reversed, which clearly is not the same. Thus it would not qualify as a palindrome if we are very strict, but when talking about palindromes in languages, capitalization, spaces, and punctuation symbols are almost always ignored, or even changed. </p><p>In some situations, even different symbols are considered equal when finding palindromes. In DNA, a double helix is formed by two paired strands of nucleotides that run in opposite directions, and the nucleotides always pair in the same way (Adenine (A) with Thymine (T); Cytosine (C) with Guanine (G)). A (single-stranded) nucleotide sequence is said to be a palindrome if it is equal to its reverse complement. For example, the DNA sequence ACCTAGGT is palindromic because its nucleotide-by-nucleotide complement is TGGATCCA, and reversing the order of the nucleotides in the complement gives the original sequence. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-76038287104952076582012-03-30T04:47:00.000-07:002012-07-03T23:12:06.377-07:00Why finding palindromes?<p>My former neighbor is an author. Many years ago he was writing a historical novel about the colonial past of the Netherlands. Unfortunately, his apartment was broken into, and the laptop that contained all material for the book was stolen. He had made no backups, and decided to write <a href="http://www.hansschellekens.nl/index.html">another book</a> instead. </p><p>Things break, disappear, get lost, or fall apart. Dead things, like laptops, washing machines, and light bulbs, need to be replaced or repaired when they break. Living things have some advanced repair mechanisms, with which they repair a broken leg, a bleeding finger, or a missing leaf. </p><p>To repair their parts, living things use the building plans available in their cells. These building plans are expressed in DNA. Each cell contains a complete copy of the DNA of an organism. When an organism builds a new cell, or repairs part of the organism that has been damaged, it uses DNA to create new cells. But what if the building plans themselves get damaged? Then reparation wouldn't be possible anymore, or would lead to malformed parts. </p><p><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-EfL2ByrANZ0/T3WaJ8PozzI/AAAAAAAAAGM/NEQCHzKBlU8/s1600/800px-DNA_palindrome.svg.png" imageanchor="1" style="clear:right; float:right; margin-left:1em; margin-bottom:1em"><img border="0" height="141" width="320" src="http://3.bp.blogspot.com/-EfL2ByrANZ0/T3WaJ8PozzI/AAAAAAAAAGM/NEQCHzKBlU8/s320/800px-DNA_palindrome.svg.png" /></a></div>The DNA mechanism has several means to repair broken DNA. Most of these means consist of maintaining copies of the DNA, using multiple copies of DNA to construct a new copy, or by repeating pieces of DNA inside the DNA itself. But one of the most intriguing mechanisms DNA uses is storing information in <i>palindromes</i>. If DNA in one half of a palindrome gets broken, the other half of the palindrome tells how to repair the broken DNA. </p><p>A palindrome is a number, verse, word or a phrase that is the same whether you read it backwards or forwards, for example the word "refer". </p><p>When I started as a PhD student in Computer Science in 1988, my professor gave the following problem assignment to me. Construct an algorithm that, when given a string, finds the longest palindrome substring occurring in the string. For example, given the string "yabadabadoo", the algorithm should return the substring "abadaba". The algorithm should not only <i>return</i> the longest palindrome substring, it should also return it <i>fast</i>. The requirement was that the algorithm should only use a number of steps comparable to the length of the argument string. This was a difficult problem, which gave me severe headaches in the months to follow. I spent four months on the problem, and solved it. The one comment I got from my professor was "The best thing about this problem and its solution is that they have absolutely no practical relevance." </p><p>He was wrong. </p><p><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-eacE23YqLUQ/T3Wc5oVB2PI/AAAAAAAAAGk/vbc69ey-DoU/s1600/MozartSpiegel.png" imageanchor="1" style="clear:left; float:left;margin-right:1em; margin-bottom:1em"><img border="0" height="200" width="154" src="http://3.bp.blogspot.com/-eacE23YqLUQ/T3Wc5oVB2PI/AAAAAAAAAGk/vbc69ey-DoU/s200/MozartSpiegel.png" /></a></div>Palindromes play an important role in the world around us. Palindromes appear in poetry, music, math, and in genomes: human life depends on palindromes! It is intriguing that a concept generally considered nice but useless is fundamental to our existence. </p><p>This blog discusses where palindromes come from, where they appear around us, what role they play in human life, and how you can use a computer to find them. The latter focus is inspired by my work during my PhD, and is essential for finding palindromes in DNA. The resulting program for finding palindromes can also be used to find palindromes in other sizable sources. For example, what is the longest palindrome occurring in the bible? Of course, this depends on the version of the bible, but suppose we take the King James Bible. I will reveal the longest palindrome in the bible in a later blog post, in which I describe an algorithm for finding palindromes. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0tag:blogger.com,1999:blog-5279002431366685031.post-21117130102859885892012-03-25T12:38:00.001-07:002012-04-09T23:10:17.823-07:00Introduction<p>You have probably come here via my blog post on <a href="http://johanjeuring.blogspot.com/2007/08/finding-palindromes.html">Finding palindromes</a>, which I wrote on the <a href="http://johanjeuring.blogspot.com/">blog</a> we used to accompany the ICFP programming contest in 2007. The blog post on finding palindromes receives thousands of hits per month, and readers frequently ask me to give more details of the algorithm described there, or of the programming language used. Since I don't want to use the ICFP programming contest blog for these purposes, I have started this new blog. In the coming months, I will write about palindromes on this blog: where they come from, where they appear, and how you can find them. </p><p>Any comments about my expositions are welcome. </p>Johan Jeuringhttp://www.blogger.com/profile/14035760811349164958noreply@blogger.com0