dynamic.c

GESPI 2.0动态系统模拟工具

#include "copyleft.h"

/*
    GEPASI - a simulator of metabolic pathways and other dynamical systems
    Copyright (C) 1989, 1992, 1993  Pedro Mendes
*/

/*************************************/
/*                                   */
/*        integration shells         */
/*                                   */
/*        Zortech C/C++ 3.0 r4       */
/*          MICROSOFT C 6.00         */
/*          Visual C/C++ 1.0         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/


#ifndef MSDOS
#define COMMON
#endif

#ifdef _WINDOWS
#define COMMON
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#ifdef _WINDOWS
#include <io.h>
#include <fcntl.h>
#endif

#ifdef _ZTC
#define MEM_DEBUG 1
#include "mem.h"
#else
#define mem_malloc malloc
#define mem_free free
#define mem_realloc realloc
#endif

#include "lsoda.h"
#include "globals.h"                                /* global symbols        */
#include "globvar.h"                                /* global variables      */
#include "rates.h"
#include "datab.h"
#include "pmu.h"                                    /* several utilities     */
#include "heapchk.h"

#define D_OK		0
#define D_OUTMEM	1
#define D_NOPEN		2
#define D_NOACC		3
#define D_BADDAT	4
#define D_NSPERR	5

/***************************************************/
/* fint() calculates the rates of each metabolite  */
/* for a certain set of metabolite concentrations  */
/* [like rates() but adapted to work with lsoda()] */
/***************************************************/

static void fint (int neq, double t, double *y, double *ydot )
{
 register int i, j;

 for(i=0;i<indmet;i++)
 {
  ydot[i+1] = (double) 0;
  for( j=0; j<nsteps; j++ )
  {
   if( stoiu[i][j] != (float) 0 )
      ydot[i+1] += stoiu[i][j] * rateq[j](&y[1], j);
  }
 }
 for( i=indmet; i<nmetab; i++)
 {
  ydot[i+1] = (double) 0;
  for( j=0; j<indmet; j++)
   ydot[i+1] -= ld[i][j] * ydot[j+1];
 }
}

/************************************/
/* integrate the system of ODEs and */
/* put data points in a .dyn file   */
/************************************/

