getall.c

GESPI 2.0动态系统模拟工具

  if( Buff != NULL ) Buff++;
  if ( (Buff == NULL) || (sscanf( Buff,"%d %d %d %d %d %d %d %d %d",
	&options.dat, &options.datsep, &options.datwidth,
	&options.dattit, &options.datmca, &options.datss,
	&totsel, &options.append, &options.quotes ) < 9 ) )
   return IDS_ERR_BAD_OPTIONS;
  /* read the dat filename in one line                                                                  */
  Buff = strchr( Buff, '\n' );
  if( Buff != NULL )
   Buff++;
  if ( (Buff == NULL) || (sscanf( Buff,"%s", &options.datname[0] ) != 1 ) )
   return IDS_ERR_BAD_OPTIONS;
  paux = strchr( Buff, '\n' );
  if( paux == NULL) return IDS_ERR_BAD_OPTIONS;
  /* anything selected for output?	*/

  if( totsel )
  {
   /* allocate enough memory for output pointers and titles	*/
   poutelem = (PDBL *) realloc( poutelem, totsel * sizeof( PDBL ) );
   outtit = (char *) realloc( outtit, totsel * (options.datwidth+1) * sizeof( char ) );
   if( (poutelem == NULL) || (outtit == NULL) )
    return IDS_ERR_NOEXEC;
   /* read each item in one line (index and title)                                              */
   if( ver_no < (float) 2.05 ) adj = totmet;
   else adj = 0;
   for( i=0, ptr=outtit; i<totsel; i++ )
   {
    Buff = paux+1;
    if ( (Buff == NULL) || (sscanf( Buff,"%d", &j) < 1 ) )
     return IDS_ERR_BAD_OPTIONS;

    poutelem[i] = outpel[j+adj];

    Buff = strchr( Buff, ' ' );
    if( Buff != NULL ) for( ;*Buff == ' '; Buff++);
    paux = strchr( Buff, '\n' );
    if( paux == NULL ) return IDS_ERR_BAD_OPTIONS;
    /* take care of CR if one existed before LF                                         */
    if( *(paux-1)=='\r' ) *(paux-1) = '\0';
    *paux = '\0';
    if( strlen( Buff ) > (size_t) options.datwidth ) Buff[options.datwidth] = '\0';
    strcpy( ptr, Buff );
    ptr += strlen( ptr )+1;
   }
  }
  /* read the scan options in one line                                                                  */
  Buff = paux+1;
  if ( (Buff == NULL) || (sscanf( Buff, "%d %d %d %d %d", &options.scan,
			  &options.scandens, &options.scanlog,
			  &totscan, &nlinks) < 5 ) )
   return IDS_ERR_BAD_OPTIONS;

  /* read each scan item in one line                 */
  for( i=0; i<totscan; i++ )
  {
   Buff = strchr( Buff, '\n' );
   if( Buff != NULL )
    Buff++;
   if ( (Buff == NULL) || (sscanf( Buff,"%d %le %le %d %d", &(scindex[i]),
				   &low, &high, &dens, &log ) < 5 ) )
    return IDS_ERR_BAD_OPTIONS;
   sparam[scindex[i]].low  = low;
   sparam[scindex[i]].high = high;
   if( log )
    sparam[scindex[i]].ampl = log10( high ) - log10( low );
   else
    sparam[scindex[i]].ampl = high - low;
   sparam[scindex[i]].dens = dens;
   sparam[scindex[i]].log  = log;
  }

  /* read each link in one line                 */
  for( i=0; i<nlinks; i++ )
  {
   Buff = strchr( Buff, '\n' );
   if( Buff != NULL )
    Buff++;
   if ( (Buff == NULL) || (sscanf( Buff,"%d %d %d %le", &(lindex[i]),
				   &lkt, &op, &fac ) < 4 ) )
    return IDS_ERR_BAD_OPTIONS;
   sparam[lindex[i]].linkedto  = sparam[lkt].var;
   sparam[lindex[i]].operation = op;
   sparam[lindex[i]].factor = fac;
  }

  return 0;
}


/*
	read a simulation file
*/

