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

at first, i would like to thank all the people here. i'm working with mathematica for half a year now and always found great advice and tricks.

now, i came up with a problem i couldn't find a solution for. from a file, i import a string with around 23 mio. characters which is reasonably fast (it acutually is a dna sequence). i then want to split this string into its characters and convert each of them to an integer number.

my first, naive approach was the following:

MaxMemoryUsed[]/1024^2//N
AbsoluteTiming[list=Characters[a]/.{"A"->1,"T"->2,"G"->3,"C"->4};]
ByteCount[list]/1024^2//N
MaxMemoryUsed[]/1024^2//N

which creates the following output:

134.286
{13.301206,Null}
702.261
1976.42

as you can see, the timing is quite good but the memory usage during the process is horrible. my next approach was reading from a stream and convert each character to an integer instantly:

MaxMemoryUsed[]/1024^2//N
string=StringReplace[string,{"A"->"1","T"->"2","G"->"3","C"->"4"}];
str=StringToStream[string];
AbsoluteTiming[list=Table[ToExpression[Read[str,Character]],{i,1,StringLength[string]}];]
Close[str];
ByteCount[list]/1024^2//N
MaxMemoryUsed[]/1024^2//N

which results in

134.281
{176.502736,Null}
702.252
813.208

as yo can see, this approach needs way less memory but takes considerably longer. my question now is, if there is another memory efficient way to do this which is faster than my second try?

thanks a lot!

share|improve this question
    
Is the string stored in a file as one long line, or broken up in some way ? –  image_doctor Nov 14 '12 at 10:35
    
Can't test right now, but givent hat the slow part in your second approach is turning the string of numbers into a list, what about something like ToExpression@StringCases[string, DigitCharacter], or ToCharacterCode@string - 48 –  Rojo Nov 14 '12 at 10:50
    
to wrap this up, the string is broken up into many lines of constant length –  mincos Nov 14 '12 at 16:36
add comment

5 Answers 5

up vote 9 down vote accepted
string = StringJoin @@ RandomChoice[{"A", "T", "G", "C"}, 30000000];

This seems to take little memory

Boole[# === 65] + 2 Boole[# === 84] + Boole[# === 71] 3 + 
   Boole[# === 67] 4 & /@ ToCharacterCode[string]

but twice as long as

ToCharacterCode[
  StringReplace[
   string, {"A" -> "1", "T" -> "2", "G" -> "3", "C" -> "4"}]] - 48

Compiling the first solution makes it faster than the second

Compile[{{x, _Integer}}, \[Piecewise] {
    {1, x === 65},
    {2, x === 84},
    {3, x === 71},
    {4, x === 67}
   }, RuntimeAttributes -> Listable]@ToCharacterCode[string]
share|improve this answer
    
Interestingly in version 7 the first method takes 62 seconds and ~2GB of memory, while the second takes only 1.7 seconds and ~1GB (no more than required for string). –  Mr.Wizard Nov 14 '12 at 11:34
    
@Mr.Wizard string requires as many bytes as characters. Are you using a billion character string? –  Rojo Nov 14 '12 at 11:36
    
I mean the memory required for the creation of string using the code you provided`. (What do you get for MaxMemoryUsed[] right after?) –  Mr.Wizard Nov 14 '12 at 11:38
    
@Mr.Wizard, 650MB, but I forgot I was actually restarting the kernel and loading the string from the front end every time, so after that the memory in use is the same as the max memory used, which is 49MB –  Rojo Nov 14 '12 at 11:40
    
Ah, that would be wise to find the actual memory used, wouldn't it? :) (I mean what you're doing would be wise.) –  Mr.Wizard Nov 14 '12 at 11:40
show 5 more comments

OK, we don't need command-line tools for this. We can do better, faster and with less memory.

Lets create a test file with 30 million letters:

Export["tmp/sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 30000000]]

The trick is to make the conversion parallel on all data and to keep packet-arrays. First we define a parallel, compiled function which does the conversion for one character:

repl = Compile[{{n, _Integer, 0}},
  If[n == 65, 1,
   If[n == 67, 2,
    If[n == 71, 3,
     If[n == 84, 4, -1]
     ]
    ]
   ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
  Parallelization -> True
  ]

Yes, those Ifs look pretty awful, but they get compiled down while Switch or Which doesn't. See, that we don't read characters, we read the ASCII codes instead which is a nicely packed integer array:

MaxMemoryUsed[]/1024^2 // N
AbsoluteTiming[dna = repl[Developer`ToPackedArray[ReadList["tmp/sequence.txt", Byte]]];]
MaxMemoryUsed[]/1024^2 // N

Here this takes 2 seconds and less than 500 MB of memory.

Timings

I timed all methods on my machine. I really don't know how image_doctor got his method so fast. It was the slowest here (native Linux, SSD, Intel i7 extreme, 32 GB Ram). I didn't preload the data. If a method used a string, I used Import.

Mathematica graphics

Currently fastest method

I thought, that it's Rojos Piecewise construct was faster, but as it turns out this was cosmic radiation since the ugly nested Ifs and the Piecewise solution produce exact the same compiled code. The full test code consists of a sequence of only 3 million digits

