prog$ = "RNfb5p  created 27aug'03  Simulate rxn nets"
'Output: Sctout contains conc vs time of selected species over Scan2
'RNfb.1Q-14sep96-Does conc vs t.
'RNfb.1r-16sep96-add FN Scan2
'1t2-make graph output numbers floats. Compile with 10 digit precision.
'1v-add limits to dr in DoCell 13feb97  1w-delete routine for endpoint species
'1x-skip iterations in DoCell if k!(r%)<>0  28feb97
'2a-plot a Scan2 set on same graph
'2c-center operation on Scan2 mode and kinetic data output to RNtout (#3)
'2n-bugs fixed in Scan2 to display curves in proper order
'2p-link rate constants if e%()=-1
'2q-add ability to set time delays for rate consts.
'   Format: 'set'CR rxnCR tDelayk!(rxn) - works!
'2r-add scan tDelayk!()
'2T-workhorse for 1997 Oncogene paper
'2u-15nov97-add ability to make a species that is the sum of species,
'   using rxnID$() & "ADD"
'2v-17nov97-add ratio species using "DIV"
'2V5-30jun03-fix UpdateLinkedConstants
'5d-27july03-RndScan1-Simulate2  
'5m-add mode3%=-2 13aug'03     5n-add mode3%=-1  15aug'03
'5p-add kscan1   5p1-remove the bar graph calculation
'5R-add hysteresis (reversibility) check - used for HIF simulations.
    mode3% options are used to find opimum k-sets for switch-like response or
    to survey k-space.
 
' -------------------- Globals --------------------------
RxnFi$="RxnFi$ not yet set"
mode1%=0
mode2%=0
mode3%=0                                          'for RndScan1

RxnFi$=""
nrxns%=12                                         'max num of reactions
cells%=2                                          'max num of cells
nspecies%=12                                      'max num of molecular species
nplots%=1
hscale!=1
gheight1%=0
tmax! =100
dt1! =.1

mode1$=""
outputspecies%=1                                  'for RndScan
varspecies%=1
x1varsp!=0
x2varsp!=0
lowerlimit!=0                                     
rndfactor!=1
rndfactorlimit!=0
sim%=0
basisSet%=0

kdist!=0

DEF FN showTopics USING gShowTopics&

v0%= 200                                          'Position of horizontal axis
hmax%= 1100

black%=1
nsteps%=100
nsubsteps%=10
rstep%=0
stepslope1!=0
stepslope2!=0

'total steps = nsteps%*nsubsteps%
'total time  = nsteps%*nsubsteps%*dt = nsteps%*dt1
dt!=dt1!/nsubsteps%                               'dt1 = dt * nsubsteps%
tout%=1                                           'flag to record kinetic curves in RNtout.
_gtop=10                                          '10  20
_gwidth=850                                       '640  950
_gheight=750                                      '850   575  394
glines!=(_gheight-_gtop)/12.5

xx!=0                                             'x scan parameter
yy!=0                                             'y scan parameter
iplot%=0                                          'index of plots outputed to RNtout

itemtype1%=0
itemtype2%=0
scanitem1%=0
scanitem2%=0

_DIM1=12
_DIM2=400
_DIM3=200
DIM k!(_dim3)                                     'rate constants
DIM k0!(_dim3)                                    'initial rate constants
DIM kBasis!(_dim3)
DIM kRnd!(_dim3)
DIM x!(_dim3)                                     'concentrations (molecules per cell)
DIM x0!(_dim3)                                    'initial concentrations
DIM plotsp%(_dim1)
DIM xt!(_dim1,_dim2)                              'conc (selected species , time)
DIM xs!(_dim1)                                    'list of selected species
DIM xmax!(_dim1)                                  'max conc of selected species

DIM i1%(_dim3),i2%(_dim3),i3%(_dim3)              'input species
DIM o1%(_dim3),o2%(_dim3),o3%(_dim3)              'output species
DIM e%(_dim3)                                     'enzyme species
DIM rxnID$(_dim3)                                 
DIM vscale!(_dim3)
DIM d!(_dim3)                                     'Change in conc (rxn)    (molec/cell)
DIM x$(_dim3)
DIM tDelayk!(_dim3)                               't when k!() stops being 0

END GLOBALS

RANDOMIZE

' ----------- Functions -------------

LOCAL FN ReadRxns
  PRINT #3,"Reactions:"
  nrxns%=0
  DO
    nrxns%=nrxns%+1
    INPUT#1,j%
    INPUT#1,i1%(j%),i2%(j%),i3%(j%),o1%(j%),o2%(j%),o3%(j%),e%(j%),k0!(j%),q$
    rxnID$(j%)=q$
    LONG IF e%(j%)<0                              '2p
      k!(j%) = k0!(j%) * k!(-(e%(j% )))
    XELSE
      k!(j%)=k0!(j%)
    END IF                                        
    PRINT j%;i1%(j%);i2%(j%);i3%(j%);o1%(j%);o2%(j%);o3%(j%);e%(j%);k!(j%),q$
    PRINT #3,j%;",";i1%(j%);",";i2%(j%);",";i3%(j%);",";o1%(j%);",";o2%(j%);",";o3%(j%);",";e%(j%);",";k!(j%);",";q$
    rxnID$(j%)=q$
  UNTIL q$="end"
  nrxns%=nrxns%-1
  PRINT "nrxns%=";nrxns%
  PRINT #3,"nrxns%=";nrxns%
