2068 lines
60 KiB
C++
2068 lines
60 KiB
C++
#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> ¶, 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;
|
|
}
|
|
|
|
}
|