Message Boards Message Boards

Compiled function uses too much memory

Posted 4 years ago

This algorithm calculates the probability that S of n biological mtMRCA lineages survive after v generations. In the beginning (generation v=1), there are n mothers that have each randomly between 0 and d daughters. I have a fast analytical solution. But I need to check it with this algorithm. The uncompiled function is slow of course, but it does not need much memory. Unfortunately, memory usage of the compiled version is exponential for d > 2, n > 4 and V > 7 for a reason I don't understand. Can anyone help?

ClearAll["Global`*"]
Clear[distribution];
distribution = 
Compile[{{d, _Integer}, {n, _Integer}, {S, _Integer}, {z, _Integer}, \
{V, _Integer}},
  Module[{w, r, v, s, c, K},
   w = Table[0, {V}]; r = Table[0, {n}];
   Do[
    v = 1; s = 0;
    Do[
     c = RandomInteger[{0, d}];
     If[c != 0, s++; r[[k]] = c],
     {k, 1, n}];
    If[s == S, w[[v]]++];
    While[s != 0,
     v++; If[v == V, Break[]];
     K = s; s = 0;
     Do[
      c = Total[RandomInteger[{0, d}, r[[k]]]];
      If[c != 0, s++; r[[k]] = c],
      {k, 1, K}];
     If[s == S, w[[v]]++]],
    {z}]; Return[Table[{v, w[[v]]/z // N}, {v, 1, V - 1}]]],
  CompilationTarget -> "C"]
POSTED BY: Ulrich Utiger
10 Replies

Ulrich,

Glad this helps. I look forward to the article (although it is a bit outside of my area of expertise).

Also, I would change the title of this post to make it more useful to others.

Regards,

Neil

POSTED BY: Neil Singer
Posted 4 years ago

Thank you very much Neil for these suggestions.

Meanwhile, I also realized that for the end-result (the conditional probability that more than 1 lineages were present at the time of the mitochondrial Eve about 200'000 years ago, knowing that at present there is only 1) I don't need to calculate 200'000/20=10'000 generations because it's converging rapidly only after 20 generations or so. In other words, calculating more than 20 generations is not necessary because the result doesn't change significantly.

What I still need is to fit the code with paleoanthropological data, especially the probability that a mother has 0 daughters. This is not the same than in modern societies. My guess is that this is only conditioned by sterility. This is very important to know because if its probability is zero, lineages do not die out at all... I'll let you know when I finish my article.

POSTED BY: Ulrich Utiger

Ulrich,

Martijn gave you some great advice on speeding up your loops. I also suggest that you use SeedRandom[]. If you use SeedRandom you can get the same result every time. This helps with timing the routine and verifying that you do not inadvertently add a bug when improving your code.

For example I can do this and get consistent timing and always the same results:

SeedRandom[3]; AbsoluteTiming[
 distributionF1C[3, 5, 3, 10000, 10]]

{0.082987, {{1., 0.261}, {2., 0.329}, {3., 0.3443}, {4., 
   0.3518}, {5., 0.3516}, {6., 0.3493}, {7., 0.349}, {8., 
   0.3476}, {9., 0.3475}}}

You also asked

I could use the code

r = RandomInteger[{0, d}, n]; s = n - Count[r, 0]; only for the first generation because the lineages that do not die out must be preserved. I don't see how I could use it for generations V>1

The best way to avoid looping is to think in terms of lists. For example, you currently represent a generation with a list:

SeedRandom[3]; first = RandomInteger[3, 4]

Which is {1, 3, 2, 2}

You can then generate a new list with the next generation

Map[RandomInteger[3, #] &, first]

to get {{2}, {0, 0, 2}, {2, 0}, {0, 0}} but you only care about the total in each lineage so use this instead

Map[Total[RandomInteger[3, #]] &, first]

to get your next generation, {2, 2, 2, 0}

Now you can do this for V generations and get a list of every generation on the way using NestList

SeedRandom[3]; NestList[
 Function[x, Map[Total[RandomInteger[3, #]] &, x]], 
 RandomInteger[3, 4], 8 - 2]

to get every generation (7 of them) total for each lineage:

{{1, 3, 2, 2}, {2, 2, 2, 0}, {3, 5, 4, 0}, {7, 7, 2, 0}, {13, 7, 5, 0}, {23, 14, 6, 0}, {43, 23, 14, 0}}

I used V-2 because your code only did 7 generations for entering V=8 and NestList includes the first one automatically so V=1 gives 2 generations.

Now I just need to count up the totals that equal S and run a bunch of samples and divide by the number of samples so here is the final code:

mydist[d_, n_, S_, z_, V_] :=
 Total[Table[
    Map[Boole[(n - Count[#, 0]) == S] &, 
     NestList[Function[x, Map[Total[RandomInteger[d, #]] &, x]], 
      RandomInteger[d, n], V - 2]], z]]/N[z]

To test it I get the same answer as the compiled loop above in 10,000 samples but its 17% faster than the compiled loop version.

SeedRandom[3]; AbsoluteTiming[mydist[3, 5, 3, 10000, 10]]

{0.070848, {0.261, 0.329, 0.3443, 0.3518, 0.3516, 0.3493, 
  0.349, 0.3476, 0.3475}}

I tried to compile this loop and it did not speed up at all -- in fact its slightly slower! (I am not sure why but maybe MMA is already fast). Maybe Martijn or someone else might see if I did something wrong or if this code can be improved for compiling. Here is that code

mydistc := 
 Compile[{{d, _Integer}, {n, _Integer}, {S, _Integer}, {z, _Integer}, \
{V, _Integer}}, 
  Total[Table[
     Map[Boole[(n - Count[#, 0]) == S] &, 
      NestList[Function[x, Map[Total[RandomInteger[d, #]] &, x]], 
       RandomInteger[d, n], V - 2]], z]]/N[z]]

I hope this helps,

Regards,

Neil

POSTED BY: Neil Singer
Posted 4 years ago

Yes, I have VisualStudio installed. And to my surprise your code compiled flawlessly... Normally, when I did this in the past, it always refused some code that I needed to render C-compatible.

The next step is now to find a probability distribution for the amount of daughters because in reality this is not a uniform probability as rendered with RandomInteger. A mother, for instance, that has 0 daughters is much less probable than if she has 3, in any case this is true for hunter-gatherers. I tested it with RandomVariate and the binomial distribution, which yields quite a different end-result than with RandomInteger.

POSTED BY: Ulrich Utiger
Posted 4 years ago

Dear Martijn,

Thank you very much for having taken of your time to improve my code. This comes just at the right time because I realized that my analytical solution is only fast for the first 10 generations or so, the complexity increasing exponentially for more generations...

I used r[x] because I didn't know that RandomInteger can also be mapped with a pure function.

I will still try to C-compile it, even though I remember that this was always a very difficult task, and evaluate it parallely.

Very glad that it is about ten times faster now.

Regards, Ulrich

POSTED BY: Ulrich Utiger

As long as you have a C compiler installed on your system this is not difficult.

Just check if one is availible using

Needs["CCompilerDriver`"]
DefaultCCompiler[]
CCompilers[Full]

I quickly checked and compiling with CompilationTarget -> "C" of the compiled function in my previous post makes it faster.

(d3 = distributionF1C[3, 5, 3, 10000, 10]) // RepeatedTiming
(d3 = distributionF1C2[3, 5, 3, 10000, 10]) // RepeatedTiming

Out[17]= {0.12, {{1., 0.2657}, {2., 0.3424}, {3., 0.3544}, {4., 
   0.3604}, {5., 0.3613}, {6., 0.3611}, {7., 0.3586}, {8., 
   0.3576}, {9., 0.3574}}}

Out[18]= {0.083, {{1., 0.2731}, {2., 0.3305}, {3., 0.3393}, {4., 
   0.3391}, {5., 0.3397}, {6., 0.3402}, {7., 0.3395}, {8., 
   0.3399}, {9., 0.3393}}}
POSTED BY: Martijn Froeling
Posted 4 years ago

My code above has an error somewhere, because it does not correspond exactly to the analytical solution.

I could use the code

r = RandomInteger[{0, d}, n];
        s = n - Count[r, 0];

only for the first generation because the lineages that do not die out must be preserved. I don't see how I could use it for generations v>1.

Finally, I have done it with this code, which corresponds exactly to the analytical solution and is also somewhat faster :

distribution[d_, n_, S_, z_, V_] := Module[{w, r, v, g, s},
   w = ConstantArray[0, V]; r[x_] := RandomInteger[d, x];
   Do[
    v = 1;
    g = Select[RandomInteger[d, n], # > 0 &];
    s = Length[g];
    If[s == S, w[[v]]++];
    While[s != 0,
     v++; If[v == V, Break[]];
     g = Select[Total[Map[r, g], {2}], # > 0 &];
     s = Length[g];
     If[s == S, w[[v]]++]],
    {z}];
   Return[Table[{v, w[[v]]/z // N}, {v, 1, V - 1}]]];

The core of it is

g = Select[Total[Map[r, g], {2}], # > 0 &];

For instance, g={2,1,3} from the previous generation means that the 1st lineage has 2 daughters, the 2nd 1 and the 3rd 3. Map then transfers 2, 1 and 3 to the random function r defined above. For d=5 for instance, the result could be {{1,3},{0},{5,0,2}}, which means that the 1st daughter from the 1st lineage has 1 daughter and the 2nd 3. Idem for the other lineages. So the 1st lineage has 4 daughters in the next generation. This is calculated with Total[{{1,3},{0},{5,0,2}},{2}]={4,0,7}. Select removes the lineages that have no daughters and the result is {4,7}. The length of this list is the number s of lineages that survived.

If someone has an idea how it could still be improved, you are welcome.

POSTED BY: Ulrich Utiger

Select is generally slow, for small lists DeleteCases is faster for larger lists Pick will win.

time = Transpose[1000000 Table[
     b = RandomInteger[5, i];
     t1 = First@RepeatedTiming@Select[b, # > 0 &];
     t2 = First@RepeatedTiming@DeleteCases[b, 0];
     t3 = First@RepeatedTiming@Pick[b, Unitize[b], 1];
     {t1, t2, t3}
     , {i, 1, 20, 1}]];
ListLinePlot[time, PlotRange -> Full, 
 PlotLegends -> {"Select", "DeleteCases", "Pick"}, 
 AxesLabel -> {"N", 
   "time [\!\(\*SuperscriptBox[\(10\), \(-3\)]\) ms]"}]

time = Transpose[1000 Table[
     b = RandomInteger[5, Round[10^i]];
     t1 = First@RepeatedTiming@Select[b, # > 0 &];
     t2 = First@RepeatedTiming@DeleteCases[b, 0];
     t3 = First@RepeatedTiming@Pick[b, Unitize[b], 1];
     {t1, t2, t3}
     , {i, 1, 5, 1}]];
ListLinePlot[time, PlotRange -> Full, 
 PlotLegends -> {"Select", "DeleteCases", "Pick"}, 
 AxesLabel -> {"\!\(\*SuperscriptBox[\(10\), \(N\)]\)", "time [ms]"}]

enter image description here

Also i think the function call to r[x_] is not needed but not sure how much difference that makes. But compile still will gain you a lot.

distributionF1[daugthers_, mothers_, lin_, itt_, gen_] := 
 Block[{w, g, s, v},
  w = ConstantArray[0, gen - 1];

  Do[
   v = 1;
   g = DeleteCases[RandomInteger[daugthers, mothers], 0];
   s = Length[g];
   If[s == lin, w[[v]]++];

   While[s != 0,
    v++;
    If[v == gen, Break[]];
    g = DeleteCases[Total[Map[RandomInteger[daugthers, #] &, g], {2}],
       0];
    s = Length[g];
    If[s == lin, w[[v]]++];
    ];

   , {itt}];

  Thread[{Range[gen - 1], N[w/itt]}]
  ]

distributionF1C = 
  Compile[{{daugthers, _Integer}, {mothers, _Integer}, {lin, _Integer}, {itt, _Integer}, {gen, _Integer}},
   Block[{w, g, s, v},
    w = ConstantArray[0, gen - 1];
    Do[
     v = 1;
     g = DeleteCases[RandomInteger[daugthers, mothers], 0];
     s = Length[g];
     If[s == lin, w[[v]]++];

     While[s != 0,
      v++;
      If[v == gen, Break[]];
      g = DeleteCases[Map[Total[RandomInteger[daugthers, #]] &, g], 0];
      s = Length[g];
      If[s == lin, w[[v]]++];
      ];
     , {itt}];

    Thread[{Range[gen - 1], N[w/itt]}]
    ]
   ];

The gain is small with DeleteCases, but with compile its faster.

In[963]:= (d1 = distribution[3, 5, 3, 10000, 10]) // RepeatedTiming
(d2 = distributionF1[3, 5, 3, 10000, 10]) // RepeatedTiming
(d3 = distributionF1C[3, 5, 3, 10000, 10]) // RepeatedTiming

ListLinePlot[{d1, d2, d3}, PlotRange -> Full]

Out[963]= {1.09, {{1, 0.2633}, {2, 0.326}, {3, 0.3396}, {4, 
   0.3459}, {5, 0.3469}, {6, 0.346}, {7, 0.345}, {8, 0.3437}, {9, 
   0.3433}}}

Out[964]= {0.86, {{1, 0.2629}, {2, 0.3235}, {3, 0.3452}, {4, 
   0.3438}, {5, 0.345}, {6, 0.3447}, {7, 0.344}, {8, 0.3426}, {9, 
   0.3422}}}

Out[965]= {0.12, {{1., 0.2633}, {2., 0.3296}, {3., 0.3428}, {4., 
   0.3451}, {5., 0.3442}, {6., 0.3409}, {7., 0.3414}, {8., 
   0.342}, {9., 0.3428}}}

enter image description here

POSTED BY: Martijn Froeling
Posted 4 years ago

Hey Neil, thank you very much for your suggestions. I will try out your code and then get back.

POSTED BY: Ulrich Utiger

Ulrich,

Looping constructs like "Do" and "While" are not generally a good idea in MMA. Also, Table[0,{n}] is slow for large n. Use

ConstantArray[0, n]

To convince yourself of this do

AbsoluteTiming[Table[0, {10000000}];]
AbsoluteTiming[ConstantArray[0, 10000000];]

on my machine the ConstantArray is 6 times faster.

Now for looping.

r = ConstantArray[0, n];

 Do[c = RandomInteger[{0, d}];
  If[c != 0, s++; r[[k]] = c], {k, 1, n}];

is really slow (and I guess that it will use considerably more memory). Since (as far as I can tell) you are generating a list of random integers and counting up the non zero values, You can do the same operation very quickly using lists with

r = RandomInteger[{0, d}, n];
s = n - Count[r, 0];

By making just those two changes, I computed the average of a bunch of runs:

Mean[
 Table[AbsoluteTiming[distribution[5, 5, 3, 4, 7];], {100}]]

vs

Mean[
 Table[AbsoluteTiming[distributionF2[5, 5, 3, 4, 7];], {100}]]

The speedup is 3 orders of magnitude!

Here is the code I used. But I would not stop here. You should get rid of the rest of the looping and the code will run so fast that compiling may not be necessary.

distributionF2 = 
 Compile[{{d, _Integer}, {n, _Integer}, {S, _Integer}, {z, _Integer}, \
{V, _Integer}}, 
  Module[{w, r, v, s, c, K}, w = ConstantArray[0, V]; Do[v = 1;
    r = RandomInteger[{0, d}, n];
    s = n - Count[r, 0];
    If[s == S, w[[v]]++];
    While[s != 0, v++; If[v == V, Break[]];
     K = s; s = 0;
     Do[c = Total[RandomInteger[{0, d}, r[[k]]]];
      If[c != 0, s++; r[[k]] = c], {k, 1, K}];
     If[s == S, w[[v]]++]], {z}]; 
   Return[Table[{v, w[[v]]/z // N}, {v, 1, V - 1}]]], 
  CompilationTarget -> "C"]

In fact with this code I can run this (which is impossible with your previous code):

In[9]:= Mean[
 Table[AbsoluteTiming[distributionF2[10, 8, 3, 4, 10];], {100}]]

Out[9]= {0.0884735, Null}

I have not taken the time to understand what your other loops do but I suspect you can make similar changes.

Regards,

Neil

POSTED BY: Neil Singer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract