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