optim_private_trust.c

ASUFIT－Matlab－全局拟合程序

     * fc = fb;
     */
    mlfAssign(&fc, fb);
    /*
     * % Main loop, exit from middle of the loop
     * while fb ~= 0
     */
    while (mlfTobool(mlfNe(fb, mlfScalar(0.0)))) {
        /*
         * % Insure that b is the best result so far, a is the previous
         * % value of b, and c is on the opposite of the zero from b.
         * if (fb > 0) == (fc > 0)
         */
        if (mlfTobool(
              mlfEq(mlfGt(fb, mlfScalar(0.0)), mlfGt(fc, mlfScalar(0.0))))) {
            /*
             * c = a;  fc = fa;
             */
            mlfAssign(c, a);
            mlfAssign(&fc, fa);
            /*
             * d = b - a;  e = d;
             */
            mlfAssign(&d, mlfMinus(b, a));
            mlfAssign(&e, d);
        /*
         * end
         */
        }
        /*
         * if abs(fc) < abs(fb)
         */
        if (mlfTobool(mlfLt(mlfAbs(fc), mlfAbs(fb)))) {
            /*
             * a = b;    b = c;    c = a;
             */
            mlfAssign(&a, b);
            mlfAssign(&b, *c);
            mlfAssign(c, a);
            /*
             * fa = fb;  fb = fc;  fc = fa;
             */
            mlfAssign(&fa, fb);
            mlfAssign(&fb, fc);
            mlfAssign(&fc, fa);
        /*
         * end
         */
        }
        /*
         * 
         * % Convergence test and possible exit
         * %
         * if itfun > itbnd, break; end
         */
        if (mlfTobool(mlfGt(*itfun, itbnd))) {
            break;
        }
        /*
         * m = 0.5*(c - b);
         */
        mlfAssign(&m, mlfMtimes(mlfScalar(0.5), mlfMinus(*c, b)));
        /*
         * toler = 2.0*tol*max(abs(b),1.0);
         */
        mlfAssign(
          &toler,
          mlfMtimes(
            mlfMtimes(mlfScalar(2.0), tol),
            mlfMax(NULL, mlfAbs(b), mlfScalar(1.0), NULL)));
        /*
         * if (abs(m) <= toler) + (fb == 0.0), break, end
         */
        if (mlfTobool(
              mlfPlus(mlfLe(mlfAbs(m), toler), mlfEq(fb, mlfScalar(0.0))))) {
            break;
        }
        /*
         * 
         * % Choose bisection or interpolation
         * if (abs(e) < toler) + (abs(fa) <= abs(fb))
         */
        if (mlfTobool(
              mlfPlus(
                mlfLt(mlfAbs(e), toler), mlfLe(mlfAbs(fa), mlfAbs(fb))))) {
            /*
             * % Bisection
             * d = m;  e = m;
             */
            mlfAssign(&d, m);
            mlfAssign(&e, m);
        /*
         * else
         */
        } else {
            /*
             * % Interpolation
             * s = fb/fa;
             */
            mlfAssign(&s, mlfMrdivide(fb, fa));
            /*
             * if (a == c)
             */
            if (mlfTobool(mlfEq(a, *c))) {
                /*
                 * % Linear interpolation
                 * p = 2.0*m*s;
                 */
                mlfAssign(&p, mlfMtimes(mlfMtimes(mlfScalar(2.0), m), s));
                /*
                 * q = 1.0 - s;
                 */
                mlfAssign(&q, mlfMinus(mlfScalar(1.0), s));
            /*
             * else
             */
            } else {
                /*
                 * % Inverse quadratic interpolation
                 * q = fa/fc;
                 */
                mlfAssign(&q, mlfMrdivide(fa, fc));
                /*
                 * r = fb/fc;
                 */
                mlfAssign(&r, mlfMrdivide(fb, fc));
                /*
                 * p = s*(2.0*m*q*(q - r) - (b - a)*(r - 1.0));
                 */
                mlfAssign(
                  &p,
                  mlfMtimes(
                    s,
                    mlfMinus(
                      mlfMtimes(
                        mlfMtimes(mlfMtimes(mlfScalar(2.0), m), q),
                        mlfMinus(q, r)),
                      mlfMtimes(mlfMinus(b, a), mlfMinus(r, mlfScalar(1.0))))));
                /*
                 * q = (q - 1.0)*(r - 1.0)*(s - 1.0);
                 */
                mlfAssign(
                  &q,
                  mlfMtimes(
                    mlfMtimes(
                      mlfMinus(q, mlfScalar(1.0)), mlfMinus(r, mlfScalar(1.0))),
                    mlfMinus(s, mlfScalar(1.0))));
            /*
             * end;
             */
            }
            /*
             * if p > 0, q = -q; else p = -p; end;
             */
            if (mlfTobool(mlfGt(p, mlfScalar(0.0)))) {
                mlfAssign(&q, mlfUminus(q));
            } else {
                mlfAssign(&p, mlfUminus(p));
            }
            /*
             * % Is interpolated point acceptable
             * if (2.0*p < 3.0*m*q - abs(toler*q)) * (p < abs(0.5*e*q))
             */
            if (mlfTobool(
                  mlfMtimes(
                    mlfLt(
                      mlfMtimes(mlfScalar(2.0), p),
                      mlfMinus(
                        mlfMtimes(mlfMtimes(mlfScalar(3.0), m), q),
                        mlfAbs(mlfMtimes(toler, q)))),
                    mlfLt(
                      p,
                      mlfAbs(mlfMtimes(mlfMtimes(mlfScalar(0.5), e), q)))))) {
                /*
                 * e = d;  d = p/q;
                 */
                mlfAssign(&e, d);
                mlfAssign(&d, mlfMrdivide(p, q));
            /*
             * else
             */
            } else {
                /*
                 * d = m;  e = m;
                 */
                mlfAssign(&d, m);
                mlfAssign(&e, m);
            /*
             * end;
             */
            }
        /*
         * end % Interpolation
         */
        }
        /*
         * 
         * % Next point
         * a = b;
         */
        mlfAssign(&a, b);
        /*
         * fa = fb;
         */
        mlfAssign(&fa, fb);
        /*
         * if abs(d) > toler, b = b + d;
         */
        if (mlfTobool(mlfGt(mlfAbs(d), toler))) {
            mlfAssign(&b, mlfPlus(b, d));
        /*
         * else if b > c, b = b - toler;
         */
        } else {
            if (mlfTobool(mlfGt(b, *c))) {
                mlfAssign(&b, mlfMinus(b, toler));
            /*
             * else b = b + toler;
             */
            } else {
                mlfAssign(&b, mlfPlus(b, toler));
            /*
             * end
             */
            }
        /*
         * end
         */
        }
        /*
         * fb = feval(FunFcn,b,eigval,alpha,delta);
         */
        mlfAssign(
          &fb,
          mlfFeval(
            mclValueVarargout(),
            mclFevalLookup(FunFcn, 2, local_function_table_),
            b,
            eigval,
            alpha,
            delta,
            NULL));
        /*
         * itfun = itfun + 1;
         */
        mlfAssign(itfun, mlfPlus(*itfun, mlfScalar(1.0)));
        /*
         * if trace
         */
        if (mlfTobool(trace)) {
            /*
             * %home
             * step = [b fb]
             */
            mlfAssign(&step, mlfHorzcat(b, fb, NULL));
            mclPrintArray(step, "step");
        /*
         * end
         */
        }
    /*
     * end % Main loop
     */
    }
    mclValidateOutputs("trust/rfzero", 3, nargout_, &b, c, itfun);
    mxDestroyArray(a);
    mxDestroyArray(d);
    mxDestroyArray(dx);
    mxDestroyArray(e);
    mxDestroyArray(fa);
    mxDestroyArray(fb);
    mxDestroyArray(fc);
    mxDestroyArray(init);
    mxDestroyArray(m);
    mxDestroyArray(nargin_);
    mxDestroyArray(p);
    mxDestroyArray(q);
    mxDestroyArray(r);
    mxDestroyArray(s);
    mxDestroyArray(sign);
    mxDestroyArray(step);
    mxDestroyArray(tol);
    mxDestroyArray(toler);
    mxDestroyArray(trace);
    /*
     * 
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * 
     * 
     * 
     * 
     */
    return b;
}

