scipy output treats "I" state variable as imaginary

[sdwfrost]

Running the following XML, sir_ode.vf

<?xml version="1.0"?>
    <VectorField Name="sir_ode">
    <Parameter Name="beta" DefaultValue="0.1" Description="Transmission parameter"/>
    <Parameter Name="mu" DefaultValue="0.05" Description="Recovery rate"/>
    <StateVariable Name="S"  Formula="-beta*S*I" DefaultInitialCondition="0.99"/>
    <StateVariable Name="I" Formula="beta*S*I-mu*I" DefaultInitialCondition="0.01"/>
    <StateVariable Name="R" Formula="mu*I" DefaultInitialCondition="0.0"/>
</VectorField>

using

vfgen scipy:func=yes sir_ode.vf

interprets the state variable I as imaginary:

#
# sir_ode.py
#
# Python file for the vector field named: sir_ode
#

"""
This module implements the vector field name "sir_ode" as
the function vectorfield().  The Jacobian of the vector field
is computed by jacobian().  These functions can be used with
the SciPy odeint function.

For example:

    from scipy.integrate import odeint
    import sir_ode

    params = [beta,mu]   # Assume the parameters have been set elsewhere
    t = [i/10.0 for i in range(0, 101)]
    ic = [1.0,0.0,0.0]
    sol = odeint(sir_ode.vectorfield, ic, t, args=(params,), Dfun=sir_ode.jacobian)

This file was generated by the program VFGEN, version: 2.6.0.dev0
Generated on  9-Aug-2018 at 11:27

"""

from math import *
import numpy

#
# The vector field.
#

def vectorfield(y_,t_,p_):
    """
    The vector field function for the vector field "sir_ode"
    Arguments:
        y_ :  vector of the state variables:
                  y_[0] is S
                  y_[1] is I
                  y_[2] is R
        t_ :  time
        p_ :  vector of the parameters
                  p_[0] is beta
                  p_[1] is mu
    """
    S          = y_[0]
    I          = y_[1]
    R          = y_[2]

    beta       = p_[0]
    mu         = p_[1]

    f_ = numpy.zeros((3,))
    f_[0] = std::complex<double>(0.0,-1.0)*S*beta
    f_[1] =  std::complex<double>(0.0,-1.0)*mu+std::complex<double>(0.0,1.0)*S*beta
    f_[2] = std::complex<double>(0.0,1.0)*mu

    return f_

#
#  The Jacobian.
#

def jacobian(y_, t_, p_):
    """
    The Jacobian of the vector field "sir_ode"
    Arguments:
        y_ :  vector of the state variables:
                  y_[0] is S
                  y_[1] is I
                  y_[2] is R
        t_ :  time
        p_ :  vector of the parameters
                  p_[0] is beta
                  p_[1] is mu
    """
    S          = y_[0]
    I          = y_[1]
    R          = y_[2]
    beta       = p_[0]
    mu         = p_[1]

    # Create the Jacobian matrix:
    jac_ = numpy.zeros((3,3))
    jac_[0,0] = std::complex<double>(0.0,-1.0)*beta
    jac_[1,0] = std::complex<double>(0.0,1.0)*beta
    return jac_

This doesn't happen e.g. with R:

#
# sir_ode.R
#
# R vector field functions for: sir_ode
#
# This file was generated by the program VFGEN, version: 2.6.0.dev0
# Generated on  9-Aug-2018 at 11:28
#

#
# sir_ode(t, state, parameters)
#
# The vector field function
#
sir_ode <- function(t, state, parameters) {
    S          <- state[1]
    I          <- state[2]
    R          <- state[3]
    beta       <- parameters[1]
    mu         <- parameters[2]

    vf_ <- vector(len = 3)
    vf_[1] = -I*S*beta;
    vf_[2] = -I*mu+I*S*beta;
    vf_[3] = I*mu;
    return(list(vf_))
}

#
# sir_ode_jac(t, state, parameters)
#
# The jacobian function
#
sir_ode_jac <- function(t, state, parameters) {
    S          <- state[1]
    I          <- state[2]
    R          <- state[3]
    beta       <- parameters[1]
    mu         <- parameters[2]
    jac_ = matrix(nrow = 3, ncol = 3)
    jac_[1,1] = -I*beta
    jac_[1,2] = 0
    jac_[1,3] = 0
    jac_[2,1] = I*beta
    jac_[2,2] = 0
    jac_[2,3] = 0
    jac_[3,1] = 0
    jac_[3,2] = 0
    jac_[3,3] = 0
    return(jac_)
}

[WarrenWeckesser]

Thanks for reporting the issue. I have run into this before, and I lazily avoided the problem by changing the name of the variable to In. See, for example, SIRdelay.vf. That's not a very satisfactory work-around! I'm traveling for a couple weeks, but I'll take another look at this when I get back.

[WarrenWeckesser]

After another very long hiatus from work on vfgen, I'm back to making some updates. For this problem, I'm taking the easy way out and simply not allowing the name 'I'; see https://github.com/WarrenWeckesser/vfgen/commit/10316399e34d2c9ac5b433c23a5fe6cb64a835de.