tBuLi / symfit

Symbolic Fitting; fitting as it should be.
http://symfit.readthedocs.org
MIT License
233 stars 17 forks source link

Add the example "ODEModels as subproblems using CallableNumericalModel" (Issue #204) #248

Closed antonykamp closed 5 years ago

antonykamp commented 5 years ago

I renamed the first example about the CallableNumerical Model and added a second example about ODEModels as subproblems using CallableNumericalModel.

This pull request fixes #204

tBuLi commented 5 years ago

Thanks for making this PR! The notebooks are a bit hard to review so I'll do that here:

There have been some changes to the (CallableNumericalModel) syntax which will make it even better, I'm sorry I should've mentioned this in #204.

Change

ode_model = ODEModel({D(y, x): k * y}, initial={x: 0.0, y: 1.0})
model = CallableNumericalModel({z: fun(ode_model)},independent_vars=[x], params=[k,k2])

To

w = Variable('w')
ode_model = ODEModel({D(y, x): k * y}, initial={x: 0.0, y: 1.0})
model_dict = {
    z: 2 * w + k2,
    w: ode_model(x=x, k=k)
}
model = CallableNumericalModel(model_dict, connectivity_mapping={z: {w, k2}, w: {x, k}})

This is a new syntax for initiating CallableNumericalModel which wasn't possible at the time this issue first arose. This means the wrapping function can be removed entirely. And the independent_vars and params keywords to CallableNumericalModel have been deprecated in favor of connectivity_mapping, I think from this example you can see why :).

Additionally, remove the extra space in your prints 'k = ', print will add a space automatically.

Finally, naming-wise I would prefer to name k2 something else, like b. k to me sounds like a rate constant, whereas b is traditionally the constant term in a x + b.

I apologize for not pointing this out earlier.

antonykamp commented 5 years ago

Everything is fine :D It's not a big problem to change the code a bit. Unfortunately, I didn't find information about the CallableNumericalModel, but this is another issue :D

I changed the printing thing and the naming of the parameters, but with I got an error by creating the second entry of model_dict: cannot determine truth value of Relational. Do you have an idea?

tBuLi commented 5 years ago

Hmm weird, using the exact code I gave above? Perhaps you can post your new code and the full traceback?

Edit: and yes, you are right, the CallableNumericalModel docs seem to be outdated :(.

antonykamp commented 5 years ago

The code:

from symfit import variables, Parameter, Fit, D, ODEModel, CallableNumericalModel
import numpy as np
import matplotlib.pyplot as plt

# Create data
x_data = np.linspace(0.0, 10.0, 1000)
a_expected = 0.6
b_expected = 10.0
y_data = 2 * np.exp(a_expected * x_data) + b_expected

# Initialise variables and parameters
w, x, y, z = variables('w, x, y, z')
a = Parameter('a', 0.0)
b = Parameter('b', 0.0)

# Define model
ode_model = ODEModel({D(y, x): a * y}, initial={x: 0.0, y: 1.0})
model_dict = {
    z: 2 * w + b,
    w: ode_model(x=x, a=a),
}
model = CallableNumericalModel(model_dict, connectivity_mapping={z: {w, b}, w: {x, a}})

# Apply model
fit = Fit(model,x=x_data, z=y_data)
fit_result = fit.execute()

And the Traceback:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-dbc1f7a64940> in <module>
     18 model_dict = {
     19     z: 2 * w + b,
---> 20     w: ode_model(x=x, a=a),
     21 }
     22 model = CallableNumericalModel(model_dict, connectivity_mapping={z: {w, b}, w: {x, a}})

~\Anaconda3\lib\site-packages\symfit\core\fit.py in __call__(self, *args, **kwargs)
   1699         """
   1700         Ans = variabletuple('Ans', self)
-> 1701         ans = Ans(*self.eval_components(*args, **kwargs))
   1702         return ans
   1703 

~\Anaconda3\lib\site-packages\symfit\core\fit.py in eval_components(self, *args, **kwargs)
   1646         else:
   1647             t_bigger = np.concatenate(
-> 1648                 (np.array([t_initial]), t_like[t_like > t_initial])
   1649             )
   1650             t_smaller = np.concatenate(

~\Anaconda3\lib\site-packages\sympy\core\relational.py in __nonzero__(self)
    227 
    228     def __nonzero__(self):
--> 229         raise TypeError("cannot determine truth value of Relational")
    230 
    231     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational
tBuLi commented 5 years ago

Ah yes, I forgot some important details. The line for w should be

w: lambda x, a: ode_model(x=x, a=a).y,

since models return an Ans object which is essentially a namedtuple with all the evaluated components. So we have to select the y component of the answer.

It would actually be nice if we recognize when a component returns an Ans so we could write

model_dict = {
    z: 2 * w + b,
    w: ode_model,
}

instead, since the lamba just to do the unpacking feels a bit redundant. But that should be a separate issue.

tBuLi commented 5 years ago

Thinking about it some more, w is redundant, we can safely use y to represent the output of the ODE. That is much clearer:

from symfit import variables, Parameter, Fit, D, ODEModel, CallableNumericalModel
import numpy as np
import matplotlib.pyplot as plt

# Create data
x_data = np.linspace(0.0, 10.0, 1000)
a_expected = 0.6
b_expected = 10.0
z_data = 2 * np.exp(a_expected * x_data) + b_expected

# Initialise variables and parameters
x, y, z = variables('x, y, z')
a = Parameter('a', 0.0)
b = Parameter('b', 0.0)

# Define model
ode_model = ODEModel({D(y, x): a * y}, initial={x: 0.0, y: 1.0})
model_dict = {
    z: 2 * y + b,
    y: lambda x, a: ode_model(x=x, a=a).y,
}
model = CallableNumericalModel(model_dict, connectivity_mapping={z: {y, b}, y: {x, a}})

# Apply model
fit = Fit(model, x=x_data, z=z_data)
fit_result = fit.execute()
antonykamp commented 5 years ago

I added the new syntax now. I think it is time to look into the real code of symfit :D

@tBuLi could you please look over the description of the example :)

tBuLi commented 5 years ago

The code looks very good right now. The text could use an additional read-through, because there is for example a space missing somewhere after a period. After you give it another read through to see if the text flows well, I can merge it. And let me know in #246 what you would like to work on in the core specifically :).

antonykamp commented 5 years ago

All right. I tried to improve the text to make the example clearer and the text more fluid.