gaussw.c

GESPI 2.0动态系统模拟工具

#include "copyleft.h"

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

/*************************************/
/*                                   */
/*         GWSIM - Simulation        */
/*        MS-WINDOWS front end       */
/*                                   */
/*          Gauss reduction          */
/*            dialog box             */
/*                                   */
/*           QuickC/WIN 1.0          */
/*                                   */
/*   (include here compilers that    */
/*   compiled GWSIM successfully)    */
/*                                   */
/*************************************/

#include <windows.h>
#include <math.h>
#include <string.h>
#include "globals.h"
#include "gep2.h"
#include "strtbl.h"

#define ALMOST_ZERO 1.0e-7

GLOBALHANDLE	hRSto;					/* handle to memory block w/ rstoi		*/
GLOBALHANDLE	hSto;					/* handle to memory block w/ stoi		*/
GLOBALHANDLE	hMl;					/* handle to memory block w/ ml			*/
GLOBALHANDLE	hLm;					/* handle to memory block w/ lm			*/
GLOBALHANDLE	hLD;					/* handle to memory block w/ ld			*/
GLOBALHANDLE	hMetn;					/* handle to memory block w/ metn		*/
GLOBALHANDLE	hStepn;					/* handle to memory block w/ stepn		*/
GLOBALHANDLE	hStepst;				/* handle to memory block w/ stepst		*/
GLOBALHANDLE	hImet;					/* handle to memory block w/ imet		*/
GLOBALHANDLE	hRs;					/* handle to memory block w/ rs			*/
GLOBALHANDLE	hLoo;					/* handle to memory block w/ loo		*/
float			huge *rstoi;			/* pointer to work with rstoi array		*/
float 			huge *ml;				/* pointer to work with ml array		*/
float 			huge *lm;				/* pointer to work with lm array		*/
float 			huge *ld;				/* pointer to work with ld array		*/
int				huge *imet;				/* pointer to work with imet vector		*/
int				huge *loo;				/* pointer to work with loo array		*/
int		(huge *rs)[MAX_STEP][MAX_MOL];		/* pointer to work with rs array		*/
char	(huge *metn)[NAME_L];			/* pointer to work with metname array	*/
char	(huge *stepn)[NAME_L];			/* pointer to work with metname array	*/
char	(huge *stepst)[256];			/* pointer to work with stepst array	*/
int ur[MAX_MET];                        /* permut. on metabolites  u -> g 		*/
int uc[MAX_MET];                        /* permutations on steps   u -> g		*/

void rowsw( int r1, int r2, int c );
void colsw( int c1, int c2, int r );
void rowsw_m( int r1, int r2, int c );
void colsw_m( int c1, int c2, int r );
int emptyrow( int rw, int nsteps);
void invml11( void );
void calc_ld( void );
void init_moiety( void );
void lindep( void );
int gauss( void );
void more_ext( void);
void int_ord( void );
void ext_ord( void );
void initreds( void );
void virt_step( void );
int reduce( int swapnames );

#pragma alloc_text( CODE14, rowsw, colsw, rowsw_m, colsw_m, emptyrow, invml11, calc_ld, init_moiety, lindep, gauss, int_ord, ext_ord, more_ext, initreds, virt_step, reduce )

void rowsw( int r1, int r2, int c )
{
 register int j;
 float dummy;
 double dumb;
 int dum;

 for(j=0;j<c;j++)
 {
  dum         = stoi[r1*MAX_MET + j];                      /* switch rows in stoi and */
  stoi[r1*MAX_MET + j] = stoi[r2*MAX_MET + j];
  stoi[r2*MAX_MET + j] = dum;
  dummy        = rstoi[r1*MAX_MET + j];                    /* switch rows in rstoi &  */
  rstoi[r1*MAX_MET + j] = rstoi[r2*MAX_MET + j];
  rstoi[r2*MAX_MET + j] = dummy;
 }
 j      = ur[r1];                                 /* switch rows in ur and   */
 ur[r1] = ur[r2];
 ur[r2] = j;
}

void colsw( int c1, int c2, int r )
{
 register int j;
 float dummy;
 int dum;

 for( j=0; j<r; j++ )
 {
  dum         = stoi[j*MAX_MET + c1];                      /* switch cols in stoi and */
  stoi[j*MAX_MET + c1] = stoi[j*MAX_MET + c2];
  stoi[j*MAX_MET + c2] = dum;
  dummy        = rstoi[j*MAX_MET + c1];                    /* switch cols in rstoi &  */
  rstoi[j*MAX_MET + c1] = rstoi[j*MAX_MET + c2];
  rstoi[j*MAX_MET + c2] = dummy;
 }
 j      = uc[c1];                                 /* switch cols in uc too   */
 uc[c1] = uc[c2];
 uc[c2] = j;
}

void rowsw_m( int r1, int r2, int c )             /* switch rows in ml       */
{
 register int j;
 float dummy;

 for(j=0;j<c;j++)
 {
  dummy     = ml[r1*MAX_MET + j];
  ml[r1*MAX_MET + j] = ml[r2*MAX_MET + j];
  ml[r2*MAX_MET + j] = dummy;
 }
}

void colsw_m( int c1, int c2, int r )             /* switch cols in ml       */
{
 register int j;
 float dummy;

 for(j=0;j<r;j++)
 {
  dummy     = ml[j*MAX_MET + c1];
  ml[j*MAX_MET + c1] = ml[j*MAX_MET + c2];
  ml[j*MAX_MET + c2] = dummy;
 }
}

