Scaling ODEModel initial value parameters

[JohnGoertz]

Since 0.5.2 (and this merge) symfit can now estimate the initial values of ODEModels. However, there's no clear way to provide the scale of those initial values. For example, the following code fits a logistic curve to some data. Important information is the growth rate r and the initial signal s1_0. This works fine if I give a good guess for s1_0, but I know that the local minima in the model are evenly distributed in the "log" space of the initial values. If I brute force the residuals, I find the global minima occurs at r = 0.8 and log10(s1_0) = -7.5, but there are other local minima at log10(s1_0) ~ -5.5, -6.5, -8.5, -9.5 with r = 0.6, 0.7, 0.9, and 1.0, respectively.

Is there a way around this? Perhaps sympy could be used more explicitly, setting s1_0 as a relation to an intermediate value that lives in log-space?

import numpy as np
import symfit as sf

signal1 = np.array([ 5.86039495e-03,  3.61153194e-03,  2.24330124e-03, -2.92966086e-03,
       -1.34148231e-03,  7.86443626e-04,  1.37743090e-03,  1.63747582e-05,
       -2.96159949e-03, -4.91545334e-03,  2.25744104e-03,  8.02806834e-03,
        6.86979013e-03, -3.99727245e-03, -5.43338189e-03, -4.08880449e-03,
        2.11456048e-03,  1.97120626e-02,  5.00087356e-02,  1.14839342e-01,
        2.26256041e-01,  3.75486015e-01,  5.36502495e-01,  6.94936368e-01,
        8.08767694e-01,  8.73566838e-01,  9.04374323e-01,  9.18962105e-01,
        9.33002509e-01,  9.26539485e-01,  9.22790406e-01,  9.27219537e-01,
        9.30125574e-01,  9.37782682e-01,  9.31860668e-01,  9.25248598e-01,
        9.27760548e-01,  9.21753351e-01,  9.28642494e-01,  9.33063441e-01,
        9.26756334e-01,  9.24954549e-01,  9.24742853e-01,  9.21482340e-01,
        9.26562524e-01,  9.23340814e-01,  9.32431892e-01,  9.27025424e-01,
        9.18785776e-01,  9.17351300e-01,  9.22062558e-01,  9.27517326e-01,
        9.25602973e-01,  9.27516517e-01,  9.26961663e-01,  9.31088965e-01,
        9.24303093e-01,  9.17609376e-01,  9.23493194e-01,  9.12655604e-01])
cycles = np.arange(len(signal1))

s1,c = sf.variables('s1,c')
r,K1,s1_0 = sf.parameters('r,K1,s1_0')
lg2_e = np.log2(np.exp(1))

s1_0.min = 1e-8
s1_0.value = 10**-7.5
s1_0.max = 1e-7
K1.min = 0.1
K1.value = 1
K1.max = 2
r.min = 0.1/lg2_e
r.value = 1/lg2_e
r.max = 2/lg2_e

logistic_eqn = {
    sf.D(s1,c): r*s1*(1-s1/K1)
}

logistic_model = sf.ODEModel(logistic_eqn, initial = {c:0.0,s1:s1_0})
logistic_fit = sf.Fit(logistic_model,c=cycles,s1=signal1)
fit_result = logistic_fit.execute()
#print(fit_result)
fit = logistic_model(c=cycles, **fit_result.params)

plt.plot(cycles, signal1, 'oC0')
plt.plot(cycles, fit.s1,'C0')
try:
    print(fit_result)
except:
    for param in [r,K1,s1_0]:
        print(f'{param}={fit_result.value(param):.2f}')

Here's the code for brute-forcing the residuals. Obviously, this is a very hard problem, particularly for a local minimizer...

k = 0.9190785714863928

logistic_eqn = {
    sf.D(s1,c): r*s1*(1-s1/K1)
}
inits = np.arange(-12,-2,0.1)
rates = np.arange(0.1,2,0.1)

rss = np.zeros([inits.size,rates.size])
for (i,j),_ in np.ndenumerate(rss):
    logistic_model = sf.ODEModel(logistic_eqn, initial = {c:0.0,s1:10**inits[i]})
    fit = logistic_model(c=cycles, K1=k,r=rates[j])
    residual = signal1-fit.s1
    rss[i,j] = np.sum(residual**2)

And plotting the result:

plt.subplots(1,1,figsize = [12,8])
plt.contourf(rates,inits,np.log10(rss),cmap='Blues_r')
plt.ylabel('log10 Initial')
plt.xlabel('Growth Rate (Base $e$)')
cbar = plt.colorbar()
cbar.ax.set_ylabel('Log10 RSS')

[pckroon]

To copy my comment I made on this in #279 so we have the complete discussion in one place: Very interesting point! Since it's an ODEModel we can't get the analytical Jacobian, and finding it with numerical differentiation in every step is too expensive. Some minimizers however (BFGS springs to mind) approximate it as they go along. Those should be able to pick up on the underlying structure. (Very few minimizers are truly naive). For normal models you can work around this by introducing an intermediate variable, but I don't expect that will work here due to how this is implemented. Instead of passing just a Parameter as initial value, you kind of want to be able to pass an expression that evaluates to a single value.

[tBuLi]

My thinking here is that indeed you should just leverage the symbolical power here and define an extra expression in your model dict. However, that is not currently supported yet because ODEModels are a special snowflake.

However, there is a pretty good workaround available, I recommend that you have a look at this example and adopt it to your problem. This way you should be able to define an extra component which lives in log space.

Unfortunately I don't have time right now to get into the details but let me know if you're having problems then I'll write down what I mean in more detail.

[pckroon]

Workaround aside, I still see a value in being able to provide expressions as initial values. (initial={t0: a+b}, for example). The questions are: 1) how does odeint deal with an array as initial value? 2) do we want to restrict those expressions to not allow Variables?