search

lisphilar / covid19-sir

CovsirPhy: Python library for COVID-19 analysis with phase-dependent SIR-derived ODE models.
https://lisphilar.github.io/covid19-sir/
Apache License 2.0
109 stars 44 forks source link

[Question/Docs] Error when editing _ODESolver #1107

Closed shik-design closed 1 year ago

shik-design commented 2 years ago

Summary of question

I tried to integrate the models particularly SIR using the Fractional ODE package but is giving bug

(Optional) What tried

!/usr/bin/env python

-- coding: utf-8 --

import numpy as np import pandas as pd from scipy.integrate import solve_ivp import fodeint as fd import scipy.special as special from shiks.util.term import Term from shiks.ode.mbase import ModelBase

class _ODESolver(Term): """ Solve initial value problems for a SIR-derived ODE model.

Args:
    model (shiks.ModelBase): SIR-derived ODE model
    kwargs: values of non-dimensional model parameters, including rho and sigma

Note:
    We can check non-dimensional model parameters with model.PARAMETERS class variable.
    All non-dimensional parameters must be specified with keyword arguments.
"""

def __init__(self, model, **kwargs):
    self._model = self._ensure_subclass(model, ModelBase, name="model")
    param_dict = {k: float(v) for (k, v) in kwargs.items() if isinstance(v, (float, int))}
    self._param_dict = self._ensure_kwargs(model.PARAMETERS, float, **param_dict)

def run(self, step_n, **kwargs):
    """
    Solve an initial value problem.

    Args:
        step_n (int): the number of steps
        kwargs: initial values of dimensional variables, including Susceptible

    Returns:
        pandas.DataFrame: numerical solution
            Index
                reset index: time steps
            Columns
                (int): dimensional variables of the model

    Note:
        We can check dimensional variables with model.VARIABLES class variable.
        All dimensional variables must be specified with keyword arguments.
        Total value of initial values will be regarded as total population.
    """
    # Check arguments
    step_n = self._ensure_natural_int(step_n, name="number")
    kwargs = {param: int(value) for (param, value) in kwargs.items()}
    y0_dict = self._ensure_kwargs(self._model.VARIABLES, int, **kwargs)
    # Calculate population
    population = sum(y0_dict.values())
    # Solve problem
    return self._run(step_n=step_n, y0_dict=y0_dict, population=population)

def _run(self, step_n, y0_dict, population):
    """
    Solve an initial value problem for a SIR-derived ODE model.

    Args:
        step_n (int): the number of steps
        y0_dict (dict[str, int]): initial values of dimensional variables, including Susceptible
        population (int): total population

    Returns:
        pandas.DataFrame: numerical solution
            Index
                reset index: time steps
            Columns
                (int): dimensional variables of the model
    """
    tstart, dt, tend, alpha = 0, 1, step_n, 0.5
    variables = self._model.VARIABLES[:]
    initials = [y0_dict[var] for var in variables]
    """
    sol = solve_ivp(
        fun=self._model(population=population, **self._param_dict),
       t_span=[tstart, tend],
        y0=np.array(initials, dtype=np.int64),
        t_eval=np.arange(tstart, tend + dt, dt),
        dense_output=False
    )
    """
    sol = fd.fodeint(
        alpha=alpha
        deriv=self._model(population=population, **self._param_dict),
        x_initial=np.array(initials, dtype=np.int64),
        t=np.arange(tstart, tend + dt, dt),
    )
    y_df = pd.DataFrame(data=sol["y"].T.copy(), columns=variables)
    return y_df.round().astype(np.int64)
import covsirphy as cs

(Optional) Environment

lisphilar commented 2 years ago

Do you mean you editted the source code file to use fodent instead of scipy and got a error message? Kindly share the error message.

shik-design commented 2 years ago

Yes, I edited the source code file to investigate with fodeint instead of from scipy.integrate import solve_ivp

lisphilar commented 2 years ago

Thank you for your confirmation. What was the error message?

shik-design commented 2 years ago 
from sys import path as syspath
from os import path as ospath
syspath.append(ospath.join(ospath.expanduser("~"), 'C:/Users/USER/shiks'))
syspath.append(ospath.join(ospath.expanduser("~"), 'shiks'))
import pandas as pd
import shiks as cs

Traceback (most recent call last):

