home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / DYNAMIC.C < prev    next >
C/C++ Source or Header  |  1993-08-23  |  14KB  |  436 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992, 1993  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*        integration shells         */
  11. /*                                   */
  12. /*        Zortech C/C++ 3.0 r4       */
  13. /*          MICROSOFT C 6.00         */
  14. /*          Visual C/C++ 1.0         */
  15. /*           QuickC/WIN 1.0          */
  16. /*             ULTRIX cc             */
  17. /*              GNU gcc              */
  18. /*                                   */
  19. /*   (include here compilers that    */
  20. /*   compiled GEPASI successfully)   */
  21. /*                                   */
  22. /*************************************/
  23.  
  24.  
  25. #ifndef MSDOS
  26. #define COMMON
  27. #endif
  28.  
  29. #ifdef _WINDOWS
  30. #define COMMON
  31. #endif
  32.  
  33. #include <stdio.h>
  34. #include <stdlib.h>
  35. #include <math.h>
  36. #include <string.h>
  37.  
  38. #ifdef _WINDOWS
  39. #include <io.h>
  40. #include <fcntl.h>
  41. #endif
  42.  
  43. #ifdef _ZTC
  44. #define MEM_DEBUG 1
  45. #include "mem.h"
  46. #else
  47. #define mem_malloc malloc
  48. #define mem_free free
  49. #define mem_realloc realloc
  50. #endif
  51.  
  52. #include "lsoda.h"
  53. #include "globals.h"                                /* global symbols        */
  54. #include "globvar.h"                                /* global variables      */
  55. #include "rates.h"
  56. #include "datab.h"
  57. #include "pmu.h"                                    /* several utilities     */
  58. #include "heapchk.h"
  59.  
  60. #define D_OK        0
  61. #define D_OUTMEM    1
  62. #define D_NOPEN        2
  63. #define D_NOACC        3
  64. #define D_BADDAT    4
  65. #define D_NSPERR    5
  66.  
  67. /***************************************************/
  68. /* fint() calculates the rates of each metabolite  */
  69. /* for a certain set of metabolite concentrations  */
  70. /* [like rates() but adapted to work with lsoda()] */
  71. /***************************************************/
  72.  
  73. static void fint (int neq, double t, double *y, double *ydot )
  74. {
  75.  register int i, j;
  76.  
  77.  for(i=0;i<indmet;i++)
  78.  {
  79.   ydot[i+1] = (double) 0;
  80.   for( j=0; j<nsteps; j++ )
  81.   {
  82.    if( stoiu[i][j] != (float) 0 )
  83.       ydot[i+1] += stoiu[i][j] * rateq[j](&y[1], j);
  84.   }
  85.  }
  86.  for( i=indmet; i<nmetab; i++)
  87.  {
  88.   ydot[i+1] = (double) 0;
  89.   for( j=0; j<indmet; j++)
  90.    ydot[i+1] -= ld[i][j] * ydot[j+1];
  91.  }
  92. }
  93.  
  94. /************************************/
  95. /* integrate the system of ODEs and */
  96. /* put data points in a .dyn file   */
  97. /************************************/
  98.  
  99. int dynamic( void )
  100. {
  101.  FILE *ch1;
  102.  register int i, j;
  103.  double rwork1, rwork5, rwork6, rwork7;
  104.  double abstol[MAX_MET+1], rtol[MAX_MET+1], t, incr, oldtime;
  105.  int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
  106.  int itol, itask, istate, iopt, jt;
  107.  double *xint;
  108.  char delim;
  109.  
  110. /* put xint in the heap   */
  111.  if (
  112.      ( xint = (double *) mem_malloc( (size_t) (MAX_MET + 1) * sizeof( double ) ) )
  113.      == NULL
  114.     )
  115.  {
  116.   printf("\n  * dynamic.c - fatal error: out of memory\n");
  117.   return D_OUTMEM;
  118.  }
  119.  
  120.  iwork1= iwork2= iwork5=                          /* initialize lsoda       */
  121.  iwork7= iwork8= iwork9= 0;                       /* control variables      */
  122.  rwork1= rwork5= rwork6= rwork7= (double) 0;
  123.  iwork6 = 1000;                                   /* allow lsoda 1000 steps */
  124.  itol = 1;                                        /* scalar error control   */
  125.  rtol[0] = abstol[0] = (double) 0;
  126.  for( i=1; i<=nmetab; i++)
  127.  {
  128.   rtol[i] = options.reltol;                       /* relative tolerance    */
  129.   abstol[i] = options.abstol;                     /* absolute tolerance    */
  130.  }
  131.  iwork8 = options.adams;
  132.  iwork9 = options.bdf;
  133.  itask = 1;
  134.  istate = 1;
  135.  iopt = 1;
  136.  jt = 2;
  137.  t = actime = (double) 0;
  138.  incr = endtime/options.pfo;                        /* size of time increments*/
  139.  
  140.  if( options.datsep == 0 ) delim = ' ';
  141.  else if( options.datsep == 1 ) delim = ',';
  142.  else delim = '\t';
  143.  
  144.  ch1 = fopen( dyname, "at" );                     /* open output for .dyn   */
  145.  if ( ch1 == NULL )                               /* fopen failed ?         */
  146.  {
  147.   printf("\n  * dynamic() - fatal error: couldn't open file %s\n",
  148.          dyname );
  149.   mem_free( xint );
  150.   return D_NOPEN;
  151.  }
  152.  if ( options.dattit == 0  )
  153.  {                                                /* output header          */
  154.   if( options.quotes ) fprintf(ch1,"\"" );
  155.   fprintf(ch1,"%-*s", options.datwidth, "time" );
  156.   if( options.quotes ) fprintf(ch1,"\"" );
  157.   for(j=0;j<nmetab;j++)
  158.   {
  159.    fprintf(ch1,"%c", delim );
  160.    if( options.quotes ) fprintf(ch1,"\"" );
  161.    fprintf(ch1,"[%.*s]%*s", options.datwidth-2, metname[j], options.datwidth-strlen(metname[j])-2, "" );
  162.    if( options.quotes ) fprintf(ch1,"\"" );
  163.   }
  164.   for(j=0;j<nsteps;j++)
  165.   {
  166.    fprintf(ch1,"%c", delim );
  167.    if( options.quotes ) fprintf(ch1,"\"" );
  168.    fprintf(ch1,"J(%.*s)%*s", options.datwidth-3, stepname[j], options.datwidth-strlen(stepname[j])-3, "" );
  169.    if( options.quotes ) fprintf(ch1,"\"" );
  170.   }
  171.   fprintf(ch1,"\n");
  172.  }
  173.  
  174.  for( i=1; i<=totmet; i++ ) xint[i] = x[i-1];     /* setup initial values   */
  175.  
  176.  get_cursor(&cx, &cy);                            /* store cursor position  */
  177.  
  178.  xsetf( options.debug );                          /* switch lsoda printing  */
  179.  
  180.  for(i=0;t<=endtime;i++)                          /* integration main loop  */
  181.  {
  182.   fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, t);                        /* output current state   */
  183.   for(j=0;j<nmetab;j++)
  184.     fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, x[j]);
  185.   for(j=0;j<nsteps;j++)
  186.     fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, rateq[j]( x, j) );
  187.   fprintf(ch1,"\n");
  188.  
  189.   if ( !options.debug )
  190.   {
  191.    set_cursor(cx, cy);                            /* put cursor in position */
  192. #ifndef COMMON
  193.    fprintf(stderr, "%-.2le s" ,t);                /* print current time     */
  194. #endif
  195.   }
  196.  
  197.   do
  198.   {
  199.    if (!i) istate = 1;
  200.    else if (istate<0) istate = 2;
  201.  
  202. #ifdef _WINDOWS
  203.    fclose( ch1 );
  204.    _wyield();
  205.    ch1 = fopen( dyname, "at" );                   /* open output for .dyn   */
  206. #endif
  207.  
  208.                                                   /* call the integrator    */
  209.    lsoda( fint, nmetab, xint, &t, actime, itol, rtol, abstol, itask, &istate,
  210.           iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
  211.           rwork1, rwork5, rwork6, rwork7 );
  212.    switch (istate)
  213.    {
  214.     case  1:
  215.     case  2: oldtime = t;
  216.              if ( options.debug )
  217.               printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
  218.                      nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
  219.                      (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
  220.              actime += incr;                      /* increment time         */
  221.              break;
  222.     case -1:
  223.              actime = tn;
  224.              t = oldtime;
  225.              istate = 3;
  226.              iwork6 = 2000;
  227.              break;
  228.     case -2:
  229.             printf( "\n  * dynamic() - fatal error: too much accuracy requested for" );
  230.             printf( "\n                             this computer's precision\n");
  231.             mem_free( xint );
  232.             freevectors();
  233.             return D_NOACC;
  234.     case -3:
  235.              printf( "\n  * dynamic() - fatal error: illegal data\n");
  236.              if( options.debug )
  237.              {
  238.               printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
  239.                      nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
  240.                      (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
  241.               for(j=0;j<nmetab;j++) printf("\n%-.8le", x[j]);
  242.              }
  243.              mem_free( xint );
  244.              freevectors();
  245.              return D_BADDAT;
  246.     default:
  247.             printf( "\n  * dynamic() - fatal error: istate = %d\n", istate );
  248.             mem_free( xint );
  249.             freevectors();
  250.             return D_NSPERR;
  251.    }
  252.   }
  253.   while ( istate != 2 );
  254.  
  255.   for( j=0; j<nmetab; j++) x[j] = xint[j+1];      /* update x from xint     */
  256.  }
  257.  
  258.  /* last call to lsoda to clean up allocated memory - istate 4                */
  259.  freevectors();
  260.  
  261.  /* output current state   */
  262.  fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, t);                        /* output current state   */
  263.  for(j=0;j<nmetab;j++)
  264.    fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, x[j]);
  265.  for(j=0;j<nsteps;j++)
  266.    fprintf(ch1,"%c% -*.*le", delim, options.datwidth, options.datwidth-8, rateq[j]( x, j ) );
  267.  fprintf(ch1,"\n");
  268.  
  269.  if ( options.dattit == 1  )
  270.  {                                                /* output header          */
  271.   if( options.quotes ) fprintf(ch1,"\"" );
  272.   fprintf(ch1,"%-*s", options.datwidth, "time" );
  273.   if( options.quotes ) fprintf(ch1,"\"" );
  274.   for(j=0;j<nmetab;j++)
  275.   {
  276.    fprintf(ch1,"%c", delim );
  277.    if( options.quotes ) fprintf(ch1,"\"" );
  278.    fprintf(ch1,"[%.*s]%*s", options.datwidth-2, metname[j], options.datwidth-strlen(metname[j])-2, "" );
  279.    if( options.quotes ) fprintf(ch1,"\"" );
  280.   }
  281.   for(j=0;j<nsteps;j++)
  282.   {
  283.    fprintf(ch1,"%c", delim );
  284.    if( options.quotes ) fprintf(ch1,"\"" );
  285.    fprintf(ch1,"J(%.*s)%*s", options.datwidth-3, stepname[j], options.datwidth-strlen(stepname[j])-3, "" );
  286.    if( options.quotes ) fprintf(ch1,"\"" );
  287.   }
  288.   fprintf(ch1,"\n");
  289.  }
  290.  
  291.  if( options.scan ) fprintf(ch1,"\n");
  292.  
  293.  fclose(ch1);                                     /* close .dyn file        */
  294.  
  295. #ifndef COMMON
  296.  if ( !options.debug )
  297.  {
  298.   set_cursor(cx, cy);                             /* clear the time counter */
  299.   fprintf(stderr,"             ");
  300.   set_cursor(cx, cy);
  301.  }
  302. #endif
  303.  
  304.  mem_free( xint );                                    /* free memory for xint   */
  305.  intst = (double) nst;                               /* pass the integer vars to doubles    */
  306.  nfeval = (double) nfe;
  307.  njeval = (double) nje;
  308.  return D_OK;
  309. }
  310.  
  311.  
  312. /************************************/
  313. /* integrate the system of ODEs so  */
  314. /* as to converge to a steady-state */
  315. /************************************/
  316.  
  317. int int_to_ss( void )
  318. {
  319.  register int i, j;
  320.  double rwork1, rwork5, rwork6, rwork7;
  321.  double abstol[MAX_MET+1], rtol[MAX_MET+1], t, max_rate, oldtime;
  322.  int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
  323.  int itol, itask, istate, iopt, jt;
  324.  double  *xint;
  325.  
  326.                                                   /* put xint in the heap   */
  327.  if (
  328.      ( xint = (double  *) mem_malloc( (size_t) (MAX_MET + 1) * sizeof( double ) ) )
  329.      == NULL
  330.     )
  331.  {
  332.   printf("\n  * int_to_ss() - fatal error: out of memory\n");
  333.   return D_OUTMEM;
  334.  }
  335.  
  336.  iwork1= iwork2= iwork5=                          /* initialize lsoda       */
  337.  iwork7= iwork8= iwork9= 0;                       /* control variables      */
  338.  rwork1= rwork5= rwork6= rwork7= (double) 0;
  339.  iwork6= 1000;                                    /* allow for 1000 steps   */
  340.  itol = 1;                                        /* scalar error control   */
  341.  rtol[0] = abstol[0] = (double) 0;
  342.  for( i=1; i<=nmetab; i++)
  343.  {
  344.   rtol[i] = options.reltol;                                 /* relative tolerance    */
  345.   abstol[i] = options.abstol;                               /* absolute tolerance    */
  346.  }
  347.  iwork8 = options.adams;
  348.  iwork9 = options.bdf;
  349.  itask = 1;
  350.  istate = 1;
  351.  iopt = 1;
  352.  jt = 2;
  353.  t = (double) 0;
  354.  if ( options.dyn == 0 )
  355.  {
  356.   actime = (double) 1e-4;
  357.   for( i=1; i<=totmet; i++ ) xint[i] = x0[i-1];         /* setup initial values   */
  358.  }
  359.  else
  360.   for( i=1; i<=totmet; i++ ) xint[i] = x[i-1];          /* setup initial values   */
  361.  
  362.  xsetf( options.debug );                               /* switch lsoda printing  */
  363.  for( i = 0, max_rate = INFINITY; (max_rate>options.hrcz) && (t<=endtime); i++ )
  364.  {
  365.  
  366.   do
  367.   {
  368.    if (!i) istate = 1;
  369.    else if (istate<0) istate = 2;
  370.  
  371. #ifdef _WINDOWS
  372.    _wyield();
  373. #endif
  374.  
  375.    lsoda( fint, nmetab, xint, &t, actime, itol, rtol, abstol, itask, &istate,
  376.           iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
  377.           rwork1, rwork5, rwork6, rwork7 );
  378.    switch (istate)
  379.    {
  380.     case  1:
  381.     case  2: oldtime = t;
  382.              if ( options.debug )
  383.               printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
  384.                      nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
  385.                      (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
  386.              actime *= 10;
  387.              break;
  388.     case -1:
  389.              actime = tn;
  390.              t = oldtime;
  391.              istate = 3;
  392.              iwork6 = 2000;
  393.              break;
  394.     case -2:
  395.              printf( "\n  * int_to_ss() - fatal error: too much accuracy requested for" );
  396.              printf( "\n                               this computer's precision\n");
  397.              mem_free( xint );
  398.              freevectors();
  399.              return D_NOACC;
  400.     case -3:
  401.              printf( "\n  * int_to_ss() - fatal error: illegal data\n");
  402.              if( options.debug )
  403.              {
  404.               printf("\n%1d%s %s, hu=%.2le, t=%.2le, tcur=%.2le, s %-4d, f %-4d, j %-4d",
  405.                      nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
  406.                      (mused==1) ? "adams" : "gear ", hu, t, tn, nst, nfe, nje);
  407.               for(j=0;j<nmetab;j++) printf("\nx[%s] = %-.8le", metname[j], x[j]);
  408.              }
  409.              mem_free( xint );
  410.              freevectors();
  411.              return D_BADDAT;
  412.     default:
  413.              printf( "\n  * int_to_ss() - fatal error: istate = %d\n", istate );
  414.              mem_free( xint );
  415.              freevectors();
  416.              return D_NSPERR;
  417.    }
  418.   }
  419.   while( istate != 2 );
  420.  
  421.   for( j=0; j<nmetab; j++) x[j] = xint[j+1];      /* update x from xint     */
  422.   calcrates( x );
  423.   for( j=0, max_rate = (double) 0 ; j<nmetab; j++)
  424.    if( fabs( rate[j] ) > max_rate ) max_rate = fabs( rate[j] );
  425.  }
  426.  
  427.  /* free the vectors allocated by lsoda                */
  428.  freevectors();
  429.  
  430.  mem_free( xint );                                    /* free memory            */
  431.  intst = (double) nst;                              /* pass the integer vars to doubles    */
  432.  nfeval = (double) nfe;
  433.  njeval = (double) nje;
  434.  return D_OK;
  435. }
  436.