home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / GAUSS.C < prev    next >
C/C++ Source or Header  |  1993-08-24  |  11KB  |  352 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. /*   reduction of stoicheiometry by  */
  11. /*   the Gauss method with row and   */
  12. /*           column switches         */
  13. /*                                   */
  14. /*          MICROSOFT C 6.00         */
  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. #include <stdio.h>
  26. #include <math.h>
  27. #include <string.h>
  28. #include "globals.h"
  29. #include "globvar.h"
  30.  
  31. #define ALMOST_ZERO 1.0e-7
  32.  
  33. void rowsw( int r1, int r2, int c )
  34. {
  35.  int j;
  36.  float dummy;
  37.  double dumb;
  38.  int dum;
  39.  
  40.  for(j=0;j<c;j++)
  41.  {
  42.   dum         = stoi[r1][j];                      /* switch rows in stoi and */
  43.   stoi[r1][j] = stoi[r2][j];
  44.   stoi[r2][j] = dum;
  45.   dummy        = rstoi[r1][j];                    /* switch rows in rstoi &  */
  46.   rstoi[r1][j] = rstoi[r2][j];
  47.   rstoi[r2][j] = dummy;
  48.  }
  49.  dumb = x[r1];                                    /* switch rows in x and    */
  50.  x[r1] = x[r2];
  51.  x[r2] = dumb;
  52.  dumb = x0[r1];                                   /* switch rows in x0  and  */
  53.  x0[r1] = x0[r2];
  54.  x0[r2] = dumb;
  55.  j      = ur[r1];                                 /* switch rows in ur too   */
  56.  ur[r1] = ur[r2];
  57.  ur[r2] = j;
  58.  /*strcpy( buff, metname[r1] );
  59.  strcpy( metname[r1], metname[r2] );
  60.  strcpy( metname[r2], buff );*/
  61. }
  62.  
  63. void colsw( int c1, int c2, int r )
  64. {
  65.  int j;
  66.  float dummy;
  67.  int dum;
  68.  
  69.  for( j=0; j<r; j++ )
  70.  {
  71.   dum         = stoi[j][c1];                      /* switch cols in stoi and */
  72.   stoi[j][c1] = stoi[j][c2];
  73.   stoi[j][c2] = dum;
  74.   dummy        = rstoi[j][c1];                    /* switch cols in rstoi &  */
  75.   rstoi[j][c1] = rstoi[j][c2];
  76.   rstoi[j][c2] = dummy;
  77.  }
  78.  j      = uc[c1];                                 /* switch cols in uc too   */
  79.  uc[c1] = uc[c2];
  80.  uc[c2] = j;
  81. /* strcpy( buff, stepname[c1] );
  82.  strcpy( stepname[c1], stepname[c2] );
  83.  strcpy( stepname[c2], buff );*/
  84. }
  85.  
  86. void rowsw_m( int r1, int r2, int c )             /* switch rows in ml       */
  87. {
  88.  int j;
  89.  float dummy;
  90.  
  91.  for(j=0;j<c;j++)
  92.  {
  93.   dummy     = ml[r1][j];
  94.   ml[r1][j] = ml[r2][j];
  95.   ml[r2][j] = dummy;
  96.  }
  97. }
  98.  
  99. void colsw_m( int c1, int c2, int r )             /* switch cols in ml       */
  100. {
  101.  int j;
  102.  float dummy;
  103.  
  104.  for(j=0;j<r;j++)
  105.  {
  106.   dummy     = ml[j][c1];
  107.   ml[j][c1] = ml[j][c2];
  108.   ml[j][c2] = dummy;
  109.  }
  110. }
  111.  
  112. int emptyrow( int rw, int nsteps)
  113. {
  114.  int j, ct;                              /* ct counts entries != 0  */
  115.  
  116.  for( ct=0, j=0; j<nsteps; j++)
  117.   if ( fabs( rstoi[rw][j] ) > ALMOST_ZERO ) ct++;
  118.  return ( ct );
  119. }
  120.  
  121. void invml11( void )
  122. {                                                 /* as ml is lower triangul.*/
  123.  int i,j,k;                                       /* inverting is simply to  */
  124.  float acum;                                      /* forward-substitute with */
  125.                                                   /* the identity matrix     */
  126.  for(i=0;i<nmetab; i++) ml[i][i] = 1 ;
  127.  for( j=0; j<indmet; j++)
  128.   for( k=0; k<indmet; k++)
  129.   {
  130.    for( acum=0, i=0; i<k; i++ )
  131.     acum += ml[k][i] * lm[i][j];
  132.    lm[k][j] = ( (k==j) ? (float) 1.0 : - acum ) / ml[k][k];
  133.   }
  134. }
  135.  
  136. void calc_ld( void )
  137. {
  138.  int i,j,k;
  139.  
  140.  for( i=indmet; i<nmetab; i++ )                   /* first calculate L0      */
  141.   for( j=0; j<indmet; j++ )                       /* which is                */
  142.   {
  143.    ld[i][j] = 0;
  144.    for( k=0; k<indmet; k++ )                      /*                 -1      */
  145.     ld[i][j] += ml[i][k] * lm[k][j];              /* L0 = L21 * (L11)        */
  146.   }
  147.  for( i=indmet; i<nmetab; i++ )                   /* Now map the lin. dep.   */
  148.  {
  149.   for( j=0; j<indmet; j++ )                       /* which is                */
  150.    ld[i][j] = -ld[i][j];                          /* -L0                     */
  151.   for( j=indmet; j<nmetab; j++ )                  /* concateneted with the   */
  152.    ld[i][j] = (i == j)
  153.     ? (float) 1.0 : (float) 0.0;                  /* m0 -> m part of Id.     */
  154.  }
  155. }
  156.  
  157. void init_moiety( void )
  158. {
  159.  int i,j;
  160.  
  161.  for( i=indmet; i<nmetab; i++)
  162.  {
  163.   moiety[i] = (float) 0.0;
  164.   for( j=0; j<nmetab; j++ )
  165.    moiety[i] += ld[i][j] * x0[j];
  166.  }
  167. }
  168.  
  169. void lindep( void )
  170. {
  171.  invml11();                                       /* invert ml11 into lm     */
  172.  calc_ld();                                       /* calculate L2 * -L1'     */
  173. }
  174.  
  175. int gauss( void )
  176. {
  177.  int i, j, k, flag;
  178.  float m;
  179.  
  180.  for( k=0; k<nsteps+1; k++)
  181.  {
  182.   for( flag=0, i=k; i<nmetab; i++ )
  183.    if( fabs(rstoi[i][k]) > ALMOST_ZERO )
  184.    {
  185.     flag = 1;
  186.     break;                                        /* suitable divisor        */
  187.    }
  188.   if ( flag )
  189.   {
  190.    if ( i != k )
  191.    {
  192.     rowsw( k, i, nsteps );                        /* switch row k with i     */
  193.     rowsw_m( k, i, nmetab );
  194.    }
  195.   }
  196.   else
  197.   {
  198.    do
  199.    {
  200.     for( j=k+1; j<nsteps; j++ )
  201.      if( fabs(rstoi[k][j]) > ALMOST_ZERO )
  202.      {
  203.       flag = 1;
  204.       break;                                      /* suitable value          */
  205.      }
  206.     if ( flag && ( j != k ) )
  207.     {
  208.      colsw( j, k, nmetab );                       /* switch col j with k     */
  209.      colsw_m( j, k, nmetab );
  210.     }
  211.     if( !flag )
  212.     {
  213.      for( i=0; i<nmetab; i++ )                    /* get rid of empty rows   */
  214.      {
  215.       if ( ! emptyrow( i, nsteps ) )              /* all entries in row = 0  */
  216.       {
  217.        for( j=i+1; j<nmetab; j++)
  218.         if ( emptyrow( j, nsteps ) )              /* j = first non-empty row */
  219.         {
  220.          rowsw( i, j, nsteps );                   /* switch row i with row j */
  221.          rowsw_m( i, j, nmetab );
  222.          flag = 2;
  223.          break;
  224.         }
  225.       }
  226.      }
  227.     }
  228.     if ( !flag )
  229.     {
  230.      indmet = k;                                  /* # of independent metabs */
  231.      return(-1);                                  /* flag conservation & ret */
  232.     }
  233.    }
  234.    while( flag != 1 );
  235.   }
  236.   for( i=k+1; i<nmetab; i++)
  237.   {
  238.    if ( fabs(rstoi[k][k]) > ALMOST_ZERO )
  239.    {
  240.     m =  rstoi[i][k] / rstoi[k][k];               /* calculate the divisor   */
  241.     ml[i][k] = m;                                 /* and keep its symmetric  */
  242.     if( m != (float) 0.0 )
  243.     {
  244.      for( j=k; j<nsteps; j++ )                    /* reduce another row      */
  245.      {                                           
  246.       rstoi[i][j] = rstoi[i][j] - m * rstoi[k][j];
  247.       if( fabs(rstoi[i][j]) <= ALMOST_ZERO )      /* then make it a true zero*/
  248.        rstoi[i][j] = (float) 0.0;
  249.      }
  250.     }
  251.    }
  252.   }
  253.  }
  254.  return 0;                                        /* job done: no lin. dep.! */
  255. }
  256.  
  257. void int_ord( void )                       /* start by putting all internal  */
  258. {                                          /* metabolites as the first       */
  259.  int i,j,k;                                /* m rows in stoi                 */
  260.  
  261.  for( i=0, k=0; i<totmet; i++ )
  262.   if (intmet[i])
  263.   {
  264.    ur[k] = i;                              /* record the old number          */
  265.    for( j=0; j<nsteps; j++ )
  266.     stoi[k][j] = stoiu[i][j];
  267.    k++;
  268.   }
  269.  nmetab = k;                               /* the no.of internal metabolites */
  270. }
  271.  
  272. void ext_ord( void )                       /* now put the declared exernal   */
  273. {                                          /* metabolites in the last rows   */
  274.  int i,j,k;                                /* of stoi                        */
  275.  
  276.  k = nmetab;
  277.  for( i=0; i<totmet; i++ )
  278.   if (!intmet[i])
  279.   {
  280.    ur[k] = i;                              /* record the index       */
  281.    for( j=0; j<nsteps; j++ )
  282.     stoi[k][j] = stoiu[i][j];
  283.    k++;
  284.   }
  285.  nextmet = totmet - nmetab;                /* the no.of external metabolites */
  286.  for(i=0; i<nmetab; i++) intmet[i] = 1;    /* set intmet according with the  */
  287.  for(; i<totmet; i++) intmet[i] = 0;       /* mapping                         */
  288. }
  289.  
  290. void initreds( void )
  291. {
  292.  int i,j;
  293.  
  294.  for( i=0; i<nmetab; i++)
  295.  {
  296.   for( j=0; j<nsteps; j++)
  297.    rstoi[i][j] = (float) stoi[i][j];              /* init rstoi from stoi    */
  298.   for( j=0; j<nmetab; j++)
  299.    ml[i][j] = (float) 0;                          /* init matrix of multipl. */
  300.  }
  301. }
  302.  
  303. void more_ext( void)
  304. {
  305.  int i;
  306.  
  307.  for( i=0; i<nmetab; i++ )                        /* get rid of empty rows   */
  308.  {                                                /* which are xtrnl metabs. */
  309.   if ( ! emptyrow( i, nsteps ) )                  /* all entries in row = 0  */
  310.   {
  311.    rowsw( i, nmetab-1, nsteps );                  /* put row at the end      */
  312.    intmet[ur[nmetab]] = 0;                        /* signal this 1 is xtrnl! */
  313.    nmetab--;                                      /* signal 1 less int. met. */
  314.    nextmet++;                                     /* and 1 more extern. met. */
  315.   }
  316.  }
  317. }
  318.  
  319. void virt_step( void )
  320. {
  321.  int i, j, ct;
  322.  
  323.  for( i=0; i<nsteps; i++ )                        /* get rid of empty cols   */
  324.  {                                                /* which are 0 rate steps  */
  325.   for( ct=0, j=0; j<nmetab; j++)
  326.    if ( fabs( rstoi[j][i] ) > ALMOST_ZERO ) ct++; /* ct counts entries != 0  */
  327.   if ( !ct )                                      /* all entries in col = 0  */
  328.   {
  329.    colsw( i, nsteps-1, nmetab );                  /* put col at the end      */
  330.    nsteps--;                                      /* one less active step    */
  331.   }
  332.  }
  333. }
  334.  
  335. void reduce( void )
  336. {
  337.  int i;
  338.  
  339.  for( i=0; i<nsteps; i++) uc[i] = i;       /* setup reaction perm. vector    */
  340.  int_ord();                                /* put int. metabs. at the beggin.*/
  341.  ext_ord();                                /* put ext. metabs. at the end    */
  342.  initreds();                               /* reset rstoi and ml             */
  343.  more_ext();                               /* spot implicit external metabs. */
  344.  virt_step();                              /* spot virtual steps             */
  345.  gauss();                                  /* stoi -> rstoi by gaussian red. */
  346.  for( i=0; i<nsteps; i++ )                 /* kinetu -> kinetype (using uc)  */
  347.   kinetype[i] = kinetu[uc[i]];
  348.   
  349.  depmet = nmetab - indmet;                 /* depmet = no. of dep. metabs.   */
  350.  lindep();                                 /* work out linear dependencies   */
  351. }
  352.