spcfun.c

ASUFIT－Matlab－全局拟合程序

/*
 * MATLAB Compiler: 2.0.1
 * Date: Tue May 08 21:28:16 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 "spcfun.h"

extern mxArray * amps;
extern mxArray * chisq;
extern mxArray * data;
extern mxArray * fit;
extern mxArray * irf;
extern mxArray * numTao;
extern mxArray * numTaoFixed;
extern mxArray * spcLeftEnd;
extern mxArray * spcRightEnd;
extern mxArray * t;
extern mxArray * taoFixed;
extern mxArray * textBox;
extern mxArray * weight;

static mxArray * Mspcfun_irfshift(int nargout_,
                                  mxArray * f1,
                                  mxArray * f2,
                                  mxArray * x0);
static mxArray * mlfSpcfun_irfshift(mxArray * f1, mxArray * f2, mxArray * x0);
static void mlxSpcfun_irfshift(int nlhs,
                               mxArray * plhs[],
                               int nrhs,
                               mxArray * prhs[]);

/*
 * The function "Mspcfun" is the implementation version of the "spcfun"
 * M-function from file "C:\MATLABR11\asufit24\spcfun.m" (lines 1-70). 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 f = spcfun(params)
 */
static mxArray * Mspcfun(int nargout_, mxArray * params) {
    mxArray * f = mclGetUninitializedArray();
    mxArray * R = mclGetUninitializedArray();
    mxArray * a = mclGetUninitializedArray();
    mxArray * ans = mclInitializeAns();
    mxArray * bottom = mclGetUninitializedArray();
    mxArray * count = mclGetUninitializedArray();
    mxArray * curve = mclGetUninitializedArray();
    mxArray * dataTemp = mclGetUninitializedArray();
    mxArray * i = mclGetUninitializedArray();
    mxArray * irfTemp = mclGetUninitializedArray();
    mclForLoopIterator iterator_0;
    mclForLoopIterator iterator_1;
    mxArray * numCurves = mclGetUninitializedArray();
    mxArray * numPoints = mclGetUninitializedArray();
    mxArray * shift = mclGetUninitializedArray();
    mxArray * tao = mclGetUninitializedArray();
    mxArray * taoFit = mclGetUninitializedArray();
    mxArray * top = mclGetUninitializedArray();
    mxArray * tpch = mclGetUninitializedArray();
    mxArray * z = mclGetUninitializedArray();
    mclValidateInputs("spcfun", 1, &params);
    /*
     * 
     * global data
     * global irf
     * global weight
     * global fit
     * global amps
     * global numTao
     * global numTaoFixed
     * global taoFixed
     * global t
     * global chisq
     * global textBox
     * global spcLeftEnd
     * global spcRightEnd
     * global absLeftEnd2
     * global absRightEnd2
     * %This is the function referred to in the call to "leastsq".
     * %It assumes the form:
     * %y(x)=a11*exp(-k1*t)+a12*exp(-k2*t)+. . .+a1n*exp(-kn*t)
     * %     a21*exp(-k1*t)+a22*exp(-k2*t)+. . .+a2n*exp(-kn*t)	
     * %     .
     * %     .
     * %     .
     * %     am1*exp(-k1*t)+am2*exp(-k2*t)+. . .+amn*exp(-kn*t)
     * %and makes available the vector of k values and the matrix
     * %of a values to the script file "fit.m".
     * 
     * tpch = t(2)-t(1);
     */
    mlfAssign(
      &tpch,
      mlfMinus(
        mlfIndexRef(t, "(?)", mlfScalar(2.0)),
        mlfIndexRef(t, "(?)", mlfScalar(1.0))));
    /*
     * dataTemp = data(spcLeftEnd:spcRightEnd,:);
     */
    mlfAssign(
      &dataTemp,
      mlfIndexRef(
        data,
        "(?,?)",
        mlfColon(spcLeftEnd, spcRightEnd, NULL),
        mlfCreateColonIndex()));
    /*
     * irfTemp = irf(spcLeftEnd:spcRightEnd,:);
     */
    mlfAssign(
      &irfTemp,
      mlfIndexRef(
        irf,
        "(?,?)",
        mlfColon(spcLeftEnd, spcRightEnd, NULL),
        mlfCreateColonIndex()));
    /*
     * numPoints = size(dataTemp,1);
     */
    mlfAssign(
      &numPoints, mlfSize(mclValueVarargout(), dataTemp, mlfScalar(1.0)));
    /*
     * numCurves = size(dataTemp,2);
     */
    mlfAssign(
      &numCurves, mlfSize(mclValueVarargout(), dataTemp, mlfScalar(2.0)));
    /*
     * curve = zeros(numPoints,numTao + numTaoFixed);
     */
    mlfAssign(&curve, mlfZeros(numPoints, mlfPlus(numTao, numTaoFixed), NULL));
    /*
     * fit = zeros(size(dataTemp));
     */
    mlfAssign(
      &fit, mlfZeros(mlfSize(mclValueVarargout(), dataTemp, NULL), NULL));
    /*
     * a = zeros(numTao + numTaoFixed,numCurves);
     */
    mlfAssign(&a, mlfZeros(mlfPlus(numTao, numTaoFixed), numCurves, NULL));
    /*
     * taoFit = zeros(1,numTao + numTaoFixed);
     */
    mlfAssign(
      &taoFit, mlfZeros(mlfScalar(1.0), mlfPlus(numTao, numTaoFixed), NULL));
    /*
     * tao = params(1:numTao);
     */
    mlfAssign(
      &tao, mlfIndexRef(params, "(?)", mlfColon(mlfScalar(1.0), numTao, NULL)));
    /*
     * shift = params(numTao + 1:numTao+numCurves);
     */
    mlfAssign(
      &shift,
      mlfIndexRef(
        params,
        "(?)",
        mlfColon(
          mlfPlus(numTao, mlfScalar(1.0)), mlfPlus(numTao, numCurves), NULL)));
    /*
     * shift = shift./tpch;
     */
    mlfAssign(&shift, mlfRdivide(shift, tpch));
    /*
     * taoFit(1,1:numTaoFixed) = taoFixed;
     */
    mlfIndexAssign(
      &taoFit,
      "(?,?)",
      mlfScalar(1.0),
      mlfColon(mlfScalar(1.0), numTaoFixed, NULL),
      taoFixed);
    /*
     * taoFit(1,numTaoFixed+1:numTao + numTaoFixed) = tao;
     */
    mlfIndexAssign(
      &taoFit,
      "(?,?)",
      mlfScalar(1.0),
      mlfColon(
        mlfPlus(numTaoFixed, mlfScalar(1.0)),
        mlfPlus(numTao, numTaoFixed),
        NULL),
      tao);
    /*
     * taoFit = taoFit./tpch;
     */
    mlfAssign(&taoFit, mlfRdivide(taoFit, tpch));
    /*
     * %a=params(numTao+2:length(params))';
     * bottom = floor(shift);
     */
    mlfAssign(&bottom, mlfFloor(shift));
    /*
     * top = ceil(shift);
     */
    mlfAssign(&top, mlfCeil(shift));
    /*
     * R = exp(-1*taoFit.^(-1));
     */
    mlfAssign(
      &R,
      mlfExp(mlfMtimes(mlfScalar(-1.0), mlfPower(taoFit, mlfScalar(-1.0)))));
    /*
     * for i= 1:numCurves,
     */
    for (mclForStart(&iterator_0, mlfScalar(1.0), numCurves, NULL);
         mclForNext(&iterator_0, &i);
         ) {
        /*
         * if (top(i) < 1)
         */
        if (mlfTobool(mlfLt(mlfIndexRef(top, "(?)", i), mlfScalar(1.0)))) {
            /*
             * top(i) = 1;
             */
            mlfIndexAssign(&top, "(?)", i, mlfScalar(1.0));
        /*
         * end
         */
        }
        /*
         * curve(1,:) = irfshift(irfTemp(top(i),i),irfTemp(top(i)+1,i),shift(i));
         */
        mlfIndexAssign(
          &curve,
          "(?,?)",
          mlfScalar(1.0),
          mlfCreateColonIndex(),
          mlfSpcfun_irfshift(
            mlfIndexRef(irfTemp, "(?,?)", mlfIndexRef(top, "(?)", i), i),
            mlfIndexRef(
              irfTemp,
              "(?,?)",
              mlfPlus(mlfIndexRef(top, "(?)", i), mlfScalar(1.0)),
              i),
            mlfIndexRef(shift, "(?)", i)));
        /*
         * for count = 1:(numPoints-1),
         */
        for (mclForStart(
               &iterator_1,
               mlfScalar(1.0),
               mlfMinus(numPoints, mlfScalar(1.0)),
               NULL);
             mclForNext(&iterator_1, &count);
             ) {
            /*
             * top(i) = bottom(i) + count;
             */
            mlfIndexAssign(
              &top, "(?)", i, mlfPlus(mlfIndexRef(bottom, "(?)", i), count));
            /*
             * if (top(i) < 1) top(i) = 1; end
             */
            if (mlfTobool(mlfLt(mlfIndexRef(top, "(?)", i), mlfScalar(1.0)))) {
                mlfIndexAssign(&top, "(?)", i, mlfScalar(1.0));
            }
            /*
             * if (top(i) > numPoints-2) top(i) = numPoints-2; end
             */
            if (mlfTobool(
                  mlfGt(
                    mlfIndexRef(top, "(?)", i),
                    mlfMinus(numPoints, mlfScalar(2.0))))) {
                mlfIndexAssign(
                  &top, "(?)", i, mlfMinus(numPoints, mlfScalar(2.0)));
            }
            /*
             * curve(count+1,:) = R.*(0.5*irfshift(irfTemp(top(i),i),irfTemp(top(i)+1,i),shift(i))+...