END FN

LOCAL FN ReadSpecies
  PRINT #3, "Species:"
  k%=0                                            'counter for species to be plotted; k% -> nplots% & index for plotsp%()
  FOR j%=1 TO _dim1
    plotsp%(j%)=0
  NEXT j%
  q$=""
  nspecies%=0
  DO
    nspecies%=nspecies%+1
    INPUT#1,i%
    j%=ABS(i%)
    INPUT#1,x0!(j%),q!,q$
    x$(j%)=q$
    PRINT j%;x0!(j%);q!,x$(j%)
    PRINT #3,j%;",";x0!(j%);",";q!;",";x$(j%)
    LONG IF i%<0
      k%=k%+1
      plotsp%(k%)=j%                              'j% is species number in plot k%
    END IF
    LONG IF q!<>0
      endptSpecies% = j%
      endptConc! = q!
      PRINT "endptSpecies=";endptSpecies%;"  endptConc!=";endptConc!
    END IF
  UNTIL q$="end"
  nspecies%=nspecies%-1
  nplots%=k%
  PRINT "nspecies%=";nspecies%;"  nplots%=";nplots%
  PRINT #3,"nspecies%=";nspecies%
END FN

CLEAR LOCAL FN InitializeSpecies
  FOR i%=0 TO _DIM1
    FOR j%=0 TO _DIM2
      xt!(i%,j%)=0
    NEXT j%
  NEXT i%
  FOR j%=1 TO nspecies%
    x!(j%) = x0!(j%)                              'initial concs
    IF x!(j%)=0 THEN x!(j%)=.000001
  NEXT j%
  FOR j%=1 TO nrxns%
    IF tDelayk!(j%) > 0 THEN k!(j%)=-ABS(k!(j%))
  NEXT j%
END FN

LOCAL FN Axes
  WINDOW 2
  CLS
  TEXT _monaco,10,,_srcCopy
  PLOT 0,_gtop TO _gwidth,_gtop
  FOR i%=1 TO nplots%
    k%=plotsp%(i%)                                'k% is species number in plot i%
    h%=0
    gheight1%=(_gheight-_gtop)/nplots%
    v0%= _gtop + gheight1% * i%
    PLOT 0,v0% TO _gwidth,v0%
    iline% = 1+INT(0.5+gheight1%*(i%-1)/_gheight * glines!)
    LOCATE 50,iline%
    PRINT k%;
    PRINT USING "#.###^^^^";xmax!(i%);
    PRINT USING " #.###^^^^ ";x!(k%);
    PRINT x$(k%);
    IF tout%=1 THEN PRINT #3,k%;",";
    'IF k%=endptSpecies% THEN PRINT USING"##.###"; endptTime!
  NEXT i%
  PRINT #3, " species"
END FN

LOCAL FN Splot
  WINDOW 2
  'TEXT _monaco,10,,_srcCopy
  FOR i%=1 TO nplots%
    PLOT 0,v0%
    gwid!=_gwidth
    PRINT #3,"//"
    iplot%=iplot%+1
    PRINT #3,iplot%
    PRINT #3,"species";",";plotsp%(i%)
    PRINT #3,"xx! yy!= ,";xx!;",";yy!
    FOR j%=1 TO nsteps%
      'h%=h%+hstep%
      nst!=nsteps%
      h%=INT(j%*gwid!/nst!)
      v0%= _gtop + gheight1% * i%
      temp1!=(xt!(i%,j%) / xmax!(i%))
      v%=v0%- temp1! * gheight1%
      
      PLOT h%,v%                      
      IF tout%=1 THEN PRINT #3,xt!(i%,j%);",";    'v%;v0%;i%;j%;xt!(i%,j%);xmax!(i%);temp1!;gheight1%;
    NEXT j%
    
    k%=plotsp%(i%)
    iline% = 1+INT(0.5+gheight1%*(i%-1)/_gheight * glines!)
    LOCATE 50,iline%
    PRINT SPACE$(30);
    LOCATE 50,iline%
    PRINT k%;
    PRINT USING "#.###^^^^";xmax!(i%);
    PRINT USING " #.###^^^^ ";x!(k%);
    PRINT x$(k%);
    
    IF tout%=1 THEN PRINT #3,-99.9
  NEXT i%
END FN

