chempare Translation Listing



 FORTRAN CALCULUS (VERSION 2.0) ------------------------- 06/04/2008   16:53:44
 ####  LABEL ---+---------+---------+---------+---------+---------+---------+--
~
    1        GLOBAL ALL
    2        PROBLEM CHEMPARE !   Chemical Kinetics Parameter Estimation
            !----------------------------------------------------------------|
            ! Traditional chemical kinetics parameter estimation problem     |
            ! known in the pharmaceutical industry as pharmacokinetics.      |
            ! Unknown parameters of a set of differential equations modeling |
            ! the ingestion of chemicals in the blood (say) are estimated by |
            ! nesting the "chemical trajectories simulation" inside a non-   |
            ! linear correlator, to find the trajectories that match sampled |
            ! measurements, resulting in the solution of a multipoint        |
            ! boundary-value problem of a set of differential equations.     |
            !----------------------------------------------------------------|
    3          DIMENSION TM(6),CM(6),AM(6),CFITS(6),AFITS(6)
      
              !--------------------- Measurement Data -----------------------|
    4          DATA TM/10,20,30,40,50,60/  ! Time points for measurements
    5          DATA CM/0.419,0.563,0.629,0.666,0.689,0.708/ ! C measurements
    6          DATA AM/0.483,0.281,0.191,0.134,0.097,0.065/ ! A measurements
      
              !--------------Initial conditions of the ODEs -----------------|
    7          DATA NM,A0,B0,C0,D0,DT/6,1.0,1.03,0.0,0.0,2.0/
    8          DATA B1,B2/0.015,0.015/            ! Optimization step bounds
   10>         P1=0.01 : P2=0.05                  ! Parameter starting guesses
      
              !---------- MetaCalculus: Initialize the ODE Solution ---------|
   11>         INITIATE ATHENA; FOR REACTION;
           &    EQUATIONS DADT/A,DBDT/B,DCDT/C,DDDT/D;
           &    OF T; STEP DT; TO TF
      
              !-------------- MetaCalculus: Start the Optimizer -------------|
   12>         FIND P1,P2; IN CURFIT;
           &    BY AJAX(SET); TO MATCH CFITS,AFITS
   13        END
~
      
   14        CONTROLLER SET(AJAX) ! FIND control parameters for the HERA Solver
   15          DETAIL=1    ! Detailed iteration report for every iteration
   16          DETOUT=0    ! Send iteration reports to DETAIL stream
   17          SUMOUT=1    ! Send iteration summary to SUMMARY stream
   18          DAMP=0
   19        END
~
      
            !----------- Simulation Model of Chemical Kinetic Curves --------|
   20        MODEL CURFIT
   21          CHARFUN FCSINT*2
   22          CHARACTER*2 NI
   24>         DATA IT/0/ : NI=FCSINT(IT)
   25*         @AXES(’ITER’//NI,’ITERATION ’//NI)
   31>         A=A0 : B=B0 : C=C0 : D=D0 : T=0 : TF=0
   33>         ERROR=0 : I=1
   34*         DO WHILE (I.LE.NM)
   35            TF=TF+DT
      
                !---------- MetaCalculus: Generate Kinetic Curves -----------|
   36>           INTEGRATE REACTION; BY ATHENA
   37            @CURVES(’ITER’//NI)                         ! Plot the curves
   38            IF(TF.EQ.TM(I)) THEN
                  !--- Compute error between curves and measured data -----|
   39              AFITS(I)=AM(I)-A
   40              CFITS(I)=CM(I)-C
   41              @MEASURED(’ITER’//NI,AM(I),CM(I),TM(I)) ! Plot measurement
   42              I=I+1
   43*           ENDIF
   44          END DO
              !-------------- Show graph for this iteration -----------------|
   45          @SHOW(’ITER’//NI’,’Chemical profile of iteration’//NI)
   46          IT=IT+1
   47*       END
~
      
   48        MODEL REACTION ! Rates of reaction differential equations
   49          DCDT=P1*A*B
   50          DADT=-(DCDT+(.01+P2)*A
   51          DBDT=-(DCDT+.05*B*D)
   52          DDDT=DBDT-DADT
   53        END
~
      
      
       !===================== GRAPHICS DISPLAY ROUTINES ======================|
   54        PROCEDURE AXES(GNAME,TITLE)
            !----------------------------------------------------------------|
            ! This procedure sets up a graph of chemical kinetic curves      |
            ! corresponding to each iteration of the optimizer, showing the  |
            ! plotted measurement data. Its purpose is to show the progress  |
            ! of the parameter-estimation process iteration by iteration.    |
            !----------------------------------------------------------------|
   55         CHARACTER*(*) GNAME,TITLE
   56         @GRAFIL(GNAME,’SUM’,’IMAGE’,0,0) ! Create SUMMARY-stream graph
   57         @FONT(GNAME,’COMPLEX’,’STANDARD’,0,10) ! Character size & color
   58         @AXNAME(GNAME,’X’,’TIME (MINUTES)’,’CENT’,30,36,10)
   59         @AXNAME(GNAME,’Y’,’CONCENTRATION’,’CENT’,30,36,10)
   60         @XYPLOT(GNAME,’RECT’,0.0,60.0,0.0,10.0,0,1.0,0.0,0.1,1.0,0,0,0,0)
   61         @SETUP(GNAME,’AM’,LN_NONE,CL_RED,SY_BOX,0)
   62         @SETUP(GNAME,’CM’,LN_NONE,CL_MAGENTA,SY_DELTA,0)
   63         @SETUP(GNAME,’A ’,LN_SOL,CL_RED,SY_NONE,1)
   64         @SETUP(GNAME,’C ’,LN_DASH,CL_MAGENTA,SY_NONE,1)
   65         @SETUP(GNAME,’B ’,LN_CHNDSH,CL_GREEN,SY_NONE,1)
   66         @SETUP(GNAME,’D ’,LN_CHNDOT,CL_ORANGE,SY_NONE,1)
             !------------------- Chart Title (FORE) ------------------------|
   67         @HEAD(GNAME,’CHEMICAL PARAMETERS ESTIMATION’,50,50,10)
   68         @HEAD(GNAME,TITLE,50,100,10)
   69         @MESSAGE(GNAME,’P1= ’,20.0,.9,4) ! (Blue)
   70         @NUMBER(GNAME,P1,4,999.,999.,10,0,2) ! (Red)
   71         @MESSAGE(GNAME,’P2= ’,20.0,.85,4) ! (Blue)
   72         @NUMBER(GNAME,P2,4,999.,999.,10,0,2) ! (Red)
   73        END
~
      
   74        MODEL CURVES(GNAME)  ! Generate points on integral curves
   75          CHARACTER*(*) GNAME
   76          @CURVE(GNAME,’A ’,T,A)
   77          @CURVE(GNAME,’C ’,T,C)
   78          @CURVE(GNAME,’B ’,T,B)
   79          @CURVE(GNAME,’D ’,T,D)
   80        END
~
      
   81        MODEL MEASURED(GNAME,AP,CP,TP) ! Generate measured points
   82          CHARACTER*(*) GNAME
   83          @POINT(GNAME,’AM’,TP,AP)
   84          @POINT(GNAME,’CM’,TP,CP)
   85        END
~