optim_private_snls.c

ASUFIT－Matlab－全局拟合程序

             * dd = abs(v); D = sparse(1:n,1:n,full(sqrt(dd))); grad = D*g;
             */
            mlfAssign(&dd, mlfAbs(v));
            mlfAssign(
              &D,
              mlfSparse(
                mlfColon(mlfScalar(1.0), n, NULL),
                mlfColon(mlfScalar(1.0), n, NULL),
                mlfFull(mlfSqrt(dd)),
                NULL,
                NULL,
                NULL));
            mlfAssign(&grad, mlfMtimes(D, g));
            /*
             * vpos(it,1) = posdef; sx = zeros(n,1); theta = max(.95,1-optnrm);  
             */
            mlfIndexAssign(&vpos, "(?,?)", it, mlfScalar(1.0), posdef);
            mlfAssign(&sx, mlfZeros(n, mlfScalar(1.0), NULL));
            mlfAssign(
              &theta,
              mlfMax(
                NULL, mlfScalar(.95), mlfMinus(mlfScalar(1.0), optnrm), NULL));
            /*
             * oposdef = posdef;
             */
            mlfAssign(&oposdef, posdef);
            /*
             * [sx,snod,qp,posdef,pcgit,Z] = trdog(x,g,A,fdata,D,delta,dv,...
             */
            mlfAssign(
              &sx,
              mlfOptim_private_trdog(
                &snod,
                &qp,
                &posdef,
                &pcgit,
                &Z,
                x,
                g,
                A,
                fdata,
                D,
                delta,
                dv,
                mtxmpy,
                pcmtx,
                pcflags,
                pcgtol,
                kmax,
                theta,
                l,
                u,
                Z,
                dnewt,
                mxCreateString("jacobprecon")));
            /*
             * mtxmpy,pcmtx,pcflags,pcgtol,kmax,theta,l,u,Z,dnewt,'jacobprecon');
             * 
             * if isempty(posdef), 
             */
            if (mlfTobool(mlfIsempty(posdef))) {
                /*
                 * posdef = oposdef; 
                 */
                mlfAssign(&posdef, oposdef);
            /*
             * end
             */
            }
            /*
             * nrmsx=norm(snod); 
             */
            mlfAssign(&nrmsx, mlfNorm(snod, NULL));
            /*
             * npcg=npcg+pcgit; 
             */
            mlfAssign(&npcg, mlfPlus(npcg, pcgit));
            /*
             * newx=x+sx;  
             */
            mlfAssign(&newx, mlfPlus(x, sx));
            /*
             * vpcg(it+1,1)=pcgit;
             */
            mlfIndexAssign(
              &vpcg,
              "(?,?)",
              mlfPlus(it, mlfScalar(1.0)),
              mlfScalar(1.0),
              pcgit);
            /*
             * 
             * %       Perturb?
             * [pert,newx] = perturb(newx,l,u);
             */
            mlfAssign(
              &pert,
              mlfOptim_private_perturb(
                &newx, NULL, newx, l, u, NULL, NULL, NULL));
            /*
             * vpos(it+1,1) = posdef;
             */
            mlfIndexAssign(
              &vpos,
              "(?,?)",
              mlfPlus(it, mlfScalar(1.0)),
              mlfScalar(1.0),
              posdef);
            /*
             * 
             * xcurr(:) = newx; % reshape newx for user function
             */
            mlfIndexAssign(&xcurr, "(?)", mlfCreateColonIndex(), newx);
            /*
             * %       Evaluate F and J
             * if ~gradflag %use sparse finite-differencing
             */
            if (mlfTobool(mlfNot(gradflag))) {
                /*
                 * %newfvec = feval(fname,newx,fdata);
                 * switch funfcn{1}
                 */
                mxArray * t_ = mclUnassigned();
                mlfAssign(&t_, mlfIndexRef(funfcn, "{?}", mlfScalar(1.0)));
                /*
                 * case 'fun'
                 */
                if (mclSwitchCompare(t_, mxCreateString("fun"))) {
                    /*
                     * newfvec = feval(funfcn{3},xcurr,varargin{:});
                     */
                    mlfAssign(
                      &newfvec,
                      mlfFeval(
                        mclValueVarargout(),
                        mclFevalLookup(
                          mlfIndexRef(funfcn, "{?}", mlfScalar(3.0)), 0, NULL),
                        xcurr,
                        mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                        NULL));
                    /*
                     * if ~isempty(YDATA)
                     */
                    if (mlfTobool(mlfNot(mlfIsempty(YDATA)))) {
                        /*
                         * newfvec = newfvec - YDATA;
                         */
                        mlfAssign(&newfvec, mlfMinus(newfvec, YDATA));
                    /*
                     * end
                     */
                    }
                    /*
                     * newfvec = newfvec(:);
                     */
                    mlfAssign(
                      &newfvec,
                      mlfIndexRef(newfvec, "(?)", mlfCreateColonIndex()));
                /*
                 * otherwise
                 */
                } else {
                    /*
                     * error(['Undefined calltype in ',caller]);
                     */
                    mlfError(
                      mlfHorzcat(
                        mxCreateString("Undefined calltype in "),
                        caller,
                        NULL));
                /*
                 * end
                 */
                }
                mxDestroyArray(t_);
                /*
                 * % pass in funfcn{3} since we know no gradient
                 * [newA, findiffevals] = sfdnls(xcurr,newfvec,Jstr,group,[],funfcn{3},YDATA,varargin{:});
                 */
                mlfFeval(
                  mlfVarargout(&newA, &findiffevals, NULL),
                  mlxSfdnls,
                  xcurr,
                  newfvec,
                  Jstr,
                  group,
                  mclCreateEmptyArray(),
                  mlfIndexRef(funfcn, "{?}", mlfScalar(3.0)),
                  YDATA,
                  mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                  NULL);
            /*
             * else % use user-supplied detrmination of J or dnewt
             */
            } else {
                /*
                 * findiffevals = 0; % no finite differencing
                 */
                mlfAssign(&findiffevals, mlfScalar(0.0));
                /*
                 * %[newfvec,newA] = feval(fname,newx,fdata);
                 * switch funfcn{1}
                 */
                {
                    mxArray * t_ = mclUnassigned();
                    mlfAssign(&t_, mlfIndexRef(funfcn, "{?}", mlfScalar(1.0)));
                    /*
                     * case 'fungrad'
                     */
                    if (mclSwitchCompare(t_, mxCreateString("fungrad"))) {
                        /*
                         * [newfvec,newA] = feval(funfcn{3},xcurr,varargin{:});
                         */
                        mlfFeval(
                          mlfVarargout(&newfvec, &newA, NULL),
                          mclFevalLookup(
                            mlfIndexRef(funfcn, "{?}", mlfScalar(3.0)),
                            0,
                            NULL),
                          xcurr,
                          mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                          NULL);
                        /*
                         * numGradEvals=numGradEvals+1;
                         */
                        mlfAssign(
                          &numGradEvals, mlfPlus(numGradEvals, mlfScalar(1.0)));
                        /*
                         * if ~isempty(YDATA)
                         */
                        if (mlfTobool(mlfNot(mlfIsempty(YDATA)))) {
                            /*
                             * newfvec = newfvec - YDATA;
                             */
                            mlfAssign(&newfvec, mlfMinus(newfvec, YDATA));
                        /*
                         * end
                         */
                        }
                        /*
                         * newfvec = newfvec(:);
                         */
                        mlfAssign(
                          &newfvec,
                          mlfIndexRef(newfvec, "(?)", mlfCreateColonIndex()));
                    /*
                     * case 'fun_then_grad'
                     */
                    } else if (mclSwitchCompare(
                                 t_, mxCreateString("fun_then_grad"))) {
                        /*
                         * newfvec = feval(funfcn{3},xcurr,varargin{:});
                         */
                        mlfAssign(
                          &newfvec,
                          mlfFeval(
                            mclValueVarargout(),
                            mclFevalLookup(
                              mlfIndexRef(funfcn, "{?}", mlfScalar(3.0)),
                              0,
                              NULL),
                            xcurr,
                            mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                            NULL));
                        /*
                         * if ~isempty(YDATA)
                         */
                        if (mlfTobool(mlfNot(mlfIsempty(YDATA)))) {
                            /*
                             * newfvec = newfvec - YDATA;
                             */
                            mlfAssign(&newfvec, mlfMinus(newfvec, YDATA));
                        /*
                         * end
                         */
                        }
                        /*
                         * newfvec = newfvec(:);
                         */
                        mlfAssign(
                          &newfvec,
                          mlfIndexRef(newfvec, "(?)", mlfCreateColonIndex()));
                        /*
                         * newA = feval(funfcn{4},xcurr,varargin{:});
                         */
                        mlfAssign(
                          &newA,
                          mlfFeval(
                            mclValueVarargout(),
                            mclFevalLookup(
                              mlfIndexRef(funfcn, "{?}", mlfScalar(4.0)),
                              0,
                              NULL),
                            xcurr,
                            mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                            NULL));
                        /*
                         * numGradEvals=numGradEvals+1;
                         */
                        mlfAssign(
                          &numGradEvals, mlfPlus(numGradEvals, mlfScalar(1.0)));
                    /*
                     * otherwise
                     */
                    } else {
                        /*
                         * error(['Undefined calltype in ',caller]);
                         */
                        mlfError(
                          mlfHorzcat(
                            mxCreateString("Undefined calltype in "),
                            caller,
                            NULL));
                    /*
                     * end
                     */
                    }
                    mxDestroyArray(t_);
                }
            /*
             * end
             */
            }
            /*
             * numFunEvals = numFunEvals + 1 + findiffevals;
             */
            mlfAssign(
              &numFunEvals,
              mlfPlus(mlfPlus(numFunEvals, mlfScalar(1.0)), findiffevals));
            /*
             * 
             * [mm,pp] = size(newfvec);
             */
            mlfSize(mlfVarargout(&mm, &pp, NULL), newfvec, NULL);
            /*
             * if mm < n, 
             */
            if (mlfTobool(mlfLt(mm, n))) {
                /*
                 * error('The number of eqns must be greater than n'); end
                 */
                mlfError(
                  mxCreateString("The number of eqns must be greater than n"));
            }
            /*
             * 
             * %       Accept or reject trial point
             * newgrad = feval(mtxmpy,newfvec(:,1),newA,[],-1);  
             */
            mlfAssign(
              &newgrad,
              mlfFeval(
                mclValueVarargout(),
                mclFevalLookup(mtxmpy, 0, NULL),
                mlfIndexRef(
                  newfvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0)),
                newA,
                mclCreateEmptyArray(),
                mlfScalar(-1.0),
                NULL));
            /*
             * newval = newfvec(:,1)'*newfvec(:,1);
             */
            mlfAssign(
              &newval,
              mlfMtimes(
                mlfCtranspose(
                  mlfIndexRef(
                    newfvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0))),
                mlfIndexRef(
                  newfvec, "(?,?)", mlfCreateColonIndex(), mlfScalar(1.0))));
            /*
             * aug = .5*snod'*((dv.*abs(g)).*snod); 
             */
            mlfAssign(
              &aug,
              mlfMtimes(
                mlfMtimes(mlfScalar(.5), mlfCtranspose(snod)),
                mlfTimes(mlfTimes(dv, mlfAbs(g)), snod)));
            /*
             * ratio = (0.5*(newval-val)+aug)/qp;     % we compute val = f'*f, not val = 0.5*(f'*f)
             */
            mlfAssign(
              &ratio,
              mlfMrdivide(
                mlfPlus(mlfMtimes(mlfScalar(0.5), mlfMinus(newval, val)), aug),
                qp));
            /*
             * if (ratio >= .75) & (nrmsx >= .9*delta),
             */
            {
                mxArray * a_ = mclInitialize(mlfGe(ratio, mlfScalar(.75)));
                if (mlfTobool(a_)
                    && mlfTobool(
                         mlfAnd(
                           a_,
                           mlfGe(nrmsx, mlfMtimes(mlfScalar(.9), delta))))) {
                    mxDestroyArray(a_);
                    /*
                     * delta = min(delbnd,2*delta);
                     */
                    mlfAssign(
                      &delta,
                      mlfMin(
                        NULL, delbnd, mlfMtimes(mlfScalar(2.0), delta), NULL));
                /*
                 * elseif ratio <= .25,