#ifndef INDEPENDENTALG_H
#define INDEPENDENTALG_H

#include "gamma_alg_global.h"

namespace Independ {

// 峰曲线绘制相关函数，功能是取比1道更密集的点
//stdvec interp1(PeakInfo peak, arma::vec regBase, arma::vec u);
arma::mat pchip(arma::mat x, arma::mat y, arma::mat xx = arma::mat(0u, 0u));

/* function d = pchipslopes(x,y,del)
 * PCHIPSLOPES  Derivative values for shape-preserving Piecewise Cubic Hermite Interpolation.
 * d = pchipslopes(x,y,del) computes the first derivatives, d(k) = P'(x(k)). */
arma::rowvec pchipslopes(arma::rowvec x, arma::rowvec y, arma::rowvec del);

/* @syntax [y, s] = calFcnEval(x, p)
 * @param x independent variable at which the calibration function shall be
 * evaluated
 * @param p parameter vector, where p(1) is the function type indicator
 * @return y dependent variable
 * @return s status 0: invalid parameters, 1: valid parameters */
arma::mat calFcnEval(arma::mat x, arma::colvec para);

/* list:
 * ID   p                     equation
 *  1   x1 y1 x2 y2 ...	      interpolation
 *  2   a0 a1 ...             polynomial
 *  3   a0 a1 ...             square root polynomial
 *  4   a0 a1 ...             square root of polynomial
 *  5   A E1 E2 k n           A*exp(-(E1/x)^k)*(1-exp(-(E2/n)^n))
 *  6   a0 a1 ...             ploynomial in log(y) against log(x)
 *  7   a0 a1 ...             polynomial in log(y) against 1/x
 *  8   a0 a1 ... c           polynomial in log(y) against log(c/x)
 *  9   a b ex ey             inverse exponential 1/(a*x^(-ex)+b*x^(-ey))
 * 93   A E0 E1 k m           HAE(1-2) efficiency
 * 94   A E0 E2 E3 k n        HAE(1-3) efficiency
 * 95   A E0 E1 E2 E3 l m n   HAE(1-2-3) Efficiency
 * 96   A B n                 inverse power
 * 97   a0 a1 l1 a2 l2 ...    exponential sum
 * 98   c                     constant
 * 99   ECut E0 k		      linear with cutoff
 *
 * @syntax [s, msg, ft, para] = calParaCheck(p)
 * @param p calibration parameter vector
 * @return s status 0: invalid parameters, 1: valid parameters
 * @return msg error message string
 * @return ft function type identifier
 * @return para parameter vector */
bool calParaCheck(arma::vec p);

/* wrapper function for spline baseline contribution
 *
 * returns the spline baseline component, assuming no steps above control
 * points. I.e. the spline component has to be substracted from the
 * control points beforehand.
 *
 * @syntax y = wrapBaseSpline(x, cx, cy, cdy)
 * @param x channels to evaluate spline baseline contribution
 * @param cx control points x position
 * @param cy control points y position
 * @param cdy control points slopes
 * @return y spline baseline contribution at channels x */
arma::mat wrapBaseSpline(arma::mat x, arma::mat cx, arma::mat cy, arma::mat cdy);

/*function bpp = makeBasePP(cx, cy, cdy, pc, pfwhmc, pstep)
 * turn baseline control data points into piecewise polynomial
 *
 * The baseline is defined as the sum of peak steps plus a spline component.
 * This function returns the spline component as a piecewise polynomial, that
 * can efficiently be evaluated using ppval, saving repetition of the
 * calculations performed here.
 *
 * @syntax bpp = makeBasePP(cx, cy, cdy, pc, pfwhmc, pstep)
 * @param cx baseline control points x values (channels)
 * @param cy baseline control points y values (counts)
 * @param cdy derivative at breakpoints, NaN at normal control points
 * @param pc peak centroids
 * @param pfwhmc peak FWHMa in channel units
 * @param pstep peak step heights
 * @return bpp piecewise polynomial describing spline component of baseline*/
PiecePoly makeBasePP(arma::vec cx, arma::vec cy, arma::mat cdy, arma::mat pc, arma::mat pfwhmc, arma::mat pstep);

/* @syntax s = peakSteps(x, pC, pFWHM, pSt)
 * @param x channels at which step shall be evaluated
 * @param pC centroids of peaks
 * @param pFWHM full widths at half maximum (FWHM) of peaks
 * @param pSt step heights
 * @return s step counts */
arma::mat peakSteps(arma::mat x, arma::mat pC, arma::mat pFWHM, arma::mat pSt);

/* @syntax y = gaussInt(x, s)
 * @param x upper border of integration interval
 * @param s sigma (square root of variance) of the gaussian bell function
 * @return y integral from -inf to x of G(t,s) dt */
arma::mat gaussInt(arma::mat x, double s);


/* @syntax  y = peakBaseSpline(x,bpp,na,ct,fc,st,lt,lta,ut,uta,bw,rb,rd)
 * @param x channels at which approximate spectrum shall be evaluated
 * @param bpp spline component of the baseline, as piecewise polynomial to be
 *  evaluated with ppval
 * @param na peak net areas
 * @param ct peak centroids (channels)
 * @param fc peak fwhms in channel units
 * @param st peak step heights
 * @param lt peak tailing parameters (lower tail)
 * @param lta peak tailing alpha (lower tail)
 * @param ut peak tailing parameters (upper tail)
 * @param uta peak tailing alpha (upper tail)
 * @param bw breit-wigner width (channels)
 * @param rb recoil-induced tail slope (1/channels)
 * @param rd delta-channels for left-hand neutron-peak model (channels)
 * @return y counts of the spectrum approximation */
arma::mat peakBaseSpline(arma::mat x, PiecePoly bpp, arma::vec na, arma::vec ct, arma::vec fc, arma::vec st,
          arma::vec lt, arma::vec lta, arma::vec ut, arma::vec uta, arma::vec bw, arma::vec rb, arma::vec rd);

/* function y = peakShape(x, ct, s, lt, lta, ut, uta, g, rb, rd)
 * peakShape - aatami's peak shape
 *
 * Peak shape can be one of the following models:
 * 1 pure gaussian
 * 2 pure voigt (breit-wigner convoluted with gaussian)
 * 3 1 or 2 with left or right tailing
 * 4 neutron-induced peak model
 *
 * 1) the counts are normally distributed around the centroid ct with
 * standard deviation s.
 * 2) the initial counts have breit-wigner distribution with gamma-parameter g.
 *  detector response is a convolute with gaussian with standard deviation s.
 *  The difference between pure gaussian and voigt is added.
 * 3) below Channel C0 = ct-s/lt, the gaussian part of the peak is replaced
 *   by a subexponential tail of the form b*exp(-k*(c'-x)^lta)
 *   where 0<lta<=1 is the tailing exponent.
 *   similarly, above channel C1 = ct+s/ut, the gaussian part is replaced
 *   by a subexponential tail with paramters ut and uta.
 * 4) if rb and rd are finite, the peak has the shape of the convolution of
 *   the following function with a gaussian with standard deviation s:
 *   0                   for x < ct-rd
 *   n * ( x -(ct-rd) )  for ct - rd < x < ct
 *   n * exp(-rb*(x-ct)) for x > ct
 *   where n = 1/ (rd/2 + 1/rb) is the normalization factor.
 *
 * @syntax y = peakShape(x, ct, s, lt, lta, ut, uta, g)
 * @param x channel at which function shall be evaluated
 * @param ct centroid of gaussian peak
 * @param s standard deviation of gaussian peak
 * @param lt tail parameter for lower (left) tail
 * @param lta decrease factor for lower (left) tail
 * @param ut tail parameter for upper (right) tail
 * @param uta decrease factor for upper (right) tail
 * @param g breit-wigner width, gamma (channels)
 * @return y peak count expectation in channels x */
arma::mat peakShape(arma::mat x, double ct, double s, double lt, double lta, double ut,
              double uta, double g, double rb = qQNaN(), double rd = qQNaN());

/* @syntax ct = neutronModel(x, C, sigma, beta, dC)
 * @param x channels at which peak shape shall be evaluated
 * @param C centroid of unconvoluted peak
 * @param sigma sigma of the convolution gaussian
 * @param beta beta of the exponential tail
 * @param dC width of the linear left tail
 * @return ct densitiy function of peak in channels x */
arma::mat neutronModel(arma::mat x, double C, double sigma, double beta, double dC);

/* @syntax y = gaussianModel(x, ct, s, lt, lta, ut, uta, g)
 * @param x channels at which to evaluate peak shape
 * @param ct peak centroid
 * @param s gaussian sigma
 * @param lt left tail parameter
 * @param lta left tail exponent
 * @param ut right tail parameter
 * @param uta right tail expontent
 * @param g gamma (channels) of the breit-wigner
 * @return y peak shape at channels x */
arma::mat gaussianModel(arma::mat x, double ct, double s, double lt,
                        double lta, double ut, double uta, double g);

/* convolution of exponential function with gaussian
 * area = 1/beta;
 * @syntax ct = convExponential(x, sigma, beta)
 * @param x channels
 * @param sigma sigma of the gaussian
 * @param beta slope of the exponential
 * @return ct convolution */
arma::mat convExponential(arma::mat x, double sigma, double beta);

/* convolution of triangular function with gaussian
 * area = dC/2
 * @syntax ct = convLin(x, sigma, dC)
 * @param x channels
 * @param sigma sigma of the gaussian
 * @param dC width of the triangle */
arma::mat convLin(arma::mat x, double sigma, double dC);

/* convert tailing parameter to distance of junction
 *
 * @syntax J = lTailToJunct(s, t)
 * @param s sigma
 * @param t tailing paramter
 * @return J distance of junction point from centroid */
double lTailToJunct(double s, double t);

/* @syntax alpha_use = lCheckAlpha(alpha)
 * @param alpha tailing alpha parameter
 * @return alpha_use alpha if in range, 1, if above, and lower treshold else */
double lCheckAlpha(double alpha);

/* NORMPDF Normal probability density function (pdf).
 *   Y = NORMPDF(X,MU,SIGMA) Returns the normal pdf with mean, MU,
 *   and standard deviation, SIGMA, at the values in X.
 *
 *   The size of Y is the common size of the input arguments. A scalar input
 *   functions as a constant matrix of the same size as the other inputs.
 *
 *   Default values for MU and SIGMA are 0 and 1 respectively. */
arma::mat NORMPDF(arma::mat X, double MU = 0, double SIGMA = 1);

/* function [errorcode,out1,out2,out3,out4] = distchck(nparms,arg1,arg2,arg3,arg4)
 * DISTCHCK Checks the argument list for the probability functions. */
bool distchck(arma::mat& out1, arma::mat &out2, arma::mat &out3, arma::mat &out4, int nparms,
              arma::mat arg1, arma::mat arg2 = arma::mat(0, 0),
              arma::mat arg3 = arma::mat(0, 0), arma::mat arg4 = arma::mat(0, 0));

/* @syntax areaXray = xrayAreaCorrection(sd, gamma)
 * @param sd standard deviation of the Gaussian width of the peak (channels)
 * @param gamma Breit-Wiegner width (channels) of the X-ray line
 * @return areaXray Area of the non-Gaussian component of an X-ray peak */
double xrayAreaCorrection(double sd, double gamma);

/* @syntax [v_g, width]=voigt_gauss(x,x0,sigma,gamma)
 * @param x channels at which the function shall be evaluated
 * @param x0 Centroid of the Gaussian distribution (channels)
 * @param sigma Standard Deviation of the Gaussian distribution (in channels)
 * @param gamma Width of the Lorentzian function (typically 50 - 100 eV)
 *   in channels
 * @return v_g (Voigt - Gauss)
 * @return width Full-Width Half-Maximum of the Voigt function, in channels */
arma::mat voigt_gauss(double &width, arma::mat x,
                 double x0, double sigma, double gamma);

arma::mat calDerivEval(arma::mat x, arma::colvec p);

/* @syntax [para,beta0,free_idx,beta_idx,NPara,argpt] = lParsePara(arg,argptin)
 * @param arg cell vector of arguments of fitFunction
 * @param argptin pointer to first parameter argument (index)
 * @return para cell vector of parameters
 * @return beta0 vector of initial values of free parameters
 * @return free_idx index of beta0 into para
 * @return indizes of arguments into para */
void lParsePara(arma::vec &beta0, arma::field<arma::urowvec> &beta_idx, arma::field<arma::vec> para, arma::field<arma::uvec> free_idx);

/* @syntax [sw, wt, uw, w, rs] = lParseOpt(arg, argpt, sz, idx)
 * @param arg cell vector of arguments
 * @param argpt pointer to next argument cell
 * @param sz size of y data vector
 * @param idx index of bad data points, to be stripped from weights
 * @return sw show waitbar flag
 * @return wt waitbar text
 * @return uw use weights flag
 * @return w vector of weights
 * @return rs restricted fitting flag */
void lParseOpt(int &sw, int &uw, arma::mat &w, int &rs,
               QStringList Opt, arma::mat weights, int sz, arma::umat idx);

/* @syntax para = lBetaToPara(beta, parain, free_idx, beta_idx, NPara)
 * @param beta vector of values of free parameters
 * @param parain cell vector of parameter vectors
 * @param free_idx cell vector of indizes of free components in para
 * @param beta_idx cell vector of indizes of free components in beta
 * @param NPara length of para */
arma::field<arma::vec> lBetaToPara(arma::vec beta, arma::field<arma::vec> parain,
           arma::field<arma::uvec> free_idx, arma::field<arma::urowvec> beta_idx, int NPara);

/* @syntax [s, p1, ..., pn, X2, r, J] = fitFunction(x,y,model,ip1, f1, ...
 *  ipn, fn, Opt1, [OptVal1], ...)
 * @param x data points x values (independent variable)
 * @param y data points y values (dependent variable)
 * @param model name of function implementing the data model
 * @param ip_i initial values for i'th parameter of model (must not be string)
 * @param f_i index of parameters in i parameter vector p_n that shall be varied
 * @param Opt_i string describing an option that possibly takes an argument
 * @param OptVal_i value for i'th option, options are:
 *  'Waitbar', value: title string for waitbar (else don't show waitbar)
 *  'Weights', value: vector of weiths, else no weighted fitting
 *  'Restrict' (no value) restrict free parameter values to positive values
 * @return s status of success
 * @return p_i value of the i'th parameter set after fitting
 * @return X2 mean square relative residual
 * @return J jacobian with respect to the free parameters */
bool fitFunction(arma::field<arma::vec> &para, double &X2, arma::field<arma::uvec> free_idx, arma::vec x, arma::vec y, QString model,
                 QStringList Opt = QStringList(), QString waittext = "", arma::mat weights = arma::mat(0, 0));

arma::vec FitModel(QString model, arma::vec x, arma::field<arma::vec> para);

/* call to peakBaseSpline with step ratios instead of step heights
 *
 * Wraps up a call to peakBaseSpline. Parameters are exactly the same, except
 * that instead of peak step heights the peak step ratios are given. These
 * are transformed to Step Heights according to
 * Step Height = Step Ratio * Net Area
 *
 * Allows to automatically update the step height in net area fitting.
 *
 * @syntax y=wrapStepRatio(x,cx,cy,cdy,na,ct,fc,sr,lt,lta,ut,uta,db,rb,rd)
 * @param x channels
 * @param cx baseline control points x data
 * @param cy baseline control points y data
 * @param cdy slope at control points
 * @param na peak net areas
 * @param ct peak centroids
 * @param fc peak fwhm in channels
 * @param sr peak step ratios
 * @param lt peak lower tailing parameters
 * @param lta peak lower tailing alphas
 * @param ut peak upper tailing parameters
 * @param uta peak upper tailing alphas
 * @param db peak doppler broadening in channels
 * @param rb recoil peak model beta
 * @param rd recoil peak model delta
 * @return y peak and baseline counts in channels x */
arma::mat wrapStepRatio(arma::vec x, arma::vec cx, arma::vec cy, arma::vec cdy,
                        arma::vec na, arma::vec ct, arma::vec fc, arma::vec sr, arma::vec lt,
                        arma::vec lta, arma::vec ut, arma::vec uta, arma::vec db, arma::vec rb, arma::vec rd);

/* peakBaseLin - peak with linear baseline
 *
 *   cnt = peakBaseLin(CHAN, BL, NA, C, F) returns the counts in channels CHAN
 *      for peaks with net areas NA, centroids C and channel resolutions F
 *      above a linear baseline:
 *
 *                       n          exp(-(i-C(j))^2/(2*sigma(j)))
 *      cnt(i)=a+b*i+ sum    NA(j)* --------------------------
 *                      j=1             sqrt(2 pi sigma(j)^2)
 *
 *      where BL = [b, a];
 *      and sigma(j)=F(j)/sqrt(8log(2));
 *
 * @syntax cnt = peakBaseLin(chan, BL, NA, C, F)
 * @param chan channels at which peak plus baseline shall be evaluated
 * @param BL baseline parameters [b, a]
 * @param NA peak net areas
 * @param C peak centroids
 * @param F peak fwhms in channel units
 * @author Andreas Pelikan */
arma::mat peakBaseLin(arma::mat chan, arma::vec BL, arma::vec NA, arma::vec C, arma::vec F);

/* @syntax bc = abkcnt(BL,CT,FC)
 * @syntax bc = abkcnt(BL,CT,FC,k)
 * @param BL baseline counts in every channel
 * @param CT peak centroids
 * @param FC peak fwhms in channel units
 * @param k background count width (default 1.25)
 * @return bc total background counts */
arma::mat abkcnt(arma::rowvec BL, arma::vec CT, arma::vec FC, double k = 1.25);

/* @syntax s = vecsum(y, l, h)
 * @param y data vector
 * @param l lower integration border (not necessarily integer)
 * @param h upper integration border (not necessarily integer)
 * @return s integral from l to h of the histogramm of y; size is the common size of l and h */
arma::mat vecsum(arma::rowvec y, arma::vec l, arma::vec h);

/* Create interpolation parameters for a data set
 * Transform a dataset to a aatami calibration parameter vector and it's error
 * @syntax [s, msg, p, perr] = lMakeInterp(sn, cal, data) */
bool lMakeInterp(arma::mat &p, arma::mat &perr, arma::vec x, arma::vec y, arma::vec err);

/* =====================================================刻度拟合======================================================== */

/* @syntax [s, msg, p, perr] = calFitPara(cal, x, y, err)
 * @param cal calibration name
 * @param x independent variable data
 * @param y dependent variable data
 * @param err data error
 * @return s status of fitting
 * @return msg error message
 * @return p parameter
 * @return perr parameter error */
//bool calFitPara(arma::vec &p, arma::vec &perr, CalibrationType cal, arma::mat calData);
bool calFitPara(arma::vec &p, arma::vec &perr, CalibrationType cal, arma::vec x, arma::vec y, arma::vec err = arma::vec(),
                arma::vec p0 = arma::vec(), arma::vec perr0 = arma::vec(), arma::urowvec idx = arma::urowvec());

/* @syntax [s, msg, pIni] = lGetIniPara(cal, x, y)
 * @param cal calibration string
 * @param x independend variable data
 * @param y dependent variable data
 * @return s status 0: failed, 1: succeeded
 * @return msg error message
 * @return pIni initial parameter vector */
bool lGetIniPara(arma::colvec &pIni, CalibrationType cal, arma::colvec x, arma::colvec y);

/* @syntax [issp,s,msg,p,perr]=lHandleSpecial(cal,x,y,err,p0,pErr0,idx,np)
 * @param cal calibration string
 * @param x independent variable data
 * @param y dependent variable data
 * @param err data error
 * @param p0 initial parameter estimate
 * @param pErr0 initial parameter error estimate
 * @param idx index of paramters to be fitted
 * @param np flag if initial parameters are user-specified
 * @return issp flag 0: not a special case, 1: special case
 * @return s status 0: failed, 1: succeeded
 * @return msg error message
 * @return p parameter vector
 * @return perr parameter error vector */
bool lHandleSpecial(bool &issp, arma::colvec &p, arma::colvec &perr, CalibrationType cal, arma::colvec x, arma::colvec y,
     arma::colvec err, arma::colvec p0, arma::colvec pErr0, arma::urowvec idx, bool np = true);

/* @syntax [s, msg, p, perr, x2] = lFitStandard(x, y, err, p0, pErr0, idx)
 * @param x s data points
 * @param y y data points
 * @param err data error estimates
 * @param p0 initial parameter estimate (first entry function type)
 * @param pErr0 initial parameter error estimate
 * @param idx index of free parameters (starting with 2) */
bool lFitStandard(arma::colvec &p, arma::colvec &perr, double &x2, arma::colvec x, arma::colvec y,
                  arma::colvec err, arma::colvec p0, arma::colvec pErr0, arma::urowvec idx);

/* @syntax [perr, RMS, X2, cv] = calEstimateError(p, x, y, err)
 * @syntax [perr, RMS, X2, cv] = calEstimateError(p, x, y, err, idx)
 * @param p calibration parameter coefficients, including function type code
 * @param x data points x data
 * @param y data points y data
 * @param err data points error (not yet taken into account)
 * @param idx optional, if specified, only the indexed parameters (starting with
 * 1) will be taken into account, others assumed to be fixed (limited model).
 * @return perr parameter coefficients error estimate
 * @return RMS mean square relative error
 * @return X2 chi-square value
 * @return cv covariance matrix of the jacobian */
arma::colvec calEstimateError(arma::colvec p, arma::colvec x, arma::colvec y, arma::colvec err,
     arma::urowvec idx = arma::urowvec()/*, arma::mat &RMS, double &X2, arma::mat &cv*/);

/* @syntax [s, msg, p] = calDefIniPara(ft, x, y, e)
 * @param ft function type identifier (number of equation)
 * @param x data points x values
 * @param y data points y values
 * @param e data points errors
 * @return p initial parameters (para = [ft, p])*/
arma::vec calDefIniPara(int ft, arma::vec x, arma::vec y);

/* fit aatami style polynomial to data points
 * fits a polynomial in aatami order [a0, a1, ...] to data points xi, yi
 * @syntax p = lFitPoly(x, y)
 * @param x x data points
 * @param y y data points
 * @return p polynomial coefficients*/
arma::vec lFitPoly(arma::vec x, arma::vec y);

}
#endif // INDEPENDENTALG_H