int ReadSim( char *FName )
{
 int i, j, k, ty, fl, ch1, nRc;
 char *Buff;
 char *ptr;
 unsigned int bufsize;
 struct stat fst;
 unsigned char  nmod;
 BOOL   badtop;

 /* open the file                                                                                                       */
 if( (ch1 = open( FName, O_RDONLY )) != -1 )
 {
  /* get the file statistics                                                                            */
  fstat( ch1, &fst );

  /* set the buffer to the size of the file plus one                            */
  bufsize = (unsigned int) fst.st_size;
  if( bufsize<1024 ) bufsize=1024;

  /* allocate space for the read buffer and lock it                                     */
  ptr = (char  *) mem_malloc( bufsize+1 );
  if( ptr == NULL )
   return IDS_ERR_NOEXEC;
  Buff = ptr;

  /* read the file and close it                                                                         */
  read( ch1, Buff, (unsigned int) fst.st_size );
  close( ch1 );

  /* read the topology part                                                                                     */
  if( nRc = BufToTop( Buff, &Buff ) )
  {
   mem_free( ptr );
   return nRc;
  }

  /* set the number of user-def functions to zero                                       */
  nudf = 0;
  nrateq = MAX_TYP;
  sizetr = 4;
  treestr = (char  *) mem_realloc( treestr, sizetr * sizeof( char ) );
  if( treestr == NULL )
  {
   mem_free( ptr );
   return IDS_ERR_NOEXEC;
  }
  treeptr = treestr;

  /* reset kinetic type flags                                                                           */
  for( i=0; i<nsteps; i++ ) kfl[i] = 1;
  /* check if model contains user-defined kinetic types                         */
  /* if there are any read them                                                                         */
  for( i=0; i<nsteps; i++ )
   if( (kinetu[i] >= MAX_TYP) && (kfl[i]) )
   {
    /* look for a form-feed                                                                                     */
    Buff = strchr( Buff, '\f' );
    Buff++;
    if( (nRc=BufToTree( Buff, &Buff )) != 0)
    {
     mem_free( ptr );
     term_gepasi();
    }

    for( j=0, ty=-1; j<nudf; j++ )
     if( strcmp( tree[j].descr, tr.descr ) == 0 )
     {
      ty = j;
      break;
     }
    if( ty == -1 )
    {
     if( ( nRc=addtr( i ) ) )
     {
      mem_free( ptr );
      return nRc;
     }
     fl = 1;
    }
    else
     for( j=i,k=kinetu[i]; j<nsteps; j++ )
      if( (kinetu[j]==k) && (kfl[j]) )
      {
       kinetu[j] = MAX_TYP+ty;
       kfl[j] = 0;
      }
   }

  badtop = 0;
  if( (nsteps<=0) || (nsteps>MAX_STEP) ) badtop = IDS_ERR_STEP_EXCEED;
  if( (totmet<=0) || (totmet>MAX_MET)  ) badtop = IDS_ERR_MET_EXCEED;
  /* check if topology is complete                                                      */
  for( i=0; i<nsteps; i++ )
  {
   if( kinetu[i]==NOT ) badtop = IDS_ERR_BAD_TOP;
   /* count the number of substrates in rstr	*/
   for( j=k=0; j<MAX_MOL; j++ )
    if( (*rstr)[i][j] < 0 ) k++;
   /* and check if it is the same or less than in ktype	*/
   if( ktype[kinetu[i]].nsub > k ) badtop = IDS_ERR_BAD_KINTYPE;
   /* if the type is mass action, then it must be equal	*/
   if( ktype[kinetu[i]].nsub < k )
    if( !strcmp(ktype[kinetu[i]].descr,"mass action") ) badtop = IDS_ERR_BAD_KINTYPE;
   /* count the number of products in rstr	*/
   for( j=k=0; j<MAX_MOL; j++ )
    if( (*rstr)[i][j] > 0 ) k++;
   /* and check if it is the same or less than in ktype	*/
   if( ktype[kinetu[i]].npro > k ) badtop = IDS_ERR_BAD_KINTYPE;
   /* if the type is mass action, then it must be equal	*/
   if( ktype[kinetu[i]].npro < k )
    if( !strcmp(ktype[kinetu[i]].descr,"mass action") ) badtop = IDS_ERR_BAD_KINTYPE;
   if( ktype[kinetu[i]].nmodf > 0 )
   {
    /* count the number of assigned modfs for each reaction     */
    for( j=0, nmod=0; j<totmet; j++ )
     if( (*loop)[i][j] != 0 ) nmod++;
    if( ktype[kinetu[i]].nmodf != nmod ) badtop = IDS_ERR_BAD_TOP;;
   }
  }

  /* if topology is incomplete, signal and return                               */
  if( badtop )
  {
   mem_free( ptr );
   return badtop;
  }

  /* reduce the stoicheiometry matrix                                           */
  reduce();

  if( SetParams() == -1 )
  {
   mem_free( ptr );
   return IDS_ERR_NOEXEC;
  }

  if( SetEff() == -1 )
  {
   mem_free( ptr );
   return IDS_ERR_NOEXEC;
  }

  /* setup the output and scanning elements                             */
  if( SetOutpEl() == -1 )
  {
   mem_free( ptr );
   return IDS_ERR_NOEXEC;
  }

  nRc = BufToSim( Buff );
  mem_free( ptr );

  return nRc;
 }
 else return IDS_ERR_LOAD;
}


