optim_private_trdog.c

ASUFIT－Matlab－全局拟合程序

      mlfFull(
        mlfPlus(
          mlfMtimes(mlfCtranspose(*Z), W),
          mlfMtimes(mlfMtimes(mlfCtranspose(*Z), DG), *Z))));
    /*
     * rhs=full(Z'*grad);
     */
    mlfAssign(&rhs, mlfFull(mlfMtimes(mlfCtranspose(*Z), grad)));
    /*
     * 
     * %  Determine 2-D TR soln
     * [st,qpval,po,fcnt,lambda] = trust(rhs,MM,delta);
     */
    mlfAssign(
      &st, mlfOptim_private_trust(qpval, &po, &fcnt, &lambda, rhs, MM, delta));
    /*
     * ss = Z*st;  
     */
    mlfAssign(&ss, mlfMtimes(*Z, st));
    /*
     * s = abs(diag(D)).*ss; 
     */
    mlfAssign(&s, mlfTimes(mlfAbs(mlfDiag(D, NULL)), ss));
    /*
     * s = full(s); 
     */
    mlfAssign(&s, mlfFull(s));
    /*
     * ssave = s;
     */
    mlfAssign(&ssave, s);
    /*
     * sssave = ss; 
     */
    mlfAssign(&sssave, ss);
    /*
     * stsave = st;
     */
    mlfAssign(&stsave, st);
    /*
     * 
     * % Truncate the TR solution?
     * arg = (abs(s) > 0);
     */
    mlfAssign(&arg, mlfGt(mlfAbs(s), mlfScalar(0.0)));
    /*
     * if isnan(s)
     */
    if (mlfTobool(mlfIsnan(s))) {
        /*
         * error('Trust region step contains NaN''s.')
         */
        mlfError(mxCreateString("Trust region step contains NaN's."));
    /*
     * end
     */
    }
    /*
     * % No truncation if s is zero length
     * if isempty(find(arg))   
     */
    if (mlfTobool(mlfIsempty(mlfFind(NULL, NULL, arg)))) {
        /*
         * alpha = 1;
         */
        mlfAssign(&alpha, mlfScalar(1.0));
    /*
     * else
     */
    } else {
        /*
         * mdis = inf;
         */
        mlfAssign(&mdis, mlfInf());
        /*
         * dis = max((u(arg)-x(arg))./s(arg), (l(arg)-x(arg))./s(arg));
         */
        mlfAssign(
          &dis,
          mlfMax(
            NULL,
            mlfRdivide(
              mlfMinus(mlfIndexRef(u, "(?)", arg), mlfIndexRef(x, "(?)", arg)),
              mlfIndexRef(s, "(?)", arg)),
            mlfRdivide(
              mlfMinus(mlfIndexRef(l, "(?)", arg), mlfIndexRef(x, "(?)", arg)),
              mlfIndexRef(s, "(?)", arg)),
            NULL));
        /*
         * [mmdis,ipt] = min(dis);  
         */
        mlfAssign(&mmdis, mlfMin(&ipt, dis, NULL, NULL));
        /*
         * mdis = theta*mmdis;
         */
        mlfAssign(&mdis, mlfMtimes(theta, mmdis));
        /*
         * alpha = min(1,mdis);
         */
        mlfAssign(&alpha, mlfMin(NULL, mlfScalar(1.0), mdis, NULL));
    /*
     * end
     */
    }
    /*
     * s = alpha*s; 
     */
    mlfAssign(&s, mlfMtimes(alpha, s));
    /*
     * st = alpha*st; 
     */
    mlfAssign(&st, mlfMtimes(alpha, st));
    /*
     * ss = full(alpha*ss);
     */
    mlfAssign(&ss, mlfFull(mlfMtimes(alpha, ss)));
    /*
     * qpval1 = rhs'*st + (.5*st)'*MM*st;
     */
    mlfAssign(
      &qpval1,
      mlfPlus(
        mlfMtimes(mlfCtranspose(rhs), st),
        mlfMtimes(
          mlfMtimes(mlfCtranspose(mlfMtimes(mlfScalar(.5), st)), MM), st)));
    /*
     * if n > 1 
     */
    if (mlfTobool(mlfGt(n, mlfScalar(1.0)))) {
        /*
         * %   Evaluate along the reflected direction?
         * qpval3 = inf; 
         */
        mlfAssign(&qpval3, mlfInf());
        /*
         * ssssave = mmdis*sssave;
         */
        mlfAssign(&ssssave, mlfMtimes(mmdis, sssave));
        /*
         * if norm(ssssave) < .9*delta
         */
        if (mlfTobool(
              mlfLt(mlfNorm(ssssave, NULL), mlfMtimes(mlfScalar(.9), delta)))) {
            /*
             * r = mmdis*ssave; 
             */
            mlfAssign(&r, mlfMtimes(mmdis, ssave));
            /*
             * ns = ssave; 
             */
            mlfAssign(&ns, ssave);
            /*
             * ns(ipt) = -ns(ipt); 
             */
            mlfIndexAssign(
              &ns, "(?)", ipt, mlfUminus(mlfIndexRef(ns, "(?)", ipt)));
            /*
             * nx = x+r;
             */
            mlfAssign(&nx, mlfPlus(x, r));
            /*
             * stsave = mmdis*stsave;
             */
            mlfAssign(&stsave, mlfMtimes(mmdis, stsave));
            /*
             * qpval0 = rhs'*stsave + (.5*stsave)'*MM*stsave;
             */
            mlfAssign(
              &qpval0,
              mlfPlus(
                mlfMtimes(mlfCtranspose(rhs), stsave),
                mlfMtimes(
                  mlfMtimes(
                    mlfCtranspose(mlfMtimes(mlfScalar(.5), stsave)), MM),
                  stsave)));
            /*
             * ng = feval(mtxmpy,r,H,fdata); 
             */
            mlfAssign(
              &ng,
              mlfFeval(
                mclValueVarargout(),
                mclFevalLookup(mtxmpy, 0, NULL),
                r,
                H,
                fdata,
                NULL));
            /*
             * ng = ng + g; 
             */
            mlfAssign(&ng, mlfPlus(ng, g));
            /*
             * ngrad = D*ng;
             */
            mlfAssign(&ngrad, mlfMtimes(D, ng));
            /*
             * ngrad = ngrad + DG*ssssave;
             */
            mlfAssign(&ngrad, mlfPlus(ngrad, mlfMtimes(DG, ssssave)));
            /*
             * 
             * %      nss is the reflected direction
             * nss = sssave; 
             */
            mlfAssign(&nss, sssave);
            /*
             * nss(ipt) = -nss(ipt); 
             */
            mlfIndexAssign(
              &nss, "(?)", ipt, mlfUminus(mlfIndexRef(nss, "(?)", ipt)));
            /*
             * ZZ(:,1) = nss/norm(nss);
             */
            mlfIndexAssign(
              &ZZ,
              "(?,?)",
              mlfCreateColonIndex(),
              mlfScalar(1.0),
              mlfMrdivide(nss, mlfNorm(nss, NULL)));
            /*
             * W = DM*ZZ; 
             */
            mlfAssign(&W, mlfMtimes(DM, ZZ));
            /*
             * WW = feval(mtxmpy,W,H,fdata); 
             */
            mlfAssign(
              &WW,
              mlfFeval(
                mclValueVarargout(),
                mclFevalLookup(mtxmpy, 0, NULL),
                W,
                H,
                fdata,
                NULL));
            /*
             * W = DM*WW;
             */
            mlfAssign(&W, mlfMtimes(DM, WW));
            /*
             * MM = full(ZZ'*W + ZZ'*DG*ZZ);
             */
            mlfAssign(
              &MM,
              mlfFull(
                mlfPlus(
                  mlfMtimes(mlfCtranspose(ZZ), W),
                  mlfMtimes(mlfMtimes(mlfCtranspose(ZZ), DG), ZZ))));
            /*
             * nrhs=full(ZZ'*ngrad);
             */
            mlfAssign(&nrhs, mlfFull(mlfMtimes(mlfCtranspose(ZZ), ngrad)));
            /*
             * [nss,tau] = quad1d(nss,ssssave,delta); 
             */
            mlfAssign(&nss, mlfOptim_private_quad1d(&tau, nss, ssssave, delta));
            /*
             * nst = tau/norm(nss);
             */
            mlfAssign(&nst, mlfMrdivide(tau, mlfNorm(nss, NULL)));
            /*
             * ns = abs(diag(D)).*nss; 
             */
            mlfAssign(&ns, mlfTimes(mlfAbs(mlfDiag(D, NULL)), nss));
            /*
             * ns = full(ns);
             */
            mlfAssign(&ns, mlfFull(ns));
            /*
             * 
             * %      Truncate the reflected direction?
             * arg = (abs(ns) > 0); 
             */
            mlfAssign(&arg, mlfGt(mlfAbs(ns), mlfScalar(0.0)));
            /*
             * if isnan(ns)
             */
            if (mlfTobool(mlfIsnan(ns))) {
                /*
                 * error('Reflected trust region step contains NaN''s.')
                 */
                mlfError(
                  mxCreateString(
                    "Reflected trust region step contains NaN's."));
            /*
             * end
             */
            }
            /*
             * % No truncation if s is zero length
             * if isempty(find(arg))   
             */
            if (mlfTobool(mlfIsempty(mlfFind(NULL, NULL, arg)))) {
                /*
                 * alpha = 1;
                 */
                mlfAssign(&alpha, mlfScalar(1.0));
            /*
             * else
             */
            } else {
                /*
                 * mdis = inf;
                 */
                mlfAssign(&mdis, mlfInf());
                /*
                 * dis = max((u(arg)-nx(arg))./ns(arg), (l(arg)-nx(arg))./ns(arg));
                 */
                mlfAssign(
                  &dis,
                  mlfMax(
                    NULL,
                    mlfRdivide(
                      mlfMinus(
                        mlfIndexRef(u, "(?)", arg),
                        mlfIndexRef(nx, "(?)", arg)),
                      mlfIndexRef(ns, "(?)", arg)),
                    mlfRdivide(
                      mlfMinus(
                        mlfIndexRef(l, "(?)", arg),
                        mlfIndexRef(nx, "(?)", arg)),
                      mlfIndexRef(ns, "(?)", arg)),
                    NULL));
                /*
                 * mdis = min(dis);  
                 */
                mlfAssign(&mdis, mlfMin(NULL, dis, NULL, NULL));
                /*
                 * mdis = theta*mdis;
                 */
                mlfAssign(&mdis, mlfMtimes(theta, mdis));
                /*
                 * alpha = min(1,mdis); 
                 */
                mlfAssign(&alpha, mlfMin(NULL, mlfScalar(1.0), mdis, NULL));
            /*
             * end
             */
            }
            /*
             * ns = alpha*ns; 
             */
            mlfAssign(&ns, mlfMtimes(alpha, ns));
            /*
             * nst = alpha*nst; 
             */
            mlfAssign(&nst, mlfMtimes(alpha, nst));
            /*
             * nss = full(alpha*nss);
             */
            mlfAssign(&nss, mlfFull(mlfMtimes(alpha, nss)));
            /*
             * qpval3 = qpval0 +  nrhs'*nst + (.5*nst)'*MM*nst;
             */
            mlfAssign(
              &qpval3,
              mlfPlus(
                mlfPlus(qpval0, mlfMtimes(mlfCtranspose(nrhs), nst)),
                mlfMtimes(
                  mlfMtimes(mlfCtranspose(mlfMtimes(mlfScalar(.5), nst)), MM),
                  nst)));
        /*
         * end
         */
        }
        /*
         * 
         * %   Evaluate along gradient direction
         * ZZ(:,1) = grad/norm(grad);
         */
        mlfIndexAssign(
          &ZZ,
          "(?,?)",
          mlfCreateColonIndex(),
          mlfScalar(1.0),
          mlfMrdivide(grad, mlfNorm(grad, NULL)));
        /*
         * W = DM*ZZ; 
         */
        mlfAssign(&W, mlfMtimes(DM, ZZ));
        /*
         * WW = feval(mtxmpy,W,H,fdata); 
         */
        mlfAssign(
          &WW,
          mlfFeval(
            mclValueVarargout(),
            mclFevalLookup(mtxmpy, 0, NULL),
            W,
            H,
            fdata,
            NULL));
        /*
         * W = DM*WW;
         */
        mlfAssign(&W, mlfMtimes(DM, WW));
        /*
         * MM = full(ZZ'*W + ZZ'*DG*ZZ); 
         */
        mlfAssign(
          &MM,
          mlfFull(
            mlfPlus(
              mlfMtimes(mlfCtranspose(ZZ), W),
              mlfMtimes(mlfMtimes(mlfCtranspose(ZZ), DG), ZZ))));
        /*
         * rhs=full(ZZ'*grad);
         */
        mlfAssign(&rhs, mlfFull(mlfMtimes(mlfCtranspose(ZZ), grad)));
        /*
         * [st,qpval,po,fcnt,lambda] = trust(rhs,MM,delta);
         */
        mlfAssign(
          &st,
          mlfOptim_private_trust(qpval, &po, &fcnt, &lambda, rhs, MM, delta));
        /*
         * ssg = ZZ*st; 
         */
        mlfAssign(&ssg, mlfMtimes(ZZ, st));
        /*
         * sg = abs(diag(D)).*ssg; 
         */
        mlfAssign(&sg, mlfTimes(mlfAbs(mlfDiag(D, NULL)), ssg));
        /*
         * sg = full(sg);
         */
        mlfAssign(&sg, mlfFull(sg));
        /*
         * 
         * %   Truncate the gradient direction?
         * arg = (abs(sg) > 0); 
         */
        mlfAssign(&arg, mlfGt(mlfAbs(sg), mlfScalar(0.0)));
        /*
         * if isnan(arg)
         */
        if (mlfTobool(mlfIsnan(arg))) {
            /*
             * % No truncation if s is zero length
             * error('Gradient step contains NaN''s.')
             */
            mlfError(mxCreateString("Gradient step contains NaN's."));
        /*
         * end
         */
        }
        /*
         * if  isempty(find(arg))   
         */
        if (mlfTobool(mlfIsempty(mlfFind(NULL, NULL, arg)))) {