I am trying to make a NonlinearModelFit
where the model is the result of a numerical calculation. My code works, but it is very slow and also seems to have a memory leak. I'd appreciate advice about how to improve the code, make it faster, and stop leaking.
A simplified version of my model is:
$HistoryLength=0;
model[gnIN_?NumberQ, gpIN_?NumberQ] :=
Module[{rep, eqs, sol, t},
rep = {gn -> gnIN, gp -> gpIN, n0 -> 10, p0 -> 0.01};
eqs = {
n'[t] == -gn (n0 + n[t]) f[t],
p'[t] == -gp ((p0 + p[t])) (1 - f[t]),
n[t] == p[t] +f[t],
n[0] == 10, p[0] == 10, f[0] == 0};
sol = NDSolve[eqs /. rep, {n, p, f}, {t, 0, 10}];
First[(n[#] + p[#])/(n[0] + p[0]) /. sol] &]
[Edit: maximum time in NDSolve
changed from 50 to 10. See comment below.] There are three coupled equations (one of them not a differential equation, so it could be removed easily). The output that I want for my model, however, is not n
, p
, or f
, but rather (n+p)/(n[0]+p[0])
. I have used a pure function to pass this out of the model, but I don't know if there's a better way.
We can make some fake data by
tmp = model[1,.1];
data = Table[{i, tmp[i]}, {i, 0, 50}];
I then call
MemoryInUse[]/10^6.
nlm=NonlinearModelFit[data, {model[gn,gp][t], gn>=0 && gp>=0},{gn,gp},t];
MemoryInUse[]/10^6.
This runs successfully, giving best fit parameters {gn -> 1.04935, gp -> 0.0998232}
, but it takes several minutes -- 376 seconds [Edit: revised to 68 seconds], from AbsoluteTiming
on my 3.2GHz i5 iMac running MMA 9. (In this case, I could clearly give it excellent initial guesses, but that's a bit harder in my real problem. As an aside, giving it initial guesses of {2,0.2}
causes it to take 450 seconds and return {1.48,0.098}
as best fits, badly missing on gn
.)
On this time that I'm running, the MemoryInUse
commands return 92.9
and 101.1
. Running it again will cause steady increase in memory usage, which is why I believe there is a memory leak.
Running Names["*$*"]
returns a list with many variables like n$10849.
Names["n*$*"]//Length
Names["p*$*"]//Length
returns 289
and 336
. I assume that all these n
and p
variables are being kept for some reason related to the return of the Module
. But if I run the NonlinearModelFit
again, the total memory continues to rise, but the number of n$*
and p$*
variables does not continue to increase -- it seems to stay very close to 300, indicating that some garbage collection is occurring. Also, the total memory used by those n$*
and p$*
variables, using symbolMemoryUsage
from here, does not explain the increase in memory usage of the kernel. That is,
Total[symbolMemoryUsage /@ Names["*n$*"]]
Total[symbolMemoryUsage /@ Names["*p$*"]]
returns 13872
and 16176
, which is much less than the 8 MB increase in memory use on running the NonlinearModelFit
. In my actual case, the increased memory use can reach over 5 GB (and grind my computer to a halt) before NonlinearModelFit
finishes. In this simple case, the leak is smaller, but it is still there. Running
MemoryInUse[]/10^6.
AbsoluteTiming[Do[model[1, .1], {i, 1, 1000}];]
MemoryInUse[]/10^6.
returns 104.161
, 4.34
, and 104.261
. So the memory increase is about 100 KB with each 1000 evaluations of model
. When I add a simple counter into model
to count the number of times it is evaluated during the NonlinearModelFit
, the answer is about 85000, indicating that the 8 MB increase in memory use is mostly due to a leak from model
. Similarly, the time to evaluate model
85k times should be about 370 seconds, which is close to the total evaluation time. The slowness of my code and the memory leak are all in my model
.
I understand that there's a limit to how fast a model with NDSolve
could possibly be, but I suspect it can be better than this. Previous discussions of Block
vs. Module
have left me concerned that Block
will do things I don't understand, but I haven't tried it. Any advice for coding practice (in particular, better ways to return the function from model
so it's in the correct format to put into NonlinearModelFit
), speed, or memory will be greatly appreciated.
EDIT: In case people have been dissuaded from working on this problem by the running time of my example code, I have changed the maximum time of NDSolve
in model
from 50 to 10, which makes the NonlinearModelFit
run in 68 seconds, rather than 376. The memory leak is still there and of the same magnitude as in the original question. Any comments or advice welcome!
model[1, .1]
complains withNDSolve::ndinnt : Initial condition f0 is not a number or a rectangular array of numbers.
Are you using a specific definition forf0
that you could edit into the question? $\endgroup$