LOCAL FN RecordFinalConcs
  PRINT #5," scanitem1%;scanitem2%=";scanitem1%;scanitem2%
  IF itemtype1%=0 THEN PRINT #5," x0!(scanitem1%)=";x0!(scanitem1%)     
  IF itemtype1%=1 THEN PRINT #5," k!(scanitem1%)=";k!(scanitem1%)
  IF itemtype2%=0 THEN PRINT #5," x0!(scanitem2%)="x0!(scanitem2%)    
  IF itemtype2%=1 THEN PRINT #5," k!(scanitem2%)=";k!(scanitem2%)
  PRINT #5, "   Initial concentrations:'
  FOR j%=1 TO nspecies%
    PRINT #5,j%;",";x0!(j%);",";"0";",";x$(j%)
  NEXT j%
  PRINT #5, "   Final concentrations:"
  FOR j%=1 TO nspecies%
    PRINT #5,j%;",";x!(j%);",";"0";",";x$(j%)
  NEXT j%
  PRINT #5,
END FN

LOCAL FN DoCell
  totalsteps%=nsteps%*nsubsteps%
  t! = 0                                          'time
  tstep%=0
  rstep%=0                                        'sequential number of recorded step
  's1%,s2%,s3%,s4% -- molecular species
  'd(rxn) -- change in conc (molec/cell)
  FOR j%=1 TO nplots%
    xmax!(j%)=0
  NEXT j%
  WHILE tstep% < totalsteps%
    t!=t!+dt!                                     'Next time step
    tstep%=tstep%+1
    'Determine d(r%):
    FOR r%=1 TO nrxns%
      LONG IF rxnID$(r%)<>"ADD" AND rxnID$(r%)<>"DIV"
        LONG IF k!(r%) < 0
          IF t! > tDelayk!(r%) THEN k!(r%)=ABS(k!(r%))
        END IF
        LONG IF k!(r%) > 0
          IF r%=kscan% THEN  k!(r%) = EXP( logkmin! + (t!/tmax!) * (logkmax! - logkmin!))
          'kscan% doesn't appear elsewhere, so the above statement could be deleted
          s1% = ABS(i1%(r%))                      '1st input species for rxn r%
          d!(r%) = k!(r%)*x!(s1%)*dt!             'every rxn is at least 1st order
          s2% = ABS(i2%(r%))
          IF s2%>0 THEN  d!(r%)=d!(r%)*x!(s2%)    'if rxn is 2nd order
          s3% = ABS(i3%(r%))
          IF s3%>0 THEN d!(r%)=d!(r%)*x!(s3%)     'if rxn has 3 input species
          s4% = e%(r%)                            'enzyme species
          IF s4%>0 THEN  d!(r%)=d!(r%)*x!(s4%)    'If rxn is enzyme mediated
        END IF
      END IF
    NEXT r%
    FOR r%=1 TO nrxns%
      'Deplete input species unless flagged negative:
      LONG IF k!(r%) > 0
        ss1%=i1%(r%)                              'Input species of rnx r%
        s1%=ABS(ss1%)
        ss2%=i2%(r%)
        s2%=ABS(ss2%)
        ss3%=i3%(r%)
        s3%=ABS(ss3%)
        LONG IF rxnID$(r%)="ADD" OR rxnID$(r%)="DIV"
          sum1!=x!(s1%)+x!(s2%)+x!(s3%)
          'PRINT r%;sum1!,
        XELSE
          LONG IF ss1% >0
            IF d!(r%) > x!(s1%) THEN d!(r%) = x!(s1%) - 1e-12
          END IF
          LONG IF ss2% >0
            IF d!(r%) > x!(s2%) THEN d!(r%) = x!(s2%) - 1e-12
          END IF
          LONG IF ss3% >0
            IF d!(r%) > x!(s3%) THEN d!(r%) = x!(s3%) - 1e-12
          END IF
          
          IF ss1%>0 THEN x!(s1%)=x!(s1%)-d!(r%)   'Don't deplete if flagged negative
          IF ss2%>0 THEN x!(s2%)=x!(s2%)-d!(r%)
          IF ss3%>0 THEN x!(s3%)=x!(s3%)-d!(r%)
        END IF
        
        'Enhance output species unless its conc is flagged negative:
        
        s1%=o1%(r%)                               'Now output species of the rxn
        s2%=o2%(r%)
        s3%=o3%(r%)
        SELECT CASE rxnID$(r%)
          CASE "ADD" 
            x!(s3%)=sum1!+x!(s1%)+x!(s2%)
          CASE "DIV"
            x!(s3%)=sum1!/(x!(s1%)+x!(s2%))
          CASE ELSE
            IF s1%>0 THEN x!(s1%)=x!(s1%)+d!(r%)
            IF s2%>0 THEN x!(s2%)=x!(s2%)+d!(r%)
            IF s3%>0 THEN x!(s3%)=x!(s3%)+d!(r%)
        END SELECT
      END IF
    NEXT r%
    
    LONG IF tstep% MOD nsubsteps% = 0
      rstep%=rstep%+1                             'sequential number of step recorded
      LONG IF mode1$<>"rndscan1" AND mode1$<>"kscan1"
        LOCATE 60,4
        CLS LINE 
        PRINT rstep%
      END IF
      FOR j%=1 TO nplots%
        k%=plotsp%(j%)                            'k% is species number in plot j%
        xt!(j%, rstep%) = x!(k%)
        IF xmax!(j%)<x!(k%) THEN xmax!(j%)=x!(k%)
        'LONG IF j%=endptSpecies% AND endptConc! > 0 AND tstep% >1
        'LONG IF xt!(j%, tstep%-1) < endptConc! AND xt!(j%, tstep%) >= endptConc!
        'endptTime! = tstep%*dt!
        'PRINT "xt!(j%, tstep%-1)="; xt!(j%, tstep%-1); " xt!(j%, tstep%)=";xt!(j%, tstep%)
        'PRINT "endptTime!=";endptTime!;"  species ";j%;"  tstep%=";tstep%
        'END IF
        'END IF
      NEXT j%
    END IF
  WEND
