home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 209_01 / simpfitr.doc < prev    next >
Encoding:
Text File  |  1990-03-05  |  14.8 KB  |  533 lines
  1. SIMPFITR.DOC         VERS:- 01.00  DATE:- 09/26/86  TIME:- 10:02:49 PM
  2.  
  3. documentation for C routines for simplex fitting and
  4.     quadratic approximation of least squares surface.
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
  12.         By J.A. Rupley, Tucson, Arizona
  13.  
  14.  
  15.          NOTES ON THE SIMPLEX FITTING PROGRAM XXXXFITn
  16.          AND ITS SUPPORTING MODULES
  17.  
  18.  
  19.         A.  DESCRIPTION OF THE PROGRAM:
  20.  
  21.  
  22.                 A function explicitly relating a dependent variable to
  23.         one or more independent variables and up to ten variable
  24.         parameters is fit to a set of data points (observed values of the
  25.         dependent and independent variables).  The Nelder-Mead algorithm
  26.         is used to locate the least squares minimum by a simplex search
  27.         procedure, giving estimates of the best-fit values of the
  28.         parameters.  Standard deviations of the parameters are estimated
  29.         by use of a quadratic approximation of the least squares surface
  30.         near the minimum.  The quadratic approximation, equivalent to
  31.         minimization by use of a Taylor's series approximation, also
  32.         returns an improved estimate of the parameter.  Thus the
  33.         quadratic approximation can be used as part of the minimization
  34.         procedure.
  35.  
  36.              Alternation of  several cycles (eg, 30) of simplex
  37.         minimization with one cycle of quadratic minimization quickly
  38.         converges on the least squares minimum.  The simplex search
  39.         efficiently defines the region of parameter space that contains
  40.         the minimum; bounds are easily incorporated; because derivatives
  41.         are not used, ill-behaved functions can be handled.  The
  42.         quadratic minimization, on the other hand, rapidly moves to the
  43.         least squares minimum, but it is reliable only when the parameter
  44.         estimates have been brought near the minimum.
  45.  
  46.              One cycle of quadratic minimization takes about as much time
  47.         as 30 cycles of simplex minimization.  If the quadratic
  48.         minimization fails, as it often will in the early stages of the
  49.         search, considerable time is wasted.  To reduce this waste, the
  50.         next one or several cycles of quadratic minimization are skipped
  51.         after a failure.
  52.  
  53.  
  54.  
  55.  
  56.  
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.  
  64.  
  65.  
  66.  
  67.  
  68.  
  69.                                         1
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.  
  79.         USAGE:
  80.  
  81.         A>xxxxfitn   [d:]input_file   [[d:]output_file]
  82.  
  83.            the optional output file can be LST: or a disk file;
  84.  
  85.            the required input file has the following structure:
  86.                   one-line title, with no tabs or ctrl characters;
  87.                   lines following give control variables, the starting
  88.                            simplex, and the data array;
  89.                   the order of entry is fixed by <read_data()>;
  90.                            see sample file for structure;
  91.                   the one-word descriptors must be present;
  92.                   the data format is free-form (no set field widths);
  93.                   comments at the end of the file are not read by the
  94.                             program and any other information to be
  95.                             stored should be there.
  96.  
  97.  
  98.  
  99.  
  100.  
  101.  
  102.  
  103.  
  104.  
  105.  
  106.  
  107.  
  108.  
  109.  
  110.  
  111.  
  112.  
  113.  
  114.  
  115.  
  116.  
  117.  
  118.  
  119.  
  120.  
  121.  
  122.  
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.  
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136.                                         2
  137.  
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144.  
  145.  
  146.         B.  CONTENTS OF THE C SOURCE FILES:
  147.  
  148.  
  149.         SIMPMAIN:
  150.         main program controlling input, output & fitting        = main()
  151.  
  152.  
  153.         XXXXFITn:
  154.         routines for simplex fitting that depend on the data and function
  155.         to be fit (must be changed when function fit is changed):
  156.                 all external definitions, except for <data>, which is
  157.                         defined as dummy storage in SIMPLIB1
  158.                         or its modification
  159.                 function for calculation of dependent variable  and
  160.                         weighted sum of residuals squared       = func()
  161.                 print of <data> records                         = fdatprint()
  162.                 customizable display called by <simpfit()>      = fspecial()
  163.                 other special functions called by above functions
  164.  
  165.  
  166.         SIMPLIB0:
  167.         general routines for simplex fitting:
  168.                 simplex fitting routine                          = simpfit()
  169.                 quadratic fit for standard deviations            = simpdev()
  170.                 supporting functions
  171.  
  172.  
  173.  
  174.         SIMPLIB1:
  175.         routines for data input that depend (partially) on the data
  176.         and function to be fit:
  177.                 definition of the aggregate <data>, with a
  178.                        dummy structure declaration
  179.                 usage message displayed on command line error     = use_mess()
  180.                 opening of files for input and optional output    = file()
  181.                 input of variables and data                       = read_data()
  182.  
  183.  
  184.  
  185.  
  186.  
  187.  
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202.  
  203.                                         3
  204.  
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.         C.  PROGRAM FLOW:
  214.  
  215.  
  216.              Program flow is controlled by <main()>.
  217.  
  218.  
  219.         I.  A call to <read_data()> at the start of execution initializes
  220.         the following:
  221.  
  222.         a one-line title string, without tabs or other control
  223.         characters;
  224.  
  225.         control values that are used in testing (1) for exit from the
  226.         minimization (<exit_test>), (2) for cycles on which to print
  227.         intermediate fitting results (<prt_cycle>), and (3) for cycles on
  228.         which to pass through the quadratic approximation (<maxquad_skip>
  229.  
  230.         and <quad_test>);
  231.  
  232.         <iter>, the starting interation number, generally set at 0;
  233.         the starting simplex <p>, with <nvert> vertices and <nparm>
  234.         parameters; <nvert> = 1 + number of free parameters, <nfree>;
  235.         parameters that are the same for all vertices are treated as
  236.         fixed;
  237.  
  238.         the data set, stored in the aggregate <data> of size <ndata>,
  239.         each record of <data[ndata]> having <ndatval> members.
  240.  
  241.  
  242.              On return from <read_data()>:
  243.  
  244.         <maxiter> is set at <iter> + <prt_cycle>; <nquad_skip> and
  245.         <quad_cycle>, used in resetting <quad_test>, are set at zero and
  246.         <quad_test>, respectively;
  247.  
  248.  
  249.              The first command line argument specifies the input file
  250.         from which the above information is read.  The (optional) second
  251.         argument specifies the (optional) output file or device, on which
  252.         a selection of the crt display is saved.
  253.  
  254.              A sample input file is given in a following section.
  255.  
  256.  
  257.         II.  Fitting by use of the Nelder-Mead simplex algorithm is
  258.         carried out by call of <simpfit()>.  Control is returned to main
  259.         every <prt_cycle> number of fitting cycles, for printout of the
  260.         current simplex and its centroid.
  261.  
  262.  
  263.         III.  On certain cycles of simplex minimization, determined by
  264.         the value of <quad_test>, the current simplex is passed to the
  265.         quadratic approximation routine <simpdev()>.  Printout consists
  266.         of (1) the data array with calculated values of the dependent
  267.         variable and (2) a summary of the results of the quadratic fit,
  268.  
  269.  
  270.                                         4
  271.  
  272.  
  273.  
  274.  
  275.  
  276.  
  277.  
  278.  
  279.  
  280.         including estimates of the standard deviations of the parameters.
  281.         <Simpdev()> returns a modified simplex, with, if quadratic
  282.         minimization was successful, the new low vertex closer to the
  283.         least squares minimum.  The variance-covariance matrix is
  284.         returned in <qmat[nfree][nfree]>.  The standard deviations of the
  285.  
  286.         parameters and other information are returned in the structure
  287.         <q>.
  288.  
  289.  
  290.         IV.  Minimization consists of looping through steps II and III,
  291.         until the value of <test> is less than <exit_test>:
  292.  
  293.              <exit_test> = input value
  294.              <test>      = <rms_func>/<mean_func>
  295.              <rms_func>  = (root mean square of the deviations
  296.                                       of the least squares values
  297.                                       at the simplex vertices)
  298.              <mean_func> = (mean of the least squares values)
  299.  
  300.              On exit from the minimization, a pass through <simpdev()>
  301.         gives the final estimates of the standard deviations of the
  302.         parameters.
  303.  
  304.              The alternation of a set of simplex fitting cycles with a
  305.         pass through the quadratic approximation gives rapid convergence
  306.         to the least squares minimum.  In order to control this,
  307.         <quad_test> can be used in either of two ways:
  308.  
  309.              If in the input file <quad_test> is set < 1, then
  310.         <simpdev()> is entered at the first printing cycle for which
  311.         <test> is less than <quad_test>.  Before looping back to
  312.         <simpfit()>, <quad_test> is reduced by a factor of 10.
  313.  
  314.              If in the input file <quad_test> is set >= 1 (in which case
  315.         it should be an integer multiple of <prt_cycle>), then
  316.         <simpdev()> is entered every {<quad_cycle> * (<nquad_skip> + 1)}
  317.         number of cycles.  <nquad_skip> is initially zero; it is
  318.         incremented on each failure of the quadratic minimization, up to
  319.         the limit <maxquad_skip>; on a successful quadratic minimization,
  320.         <nquad_skip> is reset at zero.  This algorithm reduces the time
  321.         spent in unsuccessful quadratic minimization.
  322.  
  323.  
  324.              The operator can choose, from within <simpfit()>, to return
  325.         to <main()> and pass through <simpdev()> for a cycle of quadratic
  326.         approximation.
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.                                         5
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.         D.  CODING IN XXXXFITn FOR THE SPECIFIC FUNCTION TO BE FIT:
  349.  
  350.  
  351.              The routine <func()> calculates for a simplex vertex the
  352.         weighted sum of residuals for the set of parameter values passed
  353.         to it.  The result is stored in <p.val>.  Bounds on the
  354.         parameters can be implemented easily, by a test that returns an
  355.         excessively large function value (eg, 1.E38) if a parameter
  356.         exceeds a bound.  Coding should be efficient, because much of the
  357.         minimization time is spent in <func()>.
  358.  
  359.              The print statement of the routine <fdatprint()> may have to
  360.         be recoded after change in the function or data.
  361.  
  362.              The routine <fspecial()> controls the customizable display
  363.         printed at the end of the standard display for each cycle of the
  364.         simplex minimization in <simpfit()>.
  365.  
  366.              Other routines should be function and data independent,
  367.         unless special manipulation of the data is required.
  368.  
  369.              NOTE: change in the model being fit should require changes
  370.         in the coding only of the function of XXXXFITn.
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.                                         6
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.         E. COMMENTS ON THE CONSTRUCTION OF THE SIMPLEX FITTING PROGRAM:
  416.  
  417.  
  418.              About 4K of memory are reserved to allow expansion of
  419.         <main()> and<func()>, for more elaborate output, functions, etc.
  420.  
  421.  
  422.              The maximum number of parameters, NPARM, is set at 10; if
  423.         more are needed, all routines (SIMPMAIN, SIMPLIB0, etc) must be
  424.         recompiled with the new value of NPARM.
  425.  
  426.  
  427.              To make more readable the coding in <func()> of the model
  428.         equation to be fit to the data:
  429.  
  430.         (1) use mnemonic member names in declaring <struct dat> in
  431.         XXXXFITn;
  432.  
  433.         (2) declare a dummy structure, <struct pnamestruct>, that is
  434.         entirely equivalent to the structure that holds the parameter
  435.         values, <pstruct>, but that has mnemonic member names; the
  436.         mnemonic dummy structure then can be used with the <pstruct>
  437.         address passed as a parameter to <func()>.
  438.  
  439.  
  440.              The DEFINITION of the aggregate <data> and the functions
  441.         <use_mess()>, <file()>, and <read_data()> are in a separate file,
  442.         SIMPLIB1; this is to to allow expansion of the aggregate <data>
  443.         by overwriting most of SIMPLIB1; the SIMPLIB1 routines are
  444.         entered only once, at the start of execution.
  445.  
  446.              the DECLARATION of <struct dat> according to the
  447.         requirements of the model described by <func()> is given in
  448.         XXXXFITn; <func()>, etc. reference the aggregate <data> as
  449.         external to XXXXFITn, but through the structure <dat>, declared
  450.         locally in XXXXFITn with mnemonic member names suitable for use
  451.         in the coding of <func()>, etc.
  452.  
  453.  
  454.  
  455.              The intent is to generalize the read of the data file and
  456.         the allocation of data storage (in SIMPLIB1), while retaining
  457.         flexibility in the declaration of <struct dat data> in XXXXFITn;
  458.         the following comments bear upon this arrangement:
  459.  
  460.              The loading of values into the aggregate <data> is done in
  461.         the SIMPLIB1 module <read_me()> by use of:
  462.  
  463.         (1) a generalized ("dummy") structure for <data>;
  464.  
  465.         (2) a read loop that moves successive double values from the
  466.         ascii data file into the storage at and above <data[0]>, without
  467.         referencing the elements of <data> by structure member or index.
  468.  
  469.  
  470.  
  471.  
  472.                                         7
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.              The "useful" declaration of the structure for <data> is in
  484.         XXXXFITn, where it is referenced by <func()> and <fdatprint()>;
  485.         <struct dat> must be changed to accord with the requirements of
  486.         <func()> and <fdatprint()>; all members of <struct dat> MUST be
  487.         of type double.
  488.  
  489.  
  490.              Change in the model being fit should not require recoding
  491.         and recompilation of <read_me()> or of any other routines except
  492.         those of XXXXFITn; of course, change in the model requires change
  493.         of <func()>, <fdatprint()>, and of the declarations of <struct
  494.         dat>
  495.         and <struct pnamestruct> in XXXXFITn.
  496.  
  497.  
  498.              The size of the <data> aggregate is limited by (1) the size
  499.         of free memory and (2) the size of SIMPLIB1, most of which can be
  500.         overwritten by data records; for this version of the program,
  501.         SIMPLIB1 corresponds to about 600 double values, and unused
  502.         memory to about 600 double values; overwriting of the code of
  503.         SIMPLIB1 may not be allowed by some compilers.
  504.  
  505.  
  506.              For the six-member structure <dat> used in this version of
  507.         the program, the maximum number of data points is more than 100
  508.         (600 double values), expandable to more than 200 (1200 double
  509.  
  510.         values) if SIMPLIB1 is recompiled with an increase of NDATA;
  511.         NDATA is currently set at 350 double values; increase of NDATA of
  512.         course decreases the amount of memory available for expansion of
  513.         the code of <main()>, <func()>, etc.
  514.  
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.                                         8
  540.  
  541.  
  542.  
  543.