1
$\begingroup$

I need to do forward iterations of a system of 16 coupled recursion equations that track the frequency of genotypes over time in a biological model. I find that the iterations run slowly (it takes ~ 50+ seconds to run 10 generations) and I am looking for a way to speed them up because I need to run more like 20,000 generations. The code is pasted below. Each generation two populations exchange migrants, then undergo selection, mating, and recombination. I only need to track/memoize the genotype frequencies at the end of the generation, given by $x[k,gf,l,t]$, where $k$ designates the population (1 or 2), $gf$ the allele at one locus (1, 2, 3, or 4), $l$ the allele at a second locus (1 or 2), and $t$ indicated the generation.

*Currently running Mathematica 11.

ClearAll["Global`*"]

xMigration[k_, gf_, l_, t_, parameters_] := (1 - m) x[k, gf, l, t - 1, parameters] + m x[3 - k, gf, l, t - 1, parameters] /. parameters;

xMating[k_, gf1_, l1_, gf2_, l2_, t_, parameters_] := If[gf1 == 1 || gf1 == 3, If[gf2 == 1 || gf2 == 2, x[k, gf1, l1, t - 1, parameters]*xMigration[k, gf2, l2, t, parameters], 0], x[k, gf1, l1, t - 1, parameters]*xMigration[k, gf2, l2, t, parameters]]/Sum[xMigration[k, 1, L, t, parameters] + xMigration[k, 2, L, t, parameters] + KroneckerDelta[Mod[gf1, 2], 0] xMigration[k, 3, L, t, parameters] + KroneckerDelta[Mod[gf1, 2], 0] xMigration[k, 4, L, t, parameters], {L, 2}];


xSelection[k_, gf1_, l1_, gf2_, l2_, t_, parameters_] := 
  ((1 - (KroneckerDelta[l1, 3 - k] + KroneckerDelta[l2, 3 - k]) s) xMating[k, gf1, l1, gf2, l2, t, parameters])/Sum[Sum[Sum[Sum[(1 - (KroneckerDelta[L1, 3 - k] + KroneckerDelta[L2, 3 - k]) s) xMating[k, GF1, L1, GF2, L2, t, parameters], {GF1, 1, 4}], {L1, 1, 2}], {GF2, 1, 4}], {L2, 1, 2}] /. parameters

xRecombination[k_, gf1_, l1_, t_, parameters_] := Sum[If[GF1 == gf1 || GF2 == gf1,(*One of the chromosomes has to have the gametophytic factor haplotype gf1*)If[L1 == l1 || L2 == l1, (*One of the chromosomes has to have the trait allele l1*) If[GF1 == GF2 && L1 == L2, 1, (*If the chromosomes have the same genotype, all of their gametes will be this genotype*) If[GF1 == GF2 || L1 == L2, 1/2, (*If the chromosomes share an allele at the gf locus OR the trait locus, recombination doesn't matter and half of their gametes will be genotype gf l*)If[GF1 == gf1 && L1 == l1, 1/2 (1 - r),(*If one chromosome has the gamete genotype, half of nonrecombinant offspring will be be gf l*) If[GF2 == gf1 && L2 == l1, 1/2 (1 - r), (*If the other chromosome has the gamete genotype, half of nonrecombinant offspring will be gf l*)1/2 r (*If neither parents has genotype pi ti, it can only be made through recombination*)]]]], 0], 0] xSelection[k, GF1, L1, GF2, L2, t, parameters] /. parameters, {GF1, 4}, {L1, 2}, {GF2, 4}, {L2, 2}]

(*Memoization*)
x[k_, gf_, l_, t_, parameters_] := x[k, gf, l, t, parameters] = xRecombination[k, gf, l, t, parameters]

(*SL*)x[1, 1, 1, 0, parameters_] := 0.05;
(*ML*)x[1, 2, 1, 0, parameters_] := 0.95;
(*FL*)x[1, 3, 1, 0, parameters_] := 0;
(*gL*)x[1, 4, 1, 0, parameters_] := 0;
(*Sl*)x[1, 1, 2, 0, parameters_] := 0;
(*Ml*)x[1, 2, 2, 0, parameters_] := 0;
(*Fl*)x[1, 3, 2, 0, parameters_] := 0;
(*gl*)x[1, 4, 2, 0, parameters_] := 0;
(*SL*)x[2, 1, 1, 0, parameters_] := 0;
(*ML*)x[2, 2, 1, 0, parameters_] := 0;
(*FL*)x[2, 3, 1, 0, parameters_] := 0;
(*gL*)x[2, 4, 1, 0, parameters_] := 0;
(*Sl*)x[2, 1, 2, 0, parameters_] := 0;
(*Ml*)x[2, 2, 2, 0, parameters_] := 0;
(*Fl*)x[2, 3, 2, 0, parameters_] := 0;
(*gl*)x[2, 4, 2, 0, parameters_] := 1;

params = {m -> 0.01, s -> 0.5, r -> 0.01};

x[1, 1, 1, 10, params]//AbsoluteTiming

I had trouble formatting the code (it runs in one line off the edge of the screen). I can fix that if someone can recommend how to do so without ruining the ability to copy and paste into Mathematica.

In general, I want to optimize my code and see if it can compete with an equivalent model coded in 'R'. One specific concern that I have is that memoization is using too much memory after thousands of generations (it appears to slow down a little). To circumvent this I tried using nested Do loops and storing the output in a table, allowing Mathematica to only remember one generation of genotype frequencies. This did little to speed things up.

Any suggestions welcome! Let me know if I need to edit the post.

$\endgroup$
  • $\begingroup$ Another idea: have you tried NestList or RecurrenceTable? $\endgroup$ – Chris K Jan 20 '17 at 23:41
1
$\begingroup$

A quick change I made was to memoize xMigration, xMating, xSelection, and xRecombination. This resulted in a 100x speedup for t=5.

Your model is pretty complicated and the notation seems a bit excessive, so there might be more opportunities to optimize that I won't have time to look into.

$\endgroup$
  • $\begingroup$ This helped enormously. It now runs 5000 generations in ~250+ seconds. $\endgroup$ – biologyUser Jan 22 '17 at 3:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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