optim_private_snls.c

ASUFIT－Matlab－全局拟合程序

     * %if ischar(typx)
     * %   typx = eval(typx);
     * %end
     * 
     * % This will be user-settable in the future:
     * mtxmpy = optimget(options,'JacobMult','atamult') ;
     */
    mlfAssign(
      &mtxmpy,
      mlfOptimget(
        options, mxCreateString("JacobMult"), mxCreateString("atamult")));
    /*
     * pcmtx = optimget(options,'Preconditioner','aprecon') ;
     */
    mlfAssign(
      &pcmtx,
      mlfOptimget(
        options, mxCreateString("Preconditioner"), mxCreateString("aprecon")));
    /*
     * showstat = optimget(options,'showstatus','off');
     */
    mlfAssign(
      &showstat,
      mlfOptimget(
        options, mxCreateString("showstatus"), mxCreateString("off")));
    /*
     * switch showstat
     * case 'iter'
     */
    if (mclSwitchCompare(showstat, mxCreateString("iter"))) {
        /*
         * showstat = 2;
         */
        mlfAssign(&showstat, mlfScalar(2.0));
    /*
     * case {'none','off'}
     */
    } else if (mclSwitchCompare(
                 showstat,
                 mlfCellhcat(
                   mxCreateString("none"), mxCreateString("off"), NULL))) {
        /*
         * showstat = 0;
         */
        mlfAssign(&showstat, mlfScalar(0.0));
    /*
     * case 'final'
     */
    } else if (mclSwitchCompare(showstat, mxCreateString("final"))) {
        /*
         * showstat = 1;
         */
        mlfAssign(&showstat, mlfScalar(1.0));
    /*
     * case 'iterplusbounds'  % if no finite bounds, this is the same as 'iter'
     */
    } else if (mclSwitchCompare(showstat, mxCreateString("iterplusbounds"))) {
        /*
         * showstat = 3;
         */
        mlfAssign(&showstat, mlfScalar(3.0));
    /*
     * otherwise
     */
    } else {
        /*
         * showstat = 1;
         */
        mlfAssign(&showstat, mlfScalar(1.0));
    /*
     * end
     */
    }
    /*
     * 
     * pcflags = optimget(options,'PrecondBandWidth') ;
     */
    mlfAssign(
      &pcflags, mlfOptimget(options, mxCreateString("PrecondBandWidth"), NULL));
    /*
     * tol2 = optimget(options,'tolx') ;
     */
    mlfAssign(&tol2, mlfOptimget(options, mxCreateString("tolx"), NULL));
    /*
     * tol1 = optimget(options,'tolfun') ;
     */
    mlfAssign(&tol1, mlfOptimget(options, mxCreateString("tolfun"), NULL));
    /*
     * tol = tol1;
     */
    mlfAssign(&tol, tol1);
    /*
     * itb = optimget(options,'maxiter') ;
     */
    mlfAssign(&itb, mlfOptimget(options, mxCreateString("maxiter"), NULL));
    /*
     * maxfunevals = optimget(options,'maxFunEvals') ;
     */
    mlfAssign(
      &maxfunevals, mlfOptimget(options, mxCreateString("maxFunEvals"), NULL));
    /*
     * pcgtol = optimget(options,'TolPCG', 0.1) ;  % pcgtol = .1;
     */
    mlfAssign(
      &pcgtol, mlfOptimget(options, mxCreateString("TolPCG"), mlfScalar(0.1)));
    /*
     * kmax = optimget(options,'MaxPCGIter', max(1,floor(n/2))) ;
     */
    mlfAssign(
      &kmax,
      mlfOptimget(
        options,
        mxCreateString("MaxPCGIter"),
        mlfMax(
          NULL,
          mlfScalar(1.0),
          mlfFloor(mlfMrdivide(n, mlfScalar(2.0))),
          NULL)));
    /*
     * %if ischar(kmax)
     * %  kmax = eval(kmax);
     * %end
     * %if ischar(maxfunevals)
     * %   maxfunevals = eval(maxfunevals);
     * %end
     * 
     * ex = 0; posdef = 1; npcg = 0; pcgit = 0;
     */
    mlfAssign(&ex, mlfScalar(0.0));
    mlfAssign(&posdef, mlfScalar(1.0));
    mlfAssign(&npcg, mlfScalar(0.0));
    mlfAssign(&pcgit, mlfScalar(0.0));
    /*
     * if strcmp(optimget(options,'DerivativeCheck'),'on')
     */
    if (mlfTobool(
          mlfStrcmp(
            mlfOptimget(options, mxCreateString("DerivativeCheck"), NULL),
            mxCreateString("on")))) {
        /*
         * warnstr = sprintf('%s\n%s\n', ...
         */
        mlfAssign(
          &warnstr,
          mlfSprintf(
            NULL,
            mxCreateString("%s\\n%s\\n"),
            mxCreateString(
              "Trust region algorithm does not curren"
              "tly check user-supplied derivatives,"),
            mxCreateString("  ignoring OPTIONS.DerivativeCheck"),
            NULL));
        /*
         * 'Trust region algorithm does not currently check user-supplied derivatives,', ...
         * '  ignoring OPTIONS.DerivativeCheck');
         * warning(warnstr);
         */
        mclAssignAns(&ans, mlfWarning(NULL, warnstr));
    /*
     * end
     */
    }
    /*
     * 
     * %tol1 = tol; tol2 = sqrt(tol1)/10;
     * vpcg(1,1) = 0; vpos(1,1) = 1;
     */
    mlfIndexAssign(
      &vpcg, "(?,?)", mlfScalar(1.0), mlfScalar(1.0), mlfScalar(0.0));
    mlfIndexAssign(
      &vpos, "(?,?)", mlfScalar(1.0), mlfScalar(1.0), mlfScalar(1.0));
    /*
     * delta = 10;nrmsx = 1; ratio = 0; 
     */
    mlfAssign(&delta, mlfScalar(10.0));
    mlfAssign(&nrmsx, mlfScalar(1.0));
    mlfAssign(&ratio, mlfScalar(0.0));
    /*
     * degen = inf; DS = speye(n);
     */
    mlfAssign(&degen, mlfInf());
    mlfAssign(&DS, mlfSpeye(n, NULL));
    /*
     * v = zeros(n,1); dv = ones(n,1); del = 10*eps;
     */
    mlfAssign(&v, mlfZeros(n, mlfScalar(1.0), NULL));
    mlfAssign(&dv, mlfOnes(n, mlfScalar(1.0), NULL));
    mlfAssign(&del, mlfMtimes(mlfScalar(10.0), mlfEps()));
    /*
     * x = xstart; oval = inf;
     */
    mlfAssign(&x, xstart);
    mlfAssign(&oval, mlfInf());
    /*
     * g = zeros(n,1); nbnds = 1; Z = [];
     */
    mlfAssign(&g, mlfZeros(n, mlfScalar(1.0), NULL));
    mlfAssign(&nbnds, mlfScalar(1.0));
    mlfAssign(&Z, mclCreateEmptyArray());
    /*
     * if isinf(u) & isinf(l)
     */
    {
        mxArray * a_ = mclInitialize(mlfIsinf(u));
        if (mlfTobool(a_) && mlfTobool(mlfAnd(a_, mlfIsinf(l)))) {
            mxDestroyArray(a_);
            /*
             * nbnds = 0; degen = -1;
             */
            mlfAssign(&nbnds, mlfScalar(0.0));
            mlfAssign(&degen, mlfScalar(-1.0));
        } else {
            mxDestroyArray(a_);
        }
    /*
     * end
     */
    }
    /*
     * 
     * xcurr(:) = xstart;
     */
    mlfIndexAssign(&xcurr, "(?)", mlfCreateColonIndex(), xstart);
    /*
     * %   Evaluate F and J
     * if ~gradflag % use sparse finite differencing
     */
    if (mlfTobool(mlfNot(gradflag))) {
        /*
         * fvec = fval;
         */
        mlfAssign(fvec, fval);
        /*
         * % Determine coloring/grouping for sparse finite-differencing
         * p = colmmd(Jstr)'; 
         */
        mlfAssign(&p, mlfCtranspose(mlfColmmd(Jstr)));
        /*
         * p = (n+1)*ones(n,1)-p; 
         */
        mlfAssign(
          &p,
          mlfMinus(
            mlfMtimes(
              mlfPlus(n, mlfScalar(1.0)), mlfOnes(n, mlfScalar(1.0), NULL)),
            p));
        /*
         * group = color(Jstr,p);
         */
        mlfAssign(&group, mlfColor(Jstr, p));
        /*
         * xcurr(:) = x;  % reshape x for user function
         */
        mlfIndexAssign(&xcurr, "(?)", mlfCreateColonIndex(), x);
        /*
         * % pass in funfcn{3} since we know no gradient: not funfcn{4}
         * [A,findiffevals] = sfdnls(xcurr,fvec,Jstr,group,[],funfcn{3},YDATA,varargin{:});
         */
        mlfFeval(
          mlfVarargout(&A, &findiffevals, NULL),
          mlxSfdnls,
          xcurr,
          *fvec,
          Jstr,
          group,
          mclCreateEmptyArray(),
          mlfIndexRef(funfcn, "{?}", mlfScalar(3.0)),
          YDATA,
          mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
          NULL);
    /*
     * else % user-supplied computation of J or dnewt 
     */
    } else {
        /*
         * %[fvec,A] = feval(fname,x,fdata);
         * fvec = fval;
         */
        mlfAssign(fvec, fval);
        /*
         * A = JACval;
         */
        mlfAssign(&A, JACval);
        /*
         * findiffevals = 0;
         */
        mlfAssign(&findiffevals, mlfScalar(0.0));
    /*
     * end
     */
    }
    /*
     * numFunEvals = numFunEvals + 1 + findiffevals;
     */
    mlfAssign(
      &numFunEvals,
      mlfPlus(mlfPlus(numFunEvals, mlfScalar(1.0)), findiffevals));
    /*
     * 
     * 
     * delbnd = max(100*norm(xstart),1);
     */
    mlfAssign(
      &delbnd,
      mlfMax(
        NULL,
        mlfMtimes(mlfScalar(100.0), mlfNorm(xstart, NULL)),
        mlfScalar(1.0),
        NULL));
    /*
     * 
     * [mm,pp] = size(fvec);
     */
    mlfSize(mlfVarargout(&mm, &pp, NULL), *fvec, NULL);
    /*
     * if mm < n, 
     */
    if (mlfTobool(mlfLt(mm, n))) {
        /*
         * error('The number of equations must not be less than n'); end
         */
        mlfError(
          mxCreateString("The number of equations must not be less than n"));
    }
    /*
     * 
     * %   Extract the Newton direction?
     * if pp == 2, 
     */
    if (mlfTobool(mlfEq(pp, mlfScalar(2.0)))) {
        /*
         * dnewt = fvec(1:n,2); 
         */
        mlfAssign(
          &dnewt,
          mlfIndexRef(
            *fvec, "(?,?)", mlfColon(mlfScalar(1.0), n, NULL), mlfScalar(2.0)));
    /*
     * end
     */
    }
    /*
     * 
     * %   Determine gradient of the nonlinear least squares function
     * g = feval(mtxmpy,fvec(:,1),A,[],-1);
     */
    mlfAssign(
      &g,
      mlfFeval(
        mclValueVarargout(),
        mclFevalLookup(mtxmpy, 0, NULL),
        mlfIndexRef(*fvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0)),
        A,
        mclCreateEmptyArray(),
        mlfScalar(-1.0),
        NULL));
    /*
     * 
     * %   Evaluate F (initial point)
     * val = fvec(:,1)'*fvec(:,1); vval(it,1) = val;
     */
    mlfAssign(
      &val,
      mlfMtimes(
        mlfCtranspose(
          mlfIndexRef(*fvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0))),
        mlfIndexRef(*fvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0))));
    mlfIndexAssign(&vval, "(?,?)", it, mlfScalar(1.0), val);
    /*
     * 
     * %   Display
     * if showstat > 1 
     */
    if (mlfTobool(mlfGt(showstat, mlfScalar(1.0)))) {
        /*
         * figtr=display1('init',itb,tol,verb,nbnds,x,g,l,u);
         */
        mlfAssign(
          &figtr,
          mlfOptim_private_display1(
            mclValueVarargout(),
            mxCreateString("init"),
            itb,
            tol,
            verb,
            nbnds,
            x,
            g,
            l,
            u,
            NULL));
    /*
     * end
     */
    }
    /*
     * if verb > 1
     */
    if (mlfTobool(mlfGt(verb, mlfScalar(1.0)))) {
        /*
         * disp(header)
         */
        mlfDisp(header);
    /*
     * end
     */
    }