home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 209_01 / linefitr.c < prev    next >
C/C++ Source or Header  |  1990-03-04  |  6KB  |  282 lines
  1. /* LINEFITR.C    VERS:- 01.00  DATE:- 09/26/86  TIME:- 09:39:08 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. test program operation:
  15.     by fit of the linear function:     y  =  a  +  b * x
  16.     to data that can be reduced by hand.
  17.  
  18. compile with:
  19.     simpmain, simplib1, simplib0
  20.     floating point functions
  21.  
  22. the precision of the floating point functions must be high, ie 14 or more
  23.     decimal digits; eight digit precision is not enough for the 
  24.     matrix operations.
  25.  
  26.  
  27. contents of the xxxxfitn module  =  
  28. declarations and routines for simplex fitting special to function to be fit:
  29.  
  30.     function for calculation of dependent variable and
  31.         weighted sum of residuals squared        = func()
  32.  
  33.     print of <data> records                    = fdatprint()
  34.                                 = fpointprint()
  35.  
  36.     other special displays                    = fspecial()
  37.  
  38.     all external definitions, except for <data>, which is 
  39.         defined as dummy storage in SIMPLIB1 
  40.         or its modification
  41.  
  42.  
  43. By J.A. Rupley, Tucson, Arizona
  44. Coded for ECO C compiler, version 3.40
  45. */
  46.  
  47.     /* page eject */
  48.  
  49. /*
  50.      To make more readable the coding in <func()> of the model 
  51. equation to be fit to the data: 
  52.  
  53. (1) use mnemonic member names in declaring <struct dat> in 
  54. XXXXFITn; 
  55.  
  56. (2) declare a dummy structure, <struct pnamestruct>, that is 
  57. entirely equivalent to the structure that holds the parameter 
  58. values, <pstruct>, but that has mnemonic member names; the 
  59. mnemonic dummy structure then can be used with the <pstruct> 
  60. address passed as a parameter to <func()>.
  61.  
  62.  
  63.      Change in the model being fit should not require recoding 
  64. and recompilation of <read_data()> or any other routines except 
  65. those of XXXXFITn; of course, change in the model requires change 
  66. of <func()>, <fdatprint()>, and of the declarations of <struct dat> 
  67. and <struct pnamestruct> in XXXXFITn.
  68.  
  69. */
  70.  
  71. /* page eject */
  72.  
  73. #include <stdio.h>
  74. #include <ctrlcnst.h>
  75.  
  76. #define NVERT        3    /* set equal to 1 + number of parameters 
  77.                 used in <func()> */
  78.  
  79.  
  80. #define NPARM        10    /* do NOT change this define */
  81.  
  82.  
  83.  
  84.                 /* EXTERNALS AND STRUCTURES */
  85.  
  86.                 /* do NOT change any structure except <dat>*/
  87.  
  88.                 /* the structure <dat> MUST be changed 
  89.                 to accord with the requirements of 
  90.                 <func()> and <fdatprint()>; */
  91.                 /* all members of <struct dat> MUST be of 
  92.                 type double. */
  93.                 /* as written the program accomodates 
  94.                 more than 100 six-member data points */
  95.                 /* see the header above for more description
  96.                 of the declaration of <struct dat> in XXXXFITn
  97.                 and the definition of the aggregate <data> in
  98.                 SIMPLIB1 */
  99.  
  100.                 /* the structure <pnamestruct> is equivalent 
  101.                 to <pstruct>, which holds the parameter values;
  102.                 use of <pnamestruct> allows mnemonic access 
  103.                 to the contents of <pstruct> */
  104.  
  105. struct dat {            
  106.     double y ;        
  107.     double yc ;
  108.     double w ;
  109.     double x ;
  110. } ;
  111.  
  112. struct pstruct {
  113.     double val ;
  114.     double parm[NPARM] ;
  115. } ;
  116.  
  117. struct pnamestruct {
  118.     double val ;
  119.     double a ;
  120.     double b ;
  121.     double dummy[8] ;
  122. } ;
  123.  
  124. struct qstruct {
  125.     int parmndx[NPARM] ;
  126.     double q[NPARM] ;
  127.     double yplus[NPARM] ;
  128.     double yminus[NPARM] ;
  129.     double a[NPARM] ;
  130.     double bdiag[NPARM] ;
  131.     double inv_bdiag[NPARM] ;
  132.     double std_dev[NPARM] ;
  133. } ;
  134.  
  135.  
  136. struct pstruct p[NVERT] ;
  137.  
  138. struct pstruct pcent ;
  139.  
  140. struct pstruct *p_p[NVERT] ;
  141.  
  142. struct pstruct pmin ;
  143.  
  144. struct qstruct q ;
  145.  
  146.  
  147. double qmat[NPARM][NPARM] ;
  148.  
  149. double mean_func, rms_func, test, rms_data ;
  150.  
  151. double yzero, ymin, ypmin, mse ;
  152.  
  153. double quad_test, exit_test ;
  154.  
  155. int iter, maxiter, nparm, nvert, ndata, maxquad_skip, prt_cycle, ndatval ;
  156.  
  157. int nfree ;
  158.  
  159. char title[80] ;
  160.  
  161. /* page eject */
  162.  
  163. /* FUNC
  164.     CALCULATION OF LEAST SQUARES FUNCTION
  165.     CODED ACCORDING TO MODEL BEING FIT
  166.     IT SHOULD BE EFFICIENT
  167.     DURING THE FIT, TIME IS MOSTLY SPENT HERE
  168.     (THE OVERHEAD ELSEWHERE IS ONLY SEVERAL SECONDS PER CYCLE)
  169. */
  170.  
  171. int func(g)
  172. struct pstruct *g ;
  173. {
  174.     void printf() ;
  175.     register int n ;
  176.  
  177.     struct pnamestruct *pnam ;
  178.  
  179.     extern struct dat     data[] ;
  180.     extern int         nparm ;
  181.     extern int         ndata ;
  182.  
  183. /*
  184. /* example of setting of bounds, not applicable to LINEFITR */
  185.     for (n = 0; n < nparm; n++)
  186.         if (g->parm[n] <= 0) {
  187.             g->val = g->val + 1.E36 ;
  188.             printf("function error\n") ;
  189.             return (ERROR) ;
  190.         }
  191. */
  192.  
  193.     pnam = g ;
  194.     pnam->val = 0 ;
  195.     for (n = 0; n < ndata; n++) {
  196.         data[n].yc = pnam->a  +  pnam->b * data[n].x;
  197.         pnam->val = pnam->val 
  198.                 + (data[n].y - data[n].yc) 
  199.                 * (data[n].y - data[n].yc) 
  200.                 * data[n].w 
  201.                 * data[n].w ;
  202.     }
  203.     return (OK) ;
  204. }                              /* END OF FUNC                    */
  205.  
  206. /* page eject */
  207.  
  208. /* FDATPRINT
  209.     PRINT DATA AND COMPARE WITH CALCULATED VALUES
  210.     CODED ACCORDING TO MODEL AND DATA
  211. */
  212.  
  213. void fdatprint(fptr)
  214. FILE *fptr ;
  215. {
  216.     register int j ;
  217.  
  218.     void fprintf() ;
  219.  
  220.     extern int         iter, ndata, maxiter ;
  221.     extern struct dat     data[] ;
  222.     extern char        title[] ;
  223.  
  224.     fprintf(fptr, "\1\f\n%-s\n\niteration number %d\n\n", title, iter) ;
  225.     fprintf(fptr, 
  226. "        yobs   ycalc   del-y  weight       x\n") ;
  227.     for (j = 0; j < ndata; j++) 
  228.         fprintf(fptr, 
  229.             "%4d %7.3f %7.3f %7.3f %7.3f %10.5f\n", 
  230.             (j+1), data[j].y, data[j].yc, (data[j].y - data[j].yc),
  231.             data[j].w, data[j].x) ;
  232.  
  233.     if (maxiter != 0)
  234.         fpointprint(fptr) ;
  235. }                /* END OF FDATPRINT            */
  236.  
  237.  
  238.  
  239.  
  240. /* FPOINTPRINT
  241.     PRINT POINTS FOR CONSTRUCTION OF PRIMARY AND SECONDARY PLOTS
  242.     CODED ACCORDING TO MODEL AND DATA
  243. */
  244.  
  245. void fpointprint(fptr)
  246. FILE *fptr ;
  247. {
  248. }                /* END OF FPOINTPRINT            */
  249.  
  250.  
  251.  
  252.  
  253. /* FSPECIAL
  254.     DISPLAY ADDITIONAL INFORMATION DURING
  255.     TRACKING OF MINIMIZATION, IN SIMPFIT()
  256. */
  257.  
  258. void fspecial(fptr)
  259. FILE *fptr ;
  260. {
  261.     register j ;
  262.  
  263.     void fprintf() ;
  264.  
  265.     extern int         maxiter ;
  266.     extern double        ypmin, yzero, quad_test ;
  267.  
  268.     if (ypmin > 0.99E38) {
  269.         fprintf(fptr, "y-pmin error") ;
  270.         for (j = 0; j < nvert; j++)
  271.             if (pmin.parm[j] < 0)
  272.                 fprintf(fptr, " (parm(%d) < 0)", j) ;
  273.     }
  274.     else if (ypmin > yzero)
  275.         fprintf(fptr, "y-pmin > yzero") ;
  276.  
  277.     fprintf(fptr, "    next_prt= %d  next_quad= %7.2e\n", 
  278.         maxiter, quad_test) ;
  279.  
  280.     return ;
  281. }                /* END OF FSPECIAL            */
  282.