END FN                                            'END FN DoCell

LOCAL FN UpdateLinkedConstants                    '2p
  FOR i%=1 TO nrxns%
    IF e%(i%)<0 THEN k!(i%) = k!(-(e%(i%)))*k0!(i%)
    'IF e%(i%)<0 THEN PRINT #4, "i;e;k(-e);k0;k(i)=";i%;e%(i%);k!(-(e%(i%)));k0!(i%);k!(i%)
    'IF e%(i%)<0 THEN PRINT i%;k!(i%);
  NEXT i%
  PRINT
END FN

LOCAL FN UpdateLinkedBasisConstants               '2p
  FOR i%=1 TO nrxns%
    IF e%(i%)<0 THEN kBasis!(i%) = kBasis!(-(e%(i%)))*k0!(i%)
    IF e%(i%)<0 THEN PRINT #4, "i;e;k(-e);k0;k(i)=";i%;e%(i%);kBasis!(-(e%(i%)));k0!(i%);k!(i%)
  NEXT i%
END FN

LOCAL FN UpdateLinkedRndConstants                 '2p
  FOR i%=1 TO nrxns%
    IF e%(i%)<0 THEN kRnd!(i%) = kRnd!(-(e%(i%)))*k0!(i%)
    IF e%(i%)<0 THEN PRINT #4, "i;e;k(-e);k0;k(i)=";i%;e%(i%);kRnd!(-(e%(i%)));k0!(i%);k!(i%)
  NEXT i%
END FN

LOCAL FN Scan2
  WINDOW 2
  graph%=1
  FN Axes
  'WINDOW 1
  OPEN "O",2,"Scout"
  OPEN "O",4,"RnxNetUpdates"
  OPEN "O",5,"FinalConcs"
  INPUT #1,sp%
  INPUT #1,itemtype1%,scanitem1%,mode1%,x1!,x2!,xd!
  INPUT #1,itemtype2%,scanitem2%,mode2%,y1!,y2!,yd!
  PRINT #3,sp%
  PRINT #3,itemtype1%;",";scanitem1%;",";mode1%;",";x1!;",";x2!;",";xd!
  PRINT #3,itemtype2%;",";scanitem2%;",";mode2%;",";y1!;",";y2!;",";yd!
  PRINT #3,"**"
  PRINT #3, RxnFi$
  PRINT #3, "itemtype,scanitem"
  PRINT #3, itemtype1%;",";scanitem1%
  PRINT #3, itemtype2%;",";scanitem2%
  
  yy!=y1!                                         'record the yy values to be used
  PRINT #2, "y1=";y1!;
  DO                                               
    LONG IF mode2%=1
      yy!=yy!*yd!
    XELSE
      yy!=yy!+yd!
    END IF
    PRINT #2, " yy=";yy!;                              
  UNTIL yy!>=y2!
  PRINT #2,
  'Begin simulations
  xx!=x1!                                         'Set initial value of xx
  IF xx!=0 THEN xx!=.000001
  DO
    yy!=y1!                                       'Set initial value of yy
    IF yy!=0 THEN yy!=.000001
    PRINT #2, " xx=";xx!;
    DO                                            'Are xx,yy init concs (x0) or rate consts (k)? 
      IF itemtype1%=0 THEN x0!(scanitem1%)=xx!      
      IF itemtype1%=1 THEN k!(scanitem1%)=xx!
      IF itemtype1%=2 THEN tDelayk!(scanitem1%)=xx!
      'yy is init conc (x0) or rate const (k)
      IF itemtype2%=0 THEN x0!(scanitem2%)=yy!    
      IF itemtype2%=1 THEN k!(scanitem2%)=yy!
      IF itemtype2%=2 THEN tDelayk!(scanitem2%)=yy!
      
      FN UpdateLinkedConstants                    'UpdateLinkedConstants
      
      FN InitializeSpecies                        'InitializeSpecies
      
      FN DoCell                                   'DoCell
      FN RecordFinalConcs
      
      IF graph%=1 THEN FN Splot
      LOCATE 60,3
      CLS LINE
      PRINT xx!;yy!,
      'PRINT #3, "end ,";xx!;",";yy!
      PRINT #2, " sp,x(sp)=";sp%;x!(sp%);
      
      key$=INKEY$
      LONG IF key$="." 
        INPUT q$
        IF q$="." THEN STOP
      XELSE
        IF key$="g" THEN graph%=0
      XELSE
        IF key$="G" THEN graph%=1
      END IF
      
      LONG IF mode2%=1                            'increment yy
        yy!=yy!*yd!
      XELSE
        yy!=yy!+yd!
      END IF
      'FN UpdateLinkedConstants                    '2p-UpdateLinkedConstants
    UNTIL yy!>y2!                                 'last yy
    PRINT #2,
    
    LONG IF mode1%=1                              'increment xx
      xx!=xx!*xd!
    XELSE
      xx!=xx!+xd!
    END IF
  UNTIL xx!>x2!                                   'last xx
  CLOSE #2
  CLOSE #4
  CLOSE #5