Export["sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 3000000]];

which is called 100 times.

halirutan

f = Compile[{{n, _Integer, 0}}, 
  If[n == 65, 1, If[n == 67, 2, If[n == 71, 3, If[n == 84, 4, -1]]]], 
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
  Parallelization -> True];

Mean@Table[First@AbsoluteTiming[
    f[ToCharacterCode[Import["tmp/sequence.txt"]]];], {100}]

(* Out[2]= 0.181717 *)

Rojo compiled to "C" and parallel

f = Compile[{{x, _Integer}}, Piecewise[{{1, x === 65}, 
           {2, x === 84}, {3, x === 71}, {4, x === 67}}], 
       CompilationTarget -> "C", RuntimeAttributes -> {Listable},
   Parallelization -> True];

Mean@Table[First@AbsoluteTiming[f[
      ToCharacterCode[Import["tmp/sequence.txt"]]];], {100}]

(* Out[3]= 0.183226 *)

Note here, how mystical it is, that Piecewise (and Ifs with then branch (!)) get compiled correctly to "C", while Switch and Which stays in a MainEvaluate call.

share|improve this answer
    
I don't think the memory can be brought below twice the size of the packed array of integers. Unless you do a very slow in place modification. Your If version is faster than my Boole, +1. The parallelization in my computer doesn't make it faster this time. And The compilation to C I bet it doesn't if its just for one. I don't have a suitable compiler here :( –  Rojo Nov 14 '12 at 11:46
    
@Rojo The compilation to "C" is not to speed up the computation itself! It's that the vector-registers are used which is most likely not done by the Mathematica-VM. –  halirutan Nov 14 '12 at 11:51
    
Any reason why you used nested If[]s as opposed to a Switch[]? –  J. M. Nov 14 '12 at 11:57
    
@J.M. Yes, the switch doesn't compile and stays as MainEvaluate in the code. –  halirutan Nov 14 '12 at 12:03
    
The piecewise version in my edit generates the exact (except for defaulting to 0 instead of -1) same code (CompilePrint checked) as the nested If, avoiding the uglyness :D –  Rojo Nov 14 '12 at 12:16
show 8 more comments

I read 30 million characters stored in around 10,000 lines in a text file into a list of integers like this.

dna = Flatten[IntegerDigits /@ ReadList["!cat /tmp/dna.tmp | tr [ACGT] [1234]"]];
//AbsoluteTiming

{2.766828, Null}

Memory usage was 0.93GB for 29.6 mega characters.

The 0.93GB is effectively the space required for the 30 million integers. I guess this storage amount could be compressed if required as the underlying alphabet only has 4 symbols.

share|improve this answer
1  
That's the space required if they aren't packed. Otherwise they take 4 bytes each, so 120MB –  Rojo Nov 14 '12 at 11:28
    
Do you mean I should add ToPackedArray, I hadn't spotted anyone else using it? Is there some automatic packing going on somewhere ? Ah I see @halirutan has used it :) –  image_doctor Nov 14 '12 at 11:57
    
If you add ToPackedArray, you'll get an array of 120MB, but if you had to add it its because in the process it got unpacked, so somewhere along the line it used the 0.93GB. Some alternatives never leave the packed array realm. And those, use 2x120MB because they have to create a copy to turn the array into what the OP wants, and MMA always copies everything :P –  Rojo Nov 14 '12 at 12:03
    
The 0.93GB was the reported ByteCount for the variable dna. Is it automatically expanding up the report of memory usage from the dna variable that has somewhere been automatically packed? –  image_doctor Nov 14 '12 at 12:07
    
It just means that the dna variable is unpacked then. Pass it through ToPackedArray and then it should report 120MB –  Rojo Nov 14 '12 at 12:08
add comment

On my PC, using a lookup table is slightly faster than @halirutan's If or Piecewise method.

Export["tmp/sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 30000000]]

lookup = ConstantArray[0, 256];
lookup[[ToCharacterCode["ACGT"]]] = Range[4];

f = Compile[{{x, _Integer}, {lut, _Integer, 1}}, lut[[x]],
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
   Parallelization -> True];

AbsoluteTiming[f[ToCharacterCode[Import["tmp/sequence.txt"]], lookup];]

Takes 3.1198962s, compared to 3.2079074 with the Piecewise code in @halirutan's answer.

If you don't care what the function returns for characters other than A,C,G,T, this is even faster:

f = Compile[{{x, _Integer, 1}},
   UnitStep[x - 65] + UnitStep[x - 84] + UnitStep[x - 71] + UnitStep[x - 67]];
AbsoluteTiming[f[ToCharacterCode[Import["sequence.txt"]]];]

Takes 2.9723774s

(Probably because UnitStep can somehow use SIMD instructions and If and Part cannot)

share|improve this answer
    
your first version is faster for me as well. –  mincos Nov 14 '12 at 13:29
add comment

Though not quite as fast as Rojo's second method this is a little shorter and uses less memory on my machine:

ToCharacterCode @ string /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4}

You could also read your data as bytes to bypass ToCharacterCode like this:

BinaryReadList["gen.txt"] /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4}

However, his new compiled form is as fast and uses much less memory.

Comparison

First generating test data:

string = StringJoin @@ RandomChoice[{"A", "T", "G", "C"}, 30000000];

DumpSave["bigGen.txt", string];

Quit[], then:

ToCharacterCode[StringReplace[
     ReadList["bigGen.txt", String],
     {"A" -> "1", "T" -> "2", "G" -> "3", "C" -> "4"}]] - 48 // 
  Timing // First

MaxMemoryUsed[]

2.075

850227296

Quit[] again, and:

BinaryReadList["bigGen.txt"] /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4} // Timing // First

MaxMemoryUsed[] 

2.387

614962712

Now with Rojo's compiled piecewise form:

f = Compile[{{x, _Integer}},
   \[Piecewise]{{1, x === 65}, {2, x === 84}, {3, x === 71}, {4, x === 67}}];

f /@ BinaryReadList["bigGen.txt"] // Timing // First

MaxMemoryUsed[] 

2.262

254977528

This should be even better in Version 8 with the enhanced Compile.

share|improve this answer
    
Liked the BinaryReadList –  Rojo Nov 14 '12 at 12:17
add comment

Your Answer

 
discard

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

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