ABC for parameter estimation in ODE-based model

lorenzoFabbri commented 6 years ago

Summary:

I am trying to use an ABC method in order to estimate the parameters of a system of ordinary differential equations.

Description:

I'd like to use ELFI and I'm using code gathered from the tutorials. The problem is the following: one of the variables of the model is computed from a conservation equation. That is, an ODE is related to a quantity which is computed as the difference between the total amount of that quantity (which has to be estimated) and the free amount of that quantity (which is passed to the function as an initial condition). The problem is that in doing so, I am taking the difference between a prior distribution and a float:

TypeError: In executing node '_simulator_feee': unsupported operand type(s) for -: 'Prior' and 'float'.

Reproducible Steps:

import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
from pandas import read_csv
mcn_observed = np.asarray(read_csv("../data/mcn_sol.csv", sep = ",", header = None))
mcn_observed = mcn_observed.flatten()

# This is the system of ODEs. Note the conservation equations
def MCN(states, t, CCP, params):

    # Parameters to be fitted
    k_SMPF,k_ASS,k_DISS,k_DRUM1,k_IRUM1,k_I2RUM1,k_1SLP1,V_2SLP1,J_1SLP1,J_2SLP1, \
    k_1IE,V_2IE,J_1IE,J_2IE,k_DMPFRUM1,V_SRUM1,k_ARUM1,k_DRUM1P,V_WEE1,k_WEE1, \
    J_1WEE1,J_2WEE1,WEE1T,k_CDC25,V_CDC25,J_1CDC25,J_2CDC25,mu,k_1WEE1,k_2WEE1, \
    k_1CDC25,k_2CDC25,k_D1CYC,k_D2CYC,CDC25T,IET,SLP1T,alpha,k_DX = params

    k_ASS    = 100 # 1
    k_I2RUM1 = 50  # 5
    k_ARUM1  = 35  # 16
    k_DRUM1P = 250 # 17

    # Variables of the system, coming from the initial conditions
    MPF, MPFP, SLP1A, IEA, MPFRUM1, RUM1, RUM1P, WEE1, CDC25P, MASS = states

    # Conservation equations
    SLP1  = SLP1T - SLP1A
    IE    = IET - IEA
    WEE1P = WEE1T - WEE1
    CDC25 = CDC25T - CDC25P

    dMPFdt = k_SMPF*MASS - k_WEE1*MPF + k_CDC25*MPFP - \
            (k_D1CYC+k_D2CYC*SLP1A)*MPF - k_ASS*RUM1*MPF + \
            + (k_DISS+k_DRUM1+k_IRUM1)*MPFRUM1 + k_DX*CCP*MPFRUM1
    dMPFPdt = k_WEE1*MPF - k_CDC25*MPFP - (k_D1CYC+k_D2CYC*SLP1A)*MPFP
    dSLP1Adt = k_1SLP1*IEA*SLP1/(J_1SLP1+SLP1) - \
            V_2SLP1*SLP1A/(J_2SLP1+SLP1A)
    dIEAdt = k_1IE*(MPF+alpha*MPFP)*IE/(J_1IE+IE) - V_2IE*IEA/(J_2IE+IEA)
    dMPFRUM1dt = k_ASS*RUM1*MPF - \
            (k_DISS+k_DRUM1+k_IRUM1+k_DMPFRUM1)*MPFRUM1 - \
            k_DX*CCP*MPFRUM1
    dRUM1dt = V_SRUM1 - k_ASS*RUM1*MPF + (k_DISS+k_DMPFRUM1)*MPFRUM1 - \
            k_I2RUM1*MPFP*RUM1 + k_ARUM1*RUM1P - k_DRUM1*RUM1 - \
            k_DX*CCP*RUM1 + k_DX*CCP*MPFRUM1 - k_DX*CCP*RUM1
    dRUM1Pdt = k_IRUM1*MPFRUM1 + k_I2RUM1*MPFP*RUM1 - \
            k_ARUM1*RUM1P - k_DRUM1P*(MPF+alpha*MPFP)*RUM1P - \
            k_DRUM1*RUM1P + k_DX*CCP*RUM1 + k_DX*CCP*MPFRUM1 - \
            k_DX*CCP*RUM1P
    dWEE1dt = V_WEE1*WEE1P/(J_1WEE1+WEE1P) - \
            k_WEE1*(MPF+alpha*MPFP)*WEE1/(J_2WEE1+WEE1)
    dCDC25Pdt = k_CDC25*(MPF+alpha*MPFP)*CDC25/(J_1CDC25+CDC25) - \
            V_CDC25*CDC25P/(J_2CDC25+CDC25P)
    dMASSdt = mu*MASS

    system = [dMPFdt, dMPFPdt, dSLP1Adt, dIEAdt, dMPFRUM1dt, \
              dRUM1dt, dRUM1Pdt, dWEE1dt, dCDC25Pdt, dMASSdt]

    return system