int dynamic( void )
{
 FILE *ch1;
 register int i, j;
 double rwork1, rwork5, rwork6, rwork7;
 double abstol[MAX_MET+1], rtol[MAX_MET+1], t, incr, oldtime;
 int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
 int itol, itask, istate, iopt, jt;
 double *xint;
 char delim;

/* put xint in the heap   */
 if (
     ( xint = (double *) mem_malloc( (size_t) (MAX_MET + 1) * sizeof( double ) ) )
     == NULL
    )
 {
  printf("\n  * dynamic.c - fatal error: out of memory\n");
  return D_OUTMEM;
 }

 iwork1= iwork2= iwork5=                          /* initialize lsoda       */
 iwork7= iwork8= iwork9= 0;                       /* control variables      */
 rwork1= rwork5= rwork6= rwork7= (double) 0;
 iwork6 = 1000;                                   /* allow lsoda 1000 steps */
 itol = 1;                                        /* scalar error control   */
 rtol[0] = abstol[0] = (double) 0;
 for( i=1; i<=nmetab; i++)
 {
  rtol[i] = options.reltol;                       /* relative tolerance    */
  abstol[i] = options.abstol;                     /* absolute tolerance    */
 }
 iwork8 = options.adams;
 iwork9 = options.bdf;
 itask = 1;
 istate = 1;
 iopt = 1;
 jt = 2;
 t = actime = (double) 0;
 incr = endtime/options.pfo;                        /* size of time increments*/

 if( options.datsep == 0 ) delim = ' ';
 else if( options.datsep == 1 ) delim = ',';
 else delim = '\t';

 ch1 = fopen( dyname, "at" );                     /* open output for .dyn   */
 if ( ch1 == NULL )                               /* fopen failed ?         */
 {
  printf("\n  * dynamic() - fatal error: couldn't open file %s\n",
         dyname );
  mem_free( xint );
  return D_NOPEN;
 }
 if ( options.dattit == 0  )
 {                                                /* output header          */
  if( options.quotes ) fprintf(ch1,"\"" );
  fprintf(ch1,"%-*s", options.datwidth, "time" );
  if( options.quotes ) fprintf(ch1,"\"" );
  for(j=0;j<nmetab;j++)
  {
   fprintf(ch1,"%c", delim );
   if( options.quotes ) fprintf(ch1,"\"" );
   fprintf(ch1,"[%.*s]%*s", options.datwidth-2, metname[j], options.datwidth-strlen(metname[j])-2, "" );
   if( options.quotes ) fprintf(ch1,"\"" );
  }
  for(j=0;j<nsteps;j++)
  {
   fprintf(ch1,"%c", delim );
   if( options.quotes ) fprintf(ch1,"\"" );
   fprintf(ch1,"J(%.*s)%*s", options.datwidth-3, stepname[j], options.datwidth-strlen(stepname[j])-3, "" );
   if( options.quotes ) fprintf(ch1,"\"" );
  }
  fprintf(ch1,"\n");
 }

 for( i=1; i<=totmet; i++ ) xint[i] = x[i-1];     /* setup initial values   */

 get_cursor(&cx, &cy);                            /* store cursor position  */

 xsetf( options.debug );                          /* switch lsoda printing  */

 for(i=0;t<=endtime;i++)                          /* integration main loop  */
 {
  fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, t);                        /* output current state   */
  for(j=0;j<nmetab;j++)
    fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, x[j]);
  for(j=0;j<nsteps;j++)
    fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, rateq[j]( x, j) );
  fprintf(ch1,"\n");

  if ( !options.debug )
  {
   set_cursor(cx, cy);                            /* put cursor in position */
#ifndef COMMON
   fprintf(stderr, "%-.2le s" ,t);                /* print current time     */
#endif
  }

  do
  {
   if (!i) istate = 1;
   else if (istate<0) istate = 2;

#ifdef _WINDOWS
   fclose( ch1 );
   _wyield();
   ch1 = fopen( dyname, "at" );                   /* open output for .dyn   */
#endif

                                                  /* call the integrator    */
   lsoda( fint, nmetab, xint, &t, actime, itol, rtol, abstol, itask, &istate,
          iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
          rwork1, rwork5, rwork6, rwork7 );
   switch (istate)
   {
    case  1:
    case  2: oldtime = t;
             if ( options.debug )
              printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
                     nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
                     (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
             actime += incr;                      /* increment time         */
             break;
    case -1:
             actime = tn;
             t = oldtime;
             istate = 3;
             iwork6 = 2000;
             break;
    case -2:
            printf( "\n  * dynamic() - fatal error: too much accuracy requested for" );
            printf( "\n                             this computer's precision\n");
            mem_free( xint );
            freevectors();
            return D_NOACC;
    case -3:
             printf( "\n  * dynamic() - fatal error: illegal data\n");
             if( options.debug )
             {
              printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
                     nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
                     (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
              for(j=0;j<nmetab;j++) printf("\n%-.8le", x[j]);
             }
             mem_free( xint );
             freevectors();
             return D_BADDAT;
    default:
            printf( "\n  * dynamic() - fatal error: istate = %d\n", istate );
            mem_free( xint );
            freevectors();
            return D_NSPERR;
   }
  }
  while ( istate != 2 );

  for( j=0; j<nmetab; j++) x[j] = xint[j+1];      /* update x from xint     */
 }

 /* last call to lsoda to clean up allocated memory - istate 4				*/
 freevectors();

 /* output current state   */
 fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, t);                        /* output current state   */
 for(j=0;j<nmetab;j++)
   fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, x[j]);
 for(j=0;j<nsteps;j++)
   fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, rateq[j]( x, j ) );
 fprintf(ch1,"\n");

 if ( options.dattit == 1  )
 {                                                /* output header          */
  if( options.quotes ) fprintf(ch1,"\"" );
  fprintf(ch1,"%-*s", options.datwidth, "time" );
  if( options.quotes ) fprintf(ch1,"\"" );
  for(j=0;j<nmetab;j++)
  {
   fprintf(ch1,"%c", delim );
   if( options.quotes ) fprintf(ch1,"\"" );
   fprintf(ch1,"[%.*s]%*s", options.datwidth-2, metname[j], options.datwidth-strlen(metname[j])-2, "" );
   if( options.quotes ) fprintf(ch1,"\"" );
  }
  for(j=0;j<nsteps;j++)
  {
   fprintf(ch1,"%c", delim );
   if( options.quotes ) fprintf(ch1,"\"" );
   fprintf(ch1,"J(%.*s)%*s", options.datwidth-3, stepname[j], options.datwidth-strlen(stepname[j])-3, "" );
   if( options.quotes ) fprintf(ch1,"\"" );
  }
  fprintf(ch1,"\n");
 }

 if( options.scan ) fprintf(ch1,"\n");

 fclose(ch1);                                     /* close .dyn file        */

