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 >
C/C++ Source or Header  |  1990-03-05  |  19KB  |  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