cubici3.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 "cubici3.h"

/*
 * The function "Mcubici3" is the implementation version of the "cubici3"
 * M-function from file "C:\MATLABR11\toolbox\optim\cubici3.m" (lines 1-30). 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 [stepmin,fbest]=cubici3(fnew,fold,graddnew,graddold,stepsize)
 */
static mxArray * Mcubici3(mxArray * * fbest,
                          int nargout_,
                          mxArray * fnew,
                          mxArray * fold,
                          mxArray * graddnew,
                          mxArray * graddold,
                          mxArray * stepsize) {
    mxArray * stepmin = mclGetUninitializedArray();
    mxArray * abd = mclGetUninitializedArray();
    mxArray * amat = mclGetUninitializedArray();
    mxArray * bmat = mclGetUninitializedArray();
    mxArray * root = mclGetUninitializedArray();
    mxArray * x1 = mclGetUninitializedArray();
    mclValidateInputs(
      "cubici3", 5, &fnew, &fold, &graddnew, &graddold, &stepsize);
    mclCopyArray(&fnew);
    /*
     * %CUBICI3  Cubicly interpolates 2 points and gradients to find step and min.
     * %
     * %   This function uses cubic interpolation and the values of 
     * %   two points and their gradients in order estimate the minimum of a 
     * %   a function along a line.
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.10 $  $Date: 1997/11/29 01:23:03 $
     * %   Andy Grace 7-9-90.
     * if isinf(fnew), fnew=1/eps; end
     */
    if (mlfTobool(mlfIsinf(fnew))) {
        mlfAssign(&fnew, mlfMrdivide(mlfScalar(1.0), mlfEps()));
    }
    /*
     * amat=[1/3*stepsize^3 , 0.5*stepsize^2; stepsize^2     stepsize];
     */
    mlfAssign(
      &amat,
      mlfVertcat(
        mlfHorzcat(
          mlfMtimes(
            mlfMrdivide(mlfScalar(1.0), mlfScalar(3.0)),
            mlfMpower(stepsize, mlfScalar(3.0))),
          mlfMtimes(mlfScalar(0.5), mlfMpower(stepsize, mlfScalar(2.0))),
          NULL),
        mlfHorzcat(mlfMpower(stepsize, mlfScalar(2.0)), stepsize, NULL),
        NULL));
    /*
     * bmat=[fnew-graddold*stepsize-fold; graddnew-graddold];
     */
    mlfAssign(
      &bmat,
      mlfVertcat(
        mlfHorzcat(
          mlfMinus(mlfMinus(fnew, mlfMtimes(graddold, stepsize)), fold), NULL),
        mlfHorzcat(mlfMinus(graddnew, graddold), NULL),
        NULL));
    /*
     * abd=amat\bmat;
     */
    mlfAssign(&abd, mlfMldivide(amat, bmat));
    /*
     * if isequal(abd(1),0)
     */
    if (mlfTobool(
          mlfIsequal(
            mlfIndexRef(abd, "(?)", mlfScalar(1.0)), mlfScalar(0.0), NULL))) {
        /*
         * x1 = -graddold/abd(2);
         */
        mlfAssign(
          &x1,
          mlfMrdivide(
            mlfUminus(graddold), mlfIndexRef(abd, "(?)", mlfScalar(2.0))));
        /*
         * stepmin = x1;
         */
        mlfAssign(&stepmin, x1);
    /*
     * else
     */
    } else {
        /*
         * root=real(sqrt(abd(2)^2-4*abd(1)*graddold));
         */
        mlfAssign(
          &root,
          mlfReal(
            mlfSqrt(
              mlfMinus(
                mlfMpower(
                  mlfIndexRef(abd, "(?)", mlfScalar(2.0)), mlfScalar(2.0)),
                mlfMtimes(
                  mlfMtimes(
                    mlfScalar(4.0), mlfIndexRef(abd, "(?)", mlfScalar(1.0))),
                  graddold)))));
        /*
         * x1=(-abd(2)+root)/(2*abd(1));
         */
        mlfAssign(
          &x1,
          mlfMrdivide(
            mlfPlus(mlfUminus(mlfIndexRef(abd, "(?)", mlfScalar(2.0))), root),
            mlfMtimes(
              mlfScalar(2.0), mlfIndexRef(abd, "(?)", mlfScalar(1.0)))));
        /*
         * if 2*abd(1)*x1+abd(2)>0
         */
        if (mlfTobool(
              mlfGt(
                mlfPlus(
                  mlfMtimes(
                    mlfMtimes(
                      mlfScalar(2.0), mlfIndexRef(abd, "(?)", mlfScalar(1.0))),
                    x1),
                  mlfIndexRef(abd, "(?)", mlfScalar(2.0))),
                mlfScalar(0.0)))) {
            /*
             * stepmin=x1;
             */
            mlfAssign(&stepmin, x1);
        /*
         * else
         */
        } else {
            /*
             * stepmin=(-abd(2)-root)/(2*abd(1));;
             */
            mlfAssign(
              &stepmin,
              mlfMrdivide(
                mlfMinus(
                  mlfUminus(mlfIndexRef(abd, "(?)", mlfScalar(2.0))), root),
                mlfMtimes(
                  mlfScalar(2.0), mlfIndexRef(abd, "(?)", mlfScalar(1.0)))));
        /*
         * end
         */
        }
    /*
     * end
     */
    }
    /*
     * 
     * if stepmin<0,  stepmin=-stepmin; end
     */
    if (mlfTobool(mlfLt(stepmin, mlfScalar(0.0)))) {
        mlfAssign(&stepmin, mlfUminus(stepmin));
    }
    /*
     * fbest=1/3*abd(1)*stepmin^3+0.5*abd(2)*stepmin^2+graddold*stepmin+fold;
     */
    mlfAssign(
      fbest,
      mlfPlus(
        mlfPlus(
          mlfPlus(
            mlfMtimes(
              mlfMtimes(
                mlfMrdivide(mlfScalar(1.0), mlfScalar(3.0)),
                mlfIndexRef(abd, "(?)", mlfScalar(1.0))),
              mlfMpower(stepmin, mlfScalar(3.0))),
            mlfMtimes(
              mlfMtimes(
                mlfScalar(0.5), mlfIndexRef(abd, "(?)", mlfScalar(2.0))),
              mlfMpower(stepmin, mlfScalar(2.0)))),
          mlfMtimes(graddold, stepmin)),
        fold));
    mclValidateOutputs("cubici3", 2, nargout_, &stepmin, fbest);
    mxDestroyArray(abd);
    mxDestroyArray(amat);
    mxDestroyArray(bmat);
    mxDestroyArray(fnew);
    mxDestroyArray(root);
    mxDestroyArray(x1);
    return stepmin;
}

