ODE model fit - issues with the dimension of data arrays

[Stéphane Gorsse]

Hello,

I am new to using the lmfit library. I am trying to optimize parameters to a set of given data by fitting an ODE model. The code below works quite well, however, I have a problem with the dimension of my dataset. It seems that the sets of data must be one-dimensional (R,), however, to solve the ODE, the y-data must be (R,1). 

Many thanks in advance for your help.

Cheers

Stéphane

-------------------------------

from lmfit import minimize, Parameters, Parameter, report_fit, Model

from scipy.integrate import odeint

M = 3

b = 0.254e-9

alpha = 0.3

G = 83000

sigma_0 = 364

D = 20e-6

x = x_given[::50]  # x shape is (83,)

y = y_given[::50].reshape(-1,1) # y shape is (83,1)

def f(rho, x, k1, k2, zeta, n_g_s):

    """

    Derivative to be evaluated

    """

    return M*(1/(b*D)*(1-(n_g_s*(1-np.exp(-zeta/(n_g_s*b)*x)))/n_g_s)+k1/b*(rho)**(1/2)-k2*rho)

def g(x, rho_0, k1, k2, zeta, n_g_s):

    """

    Solving ODE.

    """

    rho = odeint(f, rho_0, x, args=(k1, k2, zeta, n_g_s))

    return rho

def iso(x,k1,k2,zeta,n_g_s):

    return sigma_0 + M*alpha*G*b*(g(x, rho_0, k1, k2, zeta, n_g_s))**(1/2)

def cine(x,zeta,n_g_s):

    return M*G*b/D*n_g_s*(1-np.exp(-zeta/(n_g_s*b)*x))

def flow(x,k1,k2,zeta,n_g_s):

    return iso(x,k1,k2,zeta,n_g_s) + cine(x,zeta,n_g_s)

mod = Model(flow)

pars = mod.make_params()

pars['k1'].set(value=0.01, min=0.001, max=0.1)

pars['k2'].set(value=4, min=1, max=9)

pars['zeta'].set(value=0.254e-7, min=0.254e-8, max=0.254e-6)

pars['n_g_s'].set(value=9, min=1, max=10)

result = mod.fit(y,pars,x=x)

print(result.fit_report())

plt.plot(x, y, 'bo')

plt.plot(x, result.init_fit, 'k--', label='initial fit')

plt.plot(x, result.best_fit, 'r-', label='best fit')

dely = result.eval_uncertainty(x=x, sigma=3)

plt.fill_between(x, result.best_fit-dely, result.best_fit+dely, color="#ABABAB",

                 label='3-$\sigma$ uncertainty band')

plt.legend(loc='best')

plt.show()

[Matt Newville]

On Tue, Jul 6, 2021 at 10:59 PM Stéphane Gorsse <sgo...@gmail.com> wrote:

Hello,

I am new to using the lmfit library. I am trying to optimize parameters to a set of given data by fitting an ODE model. The code below works quite well, however, I have a problem with the dimension of my dataset. It seems that the sets of data must be one-dimensional (R,), however, to solve the ODE, the y-data must be (R,1). 

Many thanks in advance for your help.

I'd probably use `np.ravel` or `np.squeeze`. 

--Matt 

[Stéphane Gorsse]

Thanks for your answer. I already tried ravel or flatten but it does not work unless I did not use it properly. Would you mind to precise your suggestion?

Stéphane

[Evan Groopman]

Hi, Stéphane.

Unfortunately, when I tried, your code was not executable. Therefore it is very difficult to reproduce the error you encountered or provide more precise feedback.

For instance, you did not provide x_given or y_given variables. When I generated dummy data in the shapes you suggest for x and y, the iso() function cannot execute because the rho_0 argument in g() is not defined locally in iso() or globally in your code. A minimal working example is necessary.

If you need x and y to have the same dimensionality, you could also try the opposite of squeeze or ravel, something like x = x[:, np.newaxis] or x.reshape(x.shape[0],-1).

- Evan

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/0bdb7ac1-a9af-4c6e-9890-ffcba3460b70n%40googlegroups.com.

[Stéphane Gorsse]

Hi Evan,

Thank you so much for your help. Attached is the set of data you can use.

Cheers

Stéphane

[Matt Newville]

On Wed, Jul 7, 2021 at 12:23 PM Stéphane Gorsse <sgo...@gmail.com> wrote:

Thanks for your answer. I already tried ravel or flatten but it does not work unless I did not use it properly. Would you mind to precise your suggestion?

Well, you're doing

       y = y_given[::50].reshape(-1,1) # y shape is (83,1)w  

So it appears that you sort of understand "silly array stride/dimension tricks".

Note that if you did. `y = y.squeeze()`, `y=y.flatten()`, and `y = y.ravel()` would both convert that back to a shape of (83, )

That's the sort of thing you'll want to use either at the return of your model function or maybe on the result from your ODE (or whatever else) calculation or whatever it is that claims to need extra dimensions.

I do not understand what "it does not work unless I did not use it properly" means.  You did not say what you tried.  The way I read "it does not work unless I did not use it properly" would suggest that it does work.  So, maybe I'm not understanding.  

But, yeah, you'll want to flatten/squeeze/ravel some arrays sot that value you return from your model function is 1-D.  How and where you do that is sort of up to you.

--Matt 

[Stéphane Gorsse]

I need a (83,1) to evaluate de ODE. If I use flatten/squeeze/ravel I got the "final = data + result.residual.reshape(data.shape)" " error cannot reshape array of size 6889 into shape (83,)"

Basically, I want to evaluate the parameters k1, k2, zeta and n_g_s by fitting a model to measured data x_given and y_given.

The model is called flow. It has 2 components iso and cine. And iso depends of an ODE.

The system of equations is shown in the attached file.

Stephane

- 

--
You received this message because you are subscribed to a topic in the Google Groups "lmfit-py" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lmfit-py/TV16gJDfuzE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbodY8xk-oWm-eBB%3D_tvVk3m3bP24Pp2rVWEn-N5hocntQ%40mail.gmail.com.

[Matt Newville]

Stephane,

On Thu, Jul 8, 2021 at 1:08 AM Stéphane Gorsse <sgo...@gmail.com> wrote:

I need a (83,1) to evaluate de ODE. If I use flatten/squeeze/ravel I got the "final = data + result.residual.reshape(data.shape)" " error cannot reshape array of size 6889 into shape (83,)"

This is all sort of hard for me to understand.  You seem to understand ODEs, fitting of data, and know a lot about using numpy indexing rules for manipulating the shapes of arrays.  

For lmfit.Model, the model function needs the data ("y" array) to be 1-dimensional.  The Model class uses a Parameters object, a collection of independent variables (of any type, actually), and a model function that will take the independent variables and the parameter models, does whatever calculation you want, and returns an array.  That array returned by the model function has to be 1-D, and of the same length as the data to be fit.

There really is not any more to it than that.  You'll have to write your model function to return values that correspond 1-to-1 with the data to be fit.  I think that I cannot give any more advice than that.

--Matt