optim_private_nlsq.c

ASUFIT－Matlab－全局拟合程序

/*
 * MATLAB Compiler: 2.0.1
 * Date: Tue May 08 21:28:22 2001
 * Arguments: "-B" "sgl" "-m" "-W" "mainhg" "-L" "C" "asufit.m" "absfun.m"
 * "absfunfree.m" "absfunshift.m" "absfunwidth.m" "dispfit.m" "dispfitdisp.m"
 * "dispfitfree.m" "dispfitwidth.m" "seqmodfree.m" "spcfun.m" "spcfun1.m"
 * "atamult.m" "aprecon.m" 
 */
#include "optim_private_nlsq.h"
#include "cubici1.h"
#include "cubici3.h"
#include "graderr.h"
#include "searchq.h"

extern mxArray * OPT_STEP;
extern mxArray * OPT_STOP;

static double __Array0_r[1] = { 0.0 };

static double __Array1_r[1] = { 0.0 };

static mxArray * Moptim_private_nlsq(mxArray * * CostFunction,
                                     mxArray * * JAC,
                                     mxArray * * EXITFLAG,
                                     mxArray * * OUTPUT,
                                     int nargout_,
                                     mxArray * funfcn,
                                     mxArray * x_,
                                     mxArray * verbosity,
                                     mxArray * options,
                                     mxArray * CostFunction_,
                                     mxArray * JAC_,
                                     mxArray * YDATA,
                                     mxArray * caller,
                                     mxArray * varargin);
static mxArray * Mnlsq_lsinit(mxArray * * fold,
                              mxArray * * how,
                              int nargout_,
                              mxArray * xnew,
                              mxArray * fnew,
                              mxArray * verbosity,
                              mxArray * levMarq);
static mxArray * mlfNlsq_lsinit(mxArray * * fold,
                                mxArray * * how,
                                mxArray * xnew,
                                mxArray * fnew,
                                mxArray * verbosity,
                                mxArray * levMarq);
static void mlxNlsq_lsinit(int nlhs,
                           mxArray * plhs[],
                           int nrhs,
                           mxArray * prhs[]);

static mlfFunctionTableEntry local_function_table_[1]
  = { { "lsinit", mlxNlsq_lsinit, 4, 3 } };

/*
 * The function "Moptim_private_nlsq" is the implementation version of the
 * "optim/private/nlsq" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\nlsq.m" (lines 1-413). It contains the
 * actual compiled code for that M-function. It is a static function and must
 * only be called from one of the interface functions, appearing below.
 */
/*
 * function [x,CostFunction,JAC,EXITFLAG,OUTPUT] = nlsq(funfcn,x,verbosity,options,CostFunction,JAC,YDATA,caller,varargin)
 */