/*
 * The function "mlfCubici3" contains the normal interface for the "cubici3"
 * M-function from file "C:\MATLABR11\toolbox\optim\cubici3.m" (lines 1-30).
 * This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
mxArray * mlfCubici3(mxArray * * fbest,
                     mxArray * fnew,
                     mxArray * fold,
                     mxArray * graddnew,
                     mxArray * graddold,
                     mxArray * stepsize) {
    int nargout = 1;
    mxArray * stepmin = mclGetUninitializedArray();
    mxArray * fbest__ = mclGetUninitializedArray();
    mlfEnterNewContext(1, 5, fbest, fnew, fold, graddnew, graddold, stepsize);
    if (fbest != NULL) {
        ++nargout;
    }
    stepmin
      = Mcubici3(&fbest__, nargout, fnew, fold, graddnew, graddold, stepsize);
    mlfRestorePreviousContext(
      1, 5, fbest, fnew, fold, graddnew, graddold, stepsize);
    if (fbest != NULL) {
        mclCopyOutputArg(fbest, fbest__);
    } else {
        mxDestroyArray(fbest__);
    }
    return mlfReturnValue(stepmin);
}

/*
 * The function "mlxCubici3" contains the feval interface for the "cubici3"
 * M-function from file "C:\MATLABR11\toolbox\optim\cubici3.m" (lines 1-30).
 * The feval function calls the implementation version of cubici3 through this
 * function. This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
void mlxCubici3(int nlhs, mxArray * plhs[], int nrhs, mxArray * prhs[]) {
    mxArray * mprhs[5];
    mxArray * mplhs[2];
    int i;
    if (nlhs > 2) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: cubici3 Line: 1 Column: "
            "0 The function \"cubici3\" was called with mor"
            "e than the declared number of outputs (2)"));
    }
    if (nrhs > 5) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: cubici3 Line: 1 Column:"
            " 0 The function \"cubici3\" was called with m"
            "ore than the declared number of inputs (5)"));
    }
    for (i = 0; i < 2; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 5 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 5; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(0, 5, mprhs[0], mprhs[1], mprhs[2], mprhs[3], mprhs[4]);
    mplhs[0]
      = Mcubici3(
          &mplhs[1], nlhs, mprhs[0], mprhs[1], mprhs[2], mprhs[3], mprhs[4]);
    mlfRestorePreviousContext(
      0, 5, mprhs[0], mprhs[1], mprhs[2], mprhs[3], mprhs[4]);
    plhs[0] = mplhs[0];
    for (i = 1; i < 2 && i < nlhs; ++i) {
        plhs[i] = mplhs[i];
    }
    for (; i < 2; ++i) {
        mxDestroyArray(mplhs[i]);
    }
}