home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 209_01 / ldhfitr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-03-05  |  8.7 KB  |  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.  
  53. equation to be fit to the data: 
  54.  
  55. (1) use mnemonic member names in declaring <struct dat> in 
  56.  
  57. XXXXFITn; 
  58.  
  59. (2) declare a dummy structure, <struct pnamestruct>, that is 
  60.  
  61. entirely equivalent to the structure that holds the parameter 
  62.  
  63. values, <pstruct>, but that has mnemonic member names; the 
  64.  
  65. mnemonic dummy structure then can be used with the <pstruct> 
  66.  
  67. address passed as a parameter to <func()>.
  68.  
  69.  
  70.      Change in the model being fit should not require recoding 
  71.  
  72. and recompilation of <read_data()> or any other routines except 
  73.  
  74. those of XXXXFITn; of course, change in the model requires change 
  75.  
  76. of <func()>, <fdatprint()>, and of the declarations of <struct dat> 
  77.  
  78. and <struct pnamestruct> in XXXXFITn.
  79.  
  80. */
  81.  
  82. /* page eject */
  83.  
  84. #include <stdio.h>
  85. #include <ctrlcnst.h>
  86.  
  87. #define NVERT        8    /* set equal to 1 + number of parameters 
  88.                 used in <func()> */
  89.  
  90.  
  91. #define NPARM        10    /* do NOT change this define */
  92.  
  93.  
  94.  
  95.                 /* EXTERNALS AND STRUCTURES */
  96.  
  97.                 /* do NOT change any structure except <dat>*/
  98.  
  99.                 /* the structure <dat> MUST be changed 
  100.                 to accord with the requirements of 
  101.                 <func()> and <fdatprint()>; */
  102.                 /* all members of <struct dat> MUST be of 
  103.                 type double. */
  104.                 /* as written the program accomodates 
  105.                 more than 100 six-member data points */
  106.                 /* see the header above for more description
  107.                 of the declaration of <struct dat> in XXXXFITn
  108.                 and the definition of the aggregate <data> in
  109.                 SIMPLIB1 */
  110.  
  111.                 /* the structure <pnamestruct> is equivalent 
  112.                 to <pstruct>, which holds the parameter values;
  113.                 use of <pnamestruct> allows mnemonic access 
  114.                 to the contents of <pstruct> */
  115.  
  116. struct dat {            
  117.     double y ;        
  118.     double yc ;
  119.     double w ;
  120.     double a ;
  121.     double b ;
  122.     double p ;
  123. } ;
  124.  
  125. struct pstruct {
  126.     double val ;
  127.     double parm[NPARM] ;
  128. } ;
  129.  
  130. struct pnamestruct {
  131.     double val ;
  132.     double vmax ;
  133.     double kma ;
  134.     double kmb ;
  135.     double kmab ;
  136.     double kmq_kmpq ;
  137.  
  138.     double kip_1 ;
  139.     double kib_k3k4 ;
  140.     double dummy[3] ;
  141. } ;
  142.  
  143. struct qstruct {
  144.     int parmndx[NPARM] ;
  145.     double q[NPARM] ;
  146.     double yplus[NPARM] ;
  147.     double yminus[NPARM] ;
  148.     double a[NPARM] ;
  149.     double bdiag[NPARM] ;
  150.     double inv_bdiag[NPARM] ;
  151.     double std_dev[NPARM] ;
  152. } ;
  153.  
  154.  
  155. struct pstruct p[NVERT] ;
  156.  
  157. struct pstruct pcent ;
  158.  
  159. struct pstruct *p_p[NVERT] ;
  160.  
  161. struct pstruct pmin ;
  162.  
  163. struct qstruct q ;
  164.  
  165.  
  166. double qmat[NPARM][NPARM] ;
  167.  
  168. double mean_func, rms_func, test, rms_data ;
  169.  
  170. double yzero, ymin, ypmin, mse ;
  171.  
  172. double quad_test, exit_test ;
  173.  
  174. int iter, maxiter, nparm, nvert, ndata, maxquad_skip, prt_cycle, ndatval ;
  175.  
  176. int nfree ;
  177.  
  178. char title[80] ;
  179.  
  180. /* page eject */
  181.  
  182. /* FUNC
  183.     CALCULATION OF LEAST SQUARES FUNCTION
  184.     CODED ACCORDING TO MODEL BEING FIT
  185.     IT SHOULD BE EFFICIENT
  186.     DURING THE FIT, TIME IS MOSTLY SPENT HERE
  187.     (THE OVERHEAD ELSEWHERE IS ONLY SEVERAL SECONDS PER CYCLE)
  188. */
  189.  
  190. int func(g)
  191. struct pstruct *g ;
  192. {
  193.     void printf() ;
  194.     register int n ;
  195.  
  196.     struct pnamestruct *pnam ;
  197.  
  198.     extern struct dat     data[] ;
  199.     extern int         nparm ;
  200.     extern int         ndata ;
  201.  
  202.     for (n = 0; n < nparm; n++)
  203.         if (g->parm[n] <= 0) {
  204.             g->val = 1.E38 ;
  205.             printf("function error\n") ;
  206.             return (ERROR) ;
  207.         }
  208.  
  209.     pnam = g ;
  210.     pnam->val = 0 ;
  211.     for (n = 0; n < ndata; n++) {
  212.         data[n].yc = 
  213.             pnam->vmax / 
  214.             (
  215.                   1 + pnam->kmq_kmpq * data[n].p 
  216.               + pnam->kma / data[n].a 
  217.               + pnam->kmb / data[n].b 
  218.                 * (1 + pnam->kmq_kmpq * data[n].p)
  219.                 * (1 + pnam->kip_1 * data[n].p)
  220.               + pnam->kmab / data[n].a / data[n].b 
  221.                 * (1 + pnam->kmq_kmpq * data[n].p) 
  222.               + pnam->kib_k3k4 * data[n].b
  223.             ) ;
  224.         pnam->val = pnam->val + (data[n].y - data[n].yc) 
  225.                 * (data[n].y - data[n].yc) 
  226.                 * data[n].w 
  227.                 * data[n].w ;
  228.     }
  229.     return (OK) ;
  230. }                              /* END OF FUNC                    */
  231.  
  232. /* page eject */
  233.  
  234. /* FDATPRINT
  235.     PRINT DATA AND COMPARE WITH CALCULATED VALUES
  236.     CODED ACCORDING TO MODEL AND DATA
  237. */
  238.  
  239. void fdatprint(fptr)
  240. FILE *fptr ;
  241. {
  242.     register int j ;
  243.  
  244.     void fprintf() ;
  245.  
  246.     extern int         iter, ndata, maxiter ;
  247.     extern struct dat     data[] ;
  248.     extern char        title[] ;
  249.  
  250.     fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter) ;
  251.     fprintf(fptr, 
  252. "        yobs   ycalc   del-y  weight       a         b         p\n") ;
  253.     for (j = 0; j < ndata; j++) 
  254.         fprintf(fptr, 
  255.             "%4d %7.3f %7.3f %7.3f %7.3f %10.5f %10.5f %10.5f\n", 
  256.             (j+1), data[j].y, data[j].yc, (data[j].y - data[j].yc),
  257.             data[j].w, data[j].a, data[j].b, data[j].p) ;
  258.  
  259.     if (maxiter != 0)
  260.         fpointprint(fptr) ;
  261. }                /* END OF FDATPRINT            */
  262.  
  263. /* page eject */
  264.  
  265. /* FPOINTPRINT
  266.     PRINT POINTS FOR CONSTRUCTION OF PRIMARY AND SECONDARY PLOTS
  267.     CODED ACCORDING TO MODEL AND DATA
  268. */
  269.  
  270. void fpointprint(fptr)
  271. FILE *fptr ;
  272. {
  273.     void fprintf() ;
  274.     extern struct pstruct    pcent ;
  275.     struct pnamestruct    *pnam ;
  276.     extern char         title[] ;
  277.     extern int        iter ;
  278.     double v,k1,k2,k3,k4,k5,vint,v1,v2,s1,s2 ;
  279.  
  280.     pnam = &pcent ;
  281.     v = pnam->vmax ;
  282.     k1 = pnam->kma ;
  283.     k2 = pnam->kmb ;
  284.     k3 = pnam->kmab ;
  285.     k4 = pnam->kmq_kmpq ;
  286.     k5 = pnam->kip_1 ;
  287.     vint = (1 - k1 * k2/k3)/v ;
  288.     v1 = (1 + k1/.035)/v ;
  289.     v2 = v1 + k4 * 150/v ;
  290.     s1 = (k2 + k3/.035)/v ;
  291.     s2 = (s1 + k2 * k5 * 150/v) * (1 + k4 * 150) ;
  292.  
  293.     fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter) ;
  294.     fprintf(fptr, 
  295. "For primary plot at fixed b = [Pyruvate]:\n") ;
  296.     fprintf(fptr, 
  297. "intersection (1/a = 1/[NADH], 1/V)      %10.2f%10.2f\n",-k2/k3,vint) ;
  298.     fprintf(fptr, 
  299. "1/V intercepts at b = (.7,.35,.22,.15)      %8.2f%8.2f%8.2f%8.2f\n\n",
  300.         (1+k2/.7)/v,(1+k2/.35)/v,(1+k2/.22)/v,(1+k2/.15)/v) ;
  301.  
  302.     fprintf(fptr, 
  303. "For primary plot at fixed a = [NADH]:\n") ;
  304.     fprintf(fptr, 
  305. "intersection (1/b = 1/[Pyruvate], 1/V)  %10.2f%10.2f\n",
  306.         -k1/k3,vint) ;
  307.     fprintf(fptr, 
  308. "1/V intercepts at a = (.035,.022,.015,.011) %8.2f%8.2f%8.2f%8.2f\n\n",
  309.         (1+k1/.035)/v,(1+k1/.022)/v,(1+k1/.015)/v,(1+k1/.011)/v) ;
  310.  
  311.     fprintf(fptr, 
  312. "For secondary plot:\n") ;
  313.     fprintf(fptr, 
  314. "1/a = 1/[NADH] and 1/b = 1/[Pyruvate] intercepts %9.2f%10.2f\n",-1/k1,-1/k2);
  315.     fprintf(fptr, 
  316. "1/V intercept     \t\t\t\t%10.2f\n\n",1/v) ;
  317.  
  318.  
  319.     fprintf(fptr, "For lactate inhibition, at a = [NADH] = 0.035\n") ;
  320.     fprintf(fptr, "1/V values at 1/b = 1/[Pyruvate] = 0 and 5.0\n");
  321.     fprintf(fptr, "\t\t\t- Lactate%7.2f%10.2f\n", v1, (v1 + 5 * s1)) ;
  322.     fprintf(fptr, "\t\t\t+ Lactate%7.2f%10.2f\n\n",v2, (v2 + 5 * s2)) ;
  323.  
  324.     fprintf(fptr, 
  325. "R(+/-,intercept), R(+/-,slope)%10.2f%10.2f\n",v2/v1,s2/s1);
  326.  
  327.     fprintf(fptr, 
  328. "\n\n\n\nValues of the kinetic constants\n") ;
  329.     fprintf(fptr, "parm(0) = V        %17.11e\n",v) ;
  330.     fprintf(fptr, "parm(1) = KmA      %17.11e\n",k1) ;
  331.     fprintf(fptr, "parm(2) = KmB      %17.11e\n",k2) ;
  332.     fprintf(fptr, "parm(3) = KmAB     %17.11e\n",k3) ;
  333.     fprintf(fptr, "parm(4) = KmQ/KmPQ %17.11e\n",k4) ;
  334.     fprintf(fptr, "parm(5) = 1/K(I,P) %17.11e\n",k5) ;
  335.  
  336. }                /* END OF FPOINTPRINT            */
  337.  
  338. /* page eject */
  339.  
  340. /* FSPECIAL
  341.     DISPLAY ADDITIONAL INFORMATION DURING
  342.     TRACKING OF MINIMIZATION, IN SIMPFIT()
  343. */
  344.  
  345. void fspecial(fptr)
  346. FILE *fptr ;
  347. {
  348.     register j ;
  349.  
  350.     void fprintf() ;
  351.  
  352.     extern int         maxiter ;
  353.     extern double        ypmin, yzero, quad_test ;
  354.  
  355.     if (ypmin > 0.99E38) {
  356.         fprintf(fptr, "y-pmin error") ;
  357.         for (j = 0; j < nvert; j++)
  358.             if (pmin.parm[j] < 0)
  359.                 fprintf(fptr, " (parm(%d) < 0)", j) ;
  360.     }
  361.     else if (ypmin > yzero)
  362.         fprintf(fptr, "y-pmin > yzero") ;
  363.  
  364.     fprintf(fptr, "    maxiter = %d  quad_test = %7.2e\n", 
  365.         maxiter, quad_test) ;
  366.  
  367.     return ;
  368. }                /* END OF FSPECIAL            */
  369.