/*
 * The function "mlfTrust_rfzero" contains the normal interface for the
 * "trust/rfzero" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 138-290). This function
 * processes any input arguments and passes them to the implementation version
 * of the function, appearing above.
 */
static mxArray * mlfTrust_rfzero(mxArray * * c,
                                 mxArray * * itfun,
                                 mxArray * FunFcn,
                                 mxArray * x,
                                 mxArray * itbnd,
                                 mxArray * eigval,
                                 mxArray * alpha,
                                 mxArray * delta,
                                 mxArray * tol,
                                 mxArray * trace) {
    int nargout = 1;
    mxArray * b = mclGetUninitializedArray();
    mxArray * c__ = mclGetUninitializedArray();
    mxArray * itfun__ = mclGetUninitializedArray();
    mlfEnterNewContext(
      2, 8, c, itfun, FunFcn, x, itbnd, eigval, alpha, delta, tol, trace);
    if (c != NULL) {
        ++nargout;
    }
    if (itfun != NULL) {
        ++nargout;
    }
    b
      = Mtrust_rfzero(
          &c__,
          &itfun__,
          nargout,
          FunFcn,
          x,
          itbnd,
          eigval,
          alpha,
          delta,
          tol,
          trace);
    mlfRestorePreviousContext(
      2, 8, c, itfun, FunFcn, x, itbnd, eigval, alpha, delta, tol, trace);
    if (c != NULL) {
        mclCopyOutputArg(c, c__);
    } else {
        mxDestroyArray(c__);
    }
    if (itfun != NULL) {
        mclCopyOutputArg(itfun, itfun__);
    } else {
        mxDestroyArray(itfun__);
    }
    return mlfReturnValue(b);
}