END FN                                            'END FN Scan2

LOCAL FN Start
  q$="fudge"
  WHILE q$<>""
    FN DoCell
    FN Axes
    FN Splot
    q$=""
    INPUT ""; q$
    IF q$="." THEN STOP
    IF q$="c" THEN q$="fudge"
  WEND
END FN

LOCAL FN NextBasisSet
  FOR i%=1 TO nrxns%
    FOR i%=1 TO nrxns%
      kBasis!(i%)=k0!(i%)
    NEXT i%
  NEXT i%
END FN                                            'END FN NextRndSet

LOCAL FN NextRndSet
  PRINT "NextRndSet:"
  FOR i%=1 TO nrxns%
    LONG IF RND(1) > 0.5
      kRnd!(i%)=k!(i%) * rndfactor!
    XELSE
      kRnd!(i%)=k!(i%) / rndfactor!
    END IF
    PRINT i%;rndfactor!;kRnd!(i%)
    INPUT;q$
  NEXT i%
END FN                                            'END FN NextRndSet

LOCAL FN Update_DoCell
  FN UpdateLinkedConstants                        
  FN InitializeSpecies                            
  FN DoCell
END FN

LOCAL FN MakeOptRxnFi
  OPEN "O",4,"OptRxnFi"
  PRINT #4,DATE$;" ";TIME$
  PRINT #4,"Prog$= ";prog$
  PRINT #4,"RxnFi$= ";RxnFi$
  PRINT #4,"mode1$= ";mode1$
  PRINT #4,"outputspecies%=";outputspecies%
  PRINT #4,"varspecies%=";varspecies%
  PRINT #4,"mode3%=";mode3%
  PRINT #4,"x1varsp!=";x1varsp!
  PRINT #4,"x2varsp!=";x2varsp!
  PRINT #4,"last lowerlimit!=";lowerlimit!
  PRINT #4,"last rndfactor!=";rndfactor!
  PRINT #4,"ntrials%=";ntrials%
  PRINT #4, "*"
  PRINT #4,dt1!;",";nsteps%;",";nsubsteps%
  FOR j%=1 TO nrxns%
    PRINT #4,j%;",";i1%(j%);",";i2%(j%);",";i3%(j%);",";
    PRINT #4,o1%(j%);",";o2%(j%);",";o3%(j%);",";
    PRINT #4,e%(j%);",";
    LONG IF e%(j%)<0
      PRINT #4,k0!(j%);",";rxnID$(j%)
    XELSE
      PRINT #4,k!(j%);",";rxnID$(j%)
    END IF
  NEXT j%
  PRINT #4,"0";",";"0";",";"0";",";"0";",";"0";",";"0";",";"0";",";"0";",";"0";",";"end"
  FOR j%=1 TO nspecies%
    PRINT #4,j%;",";x0!(j%);",";"0";",";x$(j%)
  NEXT j%
  PRINT #4,"0";",";"0";",";"0";",";"end"
  PRINT #4,"rndscan1"
  PRINT #4,outputspecies%;",";varspecies%;",";mode3%;",";
  PRINT #4,x1varsp!;",";x2varsp!;",";lowerlimit!;",";rndfactor!;",";ntrials%
  PRINT #4,"stop"
  PRINT #4," =========================="
  PRINT #4,"  final concentrations x!():
  FOR j%=1 TO nspecies%
    PRINT #4,j%;",";x!(j%);",";"0";",";x$(j%)
  NEXT j%
  CLOSE 4
END FN

LOCAL FN DoKey
  PRINT "1)change mode3%  2}change x1varsp!,x2varsp!  3)change rndfactor";
  PRINT "  4)change lowerlimit  5)make OptRxnFi  0)STOP : ";
  INPUT q%
  SELECT
    CASE q%=1
      PRINT"Mode3 options: 0)without selection" 
      PRINT " 1)maximize x1endpt!/x2endpt! (with limits)"
      PRINT " 2)maximize x1endpt!-x2endpt! (with limits)"
      PRINT " -1)randomize rndfactor= 9*(RND(32000)/32000!)+1"'make rndfactor betw 1 and 10
      PRINT " -2)randomize rndfactor= 10^(RND(32000)/32000!)
      PRINT mode3%; 
      INPUT " enter new mode3%:";mode3%
    CASE q%=2
      PRINT "  ";x1varsp!;x2varsp!;
      INPUT "  enter new x1varsp!,x2varsp! : ";x1varsp!,x2varsp!
    CASE q%=3
      PRINT rndfactor!;
      INPUT " enter new rndfactor!: ";rndfactor!
      PRINT #3, "new rndfactor!=";rndfactor!
    CASE q%=4
      PRINT lowerlimit!;
      INPUT " enter new lowerlimit!!: ";lowerlimit!
      PRINT #3, "new lowerlimit!=";lowerlimit!
    CASE q%=5
      FN MakeOptRxnFi
    CASE ELSE
      INPUT "  stop? ";q1$
      IF LEFT$(q1$,1)="y" OR q1$="." THEN STOP
  END SELECT