File "C:\Users\USER\anaconda3\envs\v-shiks\lib\site-packages\IPython\core\interactiveshell.py", line 3444, in run_code exec(code_obj, self.user_global_ns, self.user_ns)

File "C:\Users\USER\AppData\Local\Temp/ipykernel_11548/2979858712.py", line 14, in import shiks as cs

File "C:\Users\USER\shiks__init__.py", line 60, in from shiks.ode.ode_handler import ODEHandler

File "C:\Users\USER\shiks\ode\ode_handler.py", line 14, in from shiks.ode.ode_solver_multi import _MultiPhaseODESolver

File "C:\Users\USER\shiks\ode\ode_solver_multi.py", line 7, in from shiks.ode.ode_solver import _ODESolver

File "C:\Users\USER\shiks\ode\ode_solver.py", line 91 deriv=self._model(population=population, **self._param_dict), ^ SyntaxError: invalid syntax

shik-design commented 2 years ago

Kind Request

The source code files are attached for your perusal, kindly help me sort it out. I want to use it for my research.

The Source Code Files

shiks.zip

The Notebook

from sys import path as syspath
from os import path as ospath
syspath.append(ospath.join(ospath.expanduser("~"), 'C:/Users/USER/shiks'))
syspath.append(ospath.join(ospath.expanduser("~"), 'shiks'))
import pandas as pd
import shiks as cs

# Download and update datasets
data_loader = cs.DataLoader("input", update_interval=70090)
data_loader.lock(
    date="date", country="country", province="province",
    confirmed="confirmed", fatal="fatal", population="population",
)
jhu_data = data_loader.jhu()

# # Select country name and register the data
snl = cs.Scenario(country="Nigeria", province=None)
snl.register(jhu_data)

# Check records
snl.records()

# S-R trend analysis
snl.trend().summary()

# Parameter estimation of SIR model
s_I_R=snl.estimate(cs.SIR)

# History of reproduction number
_ = snl.history(target="Rt")
# History of parameters
_ = snl.history_rate()
_ = snl.history(target="rho")
lisphilar commented 2 years ago

To fix SyntaxError, please add comma to the second line of the following codes.

sol = fd.fodeint(
        alpha=alpha
        deriv=self._model(population=population, **self._param_dict),
        x_initial=np.array(initials, dtype=np.int64),
        t=np.arange(tstart, tend + dt, dt),
    )

Updated:

sol = fd.fodeint(
        alpha=alpha,
        deriv=self._model(population=population, **self._param_dict),
        x_initial=np.array(initials, dtype=np.int64),
        t=np.arange(tstart, tend + dt, dt),
    )
shik-design commented 2 years ago

Thank you for prompt response, I really appreciate. Please, exercise some patience with me. I am not an expert in python/ data science. Why trying to edit the code an error pop out.

Code Run

Parameter estimation of SIR model

s_I_R=snl.estimate(cs.SIR)

The Error

RemoteTraceback Traceback (most recent call last) RemoteTraceback: """

Traceback (most recent call last):
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\multiprocessing\pool.py", line 121, in worker
    result = (True, func(*args, **kwds))
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\multiprocessing\pool.py", line 44, in mapstar
    return list(map(*args))
  File "C:\Users\USER\shiks\ode\ode_handler.py", line 122, in _score_tau
    sim_df = solver.simulate(*info_dict.values())
  File "C:\Users\USER\shiks\ode\ode_solver_multi.py", line 93, in simulate
    solved_df = solver.run(step_n=step_n, **y0_dict)
  File "C:\Users\USER\shiks\ode\ode_solver.py", line 58, in run
    return self._run(step_n=step_n, y0_dict=y0_dict, population=population)
  File "C:\Users\USER\shiks\ode\ode_solver.py", line 93, in _run
    tspan=np.arange(tstart, tend + dt, dt),
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\site-packages\fodeint\fodeint.py", line 90, in fodeint
    (d, a, f, y0, tspan) = _check_args(a, f, y0, tspan)
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\site-packages\fodeint\fodeint.py", line 61, in _check_args
    if len(f(y0, tspan[0])) != d:
  File "C:\Users\USER\shiks\ode\sir.py", line 67, in __call__
    s,i,r =X
TypeError: cannot unpack non-iterable numpy.int32 object

"""

lisphilar commented 2 years ago

This error could be from the difference of scipy.integrate.solve_ivp (covsirphy uses) and fodeint.fodeint.

