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