#include "IndependentAlg.h"
#include "fitfunc.h"
#include "genenalfunc.h"
#include "matlab_func.h"
#include <QFile>
#include <QCoreApplication>

const long double xrayArea_g = 0.21212121212121213L;

using namespace arma;
void PrintMat2fvec(QString name, field<vec>& vecc);
void PrintColvec(QString name, colvec& dataMat);


namespace Independ {

mat calFcnEval(mat x, colvec para)
{
    mat y = zeros(size(x));
    if(x.is_empty() || para.size() < 2) return y;

    int funType = para(0u);
    vec p = para;
    p.shed_row(0u);

    if(calParaCheck(para))
    {
        switch(funType)
        {
        case 1: y = Fcn::l1interpolate(x, p); break;
        case 2: y = Fcn::l2polynomial(x, p); break;
        case 3: y = Fcn::l3sqrtpoly(x, p); break;
        case 4: y = Fcn::l4sqrtofpoly(x, p); break;
        case 5: y = Fcn::l5hteff(x, p); break;
        case 6: y = Fcn::l6loglogpoly(x, p); break;
        case 7: y = Fcn::l7logpoly(x, p); break;
        case 8: y = Fcn::l8loginvlog(x, p); break;
        case 9: y = Fcn::l9invexp(x, p); break;
        case 93: y = Fcn::l93ht2fix(x, p); break;
        case 94: y = Fcn::l94ht2free(x, p); break;
        case 95: y = Fcn::l95ht3eff(x, p); break;
        case 96: y = Fcn::l96invpow(x, p); break;
        case 97: y = Fcn::l97expsum(x, p); break;
        case 98: y = Fcn::l98constant(x, p); break;
        case 99: y = Fcn::l99cutlin(x, p); break;
        default: qDebug() << "invalid function type indicator " << funType; break;
        }
    }
    // y = reshape(y, size(x));
    return y;
}

mat calDerivEval(mat x, colvec p)
{
    mat dy;
    if(x.is_empty() || p.size() < 2) return dy;

    int funType = p(0u);
    colvec para = p.rows(1u, p.size()-1u);

    if(Independ::calParaCheck(p))
    {
        switch(funType)
        {
        case 1: dy = Deriv::l1interpolate(x, para); break;
        case 2: dy = Deriv::l2polynomial(x, para); break;
        case 3: dy = Deriv::l3sqrtpoly(x, para); break;
        case 96: dy = Deriv::l96invpow(x, para); break;
        case 97: dy = Deriv::l97expsum(x, para); break;
        case 98: dy = zeros(size(x));            break;
        case 99: dy = Deriv::l99cutlin(x, para); break;
        default: dy = Deriv::lApproxDeriv(x, p); break;
        }
    }
    return dy;
}

bool calParaCheck(vec p)
{
    bool bRet = true;
    if(p.is_empty()) return false;

    vec para = p;
    para.shed_row(0u);
    int s = para.size();
    int ft = p(0);
    switch(ft)
    {
    case 1:
        if(s % 2 != 0)
        {
            bRet = false;
            qDebug() << "parameters must be data pairs";
        }
        else if(s < 4)
        {
            bRet = false;
            qDebug() << "requires at least 2 data points";
        }
        else {
            if(any(diff(vectorise(MySpan(para, 0, s-1, 2))) <= 0))
            {
                bRet = false;
                qDebug() << "x data must be monotonic";
            }
        }
        break;

    case 2:
    case 3:
        if(s < 2)// if isempty(para)
        {
           bRet = false;
           qDebug() << "polynomial of too small degree";
        }
        else if(para(1) <= 0)
        {
           bRet = false;
           qDebug() << "first order coefficient must be positive";
        }
        break;

    case 4:
        if(s < 2)
        {
           bRet = false;
           qDebug() << "polynomial of too small degree";
        }
        else if(para(0) <= 0)
        {
           bRet = false;
           qDebug() << "constant coefficient must be positive";
        }
        else if(para(1) <= 0)
        {
           bRet = false;
           qDebug() << "first order coefficient must be positive";
        }
        break;

    case 5:
        if (s != 5)
        {
            bRet = false;
            qDebug() << "parameter vector has to have length 5";
        }
        else if (any(para <= 0))
        {
            bRet = false;
            qDebug() << "parameters have to be positive";
        }
        break;

    case 6:
        if (s < 2)
        {
            bRet = false;
            qDebug() << "at least polynomial of degree 1 required";
        }
        break;

    case 7:
    case 8:
        if (s < 3)
        {
            bRet = false;
            qDebug() << "at least polynomial of degree 2 required";
        }
        break;

    case 9:
        if (s != 4)
        {
            bRet = false;
            qDebug() << "four parameters required";
        }
        else if (any(para <= 0))
        {
            bRet = false;
            qDebug() << "parameters must be positive";
        }
        break;

    case 93:
        if (s != 5)
        {
            bRet = false;
            qDebug() << "5 parameter coefficients required";
        }
        break;

    case 94:
        if(s != 6)
        {
           bRet = false;
           qDebug() << "6 parameter coefficients required";
        }
        break;

    case 95:
        if (s != 8)
        {
            bRet = false;
            qDebug() << "8 parameter coefficients required";
        }
        break;

    case 96:
        if (s != 3)
        {
            bRet = false;
            qDebug() << "3 parameter coefficients required";
        }
        break;

    case 97:
        if (s % 2 != 1)
        {
            bRet = false;
            qDebug() << "odd number of parameters required";
        }
        break;

    case 98:
        if (s != 1)
        {
            bRet = false;
            qDebug() << "one constant reqired";
        }
        else if (!para.is_finite())
        {
            bRet = false;
            qDebug() << "constant not finite";
        }
        break;

    case 99:
        if (s != 3)
        {
            bRet = false;
            qDebug() << "three parameters required";
        }
        else if (!para.is_finite())
        {
            bRet = false;
            qDebug() << "parameters must be finite";
        }
        else if (para(2) < 0)
        {
            bRet = false;
            qDebug() << "slope must be non-negative";
        }
        break;

    default: bRet = false; break;
    }
    return bRet;
}

mat pchip(mat x, mat y, mat xx)
{
    // Check that data are acceptable and, if not, try to adjust them appropriately
    rowvec sizey;
    mat endslope;
    Matlab::chckxy(x, y, sizey, endslope, 3);

    uword n = x.size();
    rowvec t_x = x;
    rowvec h = diff(t_x);
    uword m = prod(sizey);

    // Compute slopes
    rowvec t_y = y;
    rowvec del = diff(t_y,1,1) / repmat(h,m,1);
	if (del.size() <= 0)
		return mat();

    mat slopes = zeros(size(y));
    for (uword r = 0; r<m; ++r)
    {
        //slopes(r,:) = pchipslopes(x,y(r,:),del(r,:));
        slopes.row(r) = pchipslopes(t_x, y.row(r), del.row(r));
    }

    // Compute piecewise cubic Hermite interpolant to those values and slopes
    PiecePoly v_pp = Matlab::pwch(x,y,slopes,h,del);
    v_pp.d = sizey;

    mat v;
    if (!xx.is_empty())   // if values are wanted instead, provide them
    {
        v = Matlab::ppval(xx, v_pp);
    }
    return v;
}

rowvec pchipslopes(rowvec x, rowvec y, rowvec del)
{
    //  Special case n=2, use linear interpolation.
    uword n = x.size();
    rowvec d;
    if (n == 2)
    {
        d = del(0) * ones<rowvec>(y.size());
        return d;
    }

    //  Slopes at interior points.
    //  d(k) = weighted average of del(k-1) and del(k) when they have the same sign.
    //  d(k) = 0 when del(k-1) and del(k) have opposites signs or either is zero.

    d = zeros<rowvec>(y.size());
    //k = find(sign(del(1:n-2)).*sign(del(2:n-1)) > 0);
    uvec k = find( sign(del.cols(0, n-3)) % sign(del.cols(1, n-2)) > 0 );

    rowvec h = diff(x);
    mat hs = h(k) + h(k+1u);
    mat w1 = (h(k)+hs) / (3*hs);
    mat w2 = (hs+h(k+1u)) / (3*hs);
    mat dmax = arma::max(abs(del(k)), abs(del(k+1u)));
    mat dmin = arma::min(abs(del(k)), abs(del(k+1u)));
    d(k+1u) = dmin / conj( w1 % (del(k) / dmax) + w2 % (del(k+1) / dmax) );

    //  Slopes at end points.
    //  Set d(1) and d(n) via non-centered, shape-preserving three-point formulae.
    //d(1) = ((2*h(1)+h(2))*del(1) - h(1)*del(2))/(h(1)+h(2));
    d(0) = ((2*h(0)+h(1))*del(0) - h(0)*del(1)) / (h(0)+h(1));
    if(sign(d(0)) != sign(del(0)))
    {
        d(0) = 0; //d(1) = 0;
    }
    else if (sign(del(0)) != sign(del(1)) && abs(d(0)) > abs(3*del(0)))
    {
        d(0) = 3*del(0); //d(1) = 3*del(1);
    }
    //d(n) = ((2*h(n-1)+h(n-2))*del(n-1) - h(n-1)*del(n-2))/(h(n-1)+h(n-2));
    d(n-1) = ((2*h(n-2)+h(n-3))*del(n-2) - h(n-2)*del(n-3))/(h(n-2)+h(n-3));
    if (sign(d(n-1)) != sign(del(n-2)))
    {
        d(n-1) = 0; //d(n) = 0;
    }
    else if (sign(del(n-2)) != sign(del(n-3)) && abs(d(n-1)) > abs(3*del(n-2)))
    {
        d(n-1) = 3*del(n-2); //d(n) = 3*del(n-1);
    }
    return d;
}

mat wrapBaseSpline(mat x, mat cx, mat cy, mat cdy)
{
    mat ep;
    PiecePoly pp = makeBasePP(vectorise(cx), vectorise(cy), cdy, ep, ep, ep);
    mat y = Matlab::ppval(x, pp);
    return y;
}

PiecePoly makeBasePP(vec cx, vec cy, mat cdy, mat pc, mat pfwhmc, mat pstep)
{
    PiecePoly bpp;
    if(cx.is_empty())
    {
        rowvec breaks("0,1");
        mat coefs("0");
        bpp = Matlab::mkpp(breaks, coefs);
        return bpp;
    }

    // substract steps from control point counts
    cy = cy - peakSteps(cx, pc, pfwhmc, pstep);

    // make sure first and last derivative not free
    if(cdy.size() < 1 || qIsInf(cdy(0)) || qIsInf(cdy(cdy.size()-1))) //if any(~isfinite(cdy([1,end])))
    {
        qDebug() << "Derivative of first and last control point must be given";
    }

    // get breaks in spline
    uvec idx = find_finite(cdy);

    // first spline piece
    //[tmpPP] = spline(cx(idx(1):idx(2)), [cdy(idx(1)); cy(idx(1):idx(2)); cdy(idx(2))]);
    vec t_cy;
    t_cy << cdy(idx(0)) << cdy(idx(1));
    t_cy.insert_rows(1, cy.rows(idx(0),idx(1)));
    PiecePoly tmpPP = Matlab::spline(cx.rows(idx(0),idx(1)), t_cy);
    //[breaks, coefs] = unmkpp(tmpPP);

    rowvec breaks = tmpPP.b;
    mat coefs = tmpPP.c;

    // add further pieces
    for(int i=1; i<idx.size()-1; ++i) //for i = 2 : length(idx)-1
    {
        //tmpPP = spline(cx(idx(i):idx(i+1)), [cdy(idx(i)); cy(idx(i):idx(i+1)); cdy(idx(i+1))]);
        t_cy << cdy(idx(i)) << cdy(idx(i+1));
        t_cy.insert_rows(1, cy.rows(idx(i),idx(i+1)));
        tmpPP = Matlab::spline(cx.rows(idx(i),idx(i+1)), t_cy);
        //[b, c] = unmkpp(tmpPP);
        // b(1) is same as breaks(end), i.e. cx(idx(i))
        // breaks row vector
        breaks = join_horiz( breaks, tmpPP.b.cols(1,tmpPP.b.size()-1) ); // breaks = [breaks, b(2:end)];
        // coefficients are Kx4 matrices
        coefs = join_vert(coefs, tmpPP.c); // coefs = [coefs; c];
    }

    // turn back into piecewise polynomial for use with ppval
    bpp = Matlab::mkpp(breaks, coefs);
	return bpp;
}

mat peakSteps(mat x, mat pC, mat pFWHM, mat pSt)
{
    mat s = zeros(size(x));

    mat pSigma = pFWHM / sqrt(8*log(2));

    for(int i=0; i<pC.size(); ++i) // i = 1 : length(pC)
    {
        // convolution of sharp edge with gaussian response
        s = s + pSt(i) * gaussInt(pC(i)-x, pSigma(i));
    }

    return s;
}

mat gaussInt(mat x, double s)
{
    mat y = 0.5*(1+Matlab::m_erf(x/(sqrt(2)*s)));
    return y;
}

// =====================================================start========================================================= //
mat peakBaseSpline(mat x, PiecePoly bpp, vec na, vec ct, vec fc, vec st, vec lt,
                                  vec lta, vec ut, vec uta, vec bw, vec rb, vec rd)
{
    /*if nargin < 12
        // to be removed when all callers are modified
        stk = dbstack;
        logWrite('Warning in %s, line %d: Syntax changed!', ...
        stk(end).name, stk(end).line);
        rb = repmat(NaN, size(ct));
        rd = rb;
    end*/

    // initialize count vector
    mat y = zeros(size(x));
    if(x.is_empty()) //if isempty(x)
    {
       return y;
    }

    // spline not defined, give zero spline contribution
    if(bpp.b.is_empty() || bpp.c.is_empty()) //if isempty(bpp)
    {
       bpp = Matlab::mkpp(rowvec("0, 1"), mat("0"));
    }
    /*elseif ~isstruct(bpp)
        error('bpp needs to be piecewise polynomial structure');
    end*/

    // check peak parameters match
    int NPeaks = na.size(); //NPeaks = length(na);
    /*na.reshape(NPeaks, 1);  //[na, s1] = doubleCheck(na, [NPeaks, 1]);
    ct.reshape(NPeaks, 1);  //[ct, s2] = doubleCheck(ct, [NPeaks, 1]);
    fc.reshape(NPeaks, 1);  //[fc, s3] = doubleCheck(fc, [NPeaks, 1]);
    st.reshape(NPeaks, 1);  //[st, s4] = doubleCheck(st, [NPeaks, 1]);
    lt.reshape(NPeaks, 1);  //[lt, s5] = doubleCheck(lt, [NPeaks, 1]);
    lta.reshape(NPeaks, 1); //[lta, s6] = doubleCheck(lta, [NPeaks, 1]);
    ut.reshape(NPeaks, 1);  //[ut, s7] = doubleCheck(ut, [NPeaks, 1]);
    uta.reshape(NPeaks, 1); //[uta, s8] = doubleCheck(uta, [NPeaks, 1]);
    bw.reshape(NPeaks, 1);  //[bw, s9] = doubleCheck(bw, [NPeaks, 1]);
    rb.reshape(NPeaks, 1);  //[rb, s10] = doubleCheck(rb, [NPeaks, 1]);
    rd.reshape(NPeaks, 1);  //[rd, s11] = doubleCheck(rd, [NPeaks, 1]);
    if ~all([s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11])
       error('peak data must match in length')
    end*/

    // baseline plus steps
    y = Matlab::ppval(x, bpp) + Independ::peakSteps(x, ct, fc, st);

    // fwhm to sigma
    vec sgma = fc / sqrt(8*log(2));

    // calculate peak influence width (in 1fwhm units)
    vec NAWidth = max(5, (1.5+0.5*log10(abs(na)+1)) ); // NAWidth = max(5, (1.5 + 0.5 * log10(abs(na)+1)));
    vec lowWidth = NAWidth + 3*lt/lta; // lowWidth = NAWidth + 3*lt./lta;
    vec upWidth = NAWidth + 3*ut/uta; // upWidth = NAWidth + 3*ut./uta;
    uvec sidx = find_finite(rb+rd); // sidx = find(~isnan(rb) & ~isnan(rd));

    if( any(sidx) )
    {
        //lowWidth(sidx) = lowWidth(sidx) + rd(sidx)./fc(sidx);
        lowWidth(sidx) = lowWidth(sidx) + rd(sidx) / fc(sidx);
        //upWidth(sidx) = upWidth(sidx) + 3./(fc(sidx).*rb(sidx));
        upWidth(sidx) = upWidth(sidx) + 3 / (fc(sidx) % rb(sidx));
    }

    // add steps and peaks
    for(int i=0; i<NPeaks; ++i)//for i = 1 : NPeaks
    {
        // find channels where peak has influence
        // ii = find(x >= ct(i) - lowWidth(i)*fc(i) & x <= ct(i) + upWidth(i)*fc(i) );
        uvec ii = find(x >= ct(i)-lowWidth(i)*fc(i) && x <= ct(i) + upWidth(i)*fc(i) );

        // add the peak, if it has influence to any channels
        if(!ii.is_empty()) //if ~isempty(ii)
        {
            mat tmp = peakShape(x(ii), ct(i), sgma(i), lt(i), lta(i), ut(i), uta(i), bw(i), rb(i), rd(i));
            y(ii) = y(ii) + na(i) * tmp;
        }
    }
    return y;
}

mat peakShape(mat x, double ct, double s, double lt, double lta,
                             double ut, double uta, double g, double rb, double rd)
{
    /*if nargin < 9
        logWrite('Warning: syntax changed');
        rb = NaN;
        rd = NaN;
    end*/
    mat y;
    if(qIsNaN(rb) || qIsInf(rb) || qIsNaN(rd) || qIsInf(rd)) //if all(isfinite([rb, rd]))
    {
        y = gaussianModel(x, ct, s, lt, lta, ut, uta, g);
    }
    else {
        y = neutronModel(x, ct, s, rb, rd);
    }
    return y;
}

mat neutronModel(mat x, double C, double sigma, double beta, double dC)
{
    mat cexp = convExponential(x-C, sigma, beta);
    mat clin = convLin(x-C, sigma, dC);
    mat ct = (cexp + clin) / (1/beta + dC/2);
    //ct.print("[][]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
    return ct;
}

mat gaussianModel(mat x, double ct, double s, double lt, double lta, double ut, double uta, double g)
{
    // default return value
    mat y = zeros(size(x));

    // get Junction points
    double Jl = lTailToJunct(s, lt);
    double Ju = lTailToJunct(s, ut);

    // sigma squared and normation factor
    double ss = pow(s, 2);
    double s2ps = sqrt(2*pi_)*s;

    // lower tail
    double AddLo;
    if(!qIsInf(Jl)) //if isfinite(Jl)
    {
        lta = lCheckAlpha(lta);
        double a = (2*Jl-lta*pow(Jl,2)) /(2*lta*ss);
        double k = Jl/(lta*ss);
        uvec i1 = find(x < ct - Jl);
        if(any(i1))
        {
           y(i1) = exp( a - k* pow(ct-Jl-x(i1)+1, lta) );
        }
        rowvec tmp; tmp << -Jl;
        AddLo = exp(a-k)/k - s2ps * as_scalar(Independ::gaussInt(tmp, s));
    }
    else {
        AddLo = 0;
    }

    // gaussian part
    uvec i2 = find(x >= ct - Jl && x <= ct + Ju);
    if(any(i2))
    {
        y(i2) = s2ps * NORMPDF(x(i2), ct, s);
    }

    // upper tail
    double AddHi;
    if(!qIsInf(Ju)) //if isfinite(Ju)
    {
        uta = lCheckAlpha(uta);
        double a = ( 2*Ju-uta*pow(Ju,2) ) /(2*uta*ss);
        double k = Ju/(uta*ss);
        uvec i3 = find(x > ct + Ju);
        if(any(i3))
        {
            y(i3) = exp( a - k * pow(x(i3)-ct-Ju+1, uta) );
        }
        rowvec tmp; tmp << -Ju;
        AddHi = exp(a-k)/k - s2ps * as_scalar( Independ::gaussInt(tmp, s) );
    }
    else {
        AddHi = 0;
    }

    // add difference between voigt and gaussian
    double AVoigt;
    if(g > 0)
    {
       AVoigt = s2ps * xrayAreaCorrection(s, g);
       double width;
       y = y + s2ps * voigt_gauss(width, x, ct, s, g);
    }
    else {
       AVoigt = 0;
    }

    // normalization to area 1
    y = y / (s2ps + AddLo + AddHi + AVoigt);

    return y;
}

mat convExponential(mat x, double sigma, double beta)
{
    //f1 = exp(beta.^2 .* sigma.^2 /2);
    double f1 = exp( pow(beta,2) * pow(sigma,2) / 2 );
    mat f2 = exp(-beta*x); //f2 = exp(-beta*x);
    //t1 = -(x - beta * sigma^2) / (sqrt(2)*sigma);
    mat t1 = -(x - beta * pow(sigma,2) ) / ( sqrt(2) * sigma );
    //ct = (f1.*f2.*(1-erf(t1))) / (2);
    mat ct = f1 * f2 % (1-Matlab::m_erf(t1)) / 2;
    return ct;
}

mat convLin(mat x, double sigma, double dC)
{
    mat t1 = (x + dC)/(sqrt(2)*sigma); //t1 = (x + dC)/(sqrt(2)*sigma);
    mat t2 = x/(sqrt(2)*sigma); //t2 = x/(sqrt(2)*sigma);
    //erfit = (x+dC).*(erf(t1)-erf(t2));
    mat erfit = (x+dC)%(Matlab::m_erf(t1)-Matlab::m_erf(t2));
    //t1t2 = sqrt(2/pi)*sigma*(exp(-t1.^2)-exp(-t2.^2));
    mat t1t2 = sqrt(2/pi_)*sigma*(exp(-pow(t1,2))-exp(-pow(t2,2)));
    //ct = 1/(2*dC)*(erfit + t1t2);
    mat ct = 1/(2*dC)*(erfit + t1t2);
    //double area = dC/2;
    return ct;
}

double lTailToJunct(double s, double t)
{
    double minTail = 0.01;
    double J;
    if(t < minTail)
    {
       J = qInf();
    }
    else {
       J = s/t;
    }
    return J;
}

double lCheckAlpha(double alpha)
{
    double lowTresh = 0.05;
    double alpha_use;
    if(alpha > 1) // alpha too big
    {
        alpha_use = 1;
        qDebug() << QString("alpha too big %1, setting to 1").arg(alpha);
    }
    else if(alpha >= lowTresh) // alpha in range
    {
        alpha_use = alpha;
    }
    else // alpha too small or NaN
    {
       alpha_use = lowTresh;
       qDebug() << QString("alpha too small %1, setting to %2").arg(alpha).arg(lowTresh);
    }
    return alpha_use;
}

mat NORMPDF(mat X, double MU, double SIGMA)
{
    if(X.is_empty()) // nargin < 1,
    {
        qDebug() << "Requires at least one input argument.";
    }

    //[errorcode x mu sigma] = distchck(3,x,mu,sigma);
    mat x, mu, sigma, out4;
    mu << MU; sigma << SIGMA;
    if( !distchck(x, mu, sigma, out4, 3, X, mu, sigma) )
    {
        qDebug() << "Requires non-scalar arguments to match in size.";
    }

    // Initialize Y to zero.
    mat y = zeros(size(x));
    uvec k = find(sigma > 0);
    if( any(k) )
    {
        //xn = (x(k) - mu(k)) ./ sigma(k);
        //y(k) = exp(-0.5 * xn .^2) ./ (sqrt(2*pi) .* sigma(k));
        mat xn = (x(k) - mu(k)) / sigma(k);
        y(k) = exp(-0.5 * pow(xn,2) ) / ( sqrt(2*pi_) * sigma(k) );
    }

    // Return NaN if SIGMA is negative or zero.
    uvec k1 = find(sigma <= 0);
    if( any(k1) )
    {
        //tmp = gNaN;
        //y(k1) = tmp(ones(size(k1)));
        MatAssign(y, k1, gNaN);
    }
    return y;
}

bool distchck(mat &out1, mat &out2, mat &out3, mat &out4, int nparms,
                             mat arg1, mat arg2, mat arg3, mat arg4)
{
    bool bRet = true;
    if( nparms == 1 )
    {
        out1 = arg1;
    }
    else if( nparms == 2 )
    {
        uword r1 = arg1.n_rows, c1 = arg1.n_cols; //[r1 c1] = size(arg1);
        uword r2 = arg2.n_rows, c2 = arg2.n_cols; //[r2 c2] = size(arg2);
        bool scalararg1 = (r1 * c1 == 1); //scalararg1 = (prod(size(arg1)) == 1);
        bool scalararg2 = (r2 * c2 == 1); //scalararg2 = (prod(size(arg2)) == 1);
        if(!scalararg1 && !scalararg2) //if ~scalararg1 & ~scalararg2
        {
            if(r1 != r2 || c1 != c2) return false;
        }
        if(scalararg1)
        {
            out1.set_size(r2,c2); //out1 = arg1(ones(r2,1),ones(c2,1));
            out1.fill( as_scalar(arg1) );
        }
        else { out1 = arg1; }

        if(scalararg2)
        {
            out2.set_size(r1, c1);
            out2.fill( as_scalar(arg2) ); // out2 = arg2(ones(r1,1),ones(c1,1));
        }
        else { out2 = arg2; }
    }
    else if( nparms == 3 )
    {
        uword r1 = arg1.n_rows, c1 = arg1.n_cols; //[r1 c1] = size(arg1);
        uword r2 = arg2.n_rows, c2 = arg2.n_cols; //[r2 c2] = size(arg2);
        uword r3 = arg3.n_rows, c3 = arg3.n_cols; //[r3 c3] = size(arg3);
        bool scalararg1 = (r1 * c1 == 1); //scalararg1 = (prod(size(arg1)) == 1);
        bool scalararg2 = (r2 * c2 == 1); //scalararg2 = (prod(size(arg2)) == 1);
        bool scalararg3 = (r3 * c3 == 1); //scalararg3 = (prod(size(arg3)) == 1);

        if(!scalararg1 && !scalararg2) //if ~scalararg1 & ~scalararg2
        {
            if(r1 != r2 || c1 != c2) return false;
        }
        if(!scalararg1 && !scalararg3) //if ~scalararg1 & ~scalararg3
        {
            if(r1 != r3 || c1 != c3) return false;
        }
        if(!scalararg3 && !scalararg2) //if ~scalararg3 & ~scalararg2
        {
            if(r3 != r2 || c3 != c2) return false;
        }

        uword rows = max(r1, r2, r3);
        uword columns = max(c1, c2, c3);

        if(scalararg1)
        {
            out1.set_size(rows, columns); //out1 = arg1(ones(rows,1),ones(columns,1));
            out1.fill( as_scalar(arg1) );
        }
        else out1 = arg1;

        if(scalararg2)
        {
            out2.set_size(rows, columns); //out2 = arg2(ones(rows,1),ones(columns,1));
            out2.fill( as_scalar(arg2) );
        }
        else out2 = arg2;

        if(scalararg3)
        {
            out3.set_size(rows, columns); //out3 = arg3(ones(rows,1),ones(columns,1));
            out3.fill( as_scalar(arg3) );
        }
        else out3 = arg3;

        //out4 =[];
    }
    else if( nparms == 4 )
    {
        uword r1 = arg1.n_rows, c1 = arg1.n_cols; //[r1 c1] = size(arg1);
        uword r2 = arg2.n_rows, c2 = arg2.n_cols; //[r2 c2] = size(arg2);
        uword r3 = arg3.n_rows, c3 = arg3.n_cols; //[r3 c3] = size(arg3);
        uword r4 = arg4.n_rows, c4 = arg4.n_cols; //[r4 c4] = size(arg4);
        bool scalararg1 = (r1 * c1 == 1); //scalararg1 = (prod(size(arg1)) == 1);
        bool scalararg2 = (r2 * c2 == 1); //scalararg2 = (prod(size(arg2)) == 1);
        bool scalararg3 = (r3 * c3 == 1); //scalararg3 = (prod(size(arg3)) == 1);
        bool scalararg4 = (r4 * c4 == 1); //scalararg4 = (prod(size(arg4)) == 1);

        if(!scalararg1 && !scalararg2) //if ~scalararg1 & ~scalararg2
        {
            if(r1 != r2 || c1 != c2) return false;
        }
        if(!scalararg1 && !scalararg3) //if ~scalararg1 & ~scalararg3
        {
            if(r1 != r3 || c1 != c3) return false;
        }
        if(!scalararg1 && !scalararg4) //if ~scalararg1 & ~scalararg4
        {
            if(r1 != r4 || c1 != c4) return false;
        }
        if(!scalararg3 && !scalararg2) //if ~scalararg3 & ~scalararg2
        {
            if(r3 != r2 || c3 != c2) return false;
        }
        if(!scalararg4 && !scalararg2) //if ~scalararg4 & ~scalararg2
        {
            if(r4 != r2 || c4 != c2) return false;
        }
        if(!scalararg3 && !scalararg4) //if ~scalararg3 & ~scalararg4
        {
            if(r3 != r4 || c3 != c4) return false;
        }

        uword rows = max( max(r1, r2), max(r3, r4) );
        uword columns = max( max(c1, c2), max(c3, c4) );

        if(scalararg1)
        {
           out1.set_size(rows, columns); //out1 = arg1(ones(rows,1),ones(columns,1));
           out1.fill( as_scalar(arg1) );
        }
        else out1 = arg1;

        if(scalararg2)
        {
            out2.set_size(rows, columns);
            out2.fill( as_scalar(arg2) ); //out2 = arg2(ones(rows,1),ones(columns,1));
        }
        else out2 = arg2;

        if(scalararg3)
        {
            out3.set_size(rows, columns);
            out3.fill( as_scalar(arg3) ); //out3 = arg3(ones(rows,1),ones(columns,1));
        }
        else out3 = arg3;

        if(scalararg4)
        {
            out4.set_size(rows, columns);
            out4.fill( as_scalar(arg4) ); //out4 = arg4(ones(rows,1),ones(columns,1));
        }
        else out4 = arg4;
    }
    return bRet;
}

double xrayAreaCorrection(double sd, double gamma)
{
    double areaXray = gNaN;
    rowvec sigma = zeros<rowvec>(301u);
    colvec ak = zeros<colvec>(301u);
    double g = xrayArea_g;

    //srcDir = getSetup('SourceDir');
    //load(fullfile(srcDir, 'setup', 'xrayArea.mat'), 'sigma', 'ak', 'g');
    QFile file1(QCoreApplication::applicationDirPath() + "/setup/sigma.txt");
    if(!file1.open(QIODevice::ReadOnly)) return areaXray;
    QTextStream in(&file1);
    file1.close();

    int n = 0;
    bool bOk;
    double temp;
    while(!in.atEnd() && n < 301)
    {
        temp = in.readLine().toDouble(&bOk);
        if(bOk) { sigma(n) = temp; ++n; }
    }
    if(n < 1) return areaXray;

    QFile file2(QCoreApplication::applicationDirPath() + "/setup/ak.txt");
    if(!file2.open(QIODevice::ReadOnly)) return areaXray;
    QTextStream in2(&file2);
    file2.close();

    n = 0;
    while(!in2.atEnd() && n < 301)
    {
        temp = in2.readLine().toDouble(&bOk);
        if(bOk) { ak(n) = temp; ++n; }
    }
    if(n < 1) return areaXray;

    // ind = min(find(sigma > sd));
    uvec t_ind = find(sigma > sd);
    uword ind;
    // If sd is too large (> 10), set ind to maximum value (N), i.e.,
    // X ray area = ak(N);
    // If sd is too small (< 0.4), set ind = 2, i.e.,
    // X ray area = ak(1)
    if(t_ind.is_empty()) //if isempty(ind)
    {
        ind = sigma.size() - 1u;
        sd = sigma(ind);
    }
    else
    {
        ind = min( t_ind );
        if(ind == 1)
        {
            ind = 2;
            sd = sigma(0);
        }
    }
    // Interpolation
    double area1 = ak(ind -1);
    double area2 = ak(ind);
    double area = area1 + (area2 - area1)/(sigma(ind) - sigma(ind-1))*(sd - sigma(ind-1));

    // For a fixed sigma, the X-ray area is directly proportional to gamma.
    areaXray = (gamma/g)*area;

    return areaXray;
}

mat voigt_gauss(double &width, mat x, double x0, double sigma, double gamma)
{
    x =(x-x0)/(sqrt(2)*sigma);

    double y = gamma/(2*sqrt(2)*sigma);
    double C1 = 1- 2*y/sqrt(pi_);
    double C2 = y/pi_;
    double C3 = pow(y,2);
    double C4 = -4*y/pi_;

    mat B; B << 0; // B=0;
    mat x2 = x % x; // x2=x.*x;

    // Infinite series, consider only 100 first elements
    double n2;
    for(int n=1; n<=100; ++n) //for n=1:100
    {
        n2 = n*n;
        //B = B + (1/n2)*(exp(-n2/4-x2) - 0.5*exp(-n2/4+n*x-x2) - 0.5*exp(-n2/4-n*x-x2));
        if(B.size() <= 1)
        {
            B = (1.0/n2)*(exp(-n2/4-x2) - 0.5*exp(-n2/4+n*x-x2) - 0.5*exp(-n2/4-n*x-x2));
        }
        else {
            B = B + (1.0/n2)*(exp(-n2/4-x2) - 0.5*exp(-n2/4+n*x-x2) - 0.5*exp(-n2/4-n*x-x2));
        }
    }

    // Voigt function
    mat v = exp(-x2)%(C1 + C2*x2 + C3*(1-2*x2)) + C4*B; //v = exp(-x2).*(C1 + C2*x2 + C3*(1-2*x2)) + C4*B;

    // Scaling factor A is calculted at x = 0 in such away that the Gaussian function
    // has an area of 1.
    double A = C1+C3;
    mat v_g = 1/(sigma*sqrt(2*pi_))*(v/A-exp(-x2));

    // Full-Width Half-Maximum of Breit-Wigner and Gaussian add linearly (not quadratically)
    width = 2*sqrt(2*log(2))*sigma + 0.5323*gamma;

    return v_g;
}

// =====================================================end========================================================= //

/* =====================================================拟合核心函数======================================================== */
void lParsePara(vec &beta0, field<urowvec> &beta_idx, field<vec> para, field<uvec> free_idx)
{
    uword narg = para.size();
    beta0.clear();
    beta_idx.set_size(narg); // beta_idx = {};

    // parse parameter arguments
    int lastpt = 0;
    for(uword argpt=0; argpt<narg; ++argpt)
    {
        // convert parameter vector to column
        int NFree = free_idx(argpt).size(); // NFree = length(idx);
        beta_idx(argpt) = RangeVec(lastpt, lastpt + NFree - 1); // index of free_idx{NPara} in beta
        lastpt = lastpt + NFree;

        if(NFree > 0)
        {
            beta0 = join_vert(beta0, para(argpt)( free_idx(argpt) ) );
        }
    }
}

void lParseOpt(int &sw, int &uw, mat &w, int &rs,
                              QStringList Opt, mat weights, int sz, umat idx)
{
    // default values
    sw = 0;
    //wt = "";
    uw = 0;
    w = ones(sz);
    rs = 0;

    // parse argument cell vector
    int narg = Opt.size(); // narg = length(arg);
    int argpt = 0;
    QString desc;
    while(argpt < narg)
    {
        // switch option descriptor
        desc = Opt[argpt];
        if(desc == "Waitbar")
        {
            sw = 1;
            /*wt = waitbarText;
            if ~isstr(wt)
                qDebug() << "Waitbar text must be string";
            end
            argpt = argpt + 2;*/
        }
        else if(desc == "Weights")
        {
            uw = 1;
            w = vectorise(weights);
            for(int ii=idx.size()-1; ii>=0; --ii)
            {
                w.shed_row(idx(ii));
            }
            //w(idx) = [];
            //argpt = argpt + 2;
        }
        else if(desc == "Restrict")
        {
            rs = 1;
            //argpt = argpt + 1;
        }
        else {
            qDebug() << QString("%1 is Invalid descriptor").arg(desc);
        }
        argpt = argpt + 1;
    }
}

field<vec> lBetaToPara(vec beta, field<vec> parain, field<uvec> free_idx, field<urowvec> beta_idx, int NPara)
{
    field<vec> para = parain;
    for(int i=0; i<NPara; ++i)// for i = 1 : NPara
    {
        if(!beta_idx(i).is_empty())// if any(beta_idx{i})
        {
            // para{i}(free_idx{i}) = beta(beta_idx{i});
            para(i)(free_idx(i)) = beta(beta_idx(i));
            /*for(int ii=0; ii<free_idx(i).size(); ++ii)
            {
                para(i)(free_idx(i)(ii)) = beta(beta_idx(i)(ii));
            }*/
        }
    }
    return para;
}

vec FitModel(QString model, vec x, field<vec> para)
{
    mat yfit;
    if(model == "calFcnEval")
    {
        if(!para.is_empty())
        {
            yfit = Independ::calFcnEval(x, para(0));
        }
    }
    else if(model == "wrapStepRatio")
    {
        if(para.size() > 13)
        {
            yfit = wrapStepRatio(x, para(0), para(1), para(2), para(3), para(4), para(5), para(6), para(7),
                      para(8), para(9), para(10), para(11), para(12), para(13));
        }
    }
    else if(model == "wrapBaseSpline")
    {
        if(para.size() > 2)
        {
            yfit = Independ::wrapBaseSpline(x, para(0), para(1), para(2));
        }
    }
    else if(model == "peakBaseLin")
    {
        if(para.size() > 3)
        {
            yfit = peakBaseLin(x, para(0), para(1), para(2), para(3));
        }
    }
    return vectorise(yfit);
}

mat wrapStepRatio(vec x, vec cx, vec cy, vec cdy, vec na, vec ct, vec fc,
                                 vec sr, vec lt, vec lta, vec ut, vec uta, vec db, vec rb, vec rd)
{
    //error(nargchk(15, 15, nargin))

    // change step ratio to step height and calculate peak + baseline
    vec st = na % sr;
    PiecePoly bpp = Independ::makeBasePP(cx, cy, cdy, ct, fc, st);
    mat y = Independ::peakBaseSpline(x, bpp, na, ct, fc, st, lt, lta, ut, uta, db, rb, rd);
    return y;
}

mat peakBaseLin(mat chan, vec BL, vec NA, vec C, vec F)
{
    // check input parameters
    //[BL, s] = doubleCheck(BL, [2, 1]);
    if(BL.size() != 2)
    {
        qDebug() << "Invalid BL in peakBaseLin";
    }

    uword NPeaks = NA.size();
    if(C.size() != NPeaks || F.size() != NPeaks)
    {
        qDebug() << "Peak parameters must match in size";
    }

    // sigma
    vec sg = F / sqrt(8*log(2)); //sg=F/sqrt(8*log(2));

    // baseline
    mat cnt = Matlab::polyval(BL, chan);
    // add peaks
    for(uword i=0; i<NPeaks; ++i) //for i = 1 : NPeaks
    {
        //cnt=cnt+NA(i)*normpdf(chan,C(i),sg(i));
        cnt = cnt + NA(i) * Independ::NORMPDF(chan, C(i), sg(i));
    }
    return cnt;
}
void PrintMat2xx(QString name, mat& dataMat)
{
    QFile file("C:/Result_Debug3/" + name);
	if (file.open(QIODevice::WriteOnly))
	{
		QTextStream out(&file);
		if (dataMat.is_vec())
		{
			for (int i = 0; i < dataMat.size(); ++i)
			{
				out << QString::number(dataMat(i), 'g', 10) << "\r\n";
			}
		}
		else {
			for (uword row = 0; row < dataMat.n_rows; ++row)
			{
				for (uword col = 0; col < dataMat.n_cols; ++col)
				{
					out << QString::number(dataMat(row, col), 'g', 10) << "\t";
				}
				out << "\r\n";
			}
		}
		file.close();
	}
}
bool fitFunction(field<vec> &para, double &X2, field<uvec> free_idx, vec x, vec y, QString model, QStringList Opt, QString waittext, mat weights)
{
    // default return values
    bool status = true;
    X2 = 0;
    vec r;
    mat J;

    // remove invalid data points
    uword size_y = y.size();
    uvec idx = find_nonfinite(x % y);
    if( any(idx) )
    {
        for(int i=idx.size()-1; i>=0; --i)
        {
            x.shed_row( idx(i) );
            y.shed_row( idx(i) );
        }
        size_y = size_y - idx.size();
        qDebug() << QString("Removed %d infinite data points").arg(idx.size());
    }
    if(size_y < 1) return status;

    // parse the parameter and free-index arguments
    vec beta0;
    field<urowvec> beta_idx;
    lParsePara(beta0, beta_idx, para, free_idx);

    int NPara = para.size();
    if(beta0.is_empty()) //if isempty(beta0)
    {
        //varargout = {status, para{:}, X2, r, J};
        return status;
    }

    // parse the optional arguments
    //[showwait, waittext, useweights, w, restricted] = lParseOpt(varargin, argpt, size(y), idx);
    vec w;
    int showwait, useweights, restricted;
    lParseOpt(showwait, useweights, w, restricted, Opt, weights, y.size(), idx);

    uword p = beta0.size(); // p = length(beta0); beta0 = beta0(:);

    // tolerance for relative parameter steps that must be met
    double betatol = 1.0E-4;
    // tolerance for relative square sum improvement that must be met
    double rtol = 1.0E-4;

    J = zeros<mat>(size_y, p);
    vec beta = beta0, betanew = beta;
    int maxiter = 100, iter = 0;
    // make sure that while-condition is false on first run
    double sse = 1, sseold = 0;
    double seps = sqrt( eps(1.0) );
    double srealmin = sqrt( Matlab::realmin() );
    double s10 = sqrt(10);
    vec zbeta = zeros<vec>(beta.size());
    mat eyep = eye<mat>(p, p);
    vec zerosp = zeros<vec>(p);

    if(showwait)
    {
        //hwait=waitbar(0,['Nonlinear Fitting: ',waittext]);
    }

    vec delta;
    vec yfit;
    while((norm((betanew-beta)/(beta+seps)) > betatol || abs(sseold-sse)/(sse+seps) > rtol) && iter < maxiter)
    {
        if(iter > 0) beta = betanew;

        iter = iter + 1;
        // get approx y with parameter beta and (weighted) residual

        para = lBetaToPara(beta, para, free_idx, beta_idx, NPara);
//        PrintMat2fvec("para_A_" + QString::number(iter) + ".txt", para);
        //yfit = feval(model,X,para{:});
        yfit = FitModel(model, x, para);
        //PrintMat("fitFunction_yfit.txt", yfit);

        r = (y - yfit);
        if(useweights)
        {
            r = r % w;
        }
        //PrintMat("fitFunction_r.txt", r);
        sseold = as_scalar(r.t() * r); //sseold = r'*r;

        // build approximation of (weighted) jacobian
        for(uword k=0; k<p; ++k) // for k = 1:p,
        {
            delta = zbeta;

            // for values near zero, calculate derivative with stepwidth eps
            // otherwise with sqrt(eps)*parametervalue
            if(fabs(beta(k)) < seps)
            {
                delta(k) = eps(1.0);
            }
            else {
                delta(k) = seps*beta(k);
            }

            para = lBetaToPara(beta + delta, para, free_idx, beta_idx, NPara);

            //yplus = feval(model,X,para{:});
            vec yplus = FitModel(model, x, para);
            //PrintMat("fitFunction_yplus.txt", yplus);

            if(useweights) {
                //J(:,k) = (yplus - yfit) % w / delta(k);
                J.col(k) = (yplus - yfit) % w / delta(k);
            } else {
                //J(:,k) = (yplus - yfit) / delta(k);
                J.col(k) = (yplus - yfit) / delta(k);

				vec vp = (yplus - yfit);
//				PrintMat2xx("vp.txt", vp);
            }
        }

//        PrintMat2fvec("para_B_" + QString::number(iter) + ".txt", para);
        // Levenberg-Marquardt type adjustment
        // Gauss-Newton step -> J\r
        // LM step -> inv(J'*J+constant*eye(p))*J'*r
        mat Jplus = join_vert(J, (1.0E-2)*eyep); // Jplus = [J;(1.0E-2)*eyep];
        vec rplus = join_vert(r, zerosp); // rplus = [r;zerosp];

        // step = Jplus\rplus;
        vec step;
        {
            uvec t_rplus_idx = find_finite(rplus);
            Jplus = MatRows(Jplus, t_rplus_idx);
            rplus = rplus(t_rplus_idx);
            step = solve(Jplus, rplus);

            betanew = beta + step;
        }

        // decrease step if a parameter went negative in restricted fitting
        if(restricted)
        {
            // restrict step: never reach 0, but stop at sqrt(eps)
            uvec farstepidx = find(betanew < srealmin);

            // some step went too far
            if(!farstepidx.is_empty()) // if ~isempty(farstepidx)
            {
                // all parameters too far
                if(farstepidx.size() == p) // if length(farstepidx) == p
                {
                    if(restricted == 2)
                    {
                        // stepped too far in all parameters twice in suite
                        status = 0;
                        break;
                    }
                    restricted = 2;
                }
                else {
                    restricted = 1;
                }

                double fac;
                for(uword i=0, j; i<farstepidx.size(); ++i) //for i = 1 : length(farstepidx)
                {
                    j = farstepidx(i);
                    fac = srealmin;
                    while(beta(j) <= beta(j) - fac)
                    {
                        fac = sqrt(fac);
                    }
                    step(j) = - beta(j) + fac;
                }
                betanew = beta + step;
//                PrintColvec("betaA_" + QString::number(iter) + ".txt", beta);
//                PrintColvec("stepA_" + QString::number(iter) + ".txt", step);
//                PrintColvec("betanewA_" + QString::number(iter) + ".txt", betanew);
            }
        }

        // get approx y, residual and square sum for new target
        para = lBetaToPara(betanew, para, free_idx, beta_idx, NPara);
//        PrintMat2fvec("para_C_" + QString::number(iter) + ".txt", para);
        //yfitnew = feval(model,X,para{:});
        vec yfitnew = FitModel(model, x, para);
        vec rnew = (y - yfitnew);
        if(useweights) {
            rnew = rnew % w;
        }
        sse = as_scalar(rnew.t() * rnew);

        // decrease step width if residual square sum went up (step too far)
        int iter1 = 0;
        while(sse > sseold & iter1 < 12)
        {
            step = step/s10;
            betanew = beta + step;
            para = lBetaToPara(betanew, para, free_idx, beta_idx, NPara);

            //yfitnew = feval(model,X,para{:});
            yfitnew = FitModel(model, x, para);
            rnew = (y - yfitnew);
            if(useweights) {
                rnew = rnew % w;
            }
            sse = as_scalar(rnew.t() * rnew);
            iter1 = iter1 + 1;
        }
        if(showwait) {
            //waitbar(iter/maxiter,hwait);
        }
    }

//    PrintMat2fvec("Ps_para.txt", para);
    // write beta to para
    para = lBetaToPara(beta, para, free_idx, beta_idx, NPara);

    // calculate mean squared relative distance of data points
    if(yfit.is_empty())//if ~exist('yfit', 'var')
    {
        //yfit = feval(model, X, para{:});
        yfit = FitModel(model, x, para);
    }
    r = (y - yfit) / (abs(y) + seps);
    X2 = as_scalar(r.t() * r) / size_y;

    if(showwait) {
        //waitbar(1,hwait);
        //delete(hwait);
    }

    if(iter == maxiter) {
        printf("fitFunction did NOT converge. Returning results from last iteration.");
        status = 0;
    }

    //varargout = {status, para{:}, X2, r, J};
    return status;
}

/* =====================================================求 BackgroundCount======================================================== */
mat abkcnt(rowvec BL, vec CT, vec FC, double k)
{
    //if nargin < 4, k=1.25; end
    vec L = CT - k * FC;
    vec H = CT + k * FC;
    return vecsum( BL, L, H);
}

mat vecsum(rowvec y, vec l, vec h)
{
    // sizes of input vectors
    int n = y.size();
    int m = l.size();

    if(h.size() != m)
    {
        qDebug() << "summation indizes must be of same size";
    }
    // reshape to mx1 vectors
    //l.reshape(m, 1); //[l, s1] = doubleCheck(l, [m, 1]);
    //h.reshape(m, 1); //[h, s2] = doubleCheck(h, [m, 1]);

    // return
    mat s = zeros(size(l));
    if(y.is_empty()) return s;

    //nocnt = ( l >= h | ~isfinite(l) | ~isfinite(h) );
    //uvec nocnt = OrUMat( OrUMat( l>=h, IsInf(l) ), IsInf(h) );
    //uvec nocnt = find(l>=h && (l<1 || l>n) && (h<1 || h>n) );
    //nocnt.print("nocnt = ");

    // integer and fractional part
    vec rl = round(l);

    // 0.5 added, since e.g. for integer l, only y(l)/2 accounts to the sum
    vec dl = l - rl + 0.5;

    // flag where the fractional part has influence
    //uvec bl = AndUMat( AndUMat( rl > 0, rl <= n ), dl != 0 );
    //bl.print("bl = ");

    // same for upper border
    vec rh = round(h);

    vec dh = rh - h + 0.5;

    //uvec bh = AndUMat( AndUMat( rh > 0, rh <= n), dh != 0 );
    //bh.print("bh = ");

    for(int i=0; i<m; ++i) //for i = 1 : m
    {
        if( is_finite(rl(i)) && is_finite(rh(i)) && rl(i) < rh(i) ) //if ~nocnt(i)
        {
            // sum from low to high
            //s(i) = sum( y([max(1, rl(i)) : min(n, rh(i))]) );
            if(min(n,rh(i)) < 0) throw std::string(" Row::cols(): indices out of bounds or incorrectly used");
            s(i) = sum( y.cols( max(1, rl(i))-1, min(n, rh(i))-1 ) );

            // fractional part at low border
            if( rl(i) > 0 && rl(i) <= n && dl(i) != 0 ) // if( bl(i) )
            {
                s(i) = s(i) - dl(i) * y(rl(i)-1);
            }

            // fractional part at high border
            if( rh(i) > 0 && rh(i) <= n && dh(i) != 0 ) // if( bh(i) )
            {
                s(i) = s(i) - dh(i) * y(rh(i)-1);
            }
        }
    }

    // size of l
    s.reshape( size(l) ); //s = reshape(s, sz);
    return s;
}

/* =====================================================插入拟合======================================================== */
bool lMakeInterp(mat &p, mat &perr, vec x, vec y, vec err)
{
    bool s = true;

    p << 1;
    perr.clear();

    if (x.size() != y.size())
    {
        qDebug() << "The size of x and y must be same.";
        return false;
    }
    if(err.is_empty()) err.fill(gNaN);

    if (s > 0)
    {
        uword l = x.size();
        p = ones(2*l+1, 1);
        perr = zeros(2*l, 1);
        for(uword i=0; i<l; ++i) //for i = 1 : l
        {
            p(2*i+1) = x(i); //p(2*i) = x(i);
            // perr(2*i-1) = 0;
            p(2*i+2) = y(i); //p(2*i+1) = y(i);
            perr(2*i+1) = err(i)/max(y(i), eps(1.0)); //perr(2*i) = err(i)/max(y(i), eps);
        }
    }
    //p.print("p = ");
    //perr.print("perr = ");
    return s;
}

/* =====================================================刻度拟合======================================================== */
bool calFitPara(vec &p, vec &perr, CalibrationType cal, vec x, vec y, vec err, vec p0, vec perr0, urowvec idx)
{
    bool bRet = false;

    if(x.is_empty()) return false;

    // parse optional parameters
    bool need_para, need_idx;
    if(p0.is_empty()) // if nargin == 4
    {
        need_para = true;
        need_idx = true;
    }
    else if(perr0.is_empty()) // elseif nargin == 5
    {
        need_para = false;
        need_idx = true;
    }
    else if(!idx.is_empty()) // elseif nargin == 7
    {
        need_para = false;
        need_idx = false;
    }
    else {
        qDebug() << "Wrong number of input arguments!";
        return false;
    }

    //% check if calibration type implemented, raise error if not
    //lCheckCalImplemented(cal)

    // get initial parameter depending on calibration
    colvec pIni;
    if(need_para)
    {
        bRet = lGetIniPara(pIni, cal, x, y);
        if(!bRet) return bRet;
    }
    else { pIni = p0; }
    //pIni.print("pIni = ");
    //PrintMat("calFitPara_pIni.txt", pIni);

	//pIni.print(" 1##Debug calFitPara");

    // check that parameter is reasonable
    bRet = calParaCheck(pIni);
    if(!bRet)
    {
        return bRet;//return lMakeInterp(p, perr, cal, calData);
    }

	//pIni.print(" 2##Debug calFitPara");

    // get parameter uncertainty estimate and free-index, if not specified
    urowvec freeIdx;
    colvec pErrIni;
    if(need_idx)
    {
        freeIdx = RangeVec(2, pIni.size()); // freeIdx = [2:length(pIni)];
        pErrIni.set_size(freeIdx.size());
        pErrIni.fill(gNaN);
    }
    else {
        freeIdx = idx;
        pErrIni = perr0;
    }
    //freeIdx.print("freeIdx = ");
    //PrintMat("calFitPara_pErrIni.txt", pErrIni);

    // handle special cases
    //[is_special, s, msg, p, perr] = lHandleSpecial(cal, x, y, err, pIni, pErrIni, freeIdx, need_para);
	//pIni.print(" 3##Debug calFitPara");
    bool is_special;
    bRet = lHandleSpecial(is_special, p, perr, cal, x, y, err, pIni, pErrIni, freeIdx, need_para);
    //PrintMat("calFitPara_p.txt", p);
    //PrintMat("calFitPara_perr.txt", perr);
    if(!is_special)
    {
        //[s, msg, p, perr] = lFitStandard(x(:), y(:), err(:), pIni, pErrIni, freeIdx);
        double chi;
        bRet = lFitStandard(p, perr, chi, x, y, err, pIni, pErrIni, freeIdx);
        //PrintMat("calFitPara_p.txt", p);
        //PrintMat("calFitPara_perr.txt", perr);
    }
	//pIni.print(" 4##Debug calFitPara");
    return bRet;
}

bool lGetIniPara(colvec &pIni, CalibrationType cal, colvec x, colvec y)
{
    bool bRet = true;
    colvec p;
    uword l = x.size();
    double A, B;
    switch(cal)
    {
    case Cal_Energy: // fit first or third order polynomial
        pIni << 2;
        if(l >= 6)
        {
            p = Matlab::polyfit(x, y, 3);
            pIni = join_vert(pIni, flipud(p));
        }
        else if(l >= 2)
        {
            p = Matlab::polyfit(x, y, 1);
            pIni = join_vert(pIni, flipud(p));
        }
        else if(l==1 && x(0)>0)
        {
            pIni << 2 << 0 << y(0)/x(0);
        }
        else {
            bRet = false;
            qDebug() << "to few data points";
        }
        break;

    case Cal_Resolution: // fit squareroot of first order polynomial
        if(l >= 3)
        {
            p = Matlab::polyfit(x, pow(y, 2), 1);
            pIni << 4;
            pIni = join_vert(pIni, flipud(p));
        }
        else {
            bRet = false;
            qDebug() << "to few data points";
        }
        break;

    case Cal_Efficiency:
    case Cal_Tot_efficiency:
        if(y.is_empty()) A = 0.3;
        else A = y.max();
        pIni << 94 << A << 30 << 3 << 50 << 120 << 0.8;
        break;

    case Cal_Tail: //case {'tail', 'tail_right'}
    case Cal_Tail_right:
       pIni << 99 << 0 << -3000 << .01;
       break;

    case Cal_Step_ratio: //case {'step_ratio', 'tail_alpha', 'tail_right_alpha'}
    case Cal_Tail_alpha:
    case Cal_Tail_right_alpha:
       pIni << 98 << mean(y);
       break;

    case Cal_HalflifeAnalysis: //case 'halflifeanalysis'
       if(y.is_empty()) B = 1;
       else B = y.max();
       pIni << 97 << 0 << B << 1;
       break;

    default:
        bRet = false;
        qDebug() << "invalid calibration string " << cal;
        break;
    }

    return bRet;
}
void PrintColvec1(QString name, colvec& dataMat)
{
	QFile file(QCoreApplication::applicationDirPath() + "/Result_Debug3/" + name);
	if (file.open(QIODevice::WriteOnly))
	{
		QTextStream out(&file);
		for (int i = 0; i < dataMat.size(); ++i)
		{
			out << QString::number(dataMat(i), 'g', 10) << "\r\n";
		}

		file.close();
	}
}
bool lHandleSpecial(bool &issp, colvec &p, colvec &perr,
                                   CalibrationType cal, colvec x, colvec y, colvec err,
                                   colvec p0, colvec pErr0, urowvec idx, bool np)
{
    bool bRet = false;
    issp = true; // is_special
    p.clear();
    perr.clear();
    if(np)
    {
        // no parameters specified, certain calibrations need special handling
        if(err.is_empty())
        {
            err = gNaN * ones<colvec>(x.size());
        }
        uvec dIdx;
        colvec p1, perr1, pErr1, p2, perr2, p3, perr3, perr0;
        urowvec idx1, idx2, idx3;
        double ECut, chi, chi1, chi2, chi3, A = gNaN, n, E1;
        bool s1, s2, s3;
        switch(cal) // special cases should go here
        {
        case Cal_Resolution:
//			PrintColvec1("lHandleSpecial_x.txt", x);
//			PrintColvec1("lHandleSpecial_y.txt", y);
//			PrintColvec1("lHandleSpecial_err.txt", err);
//			PrintColvec1("lHandleSpecial_p0.txt", p0);

            // fit [a, b] to data up to ECut
            ECut = 1500;
            p0 = p0.rows(0, 2); // p0 = p0(1:3);
            dIdx = find(x < ECut);
            //[s1, msg, p1, perr1, chi1] = lFitStandard(x(dIdx), y(dIdx), err(dIdx), p0, [NaN, NaN], [2, 3]);
            pErr1 << gNaN << gNaN;
            idx1 << 2 << 3;

            //p0.print("p0 = ");
            //dIdx.print("dIdx = ");
            s1 = lFitStandard(p1, perr1, chi1, x(dIdx), y(dIdx), err(dIdx), p0, pErr1, idx1);
//			PrintMat2xx("lHandleSpecial_p1.txt", p1);
//			PrintMat2xx("lHandleSpecial_perr1.txt", perr1);

            if(s1) // fix a, fit [b] and [b c] to
            {
                //[s2, garb, p2, perr2, chi2] = lFitStandard(x, y, err, p1, perr1, 3);
                idx2 << 3;
                s2 = lFitStandard(p2, perr2, chi2, x, y, err, p1, perr1, idx2);
//				PrintMat2xx("lHandleSpecial_p2.txt", p2);
//				PrintMat2xx("lHandleSpecial_perr2.txt", perr2);

                p3 = p1, perr3 = perr1;
                idx3 << 3 << 4;
                if(p3.size() < 4) p3.resize(4);
                p3(3) = 1e-6;   // p3(4) = 1e-6;
                if(perr.size() < 4) perr.resize(4);
                perr(3) = gNaN; // perr(4) = NaN;
                //[s3, garb, p3, perr3, chi3] = lFitStandard(x, y, err, p3, perr3, [3, 4]);
                s3 = lFitStandard(p3, perr3, chi3, x, y, err, p3, perr3, idx3);
//				PrintMat2xx("lHandleSpecial_p3.txt", p3);
//				PrintMat2xx("lHandleSpecial_perr3.txt", perr3);
                if(s3)
                {
                    bRet = true;
                    if(s2 && chi2 < chi3)
                    {
                        p = p2;
                        perr = perr2;
                    }
                    else {
                        p = p3;
                        perr = perr3;
                    }
                }
                else if(s2) {
                    bRet = true;
                    p = p2;
                    perr = perr2;
                }
            }
            break;

        case Cal_Step_ratio:
            // fit A+(E1/x)^n with fixed n=2.95
            if(!y.is_empty()) A = y.min() / 2;
            if( !is_finite(A) ) A = 0.0001;
            n = 2.95;
            E1 = 10;
            p0 << 96 << A << E1 << n;
            perr0 << gNaN << gNaN << gNaN;
            idx << 2 << 3;
            //[s1, msg, p1, perr1, chi1] = lFitStandard(x, y, err, p0, perr0, idx);
            s1 = lFitStandard(p1, perr1, chi1, x, y, err, p0, perr0, idx);
            if(s1)
            {
                bRet = true;
                p = p1;
                perr = perr1;
            }
            else {
                idx << 3;
                //[s, msg, p, perr] = lFitStandard(x, y, err, p0, perr0, idx);
                bRet = lFitStandard(p, perr, chi1, x, y, err, p0, perr0, idx);
            }
            break;

        case Cal_Efficiency:
        case Cal_Tot_efficiency:
            // HAE (1-2) with fixed k = 3; accept results always
            A = 0.3;
            if(!y.is_empty()) A = y.max();
            p0 << 94 << A << 30 << 3 << 50 << 120 << 0.8;
            pErr0 << gNaN << gNaN << gNaN << gNaN << gNaN << gNaN;
            idx << 2 << 3 << 5 << 6 << 7;
            //[s, msg, p, perr] = lFitStandard(x, y, err, p0, perr0, idx);
            // if failed, try a second fit (more iterations)
            bRet = lFitStandard(p, perr, chi, x, y, err, p0, pErr0, idx);
            if(!bRet)
            {
                //[s, msg, p, perr] = lFitStandard(x, y, err, p, perr, idx);
                bRet = lFitStandard(p, perr, chi, x, y, err, p, perr, idx);
            }
            break;

        default: // no special cases
            issp = false;
        }
    } else {
        issp = false; // parameters user specified, no special cases
    }
    return bRet;
}

bool lFitStandard(colvec &p, colvec &perr, double &x2, colvec x, colvec y, colvec err, colvec p0, colvec pErr0, urowvec idx)
{
    bool bRet = true;
    x2 = qInf();

    // calculate weights, if all error estimates are finite and positive
    QStringList optarg;
    colvec w;
    if(err.is_empty() || !err.is_finite() || any(err < 2*eps(1.0))) //if isempty(err) | any(~isfinite(err) | err < 2*eps)
    {
        //optarg = {};
    }
    else {
        // errors are punishments
        w = 1 / err;
        w = w / mean(w);
        optarg.append("Weights");
    }

    // fit
    //[s, p, x2, resid, J] = fitFunction(x, y, 'calFcnEval', p0, idx, optarg{:});
    field<vec> para;
    field<uvec> free_idx;
    para << p0;
//	p0.print("Standard_p0=");

    free_idx << vectorise(idx) - 1;
    bRet = Independ::fitFunction(para, x2, free_idx, x, y, "calFcnEval", optarg, "", w);
    //PrintMat("lFitStandard_resid.txt", resid);
    //PrintMat("lFitStandard_J.txt", J);
    if(!para.is_empty()) p = vectorise(para(0));

//	p.print("Standard_p=");

    if(!bRet) {
       qDebug() << "fitting didn't converge";
    }
    // assign error estimate
    perr = pErr0;
    if(perr.is_empty())
    {
        perr = calEstimateError(p, x, y, err);
    }
    else {
        urowvec idx_temp = idx-2u;
        //MatAssign(perr, idx_temp, calEstimateError(p, x, y, err, idx_temp));
        perr.resize(idx_temp.max()+1u);
        vec t_perr = calEstimateError(p, x, y, err, idx_temp);
        if(t_perr.size() == idx.size()) perr(idx_temp) = t_perr;
    }
    return bRet;
}

colvec calEstimateError(colvec p, colvec x, colvec y, colvec err, urowvec idx/*, mat &RMS, double &X2, mat &cv*/)
{
    /* The error estimates for calibration parameter coefficients are derived as follows:
     * perr(i) = RMS * JInv(i, i)
     *
     * where JInv is the inverse of J'*J, J being the jacobian matrix, and RMS is the mean square residual
     * RMS = sqrt(sum( ((yc - y) / y) ^2) / dgf)
     *
     * where yc = calFcnEval(x, p) denotes the data approximation described by the coefficients,
     * and dgf is the degree of freedom (number of data points minus number of parameter coefficients)
     *
     * The chi-square value is
     * X2 = sum( (y - yc)^2 / yc ); */

    colvec perr;
    mat cv;
    double RMS;
    double X2;
    if(p.size() < 2) return perr;

    // get default if idx missing
    if(idx.is_empty())
    {
        idx = RangeVec(0, p.size()-2); //idx = [1:length(p)-1];
    }
    colvec coeff = p.rows(1, p.size()-1);
    uword nParTot = coeff.size();
    uword nParFree = idx.size();
    uword nPts = x.size();

    colvec yc = Independ::calFcnEval(x, p);

    double seps = sqrt(eps(1.0));

    // calculate relative residual
    colvec w;
    if(err.is_empty() || !err.is_finite() || any(err <= 0))//if isempty(err) || any(~isfinite(err)) || any(err <= 0)
    {
        w = yc;
    }
    else {
        w = err;
    }
    MatAssign(w, find(w < seps), seps); // w(w<seps) = seps;

    colvec r = (y - yc) / w;

    // chi square value
    X2 = sum( (y-yc) % (y-yc) / w); // X2 = sum( (y-yc) .* (y-yc) ./ w);

    // mean residual, scaled by degrees of freedom
    RMS = norm(r)/sqrt(max(1, nPts -nParFree)); // RMS = norm(r)/sqrt(max(1, nPts -nParFree));
    mat J = zeros<mat>(nPts, nParFree); // J = zeros(nPts, nParFree);

    if(idx.is_empty())
    {
       perr.clear();
       cv.clear();
    }
    else { // calculate the jacobian
        for(uword i=0; i<nParFree; ++i)
        {
            int j = idx(i);
            double d;
            colvec delta = zeros<colvec>(nParTot+1);
            if(abs(coeff(j)) < seps) d = seps;
            else d = seps * coeff(j);

            delta(j+1) = d;
            colvec yi = Independ::calFcnEval(x, p + delta);
            J.col(i) = (yi - yc) / (d * w);
        }
        mat J2 = trans(J) * J;
        if(abs(det(J2)) < 1E-16)
        {
            J2.fill(qInf());
            qDebug() << "In function calEstimateError, J2 is singular When exec inv(J2)";
        }
        else J2 = inv(J2);

        // relative parameter error estimates
        vec tmp = J2.diag();
        perr = RMS * sqrt(J2.diag()) / max(seps, abs(coeff(idx)));

        cv = cov(J); // covariance matrix?
    }
    return perr;
}

vec calDefIniPara(int ft, vec x, vec y)
{
    vec p;
    switch (ft)
    {
        case 1: // interpolation
        {
            if (x.size() < 2) p << 0 << 0 << 1 << 1;
            else p = vectorise( join_vert(x.t(), y.t()) ); // row vector [x1, y1, x2, y2, ...]
            break;
        }
        case 2: p = lFitPoly(x, y); break; // polynomial
        case 3: p = lFitPoly(sqrt(x), y); break; // square root polynomial
        case 4: p = lFitPoly(x, pow(y,2)); break; // square root of poly
        case 5: // HT-efficiency: TBD, currently hardcoded
        {
            double A = (y.is_empty() ? 0.1 : max(y)); // A = lMax(y);
            p << A << 40 << 70 << 3 << 0.8;
            break;
        }
        case 6: p = lFitPoly(log(x), log(y)); break; // poly in log(y) against log(x)
        case 7: p = lFitPoly(1./x, log(y)%pow(x,2)); break; // polynomial in log(y) against 1/x
        case 8: // polynomial in log(y) against log(c/x)
        {
            double c = 1;
            vec pm = lFitPoly(log(c/x), log(y));
            p << c;
            p = join_vert(pm, p);
            break;
        }
        case 9: p << 1 << 1 << 1 << 1; break; // inverse exponential --TBD
        case 93: // fix shoulder efficiency
        {
            double S = (y.is_empty() ? 0.1 : max(y)); // S = lMax(y, 0.1);
            p << S << 30 << 3 << 120 << 1;
            break;
        }
        case 94: // free shoulder efficiency
        {
            double S = (y.is_empty() ? 0.1 : max(y)); // S = lMax(y, 0.1);
            p << S << 30 << 3 << 160 << 80 << 1;
            break;
        }
        case 95: // three component efficiency
        {
            double S = (y.is_empty() ? 0.1 : max(y)); // S = lMax(y, 0.1);
            p << S << 30 << 3 << 160 << 80 << 1 << 3000 << 1;
            break;
        }
        case 96: p << 1 << 60 << 1; break; // inverse power -- TBD
        case 97: p << 1e-4 << 1 << 0.1; break; // exponential sum -- TBD
        case 98: // constant
        {
            if(y.is_empty()) p << 1;
            else p << mean(y);
            break;
        }
        case 99: // liner with cutoff;
        {
            double E0 = 150;
            uvec idx = find(x > E0);
            vec pm = lFitPoly(x(idx), y(idx));
            p << E0;
            if(pm.size() > 1)
            {
                p << E0 << pm(0) << pm(1);
            }
            else p = join_vert(p, pm);
            break;
        }
        default: qDebug() << QString("Function type %1 not yet implemented").arg(ft);
    }
    return p;
}

vec lFitPoly(vec x, vec y)
{
    vec p;
    switch (x.size())
    {
        case 0: p << 0 << 1; break;
        case 1: p << as_scalar(y-x) << 1; break;
        case 2:
        case 3:
        case 4:
        case 5: p = flipud( Matlab::polyfit(x, y, 1) ); break;
        default: p = flipud( Matlab::polyfit(x, y, 3) );
    }
    return p;
}

}