home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 209_01 / ldhfitr.c < prev    next >
C/C++ Source or Header  |  1990-03-04  |  9KB  |  357 lines
  1. /* LDHFITR.C    VERS:- 01.00  DATE:- 09/26/86  TIME:- 09:39:02 PM */
  2. /*
  3. cp b:$1 -e -s0
  4. z80asm $1.bb
  5. slrnk b:$1, b:simpmain, b:simplib1/s, b:simplib0/s, b:$1/n/e
  6. */
  7. /* 
  8. Description:
  9.  
  10. VERSION OF XXXXFITn, for:
  11.  
  12. nonlinear least squares fit by simplex minimizaton (Nelder-Mead algorithm).
  13.  
  14. analysis of lactate dehydrogenase initial rate kinetic data:
  15.     two substrate -- two product mechanism,
  16.     with inhibition by first product,
  17.     and with reactant and product abortive complexes;
  18.  
  19. see LDHFITR documentation for description of rate law.
  20.  
  21. compile with:
  22.     simpmain, simplib1, simplib0
  23.     floating point functions
  24.  
  25. the precision of the floating point functions must be high, ie 14 or more
  26.     decimal digits; eight digit precision is not enough for the 
  27.     matrix operations.
  28.  
  29. contents of the xxxxfitn module  =  
  30. declarations and routines for simplex fitting special to function to be fit:
  31.  
  32.     function for calculation of dependent variable and
  33.         weighted sum of residuals squared        = func()
  34.  
  35.     print of <data> records                    = fdatprint()
  36.                                 = fpointprint()
  37.  
  38.     other special displays                    = fspecial()
  39.  
  40.     all external declarations, except for <data>, which is 
  41.         defined as dummy storage in SIMPLIB1 
  42.         or its modification
  43.  
  44. By J.A. Rupley, Tucson, Arizona
  45. Coded for ECO C compiler, version 3.40
  46. */
  47.  
  48.     /* page eject */
  49.  
  50. /*
  51.      To make more readable the coding in <func()> of the model 
  52. equation to be fit to the data: 
  53.  
  54. (1) use mnemonic member names in declaring <struct dat> in 
  55. XXXXFITn; 
  56.  
  57. (2) declare a dummy structure, <struct pnamestruct>, that is 
  58. entirely equivalent to the structure that holds the parameter 
  59. values, <pstruct>, but that has mnemonic member names; the 
  60. mnemonic dummy structure then can be used with the <pstruct> 
  61. address passed as a parameter to <func()>.
  62.  
  63.  
  64.      Change in the model being fit should not require recoding 
  65. and recompilation of <read_data()> or any other routines except 
  66. those of XXXXFITn; of course, change in the model requires change 
  67. of <func()>, <fdatprint()>, and of the declarations of <struct dat> 
  68. and <struct pnamestruct> in XXXXFITn.
  69.  
  70. */
  71.  
  72. /* page eject */
  73.  
  74. #include <stdio.h>
  75. #include <ctrlcnst.h>
  76.  
  77. #define NVERT        8    /* set equal to 1 + number of parameters 
  78.                 used in <func()> */
  79.  
  80.  
  81. #define NPARM        10    /* do NOT change this define */
  82.  
  83.  
  84.  
  85.                 /* EXTERNALS AND STRUCTURES */
  86.  
  87.                 /* do NOT change any structure except <dat>*/
  88.  
  89.                 /* the structure <dat> MUST be changed 
  90.                 to accord with the requirements of 
  91.                 <func()> and <fdatprint()>; */
  92.                 /* all members of <struct dat> MUST be of 
  93.                 type double. */
  94.                 /* as written the program accomodates 
  95.                 more than 100 six-member data points */
  96.                 /* see the header above for more description
  97.                 of the declaration of <struct dat> in XXXXFITn
  98.                 and the definition of the aggregate <data> in
  99.                 SIMPLIB1 */
  100.  
  101.                 /* the structure <pnamestruct> is equivalent 
  102.                 to <pstruct>, which holds the parameter values;
  103.                 use of <pnamestruct> allows mnemonic access 
  104.                 to the contents of <pstruct> */
  105.  
  106. struct dat {            
  107.     double y ;        
  108.     double yc ;
  109.     double w ;
  110.     double a ;
  111.     double b ;
  112.     double p ;
  113. } ;
  114.  
  115. struct pstruct {
  116.     double val ;
  117.     double parm[NPARM] ;
  118. } ;
  119.  
  120. struct pnamestruct {
  121.     double val ;
  122.     double vmax ;
  123.     double kma ;
  124.     double kmb ;
  125.     double kmab ;
  126.     double kmq_kmpq ;
  127.     double kip_1 ;
  128.     double kib_k3k4 ;
  129.     double dummy[3] ;
  130. } ;
  131.  
  132. struct qstruct {
  133.     int parmndx[NPARM] ;
  134.     double q[NPARM] ;
  135.     double yplus[NPARM] ;
  136.     double yminus[NPARM] ;
  137.     double a[NPARM] ;
  138.     double bdiag[NPARM] ;
  139.     double inv_bdiag[NPARM] ;
  140.     double std_dev[NPARM] ;
  141. } ;
  142.  
  143.  
  144. struct pstruct p[NVERT] ;
  145.  
  146. struct pstruct pcent ;
  147.  
  148. struct pstruct *p_p[NVERT] ;
  149.  
  150. struct pstruct pmin ;
  151.  
  152. struct qstruct q ;
  153.  
  154.  
  155. double qmat[NPARM][NPARM] ;
  156.  
  157. double mean_func, rms_func, test, rms_data ;
  158.  
  159. double yzero, ymin, ypmin, mse ;
  160.  
  161. double quad_test, exit_test ;
  162.  
  163. int iter, maxiter, nparm, nvert, ndata, maxquad_skip, prt_cycle, ndatval ;
  164.  
  165. int nfree ;
  166.  
  167. char title[80] ;
  168.  
  169. /* page eject */
  170.  
  171. /* FUNC
  172.     CALCULATION OF LEAST SQUARES FUNCTION
  173.     CODED ACCORDING TO MODEL BEING FIT
  174.     IT SHOULD BE EFFICIENT
  175.     DURING THE FIT, TIME IS MOSTLY SPENT HERE
  176.     (THE OVERHEAD ELSEWHERE IS ONLY SEVERAL SECONDS PER CYCLE)
  177. */
  178.  
  179. int func(g)
  180. struct pstruct *g ;
  181. {
  182.     void printf() ;
  183.     register int n ;
  184.  
  185.     struct pnamestruct *pnam ;
  186.  
  187.     extern struct dat     data[] ;
  188.     extern int         nparm ;
  189.     extern int         ndata ;
  190.  
  191.     for (n = 0; n < nparm; n++)
  192.         if (g->parm[n] <= 0) {
  193.             g->val = 1.E38 ;
  194.             printf("function error\n") ;
  195.             return (ERROR) ;
  196.         }
  197.  
  198.     pnam = g ;
  199.     pnam->val = 0 ;
  200.     for (n = 0; n < ndata; n++) {
  201.         data[n].yc = 
  202.             pnam->vmax / 
  203.             (
  204.                   1 + pnam->kmq_kmpq * data[n].p 
  205.               + pnam->kma / data[n].a 
  206.               + pnam->kmb / data[n].b 
  207.                 * (1 + pnam->kmq_kmpq * data[n].p)
  208.                 * (1 + pnam->kip_1 * data[n].p)
  209.               + pnam->kmab / data[n].a / data[n].b 
  210.                 * (1 + pnam->kmq_kmpq * data[n].p) 
  211.               + pnam->kib_k3k4 * data[n].b
  212.             ) ;
  213.         pnam->val = pnam->val + (data[n].y - data[n].yc) 
  214.                 * (data[n].y - data[n].yc) 
  215.                 * data[n].w 
  216.                 * data[n].w ;
  217.     }
  218.     return (OK) ;
  219. }                              /* END OF FUNC                    */
  220.  
  221. /* page eject */
  222.  
  223. /* FDATPRINT
  224.     PRINT DATA AND COMPARE WITH CALCULATED VALUES
  225.     CODED ACCORDING TO MODEL AND DATA
  226. */
  227.  
  228. void fdatprint(fptr)
  229. FILE *fptr ;
  230. {
  231.     register int j ;
  232.  
  233.     void fprintf() ;
  234.  
  235.     extern int         iter, ndata, maxiter ;
  236.     extern struct dat     data[] ;
  237.     extern char        title[] ;
  238.  
  239.     fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter) ;
  240.     fprintf(fptr, 
  241. "        yobs   ycalc   del-y  weight       a         b         p\n") ;
  242.     for (j = 0; j < ndata; j++) 
  243.         fprintf(fptr, 
  244.             "%4d %7.3f %7.3f %7.3f %7.3f %10.5f %10.5f %10.5f\n", 
  245.             (j+1), data[j].y, data[j].yc, (data[j].y - data[j].yc),
  246.             data[j].w, data[j].a, data[j].b, data[j].p) ;
  247.  
  248.     if (maxiter != 0)
  249.         fpointprint(fptr) ;
  250. }                /* END OF FDATPRINT            */
  251.  
  252. /* page eject */
  253.  
  254. /* FPOINTPRINT
  255.     PRINT POINTS FOR CONSTRUCTION OF PRIMARY AND SECONDARY PLOTS
  256.     CODED ACCORDING TO MODEL AND DATA
  257. */
  258.  
  259. void fpointprint(fptr)
  260. FILE *fptr ;
  261. {
  262.     void fprintf() ;
  263.     extern struct pstruct    pcent ;
  264.     struct pnamestruct    *pnam ;
  265.     extern char         title[] ;
  266.     extern int        iter ;
  267.     double v,k1,k2,k3,k4,k5,vint,v1,v2,s1,s2 ;
  268.  
  269.     pnam = &pcent ;
  270.     v = pnam->vmax ;
  271.     k1 = pnam->kma ;
  272.     k2 = pnam->kmb ;
  273.     k3 = pnam->kmab ;
  274.     k4 = pnam->kmq_kmpq ;
  275.     k5 = pnam->kip_1 ;
  276.     vint = (1 - k1 * k2/k3)/v ;
  277.     v1 = (1 + k1/.035)/v ;
  278.     v2 = v1 + k4 * 150/v ;
  279.     s1 = (k2 + k3/.035)/v ;
  280.     s2 = (s1 + k2 * k5 * 150/v) * (1 + k4 * 150) ;
  281.  
  282.     fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter) ;
  283.     fprintf(fptr, 
  284. "For primary plot at fixed b = [Pyruvate]:\n") ;
  285.     fprintf(fptr, 
  286. "intersection (1/a = 1/[NADH], 1/V)      %10.2f%10.2f\n",-k2/k3,vint) ;
  287.     fprintf(fptr, 
  288. "1/V intercepts at b = (.7,.35,.22,.15)      %8.2f%8.2f%8.2f%8.2f\n\n",
  289.         (1+k2/.7)/v,(1+k2/.35)/v,(1+k2/.22)/v,(1+k2/.15)/v) ;
  290.  
  291.     fprintf(fptr, 
  292. "For primary plot at fixed a = [NADH]:\n") ;
  293.     fprintf(fptr, 
  294. "intersection (1/b = 1/[Pyruvate], 1/V)  %10.2f%10.2f\n",
  295.         -k1/k3,vint) ;
  296.     fprintf(fptr, 
  297. "1/V intercepts at a = (.035,.022,.015,.011) %8.2f%8.2f%8.2f%8.2f\n\n",
  298.         (1+k1/.035)/v,(1+k1/.022)/v,(1+k1/.015)/v,(1+k1/.011)/v) ;
  299.  
  300.     fprintf(fptr, 
  301. "For secondary plot:\n") ;
  302.     fprintf(fptr, 
  303. "1/a = 1/[NADH] and 1/b = 1/[Pyruvate] intercepts %9.2f%10.2f\n",-1/k1,-1/k2);
  304.     fprintf(fptr, 
  305. "1/V intercept     \t\t\t\t%10.2f\n\n",1/v) ;
  306.  
  307.     fprintf(fptr, "For lactate inhibition, at a = [NADH] = 0.035\n") ;
  308.     fprintf(fptr, "1/V values at 1/b = 1/[Pyruvate] = 0 and 5.0\n");
  309.     fprintf(fptr, "\t\t\t- Lactate%7.2f%10.2f\n", v1, (v1 + 5 * s1)) ;
  310.     fprintf(fptr, "\t\t\t+ Lactate%7.2f%10.2f\n\n",v2, (v2 + 5 * s2)) ;
  311.  
  312.     fprintf(fptr, 
  313. "R(+/-,intercept), R(+/-,slope)%10.2f%10.2f\n",v2/v1,s2/s1);
  314.  
  315.     fprintf(fptr, 
  316. "\n\n\n\nValues of the kinetic constants\n") ;
  317.     fprintf(fptr, "parm(0) = V        %17.11e\n",v) ;
  318.     fprintf(fptr, "parm(1) = KmA      %17.11e\n",k1) ;
  319.     fprintf(fptr, "parm(2) = KmB      %17.11e\n",k2) ;
  320.     fprintf(fptr, "parm(3) = KmAB     %17.11e\n",k3) ;
  321.     fprintf(fptr, "parm(4) = KmQ/KmPQ %17.11e\n",k4) ;
  322.     fprintf(fptr, "parm(5) = 1/K(I,P) %17.11e\n",k5) ;
  323.  
  324. }                /* END OF FPOINTPRINT            */
  325.  
  326. /* page eject */
  327.  
  328. /* FSPECIAL
  329.     DISPLAY ADDITIONAL INFORMATION DURING
  330.     TRACKING OF MINIMIZATION, IN SIMPFIT()
  331. */
  332.  
  333. void fspecial(fptr)
  334. FILE *fptr ;
  335. {
  336.     register j ;
  337.  
  338.     void fprintf() ;
  339.  
  340.     extern int         maxiter ;
  341.     extern double        ypmin, yzero, quad_test ;
  342.  
  343.     if (ypmin > 0.99E38) {
  344.         fprintf(fptr, "y-pmin error") ;
  345.         for (j = 0; j < nvert; j++)
  346.             if (pmin.parm[j] < 0)
  347.                 fprintf(fptr, " (parm(%d) < 0)", j) ;
  348.     }
  349.     else if (ypmin > yzero)
  350.         fprintf(fptr, "y-pmin > yzero") ;
  351.  
  352.     fprintf(fptr, "    maxiter = %d  quad_test = %7.2e\n", 
  353.         maxiter, quad_test) ;
  354.  
  355.     return ;
  356. }                /* END OF FSPECIAL            */
  357.