int emptyrow( int rw, int nsteps)
{
 register int j, ct;                              /* ct counts entries != 0  */

 for( ct=0, j=0; j<nsteps; j++)
  if ( fabs( rstoi[rw*MAX_MET + j] ) > ALMOST_ZERO ) ct++;
 return ( ct );
}

void invml11( void )
{                                                 /* as ml is lower triangul.*/
 register int i,j,k;                              /* inverting is simply to  */
 float acum;                                      /* forward-substitute with */
                                                  /* the identity matrix     */
 for(i=0;i<nmetab; i++) ml[i*MAX_MET + i] = 1 ;
 for( j=0; j<indmet; j++)
  for( k=0; k<indmet; k++)
  {
   for( acum=0, i=0; i<k; i++ )
    acum += ml[k*MAX_MET + i] * lm[i*MAX_MET + j];
   lm[k*MAX_MET + j] = ( (k==j) ? (float) 1.0 : - acum ) / ml[k*MAX_MET + k];
  }
}

void calc_ld( void )
{
 register int i,j,k;

 for( i=indmet; i<nmetab; i++ )                   /* first calculate L0      */
  for( j=0; j<indmet; j++ )                       /* which is                */
  {
   ld[i*MAX_MET + j] = 0;
   for( k=0; k<indmet; k++ )                      /*                 -1      */
    ld[i*MAX_MET + j] += ml[i*MAX_MET + k] * lm[k*MAX_MET + j];              /* L0 = L21 * (L11)        */
  }
 for( i=indmet; i<nmetab; i++ )                   /* Now map the lin. dep.   */
 {
  for( j=0; j<indmet; j++ )                       /* which is                */
   ld[i*MAX_MET + j] = -ld[i*MAX_MET + j];                          /* -L0                     */
  for( j=indmet; j<nmetab; j++ )                  /* concateneted with the   */
   ld[i*MAX_MET + j] = (i == j)
    ? (float) 1.0 : (float) 0.0;                  /* m0 -> m part of Id.     */
 }
}

void init_moiety( void )
{
 register int i,j;

 for( i=indmet; i<nmetab; i++)
 {
  moiety[i] = (float) 0.0;
  for( j=0; j<nmetab; j++ )
   moiety[i] += ld[i*MAX_MET + j] * xu[j];
 }
}

void lindep( void )
{
 invml11();                                       /* invert ml11 into lm     */
 calc_ld();                                       /* calculate L2 * -L1'     */
 init_moiety();                                   /* setup the moiety conc.  */
}

int gauss( void )
{
 register int i, j;
 int k, flag;
 float m;

 for( k=0; k<nsteps+1; k++)
 {
  for( flag=0, i=k; i<nmetab; i++ )
   if( fabs(rstoi[i*MAX_MET + k]) > ALMOST_ZERO )
   {
    flag = 1;
    break;                                        /* suitable divisor        */
   }
  if ( flag )
  {
   if ( i != k )
   {
    rowsw( k, i, nsteps );                        /* switch row k with i     */
    rowsw_m( k, i, nmetab );
   }
  }
  else
  {
   do
   {
    for( j=k+1; j<nsteps; j++ )
     if( fabs(rstoi[k*MAX_MET + j]) > ALMOST_ZERO )
     {
      flag = 1;
      break;                                      /* suitable value          */
     }
    if ( flag && ( j != k ) )
    {
     colsw( j, k, nmetab );                       /* switch col j with k     */
     colsw_m( j, k, nmetab );
    }
    if( !flag )
    {
     for( i=0; i<nmetab; i++ )                    /* get rid of empty rows   */
     {
      if ( ! emptyrow( i, nsteps ) )              /* all entries in row = 0  */
      {
       for( j=i+1; j<nmetab; j++ )
        if ( emptyrow( j, nsteps ) )              /* j = first non-empty row */
        {
         rowsw( i, j, nsteps );                   /* switch row i with row j */
         rowsw_m( i, j, nmetab );
         flag = 2;
         break;
        }
      }
     }
    }
    if ( !flag )
    {
     indmet = k;                                  /* # of independent metabs */
     return(-1);                                  /* flag conservation & ret */
    }
   }
   while( flag != 1 );
  }
  for( i=k+1; i<nmetab; i++)
  {
   if ( fabs(rstoi[k*MAX_MET + k]) > ALMOST_ZERO )
   {
    m =  rstoi[i*MAX_MET + k] / rstoi[k*MAX_MET + k];               /* calculate the divisor   */
    ml[i*MAX_MET + k] = m;                                 /* and keep its symmetric  */
    if( m != (float) 0.0 )
    {
     for( j=k; j<nsteps; j++ )
     {
      rstoi[i*MAX_MET + j] = rstoi[i*MAX_MET + j] - m*rstoi[k*MAX_MET + j];  /* reduce another row      */
      if( fabs(rstoi[i*MAX_MET + j]) <= ALMOST_ZERO )      /* then make it a true zero*/
       rstoi[i*MAX_MET + j] = (float) 0.0;
     }
    }
   }
  }
 }
 return 0;                                       /* job done: no lin. dep.! */
}

void int_ord( void )                       /* start by putting all internal  */
{                                          /* metabolites as the first       */
 register int i,j,k;                       /* m rows in stoi                 */

 for( i=0, k=0; i<totmet; i++ )
  if (intmet[i])
  {
   ur[k] = i;                              /* record the old number          */
   for( j=0; j<nsteps; j++ )
    stoi[k*MAX_MET + j] = stoiu[i*MAX_MET + j];
   k++;
  }