END FN

LOCAL FN Simulate2!(varSp%,x1!,x2!,outSp%)
  sim%=sim%+1
  slopefraction!=0.95
  
  FN InitializeSpecies
  
  x!(varSp%)= x1!
  FN DoCell
  x1endpt!=x!(outSp%)
  IF x1endpt!<lowerlimit! THEN x1endpt!=lowerlimit!
  xt_end! = xt!(1, rstep%)
  xt_earlier! = xt!(1, INT(slopefraction! * rstep%))
  stepslope1!= ((xt_end! - xt_earlier!) / xmax!(1)) / (1-slopefraction!)
  PRINT "  x1=";x!(varSp%);                       'x!(outSp%);"==";xt!(1, rstep%)
  PRINT " xt_end!;xt_earlier!= ";xt_end!;xt_earlier!;
  PRINT " stepslope1=";stepslope1!
  'PRINT #3, x!(varSp%);",";x!(outSp%);",";stepslope1!;",";
  
  FN InitializeSpecies
  x!(varSp%)= x2!
  FN DoCell
  x2endpt!= x!(outSp%)
  IF x2endpt!<lowerlimit! THEN x2endpt!=lowerlimit!
  IF x2endpt!=0 THEN x2endpt!=1e-9
  xt_end! = xt!(1, rstep%)
  xt_earlier! = xt!(1, INT(slopefraction! * rstep%))
  stepslope2!= ((xt_end! - xt_earlier!) / xmax!(1)) / (1-slopefraction!)
  
  PRINT "  x2=";x!(varSp%);                       'x!(outSp%);"==";xt!(1, rstep%)
  PRINT " xt_end!;xt_earlier!= ";xt_end!;xt_earlier!;
  PRINT " stepslope2=";stepslope2!
  PRINT "x1endpt;x2endpt=";x1endpt!;x2endpt!;
  
  SELECT
    CASE mode3%<=0
      result2!=0                                  'no selection
    CASE mode3%=1
      result2!= x1endpt! / x2endpt!
    CASE mode3%=2
      result2!= x1endpt!^2 / x2endpt!
    CASE mode3%=3
      result2! = x1endpt! - x2endpt!
    CASE ELSE
      PRINT "invalid mode3%=";mode3%
      INPUT q$
      STOP
  END SELECT
  
  PRINT "result2!=";result2!;"  ratio=";x1endpt!/x2endpt!;
  PRINT #3, x1endpt!;",";stepslope1!;",";x2endpt!;",";stepslope2!;",";
  PRINT #3, x1endpt!/x2endpt!;",";x1endpt!-x2endpt!;",";result2!;",";kdist!;",";
  PRINT #3, mode3%;",";sim%;",";rndfactor!;",";
END FN = result2!

LOCAL FN RndScan1
  nplots%=1
  
  'Headings:
  CLS
  INPUT #1,outputspecies%,varspecies%,mode3%
  INPUT #1,x1varsp!,x2varsp!,lowerlimit!,rndfactor!,ntrials%
  LONG IF mode3%=-2
    BEEP:PRINT:PRINT:BEEP:PRINT "Enter rndfactorlimit: ";
    INPUT rndfactorlimit!
  END IF
  PRINT #3,DATE$;" ";TIME$
  PRINT #3,"Prog$= ";prog$
  PRINT #3,"RxnFi$= ";RxnFi$
  PRINT #3,"mode1$= ";mode1$
  PRINT #3,"outputspecies%=";outputspecies%
  PRINT #3,"varspecies%=";varspecies%
  PRINT #3,"mode3%=";mode3%
  PRINT #3,"x1varsp!=";x1varsp!
  PRINT #3,"x2varsp!=";x2varsp!
  PRINT #3,"last lowerlimit!=";lowerlimit!
  PRINT #3,"last rndfactor!=";rndfactor!
  PRINT #3,"ntrials%=";ntrials%
  PRINT #3, "-" ; "," ;
  FOR i%=1 TO nrxns%
    PRINT #3,"rxn";i%;",";
  NEXT i%
  PRINT #3,"x1endpt";",";"stepslope1";",";"x2endpt";",";"stepslope2";",";
  PRINT #3,"x1endpt/x2endpt";",";"x1endpt-x2endpt";",";"result2";",";"kdist!";",";
  PRINT #3,"mode3%";",";"sim%";",";"rndfactor!"
  
  'Get the first Basis Set:
  PRINT:PRINT "FirstBasisSet":
  FOR i%=1 TO nrxns%
    LONG IF e%(i%) >= 0 
      kBasis!(i%)= k0!(i%)
      k!(i%)= kBasis!(i%)
    END IF
  NEXT i%
  FN UpdateLinkedConstants
  
  'Simulate the first basis set:
  PRINT #3, "-" , "," ;
  FOR i%=1 TO nrxns%
    PRINT #3,k!(i%);",";
  NEXT i%
  
  criterion2A! = FN Simulate2!(varspecies%,x1varsp!,x2varsp!,outputspecies%)
  PRINT "  cr2A!=";criterion2A!
  PRINT #3, criterion2A!;",";"cr2A"               'result from the first basis set
  
  'Start randomized trials:
  FOR trial%=1 TO ntrials%
    
    'key in for changes:
    key$=INKEY$
    IF key$="." THEN FN DoKey
    
    'Get the next random set:
    PRINT "NextRndSet:  trial";trial%;"  mode3=";mode3%;
    IF mode3%=-1 THEN rndfactor!= (rndfactorlimit!-1)*(RND(32000)/32000!) + 1
    IF mode3%=-2 THEN rndfactor!= EXP(LOG(rndfactorlimit!)*(RND(32000)/32000!))
    PRINT "  rndfactor=";rndfactor!;"  sim%=";sim%
    
    kdist!=0
    FOR i%=1 TO nrxns%
      LONG IF LEFT$(rxnID$(i%),1)="@"
        k!(i%) = kBasis!(i%) * EXP(LOG(rndfactor!)*(2!*RND(32000)/32000 -1))
        kdist! = kdist! + ( LOG(k!(i%)) - LOG(kBasis!(i%)) )^2
      END IF
    NEXT i%
    kdist! = SQR(kdist!)
    FN UpdateLinkedConstants
    
    'Simulate the random set:
    PRINT #3,trial%;","; 
    FOR i%=1 TO nrxns%
      PRINT #3,k!(i%);",";
    NEXT i%
    
    criterion2B! = FN Simulate2!(varspecies%,x1varsp!,x2varsp!,outputspecies%)
    PRINT "  kdist=";kdist!;"  cr2B!=";criterion2B!
    PRINT #3,criterion2B!;",";                    'result from the random set
    
    'Compare the result from basis set (2A) with result from the random set (2B)
    '  and proceed accordingly
    
    LONG IF criterion2B! > criterion2A!           'The random set is better 
      ' and becomes the new basis set
      PRINT "  New Basis set:"
      basisSet%=basisSet%+1
      FN MakeOptRxnFi 
      FOR i%=1 TO nrxns%
        IF e%(i%) >= 0 THEN kBasis!(i%)= k!(i%)
      NEXT i%
      FN UpdateLinkedConstants
      criterion2A! = criterion2B!
      PRINT "cr2A!=";criterion2A!;basisSet%
      PRINT #3,criterion2A!; "," ; "cr2A";",";basisSet%
    XELSE
      PRINT "cr2B!=";criterion2B!
      PRINT #3,criterion2B!; "," ; "cr2B"
    END IF
  NEXT trial%
  
  BEEP
  INPUT "finished RndScan1";q$
END FN

LOCAL FN kScan1
  PRINT "entered kScan1"
  nplots%=1
  
  'Headings:
  CLS
  INPUT #1,outputspecies%,varspecies%,mode3%
  INPUT #1,x1varsp!,x2varsp!,lowerlimit!,kscandev1!,nkscandevs%
  PRINT #3,DATE$;" ";TIME$
  PRINT #3,"Prog$= ";prog$
  PRINT #3,"RxnFi$= ";RxnFi$
  PRINT #3,"mode1$= ";mode1$
  PRINT #3,"outputspecies%=";outputspecies%
  PRINT #3,"varspecies%=";varspecies%
  PRINT #3,"mode3%=";mode3%
  PRINT #3,"x1varsp!=";x1varsp!
  PRINT #3,"x2varsp!=";x2varsp!
  PRINT #3,"last lowerlimit!=";lowerlimit!
  PRINT #3,"kscandev1!=";kscandev1!
  PRINT #3,"nkscandevs%=";nkscandevs%
  
  FOR i%=1 TO nrxns%                              'set kBasis() as center for variation
    kBasis!(i%)=k0!(i%)
    k!(i%)=kBasis!(i%)
  NEXT i%
  FN UpdateLinkedConstants
  
  PRINT #3,"rxn,i1%,kscandev,";
  FOR i%=1 TO nrxns%
    PRINT #3,"rxn";i%;",";
  NEXT i%
  PRINT #3,"x1endpt";",";"stepslope1";",";"x2endpt";",";"stepslope2";",";
  PRINT #3,"x1endpt/x2endpt";",";"x1endpt-x2endpt";",";"result2";",";"kdist!";",";
  PRINT #3,"mode3%";",";"sim%";",";"rndfactor!"
  FOR i%=1 TO nrxns%
    LONG IF LEFT$(rxnID$(i%),1)="@"               'the next k-set itself
      FOR i1% = -nkscandevs% TO nkscandevs%       'do the k-set variations         
        kscandev! = kscandev1!^i1%                'the next k-set variation
        k!(i%)=kBasis!(i%)*kscandev!
        FN UpdateLinkedConstants
        PRINT#3,i%;",";i1%;",";kscandev!;",";
        FOR i2%=1 TO nrxns%
          PRINT #3,k!(i2%);",";
        NEXT i2%
        criterion2! = FN Simulate2!(varspecies%,x1varsp!,x2varsp!,outputspecies%)
        PRINT:PRINT "  rxn";i%;"   i1%,kscandev!= ";i1%;kscandev!;"  criterion2=";criterion2!
        PRINT#3,criterion2!
        key$=INKEY$
        IF key$="." THEN FN DoKey
      NEXT i1%
      k!(i%)=kBasis!(i%)                          'restore the original value of k!(i%)
    END IF
  NEXT i%
END FN                                        

LOCAL FN Triage
  WHILE 1
    PRINT "Started Triage"
    'FN Start
    'q1$=""
    'INPUT "finished 1st pass thru Start";q1$
    'LONG IF q1$="."
    'FN Start
    'INPUT#1,q$
    'EXIT FN
    'END IF
    q$=""
    INPUT#1,mode1$
    PRINT mode1$
    SELECT
      CASE mode1$="set"
        INPUT #1,j%
        INPUT #1,tDelayk!(j%)
      CASE mode1$="rxn"
        INPUT#1,j%
        INPUT#1,k!(j%)
        PRINT j%;i1%(j%);i2%(j%);i3%(j%);o1%(j%);o2%(j%);o3%(j%);e%(j%);k!(j%);
        INPUT;fudge$
        FN Start
      CASE mode1$="species"
        INPUT#1,i%
        j%=ABS(i%)
        INPUT#1,x!(j%)
        FN Start
      CASE mode1$="continue"
        FN Start
      CASE mode1$="scan species" OR mode1$="logscan species"
        OPEN "O",2,"RNout"
        INPUT #1, sp%,x1!,x2!,xd!
        PRINT q$;sp%;x1!;x2!;xd!
        INPUT; fudge$
        FOR j%=0 TO nspecies%
          PRINT #2, USING "###     &"; j%,",",
        NEXT j%
        PRINT #2, fudge%
        xx!=x1!
        IF xx!=0 THEN xx!=.000001
        DO
          LONG IF mode1$="logscan species"
            xx!=xx!*xd!
          XELSE
            xx!=xx!+xd!
          END IF
          LOCATE 1,1,1
          PRINT xx!
          FN InitializeSpecies
          x0!(sp%)=xx!
          IF x0!(sp%)=0 THEN x0!(sp%)=.000001
          x!(sp%)=x0!(sp%)
          FN DoCell
          
          PRINT #2, USING "##.### &"; x0!(sp%),",",
          FOR j%=1 TO nspecies%
            PRINT #2, USING "##.#### &"; x!(j%),",",
          NEXT j%
          PRINT #2, fudge%
        UNTIL xx!>=x2!
        CLOSE #2
      CASE mode1$="scan2"
        FN Scan2
        BEEP
        INPUT; q$
      CASE mode1$="rndscan1"
        FN RndScan1
        CLOSE #3
        BEEP
        INPUT; q$
      CASE mode1$="kscan1"
        FN kScan1
        CLOSE #3
        BEEP
        PRINT "finished kScan1 in Triage"
        INPUT; q$
      CASE ELSE                                   'stop
        CLOSE #1
        INPUT q$
        STOP
    END SELECT
  WEND
END FN


'--------------------
'-------Main--------
BREAK OFF
WINDOW -2,,(20,_gtop)-(_gwidth,_gheight),_dialogPlain,_updateVisRgn
WINDOW  1,,(20,_gtop)-(_gwidth,_gheight),_dialogPlain,_updateVisRgn
IF tout%=1 THEN OPEN "O",3,"RNtout"
WHILE 1
  RxnFi$=FILES$(_fOpen,,"Open a rxnFile",volRefNum%)
  OPEN "I",1,RxnFi$,,volRefNum%
  lin$=""
  LOCATE 0,2
  PRINT "OPENed  ";RxnFi$;"   volRefNum%=";volRefNum%
  WHILE lin$<>"*"
    INPUT#1,lin$
    PRINT lin$
    PRINT #3, lin$
  WEND
  
  INPUT#1,dt1!,nsteps%,nsubsteps%
  PRINT #3,dt1!;",";nsteps%;",";nsubsteps%
  dt!=dt1!/nsubsteps%
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  q$=""
  
  FN ReadRxns
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  INPUT; q$
  IF q$="." THEN STOP
  
  FN ReadSpecies
  PRINT "dt1!,dt!=";dt1!;dt!;"  nsteps%=";nsteps%,"  nsubsteps%=",nsubsteps%
  INPUT; q$
  IF q$="." THEN STOP
  
  PRINT "in Main, ready to start InitializeSpecies"
  FN InitializeSpecies
  PRINT"finished Initializespecies; ready for Triage"
  FN Triage
  WINDOW 1
  INPUT "END RxnNet program?";q$
  IF q$="y" THEN STOP
  CLOSE #1
  IF tout%=1 THEN CLOSE #3
  CLS
  WINDOW -2
  CLS
WEND
END