# Priors
k_SMPF     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_SMPF')
k_ASS      = elfi.Prior('uniform', 0.0, 1.0, name = 'k_ASS')
k_DISS     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DISS')
k_DRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DRUM1')
k_IRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_IRUM1')
k_I2RUM1   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_I2RUM1')
k_1SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1SLP1')
V_2SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_2SLP1')
J_1SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1SLP1')
J_2SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2SLP1')
k_1IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1IE')
V_2IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'V_2IE')
J_1IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1IE')
J_2IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2IE')
k_DMPFRUM1 = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DMPFRUM1')
V_SRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_SRUM1')
k_ARUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_ARUM1')
k_DRUM1P   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DRUM1P')
V_WEE1     = elfi.Prior('uniform', 0.0, 1.0, name = 'V_WEE1')
k_WEE1     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_WEE1')
J_1WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1WEE1')
J_2WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2WEE1')
WEE1T      = elfi.Prior('uniform', 0.0, 1.0, name = 'WEE1T')
k_CDC25    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_CDC25')
V_CDC25    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_CDC25')
J_1CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1CDC25')
J_2CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2CDC25')
mu         = elfi.Prior('uniform', 0.0, 1.0, name = 'mu')
k_1WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1WEE1')
k_2WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_2WEE1')
k_1CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1CDC25')
k_2CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_2CDC25')
k_D1CYC    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_D1CYC')
k_D2CYC    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_D2CYC')
CDC25T     = elfi.Prior('uniform', 0.0, 1.0, name = 'CDC25T')
IET        = elfi.Prior('uniform', 0.0, 1.0, name = 'IET')
SLP1T      = elfi.Prior('uniform', 0.0, 1.0, name = 'SLP1T')
alpha      = elfi.Prior('uniform', 0.0, 1.0, name = 'alpha')
k_DX       = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DX')

# Function for the simulator
def sim(params, batch_size = 1, random_state = None):

    # Initial conditions
    ivp = np.random.randint(low = 0, high = 5, size = 10)

    # Time
    stop_time = 500
    t = np.linspace(0, stop_time, stop_time*10)

    # Odeint
    res = odeint(MCN, ivp, t, args = (0, params))
    res = np.asarray(result).flatten()

    return res

params = list((k_SMPF,k_ASS,k_DISS,k_DRUM1,k_IRUM1,k_I2RUM1,k_1SLP1,V_2SLP1,J_1SLP1,J_2SLP1, \
k_1IE,V_2IE,J_1IE,J_2IE,k_DMPFRUM1,V_SRUM1,k_ARUM1,k_DRUM1P,V_WEE1,k_WEE1, \
J_1WEE1,J_2WEE1,WEE1T,k_CDC25,V_CDC25,J_1CDC25,J_2CDC25,mu,k_1WEE1,k_2WEE1, \
k_1CDC25,k_2CDC25,k_D1CYC,k_D2CYC,CDC25T,IET,SLP1T,alpha,k_DX))

elfi_sim = elfi.Simulator(sim, params, observed = mcn_observed)
elfi_sim.generate()

When I run elfi_sim.generate(), I get the error. How can I deal with this situation? Moreover, when I first tried to use ELFI for this very same problem, I did not get this error but it was complaining about the fact that res from the odeint step was not 1D.

ELFI Version: latest version

Python Version: 3.5

Operating System: Mac OS

vuolleko commented 6 years ago

Hi,

The code reads an input file and so this is not directly reproducible. But a couple of notes:

you forgot to unpack the params given to the simulator, i.e. elfi_sim = elfi.Simulator(sim, *params, observed = mcn_observed); now the priors are not set as parents to the simulator
I suggest using elfi.draw to check your ELFI model (you will need Graphviz)
a simple way to subtract a float from a prior is to define an operation node for that, i.e. modified_prior = elfi.Operation(lambda p: p - my_float, original_prior) and then set this node as the parent for the simulator
however, the previous is not needed once the parents are set properly, as ELFI will then pass values to the simulator (as arrays)

lorenzoFabbri commented 6 years ago

Thank you very much for the fast reply. I still have some problems: I changed the code for elfi.Simulator and then I had to change the function sim replacing params with the "un-packed" version (I wrote all of them). But I still get that problem with the subtraction.