void eff_rateq( void )
{
 int i,j,k,l;

 for( i=0; i<nsteps; i++)
 {
  l = 0;
  /* assign indexes to the substrates           */
  for( k=0; k<MAX_MOL; k++ )
   if( (*rstr)[i][k] < 0 )
    eff[i][l++] = -(*rstr)[i][k]-1;
  /* assign indexes to the products                     */
  for( k=0; k<MAX_MOL; k++ )
   if( (*rstr)[i][k] > 0 )
    eff[i][l++] = (*rstr)[i][k]-1;
  /* assign indexes to the modifiers            */
  for( j=0; j < (int) ktype[kinetype[i]].nmodf; j++ )
   for( k=0; k<totmet; k++ )
    if( (int) (*loop)[i][k] == (j+1) )
     eff[i][l++] = k;
 }

 for( i=0; i<nsteps; i++ )                        /* first approach: all    */
  for( j=0; j<totmet; j++ )                       /* derivatives are zero   */
   partder[i][j] = ret_zero;

 for( i=0; i<nsteps; i++ )                        /* now point to the funcs */
 {
  switch (kinetype[i])                            /* that are != zero       */
  {
   case I01:
   case I10:  rateq[i]  = i0;
	      break;
   case I11:  rateq[i]  = i11;
	      partder[i][eff[i][0]] = d1Xs1;
	      break;
   case R11:  rateq[i] = r11;
	      partder[i][eff[i][0]] = d1Xs1;
	      partder[i][eff[i][1]] = dX1p1;
	      break;
   case I12:  rateq[i] = i12;
	      partder[i][eff[i][0]] = d1Xs1;
	      break;
   case R12:  rateq[i] = r12;
	      partder[i][eff[i][0]] = d1Xs1;
	      partder[i][eff[i][1]] = dX2p1s1;
	      partder[i][eff[i][2]] = dX2p2s1;
	      break;
   case I13:  rateq[i] = i13;
	      partder[i][eff[i][0]] = d1Xs1;
	      break;
   case R13:  rateq[i] = r13;
	      partder[i][eff[i][0]] = d1Xs1;
	      partder[i][eff[i][1]] = dX3p1s1;
	      partder[i][eff[i][2]] = dX3p2s1;
	      partder[i][eff[i][3]] = dX3p3s1;
	      break;
   case I21:  rateq[i] = i21;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      break;
   case R21:  rateq[i] = r21;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      partder[i][eff[i][2]] = dX1p1;
	      break;
   case I22:  rateq[i] = i22;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      break;
   case R22:  rateq[i] = r22;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      partder[i][eff[i][2]] = dX2p1s2;
	      partder[i][eff[i][3]] = dX2p2s2;
	      break;
   case I23:  rateq[i] = i23;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      break;
   case R23:  rateq[i] = r23;
	      partder[i][eff[i][0]] = d2Xs1;
	      partder[i][eff[i][1]] = d2Xs2;
	      partder[i][eff[i][2]] = dX3p1s2;
	      partder[i][eff[i][3]] = dX3p2s2;
	      partder[i][eff[i][4]] = dX3p3s2;
	      break;
   case I31:  rateq[i] = i31;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      break;
   case R31:  rateq[i] = r31;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      partder[i][eff[i][3]] = dX1p1;
	      break;
   case I32:  rateq[i] = i32;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      break;
   case R32:  rateq[i] = r32;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      partder[i][eff[i][3]] = dX2p1s3;
	      partder[i][eff[i][4]] = dX2p2s3;
	      break;
   case I33:  rateq[i] = i33;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      break;
   case R33:  rateq[i] = r33;
	      partder[i][eff[i][0]] = d3Xs1;
	      partder[i][eff[i][1]] = d3Xs2;
	      partder[i][eff[i][2]] = d3Xs3;
	      partder[i][eff[i][3]] = dX3p1s3;
	      partder[i][eff[i][4]] = dX3p2s3;
	      partder[i][eff[i][5]] = dX3p3s3;
	      break;
   case IMM:  rateq[i] = ihmm;
	      partder[i][eff[i][0]] = dihmms;
	      break;
   case RMM:  rateq[i] = hmm;
	      partder[i][eff[i][0]] = dhmms;
	      partder[i][eff[i][1]] = dhmmp;
	      break;
   case PSI:  rateq[i] = hmmsi;
	      partder[i][eff[i][0]] = dhmmsis;
	      partder[i][eff[i][1]] = dhmmsip;
	      partder[i][eff[i][2]] = dhmmsii;
	      break;
   case PCI:  rateq[i]  = hmmci;
	      partder[i][eff[i][0]] = dhmmcis;
	      partder[i][eff[i][1]] = dhmmcip;
	      partder[i][eff[i][2]] = dhmmcii;
	      break;
   case GOM:
   case MXI:  rateq[i]  = hmmmi;
	      partder[i][eff[i][0]] = dhmmmis;
	      partder[i][eff[i][1]] = dhmmmip;
	      partder[i][eff[i][2]] = dhmmmii;
	      break;
   case PSA:  rateq[i]  = hmmsa;
	      partder[i][eff[i][0]] = dhmmsas;
	      partder[i][eff[i][1]] = dhmmsap;
	      partder[i][eff[i][2]] = dhmmsaa;
	      break;
   case PCA:  rateq[i]  = hmmca;
	      partder[i][eff[i][0]] = dhmmcas;
	      partder[i][eff[i][1]] = dhmmcap;
	      partder[i][eff[i][2]] = dhmmcaa;
	      break;
   case MXA:  rateq[i]  = hmmma;
	      partder[i][eff[i][0]] = dhmmmas;
	      partder[i][eff[i][1]] = dhmmmap;
	      partder[i][eff[i][2]] = dhmmmaa;
	      break;
   case HIL:  rateq[i]  = hill;
	      partder[i][eff[i][0]] = dhills;
	      break;
   case UBS:  rateq[i]  = unibisa;
	      partder[i][eff[i][0]] = dunibisaa;
	      partder[i][eff[i][1]] = dunibisap;
	      partder[i][eff[i][2]] = dunibisaq;
	      break;
   case UBM:  rateq[i]  = unibisb;
	      partder[i][eff[i][0]] = dunibisba;
	      partder[i][eff[i][1]] = dunibisbp;
	      break;
   case ALI:  rateq[i]  = alinh;
	      for( j=0; j<3; j++ )
		   partder[i][eff[i][j]] = deff;
			  break;
   default:   rateq[i]  = translate;
	      k = tree[kinetype[i]-MAX_TYP].nsub +
		  tree[kinetype[i]-MAX_TYP].npro +
		  tree[kinetype[i]-MAX_TYP].nmodf;
	      for( j=0; j<k; j++ )
		   partder[i][eff[i][j]] = deff;
  }
 }
}