Group Abstract Group Abstract

Message Boards Message Boards

Functional Programming for Non-linear Equations

GROUPS:
Hello All,

I am on a journey to become adept at developing functional programs with Mathematica. I have a short prototype
code that I am using to learn how to convert procedural programs into functional programs and am wondering if anyone
could provide any more advise? Specifically I would like to learn how to write this code in a much more memory friendly manner.
And also learn a better way to run the main loop of the Markov Chain Monte Carlo simulation that minimizes the distance between
the experimental radioactive decay data and my model. Does anyone have any suggestions on this or how to develop into better functional programmer.

Thanks,
Pat
  Clear["Global`*"]
  
  texp = {0.2, 2.2, 4.0, 5.0, 6.0, 8.0, 11.0, 15.0, 18.0, 26.0, 33.0,
     39.0, 45.0};
  dexp = {35., 25., 22.1, 17.9, 16.8, 13.7, 12.4, 7.5, 4.9, 4.0, 2.4,
     1.4, 1.1};
  dt = Transpose[{texp, dexp}];
  dataPlot = ListPlot[
     dt,
    PlotStyle -> {Red, AbsolutePointSize[6]},
    Frame -> {True, True, False, False}, Axes -> False,
    FrameLabel -> {"Time", "Species A"},
    LabelStyle -> {FontSize -> 20, FormatType -> Bold},
    RotateLabel -> False,
    ImageSize -> {550}];
 
 runlength = 10000;
 sets = runlength/10.;
 testing = 10.0;
 Mean[dexp];
 StandardDeviation[dexp];
 stim = 0.01;
 
 chi2[ka_?NumberQ, kb_?NumberQ] := Block[{sol, A, B, c}, sol =
     NDSolve[
       { A'[t] == -ka*A[t],
        B'[t] == ka*A[t] - kb*B[t],
        c'[t] == kb*B[t],
        A[0] == 35.,
        B[0] == 0.,
        c[0] == 0.},
       {A, B, c},
       {t, 0., 100.}][[1]];
    Apply[Plus, (dexp - (A[t] /. sol /. t -> texp))^2]];
 
 DateString[];
 RandomSeed[AbsoluteTime[%]]; (* Set Random Seed based on current time *)
 
 Clear[deltaE]
 Clear[freeE]
 Clear[ka]
 Clear[kb]
 Clear[u, v]
 bestParamSets = {};
 thrownsets = {};
 (* These arrays are for storage *)
 ka = ConstantArray[0, {runlength}];
 kb = ConstantArray[0, {runlength}];
 deltaE = ConstantArray[0, {runlength}];
 freeE = ConstantArray[0, {runlength}];
 v = ConstantArray[0, {runlength}];
 u = ConstantArray[0, {runlength}];
 decayscore = ConstantArray[0, {runlength}];
 (**************************************************)
 (******** Starting the Metropolis Algorithm ********)
 (**************************************************)
 
 range = 10.0;(* Pick different initial guesses for param. sets *)
 weight = 1./10;
 ka[[0]] = RandomReal[{0.0, 0.2}];
 kb[[0]] = RandomReal[{0.0, 0.2}];
 least = 0.0004;
 max = 1.000;
 decayscore[[0]] = chi2[ka[[0]], kb[[0]]];
 beta = 1/.01;
 
 
 (tbl = Reap[
      Table[(*Open Main For Loop*)
   ka[[i]] = ka[[i - 1]] + weight*ka[[i - 1]]*RandomReal[{-1.0, 1.0}];
   kb[[i]] = kb[[i - 1]] + weight*kb[[i - 1]]*RandomReal[{-1.0, 1.0}];
       (*************************************************)
       (*This is where the evaluation of the parameter sets begins*)
       (*************************************************)
 
     If[
      least < ka[[i]] && ka[[i]] < max &&
      least < kb[[i]] && kb[[i]] < max,
      decayscore[[i]] = chi2[ka[[i]], kb[[i]]];(*This is the function evaluated for each parameter set*)
      deltaE[[i]] = decayscore[[i]] - decayscore[[i - 1]];(*The change in energy for each step*)
      freeE[[i]] = Exp[-beta*deltaE[[i]]];(*This is a free energy associated with each step*)
        v[[i]] = Min[1, freeE[[i]]];
        u[[i]] = RandomReal[{0, 1}];
        (**************************************************)
        (******** Metropolis Selection Algorithm ********)
        (**************************************************)
      If[u[[i]] < v[[i]],(*If u<v always accept that parameter set, else accept with probability u {data}]*);
        ka[[i + 1]] = ka[];
        kb[[i + 1]] = kb[];
      Sow[{True, decayscore[], v[], ka[], kb[], i}],(*bestParamSets*)
        ka[] = ka[[i - 1]];
        kb[] = kb[[i - 1]];
       ],
        ka[] = ka[[i - 1]];
        kb[] = kb[[i - 1]]
       ],
    {i, 1,runlength - 1}]][[2, 1]];
 bestParamSets = Cases[tbl, {True, data__} -> {data}];
 ) // AbsoluteTiming  (*Close Main For Loop How long does this loop take to run?*)
POSTED BY: Pat Mac
Answer
10 months ago
Is this question ill posed? I've noticed that no one has responded.
POSTED BY: Pat Mac
Answer
10 months ago
Usually questions are more specific. Your code is very big, you did not explain in detail what it does or where the problem is. Mostly people do not have much time to go through large code - I guess this is why there is no answer yet. Could you please condense your code and make you question more specific and descriptions more detailed?
POSTED BY: Vitaliy Kaurov
Answer
10 months ago
Hello. Very roughly, there are two ways for doing that.

1: Start again from scratch. List the inputs and foresee the output(s) as a function of the inputs : Output[ in1,in2...]. The program is a function composition that expresses how the inputs are progressively transformed into the output(s). In most cases, the Mathematica definition directly mirrors some mathematical formula. In many cases, it is preferable to define annex functions (subprograms) that do simpler things, in which case, your main program Output is a composition of these annex functions.

2: Convert the procedural program into a functional one. In the broad lines, you convert a sequence of instructions a=f[in1] ; b=g[in2] ; c=h[a,b] ; d=k … into a function composition by taking the last instruction, replacing the variables by their expressions and so forth :
d=k[ h[a, b] ]
d=k[ h[ f[in1], g[in2] ] ]
hence :
Output[v1_, v2_] := k[ h[ f[ v1 ], g[ v2 ] ] ]

Also, when an instruction is a repetition (Do, For, While), try to rewrite it by using Nest (simple iteration), Fold (iteration from a value and over a list), Inner or Outer (iteration over 2 lists). You may also transform selections (If, Which) into conditional definitions.

In all cases, you should take some time to learn Mathematica programming, especially rule based programming, which is a very convenient way to express how the inputs are transformed into the outputs, thanks to patterns and pattern-matching; learn also how to build data structures SomeDataType[ data1, data2...] which gather your data that rule based programs can then process as wholes. Rémi
POSTED BY: Remi Barrere
Answer
10 months ago
Thank You!
POSTED BY: Pat Mac
Answer
10 months ago