optim_private_xtrack.c

ASUFIT－Matlab－全局拟合程序

/*
 * MATLAB Compiler: 2.0.1
 * Date: Tue May 08 21:28:23 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 "optim_private_xtrack.h"
#include "subplot.h"
#include "title.h"
#include "xlabel.h"
#include "ylabel.h"

static double __Array0_r[1] = { -1.0 };

static double __Array1_r[1] = { 0.0 };

static double __Array2_r[2] = { -1.0, 1.0 };

static double __Array3_r[4] = { .5, .23, .19, .08 };

static double __Array4_r[4] = { .5, .315, .19, .05 };

static double __Array5_r[1] = { 0.0 };

static double __Array6_r[1] = { -1.0 };

static double __Array7_r[1] = { 0.0 };

static double __Array8_r[1] = { -1.0 };

static double __Array9_r[3] = { -1.0, 0.0, 1.0 };

static double __Array10_r[4] = { .5, .05, .19, .08 };

static double __Array11_r[4] = { .5, .135, .19, .05 };

/*
 * The function "Moptim_private_xtrack" is the implementation version of the
 * "optim/private/xtrack" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\xtrack.m" (lines 1-135). 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 xtrack(x, g, l, u, init)
 */
static void Moptim_private_xtrack(mxArray * x,
                                  mxArray * g,
                                  mxArray * l,
                                  mxArray * u,
                                  mxArray * init) {
    mxArray * activeg = mclGetUninitializedArray();
    mxArray * activel = mclGetUninitializedArray();
    mxArray * activeu = mclGetUninitializedArray();
    mxArray * activex = mclGetUninitializedArray();
    mxArray * ans = mclInitializeAns();
    mxArray * arg1 = mclGetUninitializedArray();
    mxArray * arg2 = mclGetUninitializedArray();
    mxArray * arg3 = mclGetUninitializedArray();
    mxArray * arg4 = mclGetUninitializedArray();
    mxArray * argl = mclGetUninitializedArray();
    mxArray * argu = mclGetUninitializedArray();
    mxArray * currsubplot = mclGetUninitializedArray();
    mxArray * dist = mclGetUninitializedArray();
    mxArray * glog = mclGetUninitializedArray();
    mxArray * index = mclGetUninitializedArray();
    mxArray * leg = mclGetUninitializedArray();
    mxArray * lin = mclGetUninitializedArray();
    mxArray * marksize = mclGetUninitializedArray();
    mxArray * maxg = mclGetUninitializedArray();
    mxArray * midpoint = mclGetUninitializedArray();
    mxArray * n = mclGetUninitializedArray();
    mxArray * nargin_ = mclGetUninitializedArray();
    mxArray * newg = mclGetUninitializedArray();
    mxArray * newx = mclGetUninitializedArray();
    mxArray * objh = mclGetUninitializedArray();
    mxArray * onen = mclGetUninitializedArray();
    mxArray * pos = mclGetUninitializedArray();
    mxArray * units = mclGetUninitializedArray();
    mxArray * w = mclGetUninitializedArray();
    mxArray * xlog = mclGetUninitializedArray();
    mxArray * zeron = mclGetUninitializedArray();
    mlfAssign(&nargin_, mlfNargin(0, x, g, l, u, init, NULL));
    mclValidateInputs("optim/private/xtrack", 5, &x, &g, &l, &u, &init);
    mclCopyArray(&init);
    /*
     * %XTRACK display graph
     * %
     * %   xtrack(x, g, l, u)
     * %
     * %   Output display showing the components of g
     * %   and the components of x, relative to the bounds l,u
     * %
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.5 $  $Date: 1998/07/09 16:30:09 $
     * 
     * set(gcf,'doublebuffer','on');
     */
    mclAssignAns(
      &ans,
      mlfNSet(
        0,
        mlfGcf(),
        mxCreateString("doublebuffer"),
        mxCreateString("on"),
        NULL));
    /*
     * 
     * if nargin < 5
     */
    if (mlfTobool(mlfLt(nargin_, mlfScalar(5.0)))) {
        /*
         * init = '';
         */
        mlfAssign(&init, mxCreateString(""));
    /*
     * end
     */
    }
    /*
     * if any(isnan(l)) | any(isnan(u))
     */
    {
        mxArray * a_ = mclInitialize(mlfAny(mlfIsnan(l), NULL));
        if (mlfTobool(a_) || mlfTobool(mlfOr(a_, mlfAny(mlfIsnan(u), NULL)))) {
            mxDestroyArray(a_);
            /*
             * error('NaN in lower or upper bounds')
             */
            mlfError(mxCreateString("NaN in lower or upper bounds"));
        } else {
            mxDestroyArray(a_);
        }
    /*
     * end
     */
    }
    /*
     * 
     * n = length(x);  
     */
    mlfAssign(&n, mlfLength(x));
    /*
     * index = (1:n)'; 
     */
    mlfAssign(&index, mlfCtranspose(mlfColon(mlfScalar(1.0), n, NULL)));
    /*
     * onen = ones(n,1);
     */
    mlfAssign(&onen, mlfOnes(n, mlfScalar(1.0), NULL));
    /*
     * 
     * maxg = max(abs(g)); 
     */
    mlfAssign(&maxg, mlfMax(NULL, mlfAbs(g), NULL, NULL));
    /*
     * maxg = max(maxg,1); 
     */
    mlfAssign(&maxg, mlfMax(NULL, maxg, mlfScalar(1.0), NULL));
    /*
     * zeron = zeros(n,1);
     */
    mlfAssign(&zeron, mlfZeros(n, mlfScalar(1.0), NULL));
    /*
     * midpoint = zeron;
     */
    mlfAssign(&midpoint, zeron);
    /*
     * 
     * arg1 = (u < inf) & (l > -inf);
     */
    mlfAssign(&arg1, mlfAnd(mlfLt(u, mlfInf()), mlfGt(l, mlfUminus(mlfInf()))));
    /*
     * arg2 = (u < inf) & ( l == -inf);
     */
    mlfAssign(&arg2, mlfAnd(mlfLt(u, mlfInf()), mlfEq(l, mlfUminus(mlfInf()))));
    /*
     * arg3 = (u == inf) & ( l > -inf);
     */
    mlfAssign(&arg3, mlfAnd(mlfEq(u, mlfInf()), mlfGt(l, mlfUminus(mlfInf()))));
    /*
     * arg4 = (u == inf) & (l == -inf);
     */
    mlfAssign(&arg4, mlfAnd(mlfEq(u, mlfInf()), mlfEq(l, mlfUminus(mlfInf()))));
    /*
     * 
     * newx = x;
     */
    mlfAssign(&newx, x);
    /*
     * 
     * % Shift and scale
     * dist = zeron;
     */
    mlfAssign(&dist, zeron);
    /*
     * dist(arg1) = min((x(arg1)-l(arg1)) ./ max(abs(l(arg1)),onen(arg1)), ...
     */
    mlfIndexAssign(
      &dist,
      "(?)",
      arg1,
      mlfMin(
        NULL,
        mlfRdivide(
          mlfMinus(mlfIndexRef(x, "(?)", arg1), mlfIndexRef(l, "(?)", arg1)),
          mlfMax(
            NULL,
            mlfAbs(mlfIndexRef(l, "(?)", arg1)),
            mlfIndexRef(onen, "(?)", arg1),
            NULL)),
        mlfRdivide(
          mlfMinus(mlfIndexRef(u, "(?)", arg1), mlfIndexRef(x, "(?)", arg1)),
          mlfMax(
            NULL,
            mlfAbs(mlfIndexRef(u, "(?)", arg1)),
            mlfIndexRef(onen, "(?)", arg1),
            NULL)),
        NULL));
    /*
     * (u(arg1)-x(arg1)) ./ max(abs(u(arg1)),onen(arg1)));
     * dist(arg2) = (u(arg2)-x(arg2)) ./ max(abs(u(arg2)), onen(arg2));
     */
    mlfIndexAssign(
      &dist,
      "(?)",
      arg2,
      mlfRdivide(
        mlfMinus(mlfIndexRef(u, "(?)", arg2), mlfIndexRef(x, "(?)", arg2)),
        mlfMax(
          NULL,
          mlfAbs(mlfIndexRef(u, "(?)", arg2)),
          mlfIndexRef(onen, "(?)", arg2),
          NULL)));
    /*
     * dist(arg3) = (x(arg3)-l(arg3)) ./ max(abs(l(arg3)), onen(arg3)); 
     */
    mlfIndexAssign(
      &dist,
      "(?)",
      arg3,
      mlfRdivide(
        mlfMinus(mlfIndexRef(x, "(?)", arg3), mlfIndexRef(l, "(?)", arg3)),
        mlfMax(
          NULL,
          mlfAbs(mlfIndexRef(l, "(?)", arg3)),
          mlfIndexRef(onen, "(?)", arg3),
          NULL)));
    /*
     * argu = (u < inf) & (dist == (u-x) ./ max(abs(u),onen));
     */
    mlfAssign(
      &argu,
      mlfAnd(
        mlfLt(u, mlfInf()),
        mlfEq(
          dist,
          mlfRdivide(mlfMinus(u, x), mlfMax(NULL, mlfAbs(u), onen, NULL)))));
    /*
     * argl = (l > -inf) & (dist == (x-l) ./ max(abs(l),onen));
     */
    mlfAssign(
      &argl,
      mlfAnd(
        mlfGt(l, mlfUminus(mlfInf())),
        mlfEq(
          dist,
          mlfRdivide(mlfMinus(x, l), mlfMax(NULL, mlfAbs(l), onen, NULL)))));
    /*
     * dist = min(dist, onen-.001*onen);
     */
    mlfAssign(
      &dist,
      mlfMin(
        NULL, dist, mlfMinus(onen, mlfMtimes(mlfScalar(.001), onen)), NULL));
    /*
     * dist = max(dist, eps*onen);
     */
    mlfAssign(&dist, mlfMax(NULL, dist, mlfMtimes(mlfEps(), onen), NULL));
    /*
     * xlog = min(-onen ./ log(dist),onen);
     */
    mlfAssign(
      &xlog,
      mlfMin(NULL, mlfRdivide(mlfUminus(onen), mlfLog(dist)), onen, NULL));
    /*
     * newx(argl) = -1*(onen(argl)- xlog(argl));
     */
    mlfIndexAssign(
      &newx,
      "(?)",
      argl,
      mlfMtimes(
        mlfScalar(-1.0),
        mlfMinus(
          mlfIndexRef(onen, "(?)", argl), mlfIndexRef(xlog, "(?)", argl))));
    /*
     * newx(argu) = 1*(onen(argu)-xlog(argu));
     */
    mlfIndexAssign(
      &newx,
      "(?)",
      argu,
      mlfMtimes(
        mlfScalar(1.0),
        mlfMinus(
          mlfIndexRef(onen, "(?)", argu), mlfIndexRef(xlog, "(?)", argu))));
    /*
     * newx(arg4) = zeron(arg4);
     */
    mlfIndexAssign(&newx, "(?)", arg4, mlfIndexRef(zeron, "(?)", arg4));
    /*
     * 
     * % Compute active constraints
     * activel=(abs(x-l)< 1e-5*max(abs(l),onen));
     */
    mlfAssign(
      &activel,
      mlfLt(
        mlfAbs(mlfMinus(x, l)),
        mlfMtimes(mlfScalar(1e-5), mlfMax(NULL, mlfAbs(l), onen, NULL))));
    /*
     * activeu=(abs(u-x)< 1e-5*max(abs(u),onen));
     */
    mlfAssign(
      &activeu,
      mlfLt(
        mlfAbs(mlfMinus(u, x)),
        mlfMtimes(mlfScalar(1e-5), mlfMax(NULL, mlfAbs(u), onen, NULL))));
    /*
     * 
     * % Scale g
     * newg = g/(maxg + 1);
     */
    mlfAssign(&newg, mlfMrdivide(g, mlfPlus(maxg, mlfScalar(1.0))));
    /*
     * w = max(abs(newg),eps*onen);
     */
    mlfAssign(&w, mlfMax(NULL, mlfAbs(newg), mlfMtimes(mlfEps(), onen), NULL));
    /*
     * glog = -onen./log(w);
     */
    mlfAssign(&glog, mlfRdivide(mlfUminus(onen), mlfLog(w)));
    /*
     * glog = min(glog,onen);
     */
    mlfAssign(&glog, mlfMin(NULL, glog, onen, NULL));
    /*
     * newg = sign(newg).*glog;
     */
    mlfAssign(&newg, mlfTimes(mlfSign(newg), glog));
    /*
     * activeg = (abs(g) < 1e-6);
     */
    mlfAssign(&activeg, mlfLt(mlfAbs(g), mlfScalar(1e-6)));
    /*
     * 
     * switch init
     * case ''       % default case
     */
    if (mclSwitchCompare(init, mxCreateString(""))) {
        /*
         * % Update Plots
         * % Upper plot
         * subplot(2,1,1) ;
         */
        mclAssignAns(
          &ans, mlfNSubplot(0, mlfScalar(2.0), mlfScalar(1.0), mlfScalar(1.0)));
        /*
         * activex = activel | activeu ;
         */
        mlfAssign(&activex, mlfOr(activel, activeu));
        /*
         * set(findobj(gca,'tag','blueline'),'XData',index(~activex),'YData',newx(~activex));
         */
        mclAssignAns(
          &ans,
          mlfNSet(
            0,
            mlfFindobj(
              mlfGca(NULL),
              mxCreateString("tag"),
              mxCreateString("blueline"),
              NULL),
            mxCreateString("XData"),
            mlfIndexRef(index, "(?)", mlfNot(activex)),
            mxCreateString("YData"),
            mlfIndexRef(newx, "(?)", mlfNot(activex)),
            NULL));
        /*
         * set(findobj(gca,'tag','redline'),'XData',index(activex),'YData',newx(activex));
         */
        mclAssignAns(
          &ans,
          mlfNSet(
            0,
            mlfFindobj(
              mlfGca(NULL),
              mxCreateString("tag"),
              mxCreateString("redline"),
              NULL),
            mxCreateString("XData"),
            mlfIndexRef(index, "(?)", activex),
            mxCreateString("YData"),
            mlfIndexRef(newx, "(?)", activex),
            NULL));
        /*
         * 
         * % Lower Plot
         * subplot(2,1,2);
         */
        mclAssignAns(
          &ans, mlfNSubplot(0, mlfScalar(2.0), mlfScalar(1.0), mlfScalar(2.0)));
        /*
         * set(findobj(gca,'tag','blueline'),'XData',index(~activeg),'YData',newg(~activeg));
         */
        mclAssignAns(
          &ans,
          mlfNSet(
            0,
            mlfFindobj(
              mlfGca(NULL),
              mxCreateString("tag"),
              mxCreateString("blueline"),
              NULL),
            mxCreateString("XData"),
            mlfIndexRef(index, "(?)", mlfNot(activeg)),
            mxCreateString("YData"),
            mlfIndexRef(newg, "(?)", mlfNot(activeg)),
            NULL));
        /*
         * set(findobj(gca,'tag','redline'),'XData',index(activeg),'YData',newg(activeg));
         */
        mclAssignAns(
          &ans,
          mlfNSet(
            0,
            mlfFindobj(
              mlfGca(NULL),
              mxCreateString("tag"),
              mxCreateString("redline"),
              NULL),
            mxCreateString("XData"),
            mlfIndexRef(index, "(?)", activeg),
            mxCreateString("YData"),
            mlfIndexRef(newg, "(?)", activeg),
            NULL));
    /*
     * 
     * case 'init'
     */
    } else if (mclSwitchCompare(init, mxCreateString("init"))) {
        /*
         * % Calculate markersize
         * units = get(gca,'units') ; set(gca,'units','points') ;
         */
        mlfAssign(&units, mlfGet(mlfGca(NULL), mxCreateString("units"), NULL));
        mclAssignAns(
          &ans,
          mlfNSet(
            0,
            mlfGca(NULL),
            mxCreateString("units"),
            mxCreateString("points"),
            NULL));
        /*
         * pos = get(gca,'position'); 
         */
        mlfAssign(&pos, mlfGet(mlfGca(NULL), mxCreateString("position"), NULL));
        /*
         * marksize = max(1,min(15,round(3*pos(3)/n)));
         */
        mlfAssign(
          &marksize,
          mlfMax(
            NULL,
            mlfScalar(1.0),
            mlfMin(