home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / GEPOUT.C < prev    next >
C/C++ Source or Header  |  1993-08-18  |  16KB  |  596 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*           output module           */
  11. /*                                   */
  12. /*     (dynamic data is output by    */
  13. /*      the dynamic() function)      */
  14. /*                                   */
  15. /*        Zortech C/C++ 3.0 r4       */
  16. /*          MICROSOFT C 6.00         */
  17. /*          Visual C/C++ 1.0         */
  18. /*           QuickC/WIN 1.0          */
  19. /*             ULTRIX cc             */
  20. /*              GNU gcc              */
  21. /*                                   */
  22. /*   (include here compilers that    */
  23. /*   compiled GEPASI successfully)   */
  24. /*                                   */
  25. /*************************************/
  26.  
  27.  
  28. #include <stdio.h>
  29. /*#include <conio.h>*/
  30. #include <string.h>
  31. #include <math.h>
  32. #include <time.h>
  33. #include "globals.h"
  34. #include "globvar.h"
  35. #include "datab.h"
  36. #include "newton.h"
  37. #include "metconan.h"
  38. #include "lsoda.h"
  39. #include "pmu.h"
  40. #include "heapchk.h"
  41.  
  42. void report_header( FILE *ch );
  43. void out_rstoi( FILE *o );
  44. void out_rinv( FILE *o );
  45.  
  46. /* print an error message             */
  47. void out_error( int n )
  48. {
  49.  printf( "\n%s\n", errormsg[n] );
  50. }
  51.  
  52. int dat_fileovr( void )
  53. {
  54.  FILE *ch1;
  55.  
  56.  /*unlink*/
  57.  ch1 = fopen( options.datname, "wt" );             /* open output for .dat   */
  58.  if ( ch1 == NULL )                                /* fopen failed ?         */
  59.  {
  60.   if( options.debug )
  61.    printf("\n  * dat_out() - fatal error: couldn't open file %s\n", options.datname );
  62.   return -1;
  63.  }
  64.  fclose( ch1 );
  65.  return 0;
  66. }
  67.  
  68. int dat_nl( void )
  69. {
  70.  FILE *ch1;
  71.  
  72.  ch1 = fopen( options.datname, "at" );             /* open output for .dat   */
  73.  if ( ch1 == NULL )                                /* fopen failed ?         */
  74.  {
  75.   if( options.debug )
  76.    printf("\n  * dat_out() - fatal error: couldn't open file %s\n", options.datname );
  77.   return -1;
  78.  }
  79.  
  80.  fprintf(ch1,"\n");
  81.  
  82.  fclose( ch1 );
  83.  return 0;
  84. }
  85.  
  86. int dat_titles( void )
  87. {
  88.  FILE *ch1;
  89.  char *ptr;
  90.  char delim;
  91.  int i;
  92.  
  93.  switch( options.datsep )
  94.  {
  95.   case 0: delim = ' '; break;
  96.   case 1: delim = ','; break;
  97.   case 2: delim = '\t'; break;
  98.   default: delim = ' ';
  99.  }
  100.  ch1 = fopen( options.datname, "at" );             /* open output for .dat   */
  101.  if ( ch1 == NULL )                                /* fopen failed ?         */
  102.  {
  103.   if( options.debug )
  104.    printf("\n  * dat_out() - fatal error: couldn't open file %s\n", options.datname );
  105.   return -1;
  106.  }
  107.  ptr = outtit;
  108.  for(i=0;i<totsel;i++)
  109.  {
  110.   if( i ) fprintf(ch1,"%c", delim );
  111.   if( options.quotes ) fprintf(ch1,"\"" );
  112.   fprintf(ch1,"%.*s%*s", options.datwidth, ptr, options.datwidth-strlen(ptr), "" );
  113.   if( options.quotes ) fprintf(ch1,"\"" );
  114.   do{ }while( *ptr++ );
  115.  }
  116.  fprintf(ch1,"\n");
  117.  
  118.  fclose( ch1 );
  119.  return 0;
  120. }
  121.  
  122. int dat_values( void )
  123. {
  124.  FILE *ch1;
  125.  int i;
  126.  char delim;
  127.  
  128.  switch( options.datsep )
  129.  {
  130.   case 0: delim = ' '; break;
  131.   case 1: delim = ','; break;
  132.   case 2: delim = '\t'; break;
  133.   default: delim = ' ';
  134.  }
  135.  ch1 = fopen( options.datname, "at" );             /* open output for .dat   */
  136.  if ( ch1 == NULL )                                /* fopen failed ?         */
  137.  {
  138.   if( options.debug )
  139.    printf("\n  * dat_out() - fatal error: couldn't open file %s\n", options.datname );
  140.   return -1;
  141.  }
  142.  
  143.  for(i=0;i<totsel;i++)
  144.  {                                                 /* output values             */
  145.   if( i ) fprintf(ch1,"%c", delim );
  146.   fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, *(poutelem[i]));
  147.  }
  148.  fprintf(ch1,"\n");
  149.  
  150.  fclose( ch1 );
  151.  return 0;
  152. }
  153.  
  154. int dat_error( void )
  155. {
  156.  FILE *ch1;
  157. /* char delim;
  158.  
  159.  switch( options.datsep )
  160.  {
  161.   case 0: delim = ' '; break;
  162.   case 1: delim = ','; break;
  163.   case 2: delim = '\t'; break;
  164.   default: delim = ' ';
  165.  }*/
  166.  ch1 = fopen( options.datname, "at" );             /* open output for .dat   */
  167.  if ( ch1 == NULL )                                /* fopen failed ?         */
  168.  {
  169.   if( options.debug )
  170.    printf("\n  * dat_error() - fatal error: couldn't open file %s\n", options.datname );
  171.   return -1;
  172.  }
  173.  
  174. /* for(i=0;i<totsel;i++)
  175.  {                      */                          /* output values                */
  176. /*  if( i ) fprintf(ch1,"%c", delim );*/
  177.   /*fprintf(ch1,"% -*.*le", options.datwidth, options.datwidth-8, (double) 0 );*/
  178. /*  for(j=0;j<options.datwidth;j++) fprintf(ch1," ");
  179.  }*/
  180.  fprintf(ch1,"\n");
  181.  
  182.  fclose( ch1 );
  183.  return 0;
  184. }
  185.  
  186. void report_init( int p )
  187. {
  188.  FILE *ch1;
  189.  char mode[4];
  190.  int i,j,k;
  191.  char *ptr, auxstr[10];
  192.  
  193.  if( options.append ) strcpy( mode, "at" );
  194.  else strcpy( mode, "wt" );
  195.  strcpy( repfile, filename[p] );
  196.  fixext( repfile, '.', ".txt" );
  197.  ch1 = fopen( repfile, mode );
  198.  if( options.append ) fprintf( ch1, "\n" );            /* separate from previous    */
  199.  report_header( ch1 );
  200.  
  201.  fprintf( ch1, "\n%s\n", topname );                 /* pathway title            */
  202.                                                     /* output struct filen         */
  203.  if ( options.structan )
  204.  {
  205.   fprintf(ch1,"\nREDUCED STOICHEIOMETRY MATRIX\n");
  206.   out_rstoi( ch1 );
  207.   fprintf(ch1,"\nINVERSE OF THE REDUCTION MATRIX\n");
  208.   out_rinv( ch1 );
  209.  }
  210.                                                     /* info about moieties */
  211.  if ( depmet > 0 )
  212.  {
  213.   fprintf(ch1,"\nCONSERVATION RELATIONSHIPS\n");
  214.   for( i=indmet; i<nmetab; i++ )
  215.   {
  216.    for(j=0, k=0;j<totmet;j++)
  217.     if ( intmet[j] )
  218.      if (ld[i][j] != (float) 0)
  219.      {
  220.       sprintf(auxstr, "%s ", k ? "+" : " " );
  221.       if (ld[i][j] > (float) 1)
  222.        sprintf(auxstr, "%s %3.1f*", k ? "+" : " ", ld[i][j] );
  223.       else if (ld[i][j] < (float) -1)
  224.             sprintf(auxstr, "- %3.1f*", fabs(ld[i][j]));
  225.            else if (ld[i][j] == (float) -1)
  226.                  sprintf(auxstr, "- ");
  227.       fprintf(ch1,"%s[%s] ", auxstr, metname[j] );
  228.       k++;
  229.      }
  230.    fprintf(ch1,"= %.4lg %s\n", moiety[i], options.concu);
  231.   }
  232.  }
  233.  else
  234.   fprintf(ch1,"\nNO CONSERVATION RELATIONSHIPS\n");
  235.  
  236.  fprintf(ch1,"\nKINETIC PARAMETERS\n");             /* output kinetic p.        */
  237.  for( i=0; i<nsteps; i++ )
  238.  {
  239.   fprintf(ch1,"%s (%s)\n",stepname[i], ktype[kinetype[i]].descr);
  240.   ptr = ktype[kinetype[i]].constnam;
  241.   for(j=0;j<(int)ktype[kinetype[i]].nconst; j++)
  242.   {
  243.    fprintf(ch1," %s = % .4le\n", ptr, params[i][j] );
  244.    do{ } while( *(ptr++) );
  245.   }
  246.  }
  247.  
  248.  fclose(ch1);
  249. }
  250.  
  251. void report_dyn( void )
  252. {
  253.  FILE *ch1;
  254.  int i;
  255.  
  256.  ch1 = fopen( repfile, "at" );
  257.  
  258.  /* output concentrations    */
  259.  fprintf( ch1, "\nRESULTS OF INTEGRATION (after %.2le %s)\n", endtime, options.timeu );
  260.  for( i=0; i<totmet; i++ )
  261.   fprintf(ch1,"[%s] initial = % .6le %s, final = % .6le %s%s\n",
  262.               metname[i], x0[i], options.concu, x[i],
  263.               options.concu, x[i] >= 0 ? "" : " BOGUS!");
  264.  
  265.  /* output fluxes            */
  266.  for( i=0; i<nsteps; i++ )
  267.   fprintf( ch1, "J(%s) = % .6le %s/%s\n", stepname[i], flux[i], options.concu, options.timeu );
  268.  
  269.  fclose(ch1);
  270. }
  271.  
  272. void report_ss( int diverge, int ccfs )
  273. {
  274.  FILE *ch1;
  275.  int i,j,k;
  276.  double moi;
  277.  char auxstr[10];
  278.  
  279.  ch1 = fopen( repfile, "at" );
  280.  
  281.  if( diverge == N_OK )
  282.  {
  283.   /* check if steady state is an equilibrium */
  284.   if ( equilibrium() )
  285.    fprintf(ch1,"\nEQUILIBRIUM SOLUTION\n");
  286.   else
  287.    fprintf(ch1,"\nSTEADY STATE SOLUTION\n");
  288.  
  289.   /* concentrations    */
  290.   for( i=0; i<totmet; i++ )
  291.   {
  292.    fprintf(ch1,"[%s] = % .6le %s, tt = % .6le %s, rate = % .3le %s/%s",
  293.                metname[i], xss[i], options.concu, tt[i], options.timeu,
  294.                rate[i], options.concu, options.timeu);
  295.    if( xss[i] < 0 ) fprintf(ch1," BOGUS!");
  296.    fprintf(ch1,"\n");
  297.   }
  298.  
  299.   /* mass conservation    */
  300.   if ( depmet > 0 )
  301.   {
  302.    for( i=indmet; i<nmetab; i++ )
  303.    {
  304.     moi = 0;
  305.     for(j=0, k=0;j<totmet;j++)
  306.      if ( intmet[j] )
  307.       if (ld[i][j] != (float) 0)
  308.       {
  309.        moi += (double)ld[i][j] * xss[j];
  310.        sprintf(auxstr, "%s ", k ? "+" : " " );
  311.        if (ld[i][j] > (float) 1)
  312.         sprintf(auxstr, "%s %3.1f*", k ? "+" : " ", ld[i][j] );
  313.        else
  314.        {
  315.         if (ld[i][j] < (float) -1) sprintf(auxstr, "- %3.1f*", fabs(ld[i][j]));
  316.         else if (ld[i][j] == (float) -1) sprintf(auxstr, "- ");
  317.        }
  318.  
  319.        fprintf( ch1,"%s[%s] ", auxstr, metname[j] );
  320.        k++;
  321.       }
  322.     fprintf( ch1, "= %.4lg %s%s\n", moi, options.concu,
  323.                   fabs(moi-moiety[i]) < options.hrcz  ? "" : " BOGUS!" );
  324.    }
  325.   }
  326.  
  327.   /* fluxes    */
  328.   for( i=0; i<nsteps; i++ )
  329.    fprintf( ch1, "J(%s) = % .6le %s/%s\n", stepname[i], flux[i], options.concu, options.timeu );
  330.  
  331.   /* Reder's elasticities    */
  332.   if ( options.nonela )
  333.   {
  334.    fprintf(ch1,"\nABSOLUTE ELASTICITIES\n");
  335.    for( i=0; i<nsteps; i++ )
  336.    {
  337.     for(j=0;j<totmet;j++)
  338.      if ( Dxv[i][j] != (double) 0 )
  339.       fprintf(ch1,"ae(%s,[%s]) =% .4le  ",
  340.                   stepname[i],
  341.                   metname[j],
  342.                   Dxv[i][j]);
  343.     fprintf(ch1,"\n");
  344.    }
  345.   }
  346.  
  347.   /* logarithmic elasticities    */
  348.   if ( options.stdela )
  349.   {
  350.    fprintf(ch1,"\nLOGARITHMIC ELASTICITIES\n");
  351.    for( i=0; i<nsteps; i++ )
  352.    {
  353.     for(j=0;j<totmet;j++)
  354.      if ( Dxv[i][j] != (float) 0 )
  355.       if( fabs( flux[i] ) >= options.hrcz )
  356.        fprintf(ch1, " e(%s,[%s]) =% .4le  ",
  357.                     stepname[i],
  358.                     metname[j],
  359.                     Dxv[i][j] * xss[j] / flux[i] );
  360.       else fprintf(ch1, " e(%s,%s) = infinity  ", stepname[i], metname[j] );
  361.      fprintf(ch1,"\n");
  362.    }
  363.   }
  364.  
  365.   /* Reder's conc. control coeffs.    */
  366.   if ( options.noncc )
  367.   {
  368.    if ( ccfs != MCA_NOCC )
  369.    {
  370.     fprintf(ch1,"\nCONCENTRATION CONTROL COEFFICIENTS (unscaled)\n");
  371.     for( i=0; i<totmet; i++)
  372.      if (intmet[i])
  373.      {
  374.       fprintf(ch1,"\n%s\n", metname[i]);
  375.       for( j=0; j<nsteps; j++)
  376.        fprintf(ch1,"C([%s],J(%s)) = % 6.4le\n", metname[i],
  377.                    stepname[j], Gamma[i][j] );
  378.      }
  379.    }
  380.   }
  381.  
  382.   /* logarithmic conc. control coeffs.    */
  383.   if ( options.stdcc )
  384.   {
  385.    if ( ccfs != MCA_NOCC )
  386.    {
  387.     fprintf(ch1,"\nCONCENTRATION CONTROL COEFFICIENTS (scaled)\n");
  388.     for( i=0; i<totmet; i++)
  389.      if (intmet[i])
  390.      {
  391.       fprintf( ch1, "\n%s\n", metname[i] );
  392.       for( j=0; j<nsteps; j++)
  393.       {
  394.        if ( xss[i] == (double) 0 )
  395.         fprintf( ch1, "C([%s],J(%s)) = not defined\n", metname[i], stepname[j] );
  396.        else
  397.         fprintf(ch1,"C([%s],J(%s)) = % 6.4le\n", metname[i],
  398.                     stepname[j],
  399.                     Gamma[i][j] *
  400.                     flux[j] / xss[i] );
  401.       }
  402.      }
  403.    }
  404.   }
  405.  
  406.   /* Reder's flux control coefficients    */
  407.   if ( options.noncc )
  408.   {
  409.    if ( ccfs != MCA_NOCC )
  410.    {
  411.     fprintf(ch1,"\nFLUX CONTROL COEFFICIENTS (unscaled)\n");
  412.     for( i=0; i<nsteps; i++)
  413.     {
  414.      fprintf(ch1,"\n%s\n", stepname[i]);
  415.       for( j=0; j<nsteps; j++ )
  416.        fprintf(ch1,"C(J(%s),%s) = % 6.4le\n",
  417.                   stepname[i], stepname[j], C[i][j] );
  418.     }
  419.    }
  420.   }
  421.  
  422.   /* logarithmic flux control coeffs.    */
  423.   if ( options.stdcc )
  424.   {
  425.    if ( ccfs != MCA_NOCC )
  426.    {
  427.     fprintf(ch1,"\nFLUX CONTROL COEFFICIENTS (scaled)\n");
  428.     for( i=0; i<nsteps; i++)
  429.     {
  430.      fprintf(ch1,"\n%s\n",stepname[i]);
  431.      for( j=0; j<nsteps; j++)
  432.      {
  433.       if ( fabs( flux[i] ) < options.hrcz )
  434.       {
  435.        if ( fabs( flux[j] ) < options.hrcz )
  436.         fprintf(ch1,"C(J(%s),%s) = not defined\n",
  437.                    stepname[i], stepname[j] );
  438.        else
  439.         fprintf(ch1,"C(J(%s),%s) = infinity\n", stepname[i], stepname[j] );
  440.       }
  441.       else
  442.        fprintf(ch1,"C(J(%s),%s) = % 6.4le\n", stepname[i], stepname[j],
  443.                    C[i][j] *
  444.                    flux[j] / flux[i] );
  445.      }
  446.     }
  447.    }
  448.   }
  449.  }
  450.  else  /* no steady state solution    */
  451.  {
  452.   fprintf(ch1,"\n**** NO STEADY-STATE! ****\n");
  453.   switch( diverge )
  454.   {
  455.    case N_NOCONV: fprintf(ch1,"\nno convergence "); break;
  456.    case N_LMIN  : fprintf(ch1,"\nstuck in local minimum "); break;
  457.    case N_JSING : fprintf(ch1,"\nsingular jacobian ");
  458.   }
  459.   fprintf(ch1,"after 1e+10 %s\n", options.timeu);
  460.  
  461.   for( i=0; i<totmet; i++ )
  462.    fprintf(ch1,"[%s] = % .6le %s, tt = % .6le %s, rate = % .3le %s/%s%s\n",
  463.                metname[i], xss[i], options.concu, tt[i], options.timeu,
  464.                rate[i], options.concu, options.timeu, xss[i] >= 0 ? "" : " BOGUS!");
  465.   for( i=0; i<nsteps; i++ )
  466.    fprintf( ch1, "J(%s) = % .6le\n", stepname[i], flux[i] );
  467.  }
  468.  fclose(ch1);
  469. }
  470.  
  471. void report_tech( void )
  472. {
  473.  FILE *ch1;
  474.  
  475.  ch1 = fopen( repfile, "at" );
  476.  
  477.  fprintf(ch1, "\n\nTECHNICAL INFORMATION");
  478.  fprintf(ch1, "\nlast ODE integration method used: %1d%s order %s",
  479.               nqu, (nqu==1) ? "st" : ( (nqu==2) ? "nd" : ( (nqu==3) ? "rd" : "th" ) ),
  480.               (mused==1) ? "adams" : "gear ");
  481.  fprintf(ch1, "\nnumber of steps taken: %-4d", nst );
  482.  fprintf(ch1, "\nnumber of function evaluations: %-4d", nfe );
  483.  fprintf(ch1, "\nnumber of jacobian evaluations: %-4d", nje );
  484.  fprintf(ch1, "\nlast step size used: %.2le", hu );
  485.  
  486.  fclose(ch1);
  487. }
  488.  
  489. void out_rstoi( FILE *o )
  490. {
  491.  int i,j;
  492.  
  493.  fprintf( o , "\nkey for step (column) numbers:" );
  494.  for( j=0; j<nsteps; j++ )
  495.   fprintf( o, "\n%2d - %s", j, stepname[j]);
  496.  fprintf( o , "\n              " );
  497.  for( j=0; j<nsteps; j++ )
  498.   fprintf( o, "%2d  ", j);
  499.  fprintf( o , "\n" );
  500.  fprintf( o , "           -" );
  501.  for( j=0; j<nsteps; j++ )
  502.   fprintf( o, "----" );
  503.  fprintf( o, "-\n" );
  504.  for( i=0; i<indmet; i++ )
  505.  {
  506.    fprintf( o, "%-10.10s | ", metname[i] );
  507.    for( j=0; j<nsteps; j++ )
  508.     fprintf( o, " % 1.0f ", rstoi[i][j] );
  509.    fprintf( o, "\n");
  510.  }
  511. }
  512.  
  513. void out_rinv( FILE *o )
  514. {
  515.  int i,j;
  516.  
  517.  fprintf( o , "\n              " );
  518.  for( j=0; j<nmetab; j++ )
  519.   fprintf( o, "%2d  ", j);
  520.  fprintf( o, "\n" );
  521.  fprintf( o , "           -" );
  522.  for( j=0; j<nmetab; j++ )
  523.   fprintf( o, "----" );
  524.  fprintf( o, "-\n" );
  525.  for( i=0; i<nmetab; i++ )
  526.  {
  527.    fprintf( o, "%-10.10s | ", metname[i] );
  528.    for( j=0; j<nmetab; j++ )
  529.     fprintf( o, " % 1.0f ", ml[i][j] );
  530.    fprintf( o, "\n");
  531.  }
  532. }
  533.  
  534.  
  535. void report_error( void )
  536. {
  537.  FILE *ch1;
  538.  
  539.  ch1 = fopen( repfile, "at" );
  540.  fprintf( ch1, "\n\nFLOATING-POINT ERROR\n%s: (%d)\n", fpstr, fperr );
  541.  fclose( ch1 );
  542. }
  543.  
  544. /* print the header of the .txt file */
  545. void report_header( FILE *ch )
  546. {
  547.  register int i, hw;
  548.  int pre, suf;
  549.  struct tm *newtime;
  550.  time_t ltime;
  551.  char strtime[26];
  552.  
  553.  /* get seconds elapsed since 1 Jan 1970    */
  554.  time( <ime );
  555.  /* Obtain greenwich mean time: */
  556.  newtime = localtime( <ime );
  557.  strcpy( strtime, asctime( newtime ) );
  558.  strtime[24] ='\0';
  559.  
  560.  if (strlen(gepasi) > strlen(repfile)) hw = strlen(gepasi);
  561.  else hw = strlen(repfile);
  562.  if ( (int) strlen(version_no) > hw) hw = (int) strlen(version_no);
  563.  for(i=0;i<hw+4;i++) fprintf(ch,"*");
  564.  fprintf(ch,"\n");
  565.  fprintf(ch,"* ");
  566.  pre = ( hw - strlen( gepasi ) ) / 2;
  567.  suf = hw - pre - strlen( gepasi );
  568.  for(i=0;i<pre;i++) fprintf(ch," ");
  569.  fprintf(ch, gepasi);
  570.  for(i=0;i<suf;i++) fprintf(ch," ");
  571.  fprintf(ch," *\n");
  572.  fprintf(ch,"* ");
  573.  pre = ( hw - strlen( version_no ) ) / 2;
  574.  suf = hw - pre - strlen( version_no );
  575.  for(i=0;i<pre;i++) fprintf(ch," ");
  576.  fprintf(ch, version_no);
  577.  for(i=0;i<suf;i++) fprintf(ch," ");
  578.  fprintf(ch," *\n");
  579.  fprintf(ch,"* ");
  580.  pre = ( hw - strlen( strtime ) ) / 2;
  581.  suf = hw - pre - strlen( strtime );
  582.  for(i=0;i<pre;i++) fprintf(ch," ");
  583.  fprintf( ch, strtime );
  584.  for(i=0;i<suf;i++) fprintf(ch," ");
  585.  fprintf(ch," *\n");
  586.  fprintf(ch,"* ");
  587.  pre = ( hw - strlen( repfile ) ) / 2;
  588.  suf = hw - pre - strlen( repfile );
  589.  for(i=0;i<pre;i++) fprintf(ch," ");
  590.  fprintf(ch, repfile);
  591.  for(i=0;i<suf;i++) fprintf(ch," ");
  592.  fprintf(ch," *\n");
  593.   for(i=0;i<hw+4;i++) fprintf(ch,"*");
  594.  fprintf(ch,"\n");
  595. }
  596.