#ifndef COMMON
 if ( !options.debug )
 {
  set_cursor(cx, cy);                             /* clear the time counter */
  fprintf(stderr,"             ");
  set_cursor(cx, cy);
 }
#endif

 mem_free( xint );                                    /* free memory for xint   */
 intst = (double) nst;						 	  /* pass the integer vars to doubles	*/
 nfeval = (double) nfe;
 njeval = (double) nje;
 return D_OK;
}


/************************************/
/* integrate the system of ODEs so  */
/* as to converge to a steady-state */
/************************************/

int int_to_ss( void )
{
 register int i, j;
 double rwork1, rwork5, rwork6, rwork7;
 double abstol[MAX_MET+1], rtol[MAX_MET+1], t, max_rate, oldtime;
 int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
 int itol, itask, istate, iopt, jt;
 double  *xint;

                                                  /* put xint in the heap   */
 if (
     ( xint = (double  *) mem_malloc( (size_t) (MAX_MET + 1) * sizeof( double ) ) )
     == NULL
    )
 {
  printf("\n  * int_to_ss() - fatal error: out of memory\n");
  return D_OUTMEM;
 }

 iwork1= iwork2= iwork5=                          /* initialize lsoda       */
 iwork7= iwork8= iwork9= 0;                       /* control variables      */
 rwork1= rwork5= rwork6= rwork7= (double) 0;
 iwork6= 1000;                                    /* allow for 1000 steps   */
 itol = 1;                                        /* scalar error control   */
 rtol[0] = abstol[0] = (double) 0;
 for( i=1; i<=nmetab; i++)
 {
  rtol[i] = options.reltol;                                 /* relative tolerance    */
  abstol[i] = options.abstol;                               /* absolute tolerance    */
 }
 iwork8 = options.adams;
 iwork9 = options.bdf;
 itask = 1;
 istate = 1;
 iopt = 1;
 jt = 2;
 t = (double) 0;
 if ( options.dyn == 0 )
 {
  actime = (double) 1e-4;
  for( i=1; i<=totmet; i++ ) xint[i] = x0[i-1];         /* setup initial values   */
 }
 else
  for( i=1; i<=totmet; i++ ) xint[i] = x[i-1];          /* setup initial values   */

 xsetf( options.debug );                               /* switch lsoda printing  */
 for( i = 0, max_rate = INFINITY; (max_rate>options.hrcz) && (t<=endtime); i++ )
 {

  do
  {
   if (!i) istate = 1;
   else if (istate<0) istate = 2;

#ifdef _WINDOWS
   _wyield();
#endif

   lsoda( fint, nmetab, xint, &t, actime, itol, rtol, abstol, itask, &istate,
          iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
          rwork1, rwork5, rwork6, rwork7 );
   switch (istate)
   {
    case  1:
    case  2: oldtime = t;
             if ( options.debug )
              printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
                     nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
                     (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
             actime *= 10;
             break;
    case -1:
             actime = tn;
             t = oldtime;
             istate = 3;
             iwork6 = 2000;
             break;
    case -2:
             printf( "\n  * int_to_ss() - fatal error: too much accuracy requested for" );
             printf( "\n                               this computer's precision\n");
             mem_free( xint );
             freevectors();
             return D_NOACC;
    case -3:
             printf( "\n  * int_to_ss() - fatal error: illegal data\n");
             if( options.debug )
             {
              printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
                     nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
                     (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
              for(j=0;j<nmetab;j++) printf("\nx[%s] = %-.8le", metname[j], x[j]);
             }
             mem_free( xint );
             freevectors();
             return D_BADDAT;
    default:
             printf( "\n  * int_to_ss() - fatal error: istate = %d\n", istate );
             mem_free( xint );
             freevectors();
             return D_NSPERR;
   }
  }
  while( istate != 2 );

  for( j=0; j<nmetab; j++) x[j] = xint[j+1];      /* update x from xint     */
  calcrates( x );
  for( j=0, max_rate = (double) 0 ; j<nmetab; j++)
   if( fabs( rate[j] ) > max_rate ) max_rate = fabs( rate[j] );
 }

 /* free the vectors allocated by lsoda				*/
 freevectors();

 mem_free( xint );                                    /* free memory            */
 intst = (double) nst;							  /* pass the integer vars to doubles	*/
 nfeval = (double) nfe;
 njeval = (double) nje;
 return D_OK;
}