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.