static mxArray * Moptim_private_nlsq(mxArray * * CostFunction,
                                     mxArray * * JAC,
                                     mxArray * * EXITFLAG,
                                     mxArray * * OUTPUT,
                                     int nargout_,
                                     mxArray * funfcn,
                                     mxArray * x_,
                                     mxArray * verbosity,
                                     mxArray * options,
                                     mxArray * CostFunction_,
                                     mxArray * JAC_,
                                     mxArray * YDATA,
                                     mxArray * caller,
                                     mxArray * varargin) {
    mxArray * x = mclGetUninitializedArray();
    mxArray * CHG = mclGetUninitializedArray();
    mxArray * DerivativeCheck = mclGetUninitializedArray();
    mxArray * DiffMaxChange = mclGetUninitializedArray();
    mxArray * DiffMinChange = mclGetUninitializedArray();
    mxArray * EstSum = mclGetUninitializedArray();
    mxArray * FIRSTF = mclGetUninitializedArray();
    mxArray * GDOLD = mclGetUninitializedArray();
    mxArray * GradF = mclGetUninitializedArray();
    mxArray * GradFactor = mclGetUninitializedArray();
    mxArray * JACFD = mclGetUninitializedArray();
    mxArray * JACOB = mclGetUninitializedArray();
    mxArray * MATL = mclGetUninitializedArray();
    mxArray * MATX = mclGetUninitializedArray();
    mxArray * NewF = mclGetUninitializedArray();
    mxArray * OLDF = mclGetUninitializedArray();
    mxArray * OLDG = mclGetUninitializedArray();
    mxArray * OLDJ = mclGetUninitializedArray();
    mxArray * OLDX = mclGetUninitializedArray();
    mxArray * OX = mclGetUninitializedArray();
    mxArray * PCNT = mclGetUninitializedArray();
    mxArray * SD = mclGetUninitializedArray();
    mxArray * XOUT = mclGetUninitializedArray();
    mxArray * currstepsize = mclGetUninitializedArray();
    mxArray * estf = mclGetUninitializedArray();
    mxArray * fbest = mclGetUninitializedArray();
    mxArray * formatstr = mclGetUninitializedArray();
    mxArray * gcnt = mclGetUninitializedArray();
    mxArray * gdnew = mclGetUninitializedArray();
    mxArray * gradflag = mclGetUninitializedArray();
    mxArray * how = mclGetUninitializedArray();
    mxArray * iter = mclGetUninitializedArray();
    mclForLoopIterator iterator_0;
    mxArray * levMarq = mclGetUninitializedArray();
    mxArray * lineSearchType = mclGetUninitializedArray();
    mxArray * maxFunEvals = mclGetUninitializedArray();
    mxArray * maxIter = mclGetUninitializedArray();
    mxArray * newf = mclGetUninitializedArray();
    mxArray * newstep = mclGetUninitializedArray();
    mxArray * nfun = mclGetUninitializedArray();
    mxArray * nonZeroSum = mclGetUninitializedArray();
    mxArray * num = mclGetUninitializedArray();
    mxArray * numFunEvals = mclGetUninitializedArray();
    mxArray * numGradEvals = mclGetUninitializedArray();
    mxArray * numberOfVariables = mclGetUninitializedArray();
    mxArray * status = mclGetUninitializedArray();
    mxArray * steplen = mclGetUninitializedArray();
    mxArray * stepsize = mclGetUninitializedArray();
    mxArray * sumabsJAC = mclGetUninitializedArray();
    mxArray * temp = mclGetUninitializedArray();
    mxArray * tolFun = mclGetUninitializedArray();
    mxArray * tolX = mclGetUninitializedArray();
    mclValidateInputs(
      "optim/private/nlsq",
      8,
      &funfcn,
      &x_,
      &verbosity,
      &options,
      &CostFunction_,
      &JAC_,
      &YDATA,
      &caller);
    mclCopyInputArg(&x, x_);
    mclCopyInputArg(CostFunction, CostFunction_);
    mclCopyInputArg(JAC, JAC_);
    /*
     * %NLSQ Solves non-linear least squares problems.
     * %   NLSQ is the core code for solving problems of the form:
     * %   min  sum {FUN(X).^2}    where FUN and X may be vectors or matrices.   
     * %             x
     * %
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.13 $  $Date: 1998/07/30 16:31:50 $
     * %   Andy Grace 7-9-90.
     * 
     * %   The default algorithm is the Levenberg-Marquardt method with a 
     * %   mixed quadratic and cubic line search procedure.  A Gauss-Newton
     * %   method is selected by setting  OPTIONS.LevenbergMarq='on'. 
     * %
     * 
     * % ------------Initialization----------------
     * 
     * 
     * XOUT = x(:);
     */
    mlfAssign(&XOUT, mlfIndexRef(x, "(?)", mlfCreateColonIndex()));
    /*
     * % numberOfVariables must be the name of this variable
     * numberOfVariables = length(XOUT);
     */
    mlfAssign(&numberOfVariables, mlfLength(XOUT));
    /*
     * how = [];
     */
    mlfAssign(&how, mclCreateEmptyArray());
    /*
     * 
     * numFunEvals = 0;
     */
    mlfAssign(&numFunEvals, mlfScalar(0.0));
    /*
     * numGradEvals = 0;
     */
    mlfAssign(&numGradEvals, mlfScalar(0.0));
    /*
     * 
     * OUTPUT = [];
     */
    mlfAssign(OUTPUT, mclCreateEmptyArray());
    /*
     * iter = 0;
     */
    mlfAssign(&iter, mlfScalar(0.0));
    /*
     * JACOB = [];
     */
    mlfAssign(&JACOB, mclCreateEmptyArray());
    /*
     * EXITFLAG = 1;  %assume convergence
     */
    mlfAssign(EXITFLAG, mlfScalar(1.0));
    /*
     * currstepsize = 0;
     */
    mlfAssign(&currstepsize, mlfScalar(0.0));
    /*
     * 
     * % Global parameters for outside control of leastsq
     * % OPT_STOP is used for prematurely stopping the optimization
     * % OPT_STEP is set to 1 during major (non-gradient finding) iterations
     * %          set to 0 during gradient finding and 2 during line search
     * %          this can be useful for plotting etc.
     * global OPT_STOP OPT_STEP;
     * OPT_STEP = 1;
     */
    mlfAssign(&OPT_STEP, mlfScalar(1.0));
    /*
     * OPT_STOP = 0;
     */
    mlfAssign(&OPT_STOP, mlfScalar(0.0));
    /*
     * 
     * formatstr=' %5.0f       %5.0f   %13.6g %12.3g %12.3g  ';
     */
    mlfAssign(
      &formatstr,
      mxCreateString(" %5.0f       %5.0f   %13.6g %12.3g %12.3g  "));
    /*
     * 
     * OPT_STEP = 0;  % No longer a major step
     */
    mlfAssign(&OPT_STEP, mlfScalar(0.0));
    /*
     * 
     * nfun=length(CostFunction);
     */
    mlfAssign(&nfun, mlfLength(*CostFunction));
    /*
     * 
     * gradflag =  strcmp(optimget(options,'Jacobian'),'on');
     */
    mlfAssign(
      &gradflag,
      mlfStrcmp(
        mlfOptimget(options, mxCreateString("Jacobian"), NULL),
        mxCreateString("on")));
    /*
     * tolX = optimget(options,'tolx');
     */
    mlfAssign(&tolX, mlfOptimget(options, mxCreateString("tolx"), NULL));
    /*
     * % lineSearchType of 0 means quadcubic (the default)
     * lineSearchType = strcmp(optimget(options,'linesearchtype'),'cubicpoly');
     */
    mlfAssign(
      &lineSearchType,
      mlfStrcmp(
        mlfOptimget(options, mxCreateString("linesearchtype"), NULL),
        mxCreateString("cubicpoly")));
    /*
     * % levMarq=1 is the default and means use levenbergmarquardt
     * levMarq = strcmp(optimget(options,'levenbergMarq'),'on');
     */
    mlfAssign(
      &levMarq,
      mlfStrcmp(
        mlfOptimget(options, mxCreateString("levenbergMarq"), NULL),
        mxCreateString("on")));
    /*
     * 
     * tolFun = optimget(options,'tolfun');
     */
    mlfAssign(&tolFun, mlfOptimget(options, mxCreateString("tolfun"), NULL));
    /*
     * DiffMinChange = optimget(options,'diffminchange');
     */
    mlfAssign(
      &DiffMinChange,
      mlfOptimget(options, mxCreateString("diffminchange"), NULL));
    /*
     * DiffMaxChange = optimget(options,'diffmaxchange');
     */
    mlfAssign(
      &DiffMaxChange,
      mlfOptimget(options, mxCreateString("diffmaxchange"), NULL));
    /*
     * DerivativeCheck = strcmp(optimget(options,'Derivativecheck'),'on');
     */
    mlfAssign(
      &DerivativeCheck,
      mlfStrcmp(
        mlfOptimget(options, mxCreateString("Derivativecheck"), NULL),
        mxCreateString("on")));
    /*
     * maxFunEvals = optimget(options,'maxfunevals');
     */
    mlfAssign(
      &maxFunEvals, mlfOptimget(options, mxCreateString("maxfunevals"), NULL));
    /*
     * maxIter = optimget(options,'maxiter');
     */
    mlfAssign(&maxIter, mlfOptimget(options, mxCreateString("maxiter"), NULL));
    /*
     * % In case the defaults were gathered from calling: optimset('fminsearch'):
     * %if ischar(maxFunEvals)
     * %   maxFunEvals = eval(maxFunEvals);
     * %end
     * %if ischar(maxIter)
     * %   maxIter = eval(maxIter);
     * %end
     * 
     * iter = 0;
     */
    mlfAssign(&iter, mlfScalar(0.0));
    /*
     * numFunEvals = 0;
     */
    mlfAssign(&numFunEvals, mlfScalar(0.0));
    /*
     * numGradEvals = 0;
     */
    mlfAssign(&numGradEvals, mlfScalar(0.0));
    /*
     * 
     * MATX=zeros(3,1);
     */
    mlfAssign(&MATX, mlfZeros(mlfScalar(3.0), mlfScalar(1.0), NULL));
    /*
     * MATL=[CostFunction'*CostFunction;0;0];
     */
    mlfAssign(
      &MATL,
      mlfVertcat(
        mlfHorzcat(
          mlfMtimes(mlfCtranspose(*CostFunction), *CostFunction), NULL),
        mlfDoubleMatrix(1, 1, __Array0_r, NULL),
        mlfDoubleMatrix(1, 1, __Array1_r, NULL),