/*
 * The function "mlxTrust_rfzero" contains the feval interface for the
 * "trust/rfzero" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 138-290). The feval
 * function calls the implementation version of trust/rfzero through this
 * function. This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
static void mlxTrust_rfzero(int nlhs,
                            mxArray * plhs[],
                            int nrhs,
                            mxArray * prhs[]) {
    mxArray * mprhs[8];
    mxArray * mplhs[3];
    int i;
    if (nlhs > 3) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: trust/rfzero Line: 138 Colum"
            "n: 0 The function \"trust/rfzero\" was called with"
            " more than the declared number of outputs (3)"));
    }
    if (nrhs > 8) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: trust/rfzero Line: 138 Colu"
            "mn: 0 The function \"trust/rfzero\" was called wi"
            "th more than the declared number of inputs (8)"));
    }
    for (i = 0; i < 3; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 8 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 8; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(
      0,
      8,
      mprhs[0],
      mprhs[1],
      mprhs[2],
      mprhs[3],
      mprhs[4],
      mprhs[5],
      mprhs[6],
      mprhs[7]);
    mplhs[0]
      = Mtrust_rfzero(
          &mplhs[1],
          &mplhs[2],
          nlhs,
          mprhs[0],
          mprhs[1],
          mprhs[2],
          mprhs[3],
          mprhs[4],
          mprhs[5],
          mprhs[6],
          mprhs[7]);
    mlfRestorePreviousContext(
      0,
      8,
      mprhs[0],
      mprhs[1],
      mprhs[2],
      mprhs[3],
      mprhs[4],
      mprhs[5],
      mprhs[6],
      mprhs[7]);
    plhs[0] = mplhs[0];
    for (i = 1; i < 3 && i < nlhs; ++i) {
        plhs[i] = mplhs[i];
    }
    for (; i < 3; ++i) {
        mxDestroyArray(mplhs[i]);
    }
}