warrenweckesser/vfgen

Source code generator for differential equation solvers.

VFGEN

A Vector Field File Generator for differential equation solvers and other computational tools. Given a specification of an ordinary or delay differential equation, VFGEN can generate source code that defines the equations for a wide assortment of solvers and tools.

• Web page: https://warrenweckesser.github.io/vfgen
• License: GPL-2
• Author: Warren Weckesser, https://warrenweckesser.github.io

This file is in the top directory for the source code distribution of the program VFGEN. This README gives the steps necessary to build VFGEN from the source code. Before building from source, check the web page https://warrenweckesser.github.io/vfgen; a compiled binary may be available for your system.

Before you can build VFGEN, you must have installed these libraries:

• GiNaC (http://www.ginac.de)
• Mini-XML 3.0 or greater(http://www.minixml.org)

If you are using Linux, binary packages (.deb and .rpm) are available for these libraries. Check your standard repositories. Be sure you install the "development" versions of the packages. In Debian-based systems, these usually end with -dev. Note that GiNaC depends on the arbitrary precision numerical library CLN, so you will also have to install CLN.

Also be sure you have pkg-config (www.freedesktop.org/software/pkgconfig) installed. Most Linux systems have pkg-config available as a binary package.

There are three methods that you can use to build the program: CMake, configure (i.e. the autotools script), or plain old Make. These are discussed below.

CMake

CMake (http://www.cmake.org/) is a platform-independent build tool. To build VFGEN with CMake, you must have CMake installed, and you must have the GiNaC (http://www.ginac.de/) and Mini-XML (http://www.easysw.com/~mike/mxml/) libraries installed. (If you are using a debian-based Linux distribution, packages for these libraries are available.) If these libraries are installed in either /usr or /usr/local, the following commands (run in the top directory) will build VFGEN:

$ cd cmake_build
$ cmake ../src
$ make

If the libraries are installed in some other location, define the environment variables CMAKE_INCLUDE_PATH and CMAKE_LIBRARY_PATH before running cmake.

By default, the command

$ make install

will install the executable file in /usr/local/bin. To change the installation directory, set the environment variable CMAKE_INSTALL_PREFIX before running cmake. Then make install will install the executable file in $(CMAKE_INSTALL_PREFIX)/bin.

configure (i.e. autotools)

VFGEN also comes with a "configure" script. The simplest way to use this script is the following sequence of commands in the top directory:

$ ./configure
$ make
$ make install

This will install the vfgen executable in /usr/local/bin. You may change the installation directory with the --prefix option to the ./configure command, e.g.

$ ./configure --prefix=/opt
$ make
$ make install

This will put the vfgen executable in /opt/bin.

The configure script has several other options. They may be listed with the command

$ ./configure --help

Plain old Makefile

If you are using Linux (or some other Unix-like system), you have the GNU C++ compiler installed (g++), you have the program pkg-config installed, and you have installed GiNaC and Mini-XML, you can use the file Makefile.vfgen to build the program:

$ cd src
$ make -f Makefile.vfgen

GiNaC version 1.3.4, and Mini-XML version 2.3 (and later versions) provide the appropriate files to work with pkg-config; I'm not sure if older versions do.

The file Makefile.vfgen is a very simple Makefile; for example, there is no install target. Besides the default vfgen target, the only other target of interest is clean; make clean will remove all the .o files.

Building from version control source

The code under version control does not include the configure script. To create the configure script:

$ aclocal
$ automake --add-missing
$ autoconf

Some additional options may be necessary. On a Mac running OS X 10.9.5, with pkg-config installed in /usr/local, aclocal might have to run as

$ aclocal -I /usr/local/share/aclocal

vfgen's Issues

Add support for the GSL odeiv2 API to the gsl and pygsl commands.

BUG: Pi is not defined in the LSODA function code even when it is used in the formulas.

For example, here is foo.vf:

<?xml version="1.0" ?>
<VectorField
    Name="foo">
<StateVariable
    Name="x"
    Formula="-Pi*x"
    DefaultInitialCondition="10.0" />
</VectorField>

In the code generated by vfgen lsoda:demo=yes foo.vf, the Pi is used in both foo_rhs and foo_jac, but it is not defined:

c
c foo_rhs.f
c
c Vector field functions for the vector field 'foo'
c These functions are to be used with the Fortran ODE solver LSODA.
c
c This file was generated by the program VFGEN, version: 2.5.0
c Generated on 23-Jun-2017 at 22:03
c
      subroutine foo_rhs(n_,t_,y_,f_)
      implicit none
      integer n_
      double precision t_, y_, f_
      dimension y_(1), f_(1)
      double precision x

c     --- State variables ---
      x          = y_(1)
c     --- The vector field ---
      f_(1) = -Pi*x

      return
      end

c
c Jacobian of the vector field 'foo'
c
c The subroutine assumes that the full Jacobian is to be computed.
c ml_ and mu_ are ignored, and nrowpd_ is assumed to be n_.
c
      subroutine foo_jac(n_,t_,y_,ml_,mu_,jac_,nrowpd_)
      implicit none
      integer n_, ml_, mu_, nrowpd_
      double precision t_, y_, jac_
      dimension y_(1), jac_(nrowpd_,1)
      double precision x

c     --- State variables ---
      x          = y_(1)
c     --- Jacobian ---
      jac_(1, 1) = -Pi

      return
      end

Pi is defined in foo_demo.f, but it is not used there (see #1):

c
c foo_demo.f
c
c Fortran 77 program that uses LSODA to solve the differential equations
c defined in the vector field 'foo'
c
c This file was generated by the program VFGEN, version: 2.5.0
c Generated on 23-Jun-2017 at 22:03
c
      program foo

      implicit none

      external foo_rhs
      external foo_jac

      double precision atol_, rtol_, y_, t_, tout_, tfinal_, rwork_
      integer iwork_
      dimension y_(1), rwork_(36), iwork_(21)
      integer neq_, i_, j_, nsteps_
      integer itol_, iopt_, itask_, istate_, jt_, lrw_, liw_
      double precision x
      double precision Pi
      Pi = 3.1415926535897932385D0
c     --- t range ---
      t_ = 0.0D0
      tfinal_  = 10.0D0
      nsteps_ = 100
c     --- Initial conditions ---
      x = 1.0000000000000000D+01
      y_(1) = x
c     --- Solver tolerances ---
      rtol_ = 1.0D-6
      atol_ = 1.0D-8
      itol_ = 1
c     --- Other LSODA parameters ---
      neq_ = 1
      itask_ = 1
      istate_ = 1
      iopt_ = 0
      lrw_ = 36
      liw_ = 21
c
c     LSODE and LSODA take the same arguments, so either may
c     be used in the loop below.  jt_ must be set as follows:
c     jt_ =  1 for LSODA
c     jt_ = 10 for LSODE, non-stiff (Adams) method
c     jt_ = 21 for LSODE, stiff (BDF) method
c     See the documentation in the Fortran file for more details.
      jt_ = 1
c     --- Print the first point ---
      write (6,49) t_, (y_(j_), j_ = 1,1)
c     --- Call DLSODA in a loop to compute the solution ---
      do 40 i_ = 1,nsteps_
          tout_ = (i_*tfinal_)/nsteps_
          call DLSODA(foo_rhs, neq_, y_, t_, tout_,
     &           itol_, rtol_, atol_, itask_, istate_, iopt_,
     &           rwork_, lrw_, iwork_, liw_,
     &           foo_jac, jt_)
          if (istate_ .lt. 0) goto 80
40        write (6,49) t_, (y_(j_), j_ = 1,1)
49        format(1X,F10.6,1E18.10)
      stop
80    write (6,89) istate_
89    format(1X,"Error: istate=",I3)
      stop
      end

VFGEN in not getting installed

Dear,
I am trying to install VFGEN. But the download page in https://warrenweckesser.github.io/vfgen/ shows no file is there. If I go to github link and download VFGEN-MASTER then also it is not installable. Please suggest what to do.

./configure

configure:

VFGEN Configure Script

Basic environment check...

configure: error: cannot find install-sh or install.sh in "." "./.." "./../.."

=====================================

BUG: A declaration of "Pi" oftens shows up in code where it is not used.

"Pi" shouldn't be declared if it is not used in the generated code.

VfGEN fails to build against libmxml 3.1

I am getting this error message on Ubuntu Focal Fossa

/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp: In member function ‘int VectorField::ReadXML(std::string)’:
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:173:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  173 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:174:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  174 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:226:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  226 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:227:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  227 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:283:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  283 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:284:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  284 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:335:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  335 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:336:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  336 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:392:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  392 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:393:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  393 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:481:33: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  481 |         for (int i = 0; i < node->value.element.num_attrs; ++i) {
      |                                 ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
/home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:482:31: error: invalid use of incomplete type ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
  482 |             string attr = node->value.element.attrs[i].name;
      |                               ^~
In file included from /home/fccoelho/Downloads/vfgen/src/vf_xml.cpp:28:
/usr/include/mxml.h:97:16: note: forward declaration of ‘mxml_node_t’ {aka ‘struct _mxml_node_s’}
   97 | typedef struct _mxml_node_s mxml_node_t; /**** An XML node. ****/
      |                ^~~~~~~~~~~~
make[2]: *** [CMakeFiles/vfgen.dir/build.make:479: CMakeFiles/vfgen.dir/vf_xml.cpp.o] Erro 1
make[1]: *** [CMakeFiles/Makefile2:76: CMakeFiles/vfgen.dir/all] Erro 2
make: *** [Makefile:130: all] Erro 2

scipy output treats "I" state variable as imaginary

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_)
}

Generated code for dde_solver

First of all, thank you for this nice piece of software of vfgen !
However I have some troubles with the files generated by the dde_solver
command.
The *.f90 defines the *_ddes and *_history subroutines that have the
following declarations

SUBROUTINE modele1_nonneutre_ddes(t,x_,Zlags_,vf_)
! Arguments
DOUBLE PRECISION :: t
DOUBLE PRECISION, DIMENSION(NEQN) :: x_, vf_
DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Zlags_
...
SUBROUTINE modele1_nonneutre_history(t,x_)
DOUBLE PRECISION :: t
DOUBLE PRECISION, DIMENSION(NEQN) :: x_
INTENT(IN) :: t
INTENT(OUT) :: x_
...

gfortran has problem with it as DDE_SOLVER does not match any of the
possible DKL_ subroutine of the DDE_SOLVER interface (in
dde_solver_m.f90). So that it raises the following error:

modele1_nonneutre_demo.f90:64:6:
 SOL = DDE_SOLVER(NVAR, modele1_nonneutre_ddes, LAGS,modele1_nonneutre_history,TSPAN,OPTIONS=OPTS)
      1
Error: There is no specific function for the generic ‘dde_solver’ at (1)

This is confirmed by directly calling the expected DKL_* subroutine in the demo :
modele1_nonneutre_demo.f90:64:18:

 SOL = DKL_2(NVAR, modele1_nonneutre_ddes, LAGS,modele1_nonneutre_history,TSPAN,OPTIONS=OPTS)
                  1
Error: Interface mismatch in dummy procedure ‘ddes’ at (1): Shape mismatch in argument 'y'

Replacing the DIMENSION(NEQN) and DIMENSION(NEQN, NLAGS) by
DIMENSION(:) and DIMENSION(:,:), respectively, allows for the
compilation, but the program then segfaults.
Do you have any clues to solve that problem ?

suggestion to vectorize Functions

Hello,
I suggest to create a vector valued function that calculates all user defined Functions in one go, to not repeat unnecessary calculations (when all of them are needed), and most importantly: for convenience.

Specifically, I would like that in the gsl_odeiv(2) module. To use the Function names productively, one could define an enum:

enum Function_names {Name_1, Name_2}; /* etc */

and later on, in the unified output function:

/*  User defined functions */
int SYSTEM_out(double t, const double y_[], void *params, double *Function)
{
 /* assignments 
  * and 
  * expression calculations 
  */
 Function[Name_1] = Value_1;
 Function[Name_2] = Value_2;
 /* etc. */
 return GSL_SUCCESS;
}

This function could also exist in addition to the scalar functions. That would be neat, in most languages, I guess.
Thank you for all the hard work!

vfgen fails to build as it's looking for vf_help_text.cpp

The package fails to build as vf_help_text.o is listed in OBJ_FILES but there isn't a corresponding .cpp file in vfgen/src. Removing vf_help_text.o from OBJ_FILES allows compilation to continue, however, I then get the following error:

vfgen.o: In function `main':
vfgen.cpp:(.text.startup+0x3783): undefined reference to `VectorField::PrintBoostOdeint(std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > > >)'
collect2: error: ld returned 1 exit status
Makefile.vfgen:53: recipe for target 'vfgen' failed
make: *** [vfgen] Error 1

This is on Ubuntu 18.04.