import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
from pandas import read_csv
mcn_observed = np.asarray(read_csv("../data/mcn_sol.csv", sep = ",", header = None))
mcn_observed = mcn_observed.flatten()

# This is the system of ODEs. Note the conservation equations
def MCN(states, t, CCP, params):

    # Parameters to be fitted
    k_SMPF,k_ASS,k_DISS,k_DRUM1,k_IRUM1,k_I2RUM1,k_1SLP1,V_2SLP1,J_1SLP1,J_2SLP1, \
    k_1IE,V_2IE,J_1IE,J_2IE,k_DMPFRUM1,V_SRUM1,k_ARUM1,k_DRUM1P,V_WEE1,k_WEE1, \
    J_1WEE1,J_2WEE1,WEE1T,k_CDC25,V_CDC25,J_1CDC25,J_2CDC25,mu,k_1WEE1,k_2WEE1, \
    k_1CDC25,k_2CDC25,k_D1CYC,k_D2CYC,CDC25T,IET,SLP1T,alpha,k_DX = params

    k_ASS    = 100 # 1
    k_I2RUM1 = 50  # 5
    k_ARUM1  = 35  # 16
    k_DRUM1P = 250 # 17

    # Variables of the system, coming from the initial conditions
    MPF, MPFP, SLP1A, IEA, MPFRUM1, RUM1, RUM1P, WEE1, CDC25P, MASS = states

    # Conservation equations
    SLP1  = SLP1T - SLP1A
    IE    = IET - IEA
    WEE1P = WEE1T - WEE1
    CDC25 = CDC25T - CDC25P

    dMPFdt = k_SMPF*MASS - k_WEE1*MPF + k_CDC25*MPFP - \
            (k_D1CYC+k_D2CYC*SLP1A)*MPF - k_ASS*RUM1*MPF + \
            + (k_DISS+k_DRUM1+k_IRUM1)*MPFRUM1 + k_DX*CCP*MPFRUM1
    dMPFPdt = k_WEE1*MPF - k_CDC25*MPFP - (k_D1CYC+k_D2CYC*SLP1A)*MPFP
    dSLP1Adt = k_1SLP1*IEA*SLP1/(J_1SLP1+SLP1) - \
            V_2SLP1*SLP1A/(J_2SLP1+SLP1A)
    dIEAdt = k_1IE*(MPF+alpha*MPFP)*IE/(J_1IE+IE) - V_2IE*IEA/(J_2IE+IEA)
    dMPFRUM1dt = k_ASS*RUM1*MPF - \
            (k_DISS+k_DRUM1+k_IRUM1+k_DMPFRUM1)*MPFRUM1 - \
            k_DX*CCP*MPFRUM1
    dRUM1dt = V_SRUM1 - k_ASS*RUM1*MPF + (k_DISS+k_DMPFRUM1)*MPFRUM1 - \
            k_I2RUM1*MPFP*RUM1 + k_ARUM1*RUM1P - k_DRUM1*RUM1 - \
            k_DX*CCP*RUM1 + k_DX*CCP*MPFRUM1 - k_DX*CCP*RUM1
    dRUM1Pdt = k_IRUM1*MPFRUM1 + k_I2RUM1*MPFP*RUM1 - \
            k_ARUM1*RUM1P - k_DRUM1P*(MPF+alpha*MPFP)*RUM1P - \
            k_DRUM1*RUM1P + k_DX*CCP*RUM1 + k_DX*CCP*MPFRUM1 - \
            k_DX*CCP*RUM1P
    dWEE1dt = V_WEE1*WEE1P/(J_1WEE1+WEE1P) - \
            k_WEE1*(MPF+alpha*MPFP)*WEE1/(J_2WEE1+WEE1)
    dCDC25Pdt = k_CDC25*(MPF+alpha*MPFP)*CDC25/(J_1CDC25+CDC25) - \
            V_CDC25*CDC25P/(J_2CDC25+CDC25P)
    dMASSdt = mu*MASS

    system = [dMPFdt, dMPFPdt, dSLP1Adt, dIEAdt, dMPFRUM1dt, \
              dRUM1dt, dRUM1Pdt, dWEE1dt, dCDC25Pdt, dMASSdt]

    return system

