1
|
18904 Views
|
4 Replies
|
3 Total Likes
View groups...
Share

# Functional Programming for Non-linear Equations

Posted 11 years ago
 Hello All,I am on a journey to become adept at developing functional programs with Mathematica. I have a short prototypecode that I am using to learn how to convert procedural programs into functional programs and am wondering if anyonecould 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 {data}]; ) // AbsoluteTiming  (*Close Main For Loop How long does this loop take to run?*)
4 Replies
Sort By:
Posted 11 years ago
 Thank You!
Posted 11 years 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 11 years 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 11 years ago
 Is this question ill posed? I've noticed that no one has responded.