home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 209_01 / simplib0.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-03-05  |  18.6 KB  |  865 lines
  1. /* SIMPLIB0.C    VERS:- 01.00  DATE:- 09/26/86  TIME:- 09:39:24 PM */
  2. /* 
  3. Description:
  4.  
  5. functions for simplex fitting:
  6.     simplex fitting routine            = simpfit()
  7.     quadratic fit for standard deviations    = simpdev()
  8.     supporting functions
  9.  
  10. By J.A. Rupley, Tucson, Arizona
  11. Coded for ECO C compiler, version 3.40
  12. */
  13.  
  14.  
  15.  
  16.  
  17.  
  18. #include <stdio.h>
  19. #include <ctrlcnst.h>
  20.  
  21. #define NPARM        10
  22.  
  23.                 /* STRUCTURES */
  24.  
  25. struct pstruct {
  26.     double val ;
  27.     double parm[NPARM] ;
  28. }  ;
  29.  
  30. struct qstruct {
  31.     int parmndx[NPARM] ;
  32.     double q[NPARM] ;
  33.     double yplus[NPARM] ;
  34.     double yminus[NPARM] ;
  35.     double a[NPARM] ;
  36.     double bdiag[NPARM] ;
  37.     double inv_bdiag[NPARM] ;
  38.     double std_dev[NPARM] ;
  39. }  ;
  40.  
  41. /* page eject */
  42.  
  43. /* SIMPFIT
  44.     SIMPLEX MINIMIZATION BY USE OF NELDER-MEAD ALGORITHM,
  45.     FOR MODEL DEFINED IN  <FUNC()> .
  46. */
  47.  
  48.                     /* defines for simpfit */
  49. #define ILLEGAL        0
  50. #define REFLECT        1
  51. #define STARCONTRACT    2
  52. #define HIGHCONTRACT    3
  53. #define SHRINK        4
  54. #define FAILREFLECT    5
  55. #define STAREXPAND      6
  56.  
  57.  
  58. int simpfit(fptr)
  59. FILE *fptr ;
  60. {
  61.     register int j, i ;
  62.  
  63.     int opcode ;
  64.     int c ;
  65.     int test_break ;
  66.     int nlow, nhigh  ;
  67.  
  68.     struct pstruct *p_best, pbar, pstar, p2star ;
  69.     
  70.     void bsort(), fprintf() ;
  71.     int getchar() ;
  72.     double sqrt() ;
  73.  
  74.     int pvalcomp(), p_swap() ;
  75.     int centroid(), func(), ptrial() ;
  76.     void fpprint(), pcopy() ;
  77.  
  78.     extern struct pstruct     p[], pcent, *p_p[] ;
  79.     extern double       exit_test, mean_func, rms_func, test, rms_data;
  80.     extern int         iter, maxiter, nparm, nvert, ndata, nfree ;
  81.  
  82.  
  83.                 /* expansion-contraction coefficients */
  84.     double alpha ;
  85.     double beta ;
  86.     double gamma ;
  87.  
  88.                     /* descriptions of expn-cntrctn operations */
  89.     static char *opname[] = {
  90.         "illegal",
  91.         "reflect",
  92.         "reflect & contract",
  93.         "contract",
  94.         "shrink on lowest vertex",
  95.         "reflect after fail to expand",
  96.  
  97.         "reflect and expand"
  98.     }  ;
  99.     
  100.     nlow =        0 ;
  101.     nhigh =        nvert - 1 ;
  102.  
  103.  
  104.  
  105.  
  106.                 /* initialization:
  107.                     coefficients
  108.                         pointers */
  109.     alpha =    1 ;
  110.     beta  = .5 ;
  111.     gamma =  2 ;
  112.  
  113.     for (j = 0; j < nvert; j++)
  114.         p_p[j] = &p[j] ;
  115.  
  116.     fprintf(fptr, 
  117.         "\n\ntype ^X to exit loop, ^S to pause\n\n") ;
  118.  
  119.  
  120.  
  121.  
  122.                 /* MAIN LOOP OF MINIMIZATION */
  123.     while (++iter) {
  124.  
  125.  
  126.                 /* ascending sort of pointers p_p 
  127.                 to struct array p */
  128.         bsort(nvert, pvalcomp, p_swap) ;
  129.  
  130.  
  131.                 /* form reduced simplex, pbar, 
  132.                 for (nvert - 1) vertices 
  133.                 ie without high vertex */
  134.         centroid(&pbar, p_p, (nvert - 1)) ;
  135.  
  136.  
  137.                 /* form pstar, new reflection trial vertex */
  138.         ptrial(&pstar, &pbar, p_p[nhigh], (1 + alpha), -alpha) ;
  139.         func(&pstar) ;
  140.  
  141.  
  142.  
  143.  
  144.  
  145.  
  146.  
  147.                 /* NELDER-MEAD ALGORITHM = test trial vertex
  148.                  vs old high and low vertices ;
  149.                 construct new trial vertices as appropriate; 
  150.                 set pointer to best new vertex */
  151.         if (pstar.val < p_p[nlow]->val){
  152.             ptrial(&p2star, &pstar, &pbar, (1 + gamma), -gamma) ;
  153.             func(&p2star) ;
  154.             if (p2star.val < p_p[nlow]->val) {
  155.                 p_best = &p2star ;
  156.                 opcode = STAREXPAND ;
  157.             } else {
  158.                 p_best = &pstar ;
  159.                 opcode = FAILREFLECT ;
  160.             }
  161.         } else  if (pstar.val <= p_p[nhigh-1]->val) {
  162.                 p_best = &pstar ;
  163.                 opcode = REFLECT ;
  164.         } else {
  165.             if (pstar.val <= p_p[nhigh]->val) {
  166.                 pcopy(p_p[nhigh], &pstar) ;
  167.                 opcode = STARCONTRACT ;
  168.             } else
  169.                 opcode = HIGHCONTRACT ;
  170.  
  171.             ptrial(&p2star, p_p[nhigh], &pbar, beta, (1-beta)) ;
  172.             func(&p2star) ;
  173.  
  174.             if (p2star.val <= p_p[nhigh]->val)
  175.                 p_best = &p2star ;
  176.             else
  177.                 opcode = SHRINK ;
  178.         }
  179.                 /* END OF NELDER-MEAD ALGORITHM */
  180.  
  181.  
  182.  
  183.  
  184.  
  185.                 /* possible exit from loop on operator kbhit 
  186.                 clear keyboard of any extra chars typed */
  187.         test_break = FALSE ;
  188.         if (KBHIT) {
  189.             c = getchar() ;
  190.             while (KBHIT) {
  191.                 getchar() ;
  192.             }
  193.             if ((c == CTRLX) || (c == CTRLC))
  194.                 test_break = TRUE ;
  195.             else if (c == CTRLS)
  196.                 getchar() ;
  197.  
  198.         }
  199.  
  200.  
  201.                 /* print results */
  202.         fprintf(fptr, "\n%5d %20.14e %20.14e %s\n",
  203.             iter, p_p[nhigh]->val, p_best->val, opname[opcode]) ;
  204.         fpprint(fptr, p_p[nhigh]) ;
  205.         fpprint(fptr, p_best) ;
  206.  
  207.  
  208.                 /* adjust simplex 
  209.                 either loop over all vertices > lowest = [0]
  210.                 to shrink on lowest
  211.                 else copy best movement into high vertex */
  212.         if (opcode == SHRINK)
  213.             for (j = 1; j < nvert; j++) {
  214.                 for (i = 0; i < nparm; i++)
  215.                     p_p[j]->parm[i] = (p_p[nlow]->parm[i] 
  216.                             + p_p[j]->parm[i])/2 ;
  217.                 func(p_p[j]) ;
  218.             }
  219.         else
  220.             pcopy(p_p[nhigh], p_best) ;
  221.  
  222.  
  223.  
  224.                 /* calculate and print new centroid */
  225.         centroid(&pcent, p_p, nvert) ;
  226.         fpprint(fptr, &pcent) ;
  227.  
  228.  
  229.                 /* calculate and print
  230.                 mean and rms of function values */
  231.         mean_func = 0 ;
  232.         for (j = 0; j < nvert; j++)
  233.             mean_func = mean_func + p[j].val ;
  234.         mean_func = mean_func/nvert ;
  235.     
  236.         rms_func = 0 ;
  237.         for (j = 0; j < nvert; j++)
  238.             rms_func = rms_func 
  239.                  + (p[j].val - mean_func) * (p[j].val - mean_func) ;
  240.         rms_func = sqrt(rms_func/nvert) ;
  241.  
  242.         test = rms_func/mean_func ;
  243.         rms_data = sqrt(p_p[nlow]->val/(ndata - nfree));        
  244.         fprintf(fptr, "mean=%20.14e rms=%20.14e test=%20.14e\n",
  245.             mean_func, rms_func, test) ;
  246.         fprintf(fptr, 
  247. "root mean square weighted error of best fit to data = %20.14e\n\n",
  248.                 rms_data) ;
  249.         fspecial(fptr) ;
  250.  
  251.  
  252.  
  253.  
  254.                 /* exit test
  255.                 calculate centroid function value if exit */
  256.         if ((test_break == TRUE) || (iter == maxiter)) {
  257.             func(&pcent) ;
  258.             break ;
  259.         } 
  260.         if (test < exit_test) {
  261.             func(&pcent) ;
  262.             if ((pcent.val - mean_func) < (2 * rms_func))
  263.                 break ;
  264.         }
  265.     }        
  266.                 /* END OF MAIN LOOP OF MINIMIZATION */
  267.  
  268.     return (OK) ;
  269. }                              /* END OF SIMPFIT                  */
  270.  
  271. /* page eject */
  272.  
  273. /* SIMPDEV
  274.     QUADRATIC FIT TO FUNCTION SURFACE FOR ESTIMATION OF:
  275.         VARIANCE-COVARIANCE MATRIX 
  276.         STANDARD DEVIATIONS OF PARAMETERS
  277. */
  278.  
  279. int simpdev(fptr)
  280. FILE *fptr ;
  281. {
  282.     register int i, j, k, l ;
  283.     int xfree, free_count ;
  284.     int c ;
  285.     int err_count ;
  286.  
  287.     double dtemp ;
  288.     double yminusij, yplusij ;
  289.  
  290.     double tmparray[NPARM] ;
  291.  
  292.     struct pstruct ptemp ;
  293.  
  294.     void fprintf() ;
  295.     int getchar() ;
  296.     double sqrt() ;
  297.  
  298.     int func() ;
  299.     void fpprint(), pcopy() ;
  300.     void fqprint(), fmatprint() ;
  301.  
  302.     void bsort();
  303.     int pvalcomp(), p_swap() ;
  304.  
  305.     extern struct pstruct     *p_p[], p[], pcent, pmin ;
  306.     extern struct qstruct     q ;
  307.     extern double         qmat[][NPARM] ;
  308.     extern double       rms_data, yzero, ymin, ypmin, mse ;
  309.     extern int         nparm, nvert, ndata, nfree ;
  310.  
  311.  
  312.  
  313.  
  314.                 /* be sure that pcent.val is set */
  315.     if (func(&pcent) == ERROR)
  316.         return (ERROR) ;
  317.  
  318.  
  319.                 /* ascending sort of pointers p_p 
  320.                 to struct array p */
  321.     bsort(nvert, pvalcomp, p_swap) ;
  322.  
  323.                 /* if lowest vertex < centroid, then
  324.                 replace centroid by lowest vertex */
  325.     if (p_p[0]->val < pcent.val) {
  326.         pcopy(&pcent, p_p[0]);
  327.         fprintf(fptr, "\n\nlowest vertex replaces centroid");
  328.     }
  329.  
  330.  
  331.                 /* set yzero */
  332.     yzero = pcent.val ;
  333.     fprintf(fptr, 
  334.         "\n\nlowest function value = yzero = %20.14e\n\n", yzero) ;
  335.     fpprint(fptr, &pcent) ;
  336.     fprintf(fptr, "\n") ;
  337.  
  338.  
  339.                 /* calculate 
  340.                 q value = avg devn of a free parameter 
  341.                         from the centroid value ;
  342.                 store the q value in structure q :
  343.  
  344.                     q.q[free_cnt] = q value ;
  345.                 also store map of free parm to parm array in
  346.                 q.parmndx[free_cnt] = index of parm 
  347.                         in p[nvert].parm[nparm] ;
  348.                 set nfree and nvert */
  349.     free_cnt = 0 ;
  350.     for (i = 0; i < nparm; i++) {
  351.         dtemp = 0 ;
  352.         for (j = 0; j < nvert; j++)
  353.             dtemp = dtemp + ABS(p[j].parm[i] - pcent.parm[i]) ;
  354.         dtemp = dtemp/nvert ;
  355.         if (ABS(dtemp/pcent.parm[i]) < 1.e-16) {
  356.             fprintf(fptr, "parameter %d is fixed\n", i) ;
  357.             continue ;
  358.         } else {
  359.             q.q[free_cnt] = dtemp ;
  360.             q.parmndx[free_cnt] = i ;
  361.             free_cnt++ ;
  362.         }    
  363.     }
  364.     nfree = free_cnt ;
  365.     if (nvert != (nfree + 1)) {
  366.         fprintf(fptr, "\nerror in count of free parameters\n") ;
  367.         return (ERROR) ;
  368.     }
  369.     fprintf(fptr, 
  370.         "\nq values for the %d free parameters out of %d total:\n", 
  371.         nfree, nparm) ;
  372.     fqprint(fptr, q.q) ;
  373.     
  374.             
  375.  
  376.  
  377.                 /* check and if necessary adjust q values
  378.                 for each free parameter, to be sure that the
  379.                 function values are OK for the 
  380.                 centroid parameter values  (+  &  -)  q ;
  381.                 store function value for parameter value
  382.                 (+  &  -)  q   in  q.yplus  and  q.yminus */
  383.     for (i = 0; i < nfree; i++) {
  384.         k = q.parmndx[i] ;
  385.         pcopy(&ptemp, &pcent) ;
  386.         fprintf(fptr, 
  387.             "\nchecking q value %d, for parameter %d: ", i, k) ;
  388.  
  389.         while (TRUE) {
  390.             fprintf(fptr, "*") ;
  391.  
  392.             while (TRUE) {
  393.                 ptemp.parm[k] = pcent.parm[k] - q.q[i] ;
  394.                 if (func(&ptemp) == ERROR)
  395.                     q.q[i] = q.q[i] / 2 ;
  396.                 else
  397.                     break ;
  398.  
  399.             }
  400.             q.yminus[i] = ptemp.val ;
  401.  
  402.             dtemp = q.q[i] ;
  403.             while (TRUE) {
  404.                 ptemp.parm[k] = pcent.parm[k] + dtemp ;
  405.                 if (func(&ptemp) == ERROR)
  406.                     dtemp = dtemp / 2 ;
  407.                 else
  408.                     break ;
  409.             }
  410.             q.yplus[i] = ptemp.val ;
  411.  
  412.             if (dtemp == q.q[i])
  413.                 break ;
  414.             else
  415.                 q.q[i] = dtemp ;
  416.         }
  417.     }
  418.     fprintf(fptr, "\n\nadjusted q values:\n") ;
  419.     fqprint(fptr, q.q) ;
  420.     fprintf(fptr, "\nyplus values:\n") ;
  421.     fqprint(fptr, q.yplus) ;
  422.     fprintf(fptr, "\nyminus values:\n") ;
  423.     fqprint(fptr, q.yminus) ;
  424.  
  425.  
  426.  
  427.                 /* calculate  
  428.                 <a>  vector
  429.                 store in q.a */
  430.                 /* calculate  
  431.                 <B>  matrix diagonal elements ;
  432.                 store in  qmat  and  q.bdiag  */
  433.     for (i = 0; i < nfree; i++) {
  434.         q.a[i] = 0.25 * (q.yplus[i] - q.yminus[i]) ;
  435.         qmat[i][i] = q.bdiag[i] = 0.5 * (q.yplus[i] + q.yminus[i]
  436.                         - 2 * yzero) ;
  437.     }
  438.     fprintf(fptr, "\n<a> vector:\n") ;
  439.     fqprint(fptr, q.a) ;
  440.     fprintf(fptr, 
  441.         "\n<B> matrix diagonal:\n") ;
  442.     fqprint(fptr, q.bdiag) ;
  443.  
  444.  
  445.  
  446.  
  447.                 /* calculate 
  448.                 <B>  matrix off-diagonal elements ;
  449.                 keep track of any function errors ;
  450.                 store results in  qmat  */
  451.     err_count = 0 ;
  452.     for (i = 1; i < nfree; i++) {
  453.  
  454.         k = q.parmndx[i] ;
  455.         fprintf(fptr, 
  456.              "\ncalculating off-diag Bij for row %d: ", i) ;
  457.         for (j = 0; j < i; j++) {
  458.             fprintf(fptr, "*") ;
  459.             l = q.parmndx[j] ;
  460.             pcopy(&ptemp, &pcent) ;
  461.  
  462.             ptemp.parm[k] = pcent.parm[k] + q.q[i] ;
  463.             ptemp.parm[l] = pcent.parm[l] + q.q[j] ;
  464.             if (func(&ptemp) == ERROR)
  465.                 ++err_count ;
  466.             yplusij = ptemp.val ;
  467.  
  468.             ptemp.parm[k] = pcent.parm[k] - q.q[i] ;
  469.             ptemp.parm[l] = pcent.parm[l] - q.q[j] ;
  470.             if (func(&ptemp) == ERROR)
  471.                 ++err_count ;
  472.             yminusij = ptemp.val ;
  473.  
  474.             qmat[j][i] = qmat[i][j] = 0.25 * (yplusij + yminusij
  475.                 + 2 * yzero - q.yplus[i] - q.yminus[i]
  476.                         - q.yplus[j] - q.yminus[j]) ;
  477.         }
  478.     }
  479.     if (err_count > 0) {
  480.         fprintf(fptr, "\n\n%d function errors in <B> matrix calculation\n",
  481.             err_count) ;
  482.         return (ERROR) ;
  483.     }
  484.     fprintf(fptr, "\n\n<B> matrix:\n") ;
  485.     fmatprint(fptr, qmat) ;
  486.  
  487.  
  488.  
  489.  
  490.                 /* invert
  491.                 <B>  matrix  */
  492.     xfree = nfree - 1 ;
  493.     for (l = 0; l < nfree; l++) {
  494.         tmparray[xfree] = 1 / qmat[0][0] ;
  495.  
  496.         for (i = 0; i < xfree; i++)
  497.             tmparray[i] = tmparray[xfree] * qmat[0][i+1] ;
  498.  
  499.         for (i = 0; i < xfree; i++) {
  500.             qmat[i][xfree] = -tmparray[xfree] * qmat[i+1][0] ;
  501.             for (j = 0; j < xfree; j++)
  502.                 qmat[i][j] = qmat[i+1][j+1]
  503.                     - tmparray[j] * qmat[i+1][0] ;
  504.         }
  505.  
  506.         for (i = 0; i < nfree; i++)
  507.             qmat[xfree][i] = tmparray[i] ;
  508.  
  509.     }
  510.     fprintf(fptr, "\n<B> matrix inverse:\n") ;
  511.     fmatprint(fptr, qmat) ;
  512.  
  513.  
  514.  
  515.  
  516.                 /* store
  517.                 inv B matrix diagonal in q.inv_bdiag */
  518.                 /* calculate
  519.                 ymin = yzero - sumij(ai * invBij * aj) */
  520.                 /* calculate 
  521.                 pmin(k) = centroid(k) - sumj(qi * invBij * aj)
  522.                     where k = parmndx[i] if parm[k] free
  523.                     and pmin = centroid if parm[k] fixed */
  524.     ymin = yzero ;
  525.     pcopy(&pmin, &pcent) ;
  526.     for (i = 0; i < nfree; i++) {
  527.         q.inv_bdiag[i] = qmat[i][i] ;
  528.         k = q.parmndx[i] ;
  529.         for (j = 0; j < nfree; j++) {
  530.             ymin = ymin - q.a[i] * qmat[i][j] * q.a[j] ;
  531.             pmin.parm[k] = pmin.parm[k] 
  532.                     - q.q[i] * qmat[i][j] * q.a[j] ;
  533.         }
  534.     }
  535.     if (func(&pmin) == ERROR)
  536.         fprintf(fptr, "\nerror in y-pmin\n") ;
  537.     ypmin = pmin.val ;
  538.  
  539.  
  540.  
  541.  
  542.  
  543.                    /* calculate 
  544.                 mse = (func val)/degrees freedom 
  545.                 rms error */
  546.                 /* calculate 
  547.                 variance-covariance matrix
  548.                 vcij = qi * invBij * qj * mse */
  549.                 /* calculate 
  550.                 standard deviations */
  551.     mse = pcent.val / (ndata - nfree) ;
  552.     rms_data = sqrt(mse) ;
  553.     for (i = 0; i < nfree; i++) {
  554.         for (j = 0; j < nfree; j++) {
  555.             qmat[i][j] = qmat[j][i] 
  556.                    = q.q[i] * qmat[i][j] * q.q[j] * mse ;
  557.         }
  558.         q.std_dev[i] = sqrt(qmat[i][i]) ;
  559.     }
  560.     fprintf(fptr, 
  561.         "\nvariance-covariance matrix:\n") ;
  562.     fmatprint(fptr, qmat) ;
  563.     
  564.  
  565.  
  566.  
  567.  
  568.                 /* replace centroid with pmin,
  569.                 if y-pmin  <  y-centroid */
  570.                 /* reconstruct simplex, based on q values
  571.                 and centroid, in preparation for more fitting ;
  572.                 nfree vertices are based on nfree q values ;
  573.                 the remaining vertex is set at the centroid 
  574.                 (or pmin if it is lower) */
  575.     fprintf(fptr,"\n\nCalculating for reconstructed simplex...\n\n") ;
  576.     if (pmin.val < pcent.val)
  577.         pcopy(&p[nfree], &pmin) ;
  578.     else 
  579.         pcopy(&p[nfree], &pcent) ;
  580.     for (i = 0; i < nfree; i++) {
  581.         pcopy(&p[i], &pcent) ;
  582.         k = q.parmndx[i] ;
  583.         p[i].parm[k] = p[i].parm[k] + q.q[i] ;
  584.         func(&p[i]) ;
  585.     }
  586.  
  587.     return (OK) ;
  588. }                              /* END OF SIMPDEV                  */
  589.  
  590. /* page eject */
  591.  
  592. /* SUPPORTING  FUNCTIONS 
  593.     NOTE THAT <FUNC()> AND <FDATPRINT()>, <MAIN()>, AND 
  594.     THE INPUT ROUTINE <READ_DATA()> AND ITS SUPPORT ARE IN 
  595.     SEPARATE FILES
  596. */
  597.  
  598. int centroid(pdest, psrc, xvert)
  599. struct pstruct *pdest, *psrc[] ;
  600. int xvert ;
  601. {
  602.     register int i, j ;
  603.  
  604.     extern int     nparm ;
  605.  
  606.     for (i = 0; i < nparm; i++) {
  607.         pdest->parm[i] = psrc[0]->parm[i] ;
  608.         for (j = 1; j < xvert; j++)
  609.             pdest->parm[i] = pdest->parm[i] 
  610.                 + (psrc[j]->parm[i] 
  611.                     - psrc[0]->parm[i]) / xvert ;
  612.     }
  613.     return (OK) ;
  614.  
  615. }                              /* END OF CENTROID                */
  616.  
  617. /* code changes parm even if all parms equal-- WRONG CODE
  618.     for (i = 0; i < nparm; i++) {
  619.         pdest->parm[i] = 0 ;
  620.         for (j = 0; j < xvert; j++)
  621.             pdest->parm[i] = pdest->parm[i] + psrc[j]->parm[i] ;
  622.         pdest->parm[i] = pdest->parm[i] / xvert ;
  623.     } */
  624.  
  625.  
  626.  
  627. int ptrial(pnew, p1, p2, coef1, coef2)
  628. struct pstruct *pnew, *p1, *p2 ;
  629. double coef1, coef2 ;
  630. {
  631.     register int i ;
  632.  
  633.     extern int     nparm ;
  634.  
  635.     for (i = 0; i < nparm; i++)
  636.         pnew->parm[i] = p1->parm[i] + 
  637.                     coef2 * (p2->parm[i] - p1->parm[i]) ;
  638.     return (OK) ;
  639. }                /* END OF PTRIAL            */
  640.  
  641. /* code that possibly changes parm even if all parms equal-- WRONG CODE
  642.  pnew->parm[i] = (coef1 * p1->parm[i]) + (coef2 * p2->parm[i]) ; */
  643.  
  644.  
  645.  
  646.  
  647. void pcopy(pdest, psrc)
  648. struct pstruct *pdest, *psrc ;
  649. {
  650.     register int i ;
  651.  
  652.     extern int     nparm ;
  653.  
  654.     pdest->val = psrc->val ;
  655.     for (i = 0; i < nparm; i++)
  656.         pdest->parm[i] = psrc->parm[i] ;
  657.  
  658. }                /* END OF PCOPY                */
  659.  
  660.  
  661.  
  662.                 /* NOTE-- the swap routine exchanges
  663.                     pointers to structure records ;
  664.                     the comp routines compare a structure
  665.                     member for two structure records, 
  666.                     which may be swapped later     */
  667.  
  668.                 /* swap contents of array elements x and y
  669.                     which are pointers to structure
  670.                     records tested in comp routines */
  671. int p_swap(x, y)
  672. int x, y ;
  673. {
  674.     int *temp ;
  675.  
  676.     extern struct pstruct     *p_p[] ;
  677.  
  678.     temp = p_p[x] ;
  679.     p_p[x] = p_p[y] ;
  680.     p_p[y] = temp ;
  681.  
  682.     return (OK) ;
  683. }                /* END OF P_SWAP                */
  684.  
  685.  
  686.  
  687.                 /* compare float contents of a structure 
  688.                 member for records x and y        */
  689. int pvalcomp(x, y)
  690. int x,y ;
  691. {
  692.     extern struct pstruct     *p_p[] ;
  693.  
  694.     if (p_p[x]->val < p_p[y]->val)
  695.         return (-1) ;
  696.     else if (p_p[x]->val > p_p[y]->val)
  697.         return (1) ;
  698.     else 
  699.         return (0) ;
  700. }                /* END OF PVALCOMP            */
  701.  
  702.  
  703.  
  704.  
  705. char *str_in(fptr, dummy)
  706. char *dummy ;
  707. FILE *fptr ;
  708. {
  709.     register int i ;
  710.  
  711.     int isspace() ;
  712.     int getc() ;
  713.     int ungetc() ;
  714.  
  715.     while (isspace(dummy[0] = getc(fptr))) {
  716.     }
  717.  
  718.     for (i = 1; !isspace(dummy[i] = getc(fptr)); i++)
  719.         ;
  720.     dummy[i] = '\0' ;
  721.  
  722.     return (&dummy[0]) ;
  723. }                /* END OF *STR_IN            */
  724.  
  725.  
  726.  
  727. void ffitprint(fptr)
  728. FILE *fptr;
  729. {
  730.     void fprintf() ;
  731.     void fpprint(), fsprint();
  732.  
  733.     extern int         iter ;
  734.     extern double         mean_func, rms_func, test, rms_data ;
  735.     extern struct pstruct     pcent ;
  736.     extern char        title[] ;
  737.  
  738.  
  739.                 /* first all vertices of simplex */
  740.     fprintf(fptr, 
  741. "\1\f\n%-s\n\nsummary of simplex fitting results:\niteration number %d\n",
  742.          title, iter) ;
  743.     fsprint(fptr, "\nsimplex:\n") ;
  744.         
  745.                 /*then centroid */
  746.     fprintf(fptr, 
  747. "\ncentroid:\nfunction value = %20.14e\n", pcent.val) ;
  748.     fpprint(fptr, &pcent) ;
  749.     
  750.                 /* then test for exit, etc */
  751.     fprintf(fptr, 
  752. "mean=%20.14e rms=%20.14e test=%20.14e\n", mean_func, rms_func, test) ;
  753.     fprintf(fptr, 
  754. "root mean square weighted error of best fit to data = %20.14e\n\n", rms_data);
  755.  
  756.  
  757. }                /* END OF FFITPRINT            */
  758.  
  759.  
  760.  
  761. void fquadprint(fptr)
  762. FILE *fptr;
  763. {
  764.     register int i, k ;
  765.  
  766.     void fprintf() ;
  767.     void fsprint();
  768.  
  769.     extern int         iter, nfree ;
  770.     extern double         rms_data ;
  771.     extern double         yzero, ymin, ypmin, mse ; 
  772.     extern double        mean_func, rms_func, test ;
  773.  
  774.     extern struct qstruct     q ;
  775.     extern struct pstruct     pcent, pmin ;
  776.     extern char        title[] ;
  777.  
  778.     fprintf(fptr, 
  779. "\1\f\n%-s\n\nsummary of results of quadratic fit\niteration number %d\n", 
  780.         title, iter) ;
  781.     fprintf(fptr, "\nmean=%20.14e rms=%20.14e test=%20.14e\n",
  782.         mean_func, rms_func, test) ;
  783.     fprintf(fptr, 
  784. "\nmse    = %20.14e   rms weighted error = %20.14e\n", mse, rms_data) ;
  785.     fprintf(fptr, 
  786. "\nyzero  = %20.14e\nymin   = %20.14e\ny-pmin = %20.14e\n\n", 
  787.         yzero, ymin, ypmin) ;
  788.     
  789.     fprintf(fptr, 
  790. "free parm    q                centroid         pmin             std_dev\n") ;
  791.  
  792.     for (i = 0; i < nfree; i++) {
  793.         k = q.parmndx[i] ;
  794.         fprintf(fptr, "%3d %3d %17.11e %17.11e %17.11e %17.11e\n", 
  795.              k, i, q.q[i], pcent.parm[k], pmin.parm[k], q.std_dev[i]) ;
  796.     }
  797. }                /* END OF FQUADPRINT            */
  798.  
  799.  
  800.  
  801.                 /* print nfree x nfree matrix 
  802.                 calls fqprint()                */
  803. void fmatprint(fptr, g)
  804. FILE *fptr ;
  805. double g[][NPARM] ;
  806. {
  807.     register int i ;
  808.  
  809.  
  810.     extern int     nfree ;
  811.  
  812.     void fqprint() ;
  813.  
  814.     for (i = 0; i < nfree; i++)
  815.         fqprint(fptr, g[i]) ;
  816. }                /* END OF FMATPRINT            */
  817.  
  818.  
  819.  
  820. void fqprint(fptr, g)
  821. FILE *fptr ;
  822. double *g ;
  823. {
  824.     register int i ;
  825.     void fprintf() ;
  826.  
  827.     extern int     nfree ;
  828.  
  829.     for (i = 0; i < nfree; i++) {
  830.         if ((i > 0) && ((i % 4) == NULL))
  831.             fprintf(fptr, "\n  ") ;
  832.         fprintf(fptr, "%19.10e", g[i]) ;
  833.     }
  834.     fprintf(fptr, "\n") ;
  835. }                /* END OF FQPRINT            */
  836.  
  837.  
  838.  
  839. void fpprint(fptr, g)
  840. FILE *fptr ;
  841. struct pstruct *g ;
  842. {
  843.     register int i ;
  844.     void fprintf() ;
  845.  
  846.     extern int     nparm ;
  847.  
  848.     for (i = 0; i < nparm; i++) {
  849.         if (i > 0 && (i % 4) == NULL)
  850.             fprintf(fptr, "\n  ") ;
  851.         fprintf(fptr, "%19.10e", g->parm[i]) ;
  852.     }
  853.     fprintf(fptr, "\n") ;
  854. }                /* END OF FPPRINT            */
  855.  
  856.  
  857.  
  858. void fsprint(fptr, string)
  859. char *string ;
  860. FILE *fptr ;
  861. {
  862.     register int j ;
  863.  
  864.     void fprintf() ;
  865.  
  866.     void fpprint() ;
  867.  
  868.     extern int         nvert ;
  869.     extern struct pstruct     p[] ;
  870.     
  871.     fprintf(fptr, "%s", string) ;
  872.     for (j = 0; j < nvert; j++) {
  873.         fprintf(fptr, "vertex = %d   function value = %20.14e\n", 
  874.             j, p[j].val) ;
  875.         fpprint(fptr, &p[j]) ;
  876.     }
  877. }                /* END OF FSPRINT            */
  878.  
  879.