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.
NestListorRecurrenceTable? $\endgroup$