#include "copyleft.h"

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

/*************************************/
/*                                   */
/*   reduction of stoicheiometry by  */
/*   the Gauss method with row and   */
/*           column switches         */
/*                                   */
/*          MICROSOFT C 6.00         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/


#include <stdio.h>
#include <math.h>
#include <string.h>
#include "globals.h"
#include "globvar.h"

#define ALMOST_ZERO 1.0e-7

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

 for(j=0;j<c;j++)
 {
  dum         = stoi[r1][j];                      /* switch rows in stoi and */
  stoi[r1][j] = stoi[r2][j];
  stoi[r2][j] = dum;
  dummy        = rstoi[r1][j];                    /* switch rows in rstoi &  */
  rstoi[r1][j] = rstoi[r2][j];
  rstoi[r2][j] = dummy;
 }
 dumb = x[r1];                                    /* switch rows in x and    */
 x[r1] = x[r2];
 x[r2] = dumb;
 dumb = x0[r1];                                   /* switch rows in x0  and  */
 x0[r1] = x0[r2];
 x0[r2] = dumb;
 j      = ur[r1];                                 /* switch rows in ur too   */
 ur[r1] = ur[r2];
 ur[r2] = j;
 /*strcpy( buff, metname[r1] );
 strcpy( metname[r1], metname[r2] );
 strcpy( metname[r2], buff );*/
}

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

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

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

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

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

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

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

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

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

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

 for( i=indmet; i<nmetab; i++ )                   /* first calculate L0      */
  for( j=0; j<indmet; j++ )                       /* which is                */
  {
   ld[i][j] = 0;
   for( k=0; k<indmet; k++ )                      /*                 -1      */
    ld[i][j] += ml[i][k] * lm[k][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][j] = -ld[i][j];                          /* -L0                     */
  for( j=indmet; j<nmetab; j++ )                  /* concateneted with the   */
   ld[i][j] = (i == j)
    ? (float) 1.0 : (float) 0.0;                  /* m0 -> m part of Id.     */
 }
}

void init_moiety( void )
{
 int i,j;

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

void lindep( void )
{
 invml11();                                       /* invert ml11 into lm     */
 calc_ld();                                       /* calculate L2 * -L1'     */
}

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

 for( k=0; k<nsteps+1; k++)
 {
  for( flag=0, i=k; i<nmetab; i++ )
   if( fabs(rstoi[i][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][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][k]) > ALMOST_ZERO )
   {
    m =  rstoi[i][k] / rstoi[k][k];               /* calculate the divisor   */
    ml[i][k] = m;                                 /* and keep its symmetric  */
    if( m != (float) 0.0 )
    {
     for( j=k; j<nsteps; j++ )                    /* reduce another row      */
     {                                           
      rstoi[i][j] = rstoi[i][j] - m * rstoi[k][j];
      if( fabs(rstoi[i][j]) <= ALMOST_ZERO )      /* then make it a true zero*/
       rstoi[i][j] = (float) 0.0;
     }
    }
   }
  }
 }
 return 0;                                        /* job done: no lin. dep.! */
}

void int_ord( void )                       /* start by putting all internal  */
{                                          /* metabolites as the first       */
 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][j] = stoiu[i][j];
   k++;
  }
 nmetab = k;                               /* the no.of internal metabolites */
}

void ext_ord( void )                       /* now put the declared exernal   */
{                                          /* metabolites in the last rows   */
 int i,j,k;                                /* of stoi                        */

 k = nmetab;
 for( i=0; i<totmet; i++ )
  if (!intmet[i])
  {
   ur[k] = i;                              /* record the index       */
   for( j=0; j<nsteps; j++ )
    stoi[k][j] = stoiu[i][j];
   k++;
  }
 nextmet = totmet - nmetab;                /* the no.of external metabolites */
 for(i=0; i<nmetab; i++) intmet[i] = 1;    /* set intmet according with the  */
 for(; i<totmet; i++) intmet[i] = 0;       /* mapping						 */
}

void initreds( void )
{
 int i,j;

 for( i=0; i<nmetab; i++)
 {
  for( j=0; j<nsteps; j++)
   rstoi[i][j] = (float) stoi[i][j];              /* init rstoi from stoi    */
  for( j=0; j<nmetab; j++)
   ml[i][j] = (float) 0;                          /* init matrix of multipl. */
 }
}

void more_ext( void)
{
 int i;

 for( i=0; i<nmetab; i++ )                        /* get rid of empty rows   */
 {                                                /* which are xtrnl metabs. */
  if ( ! emptyrow( i, nsteps ) )                  /* all entries in row = 0  */
  {
   rowsw( i, nmetab-1, nsteps );                  /* put row at the end      */
   intmet[ur[nmetab]] = 0;                        /* signal this 1 is xtrnl! */
   nmetab--;                                      /* signal 1 less int. met. */
   nextmet++;                                     /* and 1 more extern. met. */
  }
 }
}

void virt_step( void )
{
 int i, j, ct;

 for( i=0; i<nsteps; i++ )                        /* get rid of empty cols   */
 {                                                /* which are 0 rate steps  */
  for( ct=0, j=0; j<nmetab; j++)
   if ( fabs( rstoi[j][i] ) > ALMOST_ZERO ) ct++; /* ct counts entries != 0  */
  if ( !ct )                                      /* all entries in col = 0  */
  {
   colsw( i, nsteps-1, nmetab );                  /* put col at the end      */
   nsteps--;                                      /* one less active step    */
  }
 }
}

void reduce( void )
{
 int i;

 for( i=0; i<nsteps; i++) uc[i] = i;       /* setup reaction perm. vector    */
 int_ord();                                /* put int. metabs. at the beggin.*/
 ext_ord();                                /* put ext. metabs. at the end    */
 initreds();                               /* reset rstoi and ml             */
 more_ext();                               /* spot implicit external metabs. */
 virt_step();                              /* spot virtual steps             */
 gauss();                                  /* stoi -> rstoi by gaussian red. */
 for( i=0; i<nsteps; i++ )                 /* kinetu -> kinetype (using uc)  */
  kinetype[i] = kinetu[uc[i]];
  
 depmet = nmetab - indmet;                 /* depmet = no. of dep. metabs.   */
 lindep();                                 /* work out linear dependencies   */
}