# Priors
k_SMPF     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_SMPF')
k_ASS      = elfi.Prior('uniform', 0.0, 1.0, name = 'k_ASS')
k_DISS     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DISS')
k_DRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DRUM1')
k_IRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_IRUM1')
k_I2RUM1   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_I2RUM1')
k_1SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1SLP1')
V_2SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_2SLP1')
J_1SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1SLP1')
J_2SLP1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2SLP1')
k_1IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1IE')
V_2IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'V_2IE')
J_1IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1IE')
J_2IE      = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2IE')
k_DMPFRUM1 = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DMPFRUM1')
V_SRUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_SRUM1')
k_ARUM1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_ARUM1')
k_DRUM1P   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DRUM1P')
V_WEE1     = elfi.Prior('uniform', 0.0, 1.0, name = 'V_WEE1')
k_WEE1     = elfi.Prior('uniform', 0.0, 1.0, name = 'k_WEE1')
J_1WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1WEE1')
J_2WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2WEE1')
WEE1T      = elfi.Prior('uniform', 0.0, 1.0, name = 'WEE1T')
k_CDC25    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_CDC25')
V_CDC25    = elfi.Prior('uniform', 0.0, 1.0, name = 'V_CDC25')
J_1CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'J_1CDC25')
J_2CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'J_2CDC25')
mu         = elfi.Prior('uniform', 0.0, 1.0, name = 'mu')
k_1WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1WEE1')
k_2WEE1    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_2WEE1')
k_1CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_1CDC25')
k_2CDC25   = elfi.Prior('uniform', 0.0, 1.0, name = 'k_2CDC25')
k_D1CYC    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_D1CYC')
k_D2CYC    = elfi.Prior('uniform', 0.0, 1.0, name = 'k_D2CYC')
CDC25T     = elfi.Prior('uniform', 0.0, 1.0, name = 'CDC25T')
IET        = elfi.Prior('uniform', 0.0, 1.0, name = 'IET')
SLP1T      = elfi.Prior('uniform', 0.0, 1.0, name = 'SLP1T')
alpha      = elfi.Prior('uniform', 0.0, 1.0, name = 'alpha')
k_DX       = elfi.Prior('uniform', 0.0, 1.0, name = 'k_DX')

# Function for the simulator
def sim(k_SMPF,k_ASS,k_DISS,k_DRUM1,k_IRUM1,k_I2RUM1,k_1SLP1,V_2SLP1,J_1SLP1,J_2SLP1, \
k_1IE,V_2IE,J_1IE,J_2IE,k_DMPFRUM1,V_SRUM1,k_ARUM1,k_DRUM1P,V_WEE1,k_WEE1, \
J_1WEE1,J_2WEE1,WEE1T,k_CDC25,V_CDC25,J_1CDC25,J_2CDC25,mu,k_1WEE1,k_2WEE1, \
k_1CDC25,k_2CDC25,k_D1CYC,k_D2CYC,CDC25T,IET,SLP1T,alpha,k_DX, batch_size = 1, random_state = None):

    # Initial conditions
    ivp = np.random.randint(low = 0, high = 5, size = 10)

    # Time
    stop_time = 500
    t = np.linspace(0, stop_time, stop_time*10)

    # Odeint
    res = odeint(MCN, ivp, t, args = (0, params))
    res = np.asarray(result).flatten()

    return res

params = list((k_SMPF,k_ASS,k_DISS,k_DRUM1,k_IRUM1,k_I2RUM1,k_1SLP1,V_2SLP1,J_1SLP1,J_2SLP1, \
k_1IE,V_2IE,J_1IE,J_2IE,k_DMPFRUM1,V_SRUM1,k_ARUM1,k_DRUM1P,V_WEE1,k_WEE1, \
J_1WEE1,J_2WEE1,WEE1T,k_CDC25,V_CDC25,J_1CDC25,J_2CDC25,mu,k_1WEE1,k_2WEE1, \
k_1CDC25,k_2CDC25,k_D1CYC,k_D2CYC,CDC25T,IET,SLP1T,alpha,k_DX))

elfi_sim = elfi.Simulator(sim, *params, observed = mcn_observed)
elfi_sim.generate()

# Distance function and run ABC
def SSE(x, y):
    return np.atleast_1d(np.sum(x - y)**2)

elfi_dis = elfi.Distance(SSE, elfi_sim)
elfi_dis.generate()

elfi_rej = elfi.Rejection(elfi_dis, batch_size = 1)
result = elfi_rej.sample(1000, quantile = 0.01)

mcn_sol.csv.zip

vuolleko commented 6 years ago

The input arguments to sim are unused. The params list you pass to odeint is just the list of Prior nodes, set outside the function. Try setting the signature to sim(*params, batch_size = 1, random_state = None)

vuolleko commented 5 years ago

I assume this was resolved. Closing.

elfi-dev / elfi