scipy.integrate.solve_ivp: The function must be defined as fun(t, y) https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

fodeint.fodeint: The function must be defined as fun(y, t) https://github.com/mattja/fodeint/blob/d4f06a16b8e82557ada6924a2e13f7fe0960f709/fodeint/fodeint.py#L74-L75

So, covsirphy.SIR.__call__(self, t, X) is suitable for scipy.integrate.solve_ivp, but un-suitable for fodeint.fodeint.

To use fodeint.fodeint, you may need to change the order of arguments (t, X -> X, t).

Before:

def __call__(self, t, X):
        """
        Return the list of dS/dt (tau-free) etc.
        Args:
            t (int): time steps
            X (numpy.array): values of th model variables
        Returns:
            (np.array)
        """
        n = self.population
        s, i, *_ = X
        dsdt = 0 - self.rho * s * i / n
        drdt = self.sigma * i
        didt = 0 - dsdt - drdt
        return np.array([dsdt, didt, drdt])

After:

def __call__(self, X, t):
        """
        Return the list of dS/dt (tau-free) etc.
        Args:
            X (numpy.array): values of th model variables
            t (int): time steps
        Returns:
            (np.array)
        """
        n = self.population
        s, i, *_ = X
        dsdt = 0 - self.rho * s * i / n
        drdt = self.sigma * i
        didt = 0 - dsdt - drdt
        return np.array([dsdt, didt, drdt])
shik-design commented 2 years ago

Thanks a lot. Your observation was correct. When I switched the order of arguments (t, X -> X, t), it worked. However, I got an error below:

RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\multiprocessing\pool.py", line 121, in worker
    result = (True, func(*args, **kwds))
  File "C:\Users\USER\anaconda3\envs\v-shiks\lib\multiprocessing\pool.py", line 44, in mapstar
    return list(map(*args))
  File "C:\Users\USER\shiks\ode\ode_handler.py", line 122, in _score_tau
    sim_df = solver.simulate(*info_dict.values())
  File "C:\Users\USER\shiks\ode\ode_solver_multi.py", line 93, in simulate
    solved_df = solver.run(step_n=step_n, **y0_dict)
  File "C:\Users\USER\shiks\ode\ode_solver.py", line 58, in run
    return self._run(step_n=step_n, y0_dict=y0_dict, population=population)
  File "C:\Users\USER\shiks\ode\ode_solver.py", line 95, in _run
    y_df = pd.DataFrame(data=sol["y"].T.copy(), columns=variables)
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices
"""
lisphilar commented 2 years ago

It will be required to check the output ok f sol["y"].T so that we can find the difference of the functions.

lisphilar commented 2 years ago

Source codes of fodeint indicates it returns y directly (instead of named tuple of t and y). We can try

shik-design commented 2 years ago

Thank you for your insightful responses. I replace sol["y"].T with sol. The error did not come up but it take Jupyer Notebook very long to process without end. Any suggestion for me? I use MATLAB often for my analysis but I need to switch to Python

lisphilar commented 2 years ago

Good :)

Actually, Python has a problem in speed, especially for loops. fodeint runs many loops and this may be the cause.

https://github.com/mattja/fodeint/blob/master/fodeint/fodeint.py#L129-L133

Why fodeint is required for your analysis?

shik-design commented 2 years ago

I do appreciate your efforts to help in the furtherance of my research.

My target is to investigate the Memory Effect and Non-locality posed by biological or epidemiological data.

To that extent, Fractional Differential Equations (FDE) are used to Capture the above mentioned phenomena.

That is the reason why fodeint is required. I have also want to use Machine Learning algorithm for parameter estimations, this being the reason for choosing Python.

Well, I wish you consider implementing fodeint for me or any other scheme for FDE so that CovsirPhy will attract more application/usage.

lisphilar commented 2 years ago

Thank you for your proposal, but I think CovsirPhy does not have reasons to use fodeint at this time because the main purpose is to provide effective and speedy tools for data analysis with ODE models.

As an expert in the area,

I found a paper which concluded an ODE model was better than a FDE model for fitting their own data (I'm sorry if I misread the conclusion because I do not have time to read the contents today.) https://www.researchgate.net/publication/338889256_Integer_Versus_Fractional_Order_SEIR_Deterministic_and_Stochastic_Models_of_Measles

lisphilar commented 1 year ago

@shik-design Any updates regarding this issue?