fitting.cpp

Go to the documentation of this file.

#include <vector>
#include <cmath>
#include <algorithm>

#include "atomprobe/helper/maths/misc.h"
#include "atomprobe/spectrum/fitting.h"

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf.h>

using std::vector;

//DEBUG ONLY
#include <fstream>
using std::cerr;
using std::endl;
using std::ofstream;

namespace AtomProbe
{
struct DATA_FIT
{
    const vector<double> *x,*y;
};

//Return least-squares values (y-y')^2
double lsq(const vector<double>  &y, const vector<double> &yFit)
{
//  ASSERT(y.size() ==yFit.size());
    ofstream f("debugfit.txt");

    for(auto ui=0;ui<yFit.size();ui++)
        f << y[ui] << "\t" << yFit[ui] << endl;
    

    double dySqr=0;
    for(auto ui=0u;ui<y.size();ui++)
        dySqr+=(y[ui]-yFit[ui])*(y[ui]-yFit[ui]);
    return dySqr;
}

void voigtProfile(const vector<double> &x, double sigma, double gamma, vector<double> &y);

double fitVoigtFunc(const gsl_vector *v, void *params)
{
    double sigma=v->data[0];
    double gamma=v->data[1];
    double mu=v->data[2];
    double amp=v->data[3];

    DATA_FIT *vParams=(DATA_FIT*)params;

    vector<double> yFit;
    voigtProfile(*(vParams->x),sigma,gamma,mu,amp,yFit);

    const vector<double> &yV = *(vParams->y);
    return lsq(yV,yFit);
}

double fitDoniachSunjic(const gsl_vector *v, void *params)
{
    double a=v->data[0];
    double f=v->data[1];
    double mu=v->data[2];
    double amp=v->data[3];

    vector<double> xFit;
    DATA_FIT *vParams=(DATA_FIT*)params;

    vector<double> yFit;
    doniachSunjic(*(vParams->x),a,f,mu,amp,yFit);

    //Least-squares
    const vector<double> &yV = *(vParams->y);
    return lsq(yV,yFit);
}

//Arxiv : 0711.4449, Equations (1-4)
// A "smoothed" log-gaussian like profile
double fitLikeLogGaussian( const gsl_vector *v, void *params)
{
    //FIXME :Is amp & h redudant??
    //xp : x-positioning
    double xp =v->data[0];
    double sigma=v->data[1];
    double lambda=v->data[2];
    double amp=v->data[3];
    double h=v->data[4];
    
    vector<double> xFit;
    DATA_FIT *vParams=(DATA_FIT*)params;

    vector<double> yFit;
    likeLogGaussian(*(vParams->x),xp,sigma, lambda,amp,h,yFit);

    const vector<double> &yV = *(vParams->y);
    double lsqV = lsq(yV,yFit);
    //cerr << "LSQ :" << lsqV << endl;
    return lsqV;
}

//Exponential-modified normal
double fitExpNorm(const gsl_vector *v, void *params)
{
    double K=v->data[0];
    double mu=v->data[1];
    double amp=v->data[2];
    double sigma=v->data[3];

    vector<double> xFit;
    DATA_FIT *vParams=(DATA_FIT*)params;

    vector<double> yFit;
    expNorm(*(vParams->x),K,mu,sigma, amp,yFit);
    const vector<double> &yV = *(vParams->y);

    double lsqV = lsq(yV,yFit);
//  cerr << "lsq :" << lsqV << endl; 
    return lsqV; 
}

//Values should be X-positions and y- counts
bool fitVoigt(const vector<double> &x, const vector<double> &y, double &sigma, double &gamma, double &mu, double &amp, bool autoInit)
{
    // == Initial guess == 
    //find the mean of the X values
    if(autoInit)
    {
        mu= weightedMean(x,y);
        sigma= (*std::max_element(x.begin(),x.end()) - *std::min_element(x.begin(),x.end()))/4;
        gamma=1.0;
        //Find the maximum y value
        // FIXME: More stable to use area under curve, rather than max-y
        const double VOIGT_NORM_MAX=0.3; //The maximum for the voigt function is parameter dependent.
        amp = *(std::max_element(y.begin(),y.end()))*VOIGT_NORM_MAX;
    
    }
    //===================== 

    
    gsl_vector *xZero = gsl_vector_alloc (4);   
    gsl_vector *stepSize= gsl_vector_alloc (4);
    
    xZero->data[0] = sigma;
    xZero->data[1] = gamma;
    xZero->data[2] = mu;
    xZero->data[3] = amp;

    //set step sizes
    stepSize->data[0] = sigma/10;
    stepSize->data[1] = gamma/10;
    stepSize->data[2] = mu/2;
    stepSize->data[3] = amp;


    //Use nelder mead (random initial direction) to find solution 
    const gsl_multimin_fminimizer_type *t=  gsl_multimin_fminimizer_nmsimplex2rand;
    gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc (t, 4); 

    DATA_FIT vParams;
    vParams.x = &x;
    vParams.y = &y;

    gsl_multimin_function minFunc;
    minFunc.n=4;
    minFunc.f=fitVoigtFunc;
    minFunc.params=(void*)(&vParams);

    gsl_multimin_fminimizer_set (s, &minFunc, xZero, stepSize);

    const unsigned int MAX_NM_STEPS=5000;
    unsigned int it=0,status;
    do
    {
        it++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status)
            break;

        auto size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-6);

        sigma=s->x->data[0];
        gamma=s->x->data[1];
        mu=s->x->data[2];
        amp=s->x->data[3];
        if (status == GSL_SUCCESS)
        {
            gsl_vector_free(stepSize);
            gsl_vector_free(xZero);
            gsl_multimin_fminimizer_free(s);
            return true;
        }
    } while( it< MAX_NM_STEPS && status == GSL_CONTINUE);

    gsl_vector_free(stepSize);
    gsl_vector_free(xZero);
    gsl_multimin_fminimizer_free(s);
    //Failed to successfully minimise
    return false;
}

bool fitDoniachSunjic(const vector<double> &x, const vector<double> &y, double &a, double &f, double &mu, double &amp, bool autoInit)
{
    // == Initial guess == 
    //find the mean of the X values
    if(autoInit)
    {
        mu= weightedMean(x,y);
        a= (*std::max_element(x.begin(),x.end()) - *std::min_element(x.begin(),x.end()))/12;
        f=1.0;
        //Find the maximum y value
        amp = *(std::max_element(y.begin(),y.end()))*50;
    }
    //===================== 

    
    gsl_vector *xZero = gsl_vector_alloc (4);   
    gsl_vector *stepSize= gsl_vector_alloc (4);
    
    xZero->data[0] = a;
    xZero->data[1] = f;
    xZero->data[2] = mu;
    xZero->data[3] = amp;

    //set step sizes
    gsl_vector_set_all(stepSize,1.0);

    //Use nelder mead (random initial direction) to find solution 
    const gsl_multimin_fminimizer_type *t=  gsl_multimin_fminimizer_nmsimplex2rand;
    gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc (t, 4); 

    DATA_FIT vParams;
    vParams.x = &x;
    vParams.y = &y;

    gsl_multimin_function minFunc;
    minFunc.n=4;
    minFunc.f=fitDoniachSunjic;
    minFunc.params=(void*)(&vParams);

    gsl_multimin_fminimizer_set (s, &minFunc, xZero, stepSize);

    const unsigned int MAX_NM_STEPS=5000;
    unsigned int it=0,status;
    do
    {
        it++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status)
            break;

        auto size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-6);

        a=s->x->data[0];
        f=s->x->data[1];
        mu=s->x->data[2];
        amp=s->x->data[3];
        //--
        if (status == GSL_SUCCESS)
        {
            gsl_vector_free(stepSize);
            gsl_vector_free(xZero);
            gsl_multimin_fminimizer_free(s);
            return true;
        }
    } while( it< MAX_NM_STEPS && status == GSL_CONTINUE);

    gsl_vector_free(stepSize);
    gsl_vector_free(xZero);
    gsl_multimin_fminimizer_free(s);
    //Failed to successfully minimise
    return false;
}


bool fitExpNorm(const vector<double> &x, const vector<double> &y,
        double &K, double &mu,double &sigma,double &amp,bool autoInit)
{
    if(autoInit)
    {
        //These are pretty heuristic. Weighted mean initialises X position
        // K is right-tailed, so >1. amp is adjusted for K against some sample peaks
        // and sigma is a wild guess, based on some examples
        mu=weightedMean(x,y);
        K=1.5;
        //Find the maximum y value
        amp = *(std::max_element(y.begin(),y.end()))*2.0;
        sigma=0.05;
    }
    
    const unsigned int DIMENSION=4;
    gsl_vector *xZero = gsl_vector_alloc (DIMENSION);   
    gsl_vector *stepSize= gsl_vector_alloc (DIMENSION);
    
    xZero->data[0] = K;
    xZero->data[1] = mu;
    xZero->data[2] = amp;
    xZero->data[3] = sigma;

    //set step sizes
    stepSize->data[0] = K;
    stepSize->data[1] = mu;
    stepSize->data[3] = 0.1;
    stepSize->data[2] = amp*10;
    
    //Use nelder mead (random initial direction) to find solution 
    const gsl_multimin_fminimizer_type *t=  gsl_multimin_fminimizer_nmsimplex2rand;
    gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc (t, DIMENSION); 

    DATA_FIT vParams;
    vParams.x = &x;
    vParams.y = &y;

    gsl_multimin_function minFunc;
    minFunc.n=DIMENSION;
    minFunc.f=fitExpNorm;
    minFunc.params=(void*)(&vParams);

    gsl_multimin_fminimizer_set (s, &minFunc, xZero, stepSize);

    const unsigned int MAX_NM_STEPS=100;
    unsigned int it=0,status;
    do
    {
        it++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status)
            break;

        auto size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-8);

        K=s->x->data[0];
        mu=s->x->data[1];
        amp=s->x->data[2];
        sigma=s->x->data[3];

//      cerr << "K:" << K << " mu:" << mu << " amp :" << amp << endl;
        if (status == GSL_SUCCESS)
        {
            gsl_vector_free(stepSize);
            gsl_vector_free(xZero);
            gsl_multimin_fminimizer_free(s);
            return true;
        }
    } while( it< MAX_NM_STEPS && status == GSL_CONTINUE);

    gsl_vector_free(stepSize);
    gsl_vector_free(xZero);
    gsl_multimin_fminimizer_free(s);
    //Failed to successfully minimise
    return false;
}

bool fitLikeLogGaussian(const vector<double> &x, const vector<double> &y, double &lambda, double &sigma, double &xp, double &amp, double &h, bool autoInit)
{
    // == Initial guess == 
    //find the mean of the X values
    if(autoInit)
    {
        xp= weightedMean(x,y);
        lambda= (*std::max_element(x.begin(),x.end()) - *std::min_element(x.begin(),x.end()))/12;
        sigma=lambda*2;
        //Find the maximum y value
        amp = *(std::max_element(y.begin(),y.end()));
    }
    //===================== 

    
    gsl_vector *xZero = gsl_vector_alloc (5);   
    gsl_vector *stepSize= gsl_vector_alloc (5);

    //Note that we re-arrange this here
    // so that the external interface to this
    // function is consistent with others.
    // but that the function parameter orders
    // is consistent with the paper
    xZero->data[0] = xp;
    xZero->data[1] = sigma;
    xZero->data[2] = lambda;
    xZero->data[3] = amp;
    xZero->data[4] = h;

    //set step sizes

    stepSize->data[0]=xp/10; //xp
    stepSize->data[1]=sigma/2.0; //sigma
    stepSize->data[2]=0.6;//lambda
    stepSize->data[3]=amp/2;//amp
    stepSize->data[4]=.5;//h

    //Use nelder mead (random initial direction) to find solution 
    const gsl_multimin_fminimizer_type *t=  gsl_multimin_fminimizer_nmsimplex2rand;
    gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc (t, 5); 

    DATA_FIT vParams;
    vParams.x = &x;
    vParams.y = &y;

    gsl_multimin_function minFunc;
    minFunc.n=5;
    minFunc.f=fitLikeLogGaussian;
    minFunc.params=(void*)(&vParams);

    gsl_multimin_fminimizer_set (s, &minFunc, xZero, stepSize);

    const unsigned int MAX_NM_STEPS=500;
    unsigned int it=0,status;
    do
    {
        it++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status)
            break;

        auto size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-6);

        xp=s->x->data[0];
        sigma=s->x->data[1];
        lambda=s->x->data[2];
        amp=s->x->data[3];
        h = s->x->data[4];

//      cerr << "Fit : xp" <<xp << " sigma:" << sigma << " lambda:" << lambda<<
//          " amp:"  << amp << " h :" << h << endl;
        //--
        if (status == GSL_SUCCESS)
        {
            gsl_vector_free(stepSize);
            gsl_vector_free(xZero);
            gsl_multimin_fminimizer_free(s);
            return true;
        }
    } while( it< MAX_NM_STEPS && status == GSL_CONTINUE);

    gsl_vector_free(stepSize);
    gsl_vector_free(xZero);
    gsl_multimin_fminimizer_free(s);
    //Failed to successfully minimise
    return false;
}
void voigtProfile(const vector<double> &x, double sigma, double gamma, vector<double> &y)
{
    y.resize(x.size());

    //Simple approximation. I think arxiv:1901.08366v1
    // has more detailed versions
    // My old code has    L =  1/pi() * ( 0.5*gammaV./( (x-x0).^2 + 0.5*gammaV^2));

    //Full-width-half-maximum gaussian
    const double fwhmG= 2*sqrt(log(2))*sigma;
    //Full-width-half-maximum lorentizan
    const double fwhmL= 2*gamma;
    //Equation 10,Ida et al, Extended pseudo-Voigt function for approximating the Voigt profile
    const double tch = pow(fwhmG*fwhmG*fwhmG*fwhmG*fwhmG + 
            2.69269*fwhmG*fwhmG*fwhmG*fwhmG*fwhmL+ 
            2.42843*fwhmG*fwhmG*fwhmG*fwhmL*fwhmL +
            4.47163*fwhmG*fwhmG*fwhmL*fwhmL*fwhmL +
            0.07842*fwhmG*fwhmL*fwhmL*fwhmL*fwhmL +
            fwhmL*fwhmL*fwhmL*fwhmL*fwhmL,0.2);

    const double mix = 1.36603*fwhmL/tch -0.47719*(fwhmL/tch)*(fwhmL/tch) +
            0.11116*(fwhmL/tch)*(fwhmL/tch)*(fwhmL/tch) ;
    for(auto ui=0;ui<x.size();ui++)
    {
        double L,G;
        //Lorentizan (Cauchy)
        L  = 1.0/(M_PI*gamma)/( 1.0 + x[ui]*x[ui]/(gamma*gamma));
        //Gaussian
        G  = 1.0/(sigma*sqrt(2.0*M_PI))*exp(-(x[ui]*x[ui])/(2.0*sigma*sigma));

        y[ui]=(mix)*L+(1-mix)*G;
    }

    
}

//Convenience wrapper with additional parameters
void voigtProfile(const vector<double> &x, double sigma, 
    double gamma, double mu, double amp, vector<double> &y)
{
    vector<double> shiftX;
    shiftX =x;
    for( auto &v : shiftX)
        v=v-mu;

    //Generate prfoile
    voigtProfile(shiftX,sigma,gamma,y);

    //Scale y
    for(auto &v : y)
        v*=amp;
}

void doniachSunjic(const std::vector<double> &x, double &a, double &F, double &mu,double &amp,
                std::vector<double> &y)
{
    y.resize(x.size());
    for(auto ui=0;ui<x.size();ui++)
        y[ui] = amp*( M_PI*a/2 + (1-a)*atan( (x[ui]-mu)/F))/(F + (x[ui]-mu)*(x[ui]-mu));
}

void likeLogGaussian(const std::vector<double> &x, double &xp, double &sigma, double &lambda,double &amp, double &h,
                std::vector<double> &y)
{
    double z = sqrt(h*abs(lambda) + sqrt(h*h*lambda*lambda + 1));
    double A = sqrt(2.0/M_PI) * abs(lambda)/( (z+ 1.0/z)*sigma)*exp(-sigma*sigma/2.0);



    y.resize(x.size());
    for(auto ui=0;ui<x.size();ui++)
    {
        //Log contents
        double lgvc;
        lgvc=1 + 2.0*z*lambda*(x[ui]-xp)/(z*z +1);

        if(lgvc >0)
            y[ui] = amp*A*exp(-1.0/(2.0*sigma*sigma)*log(lgvc));
        else
            y[ui]=0;
    }
}

void expNorm(const std::vector<double> &x, double &K, double &mu, double &sigma,double &amp, vector<double> &y)
{
    y.resize(x.size());
    for(auto ui=0;ui<x.size();ui++)
    {
        double cInp;
        cInp=  -((x[ui]-mu)/sigma-1.0/K)/sqrt(2.0);
        y[ui] = amp/(2.0*K)*exp(1.0/(2.0*K*K)-(x[ui]-mu)/(K*sigma))*gsl_sf_erfc(cInp);
    }
}

//Find the partial derivatives of f= (c-expNom(x,k,mu,sigma,amp))^2
// f = (c - ((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ) - (β e^(1/(2 K^2)))/(2 K))^2
//Partial/partial(m)  = 
// df/dm =2 (erfc((m - x)/σ - 1/K)/(sqrt(2) K σ) + (sqrt(2/π) (x - m) e^(-((m - x)/σ - 1/K)^2))/(K σ^2)) (c - ((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ) - (β e^(1/(2 K^2)))/(2 K))
// df/dK = 2 (((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K^2 σ) + (sqrt(2/π) (x - m) e^(-((m - x)/σ - 1/K)^2))/(K^3 σ) + (β e^(1/(2 K^2)))/(2 K^2) + (β e^(1/(2 K^2)))/(2 K^4)) (c - ((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ) - (β e^(1/(2 K^2)))/(2 K))
// df/dσ=  2 (((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ^2) - (sqrt(2/π) (m - x) (x - m) e^(-((m - x)/σ - 1/K)^2))/(K σ^3)) (c - ((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ) - (β e^(1/(2 K^2)))/(2 K))
// df/da = -(e^(1/(2 K^2)) (-(a e^(1/(2 K^2)))/(2 K) + c - ((x - m) erfc((m - x)/σ - 1/K))/(sqrt(2) K σ)))/K

void expNormLSQDerivs(const gsl_vector *v, void *params, gsl_vector *df)
{

    DATA_FIT *vParams=(DATA_FIT*)params;
    const vector<double> &y = *(vParams->y);
    const vector<double> &x = *(vParams->x);
    
    float K=v->data[0];
    float m=v->data[1];
    float amp=v->data[2];
    float sigma=v->data[3];

    double dfdK=0; // K
    double dfdm=0; //Mu
    double dfda=0;// amp
    double dfds=0; //sigma

    for(auto ui=0;ui<y.size();ui++)
    {
        float c;
        c = y[ui];
        dfdm+=2.0*(gsl_sf_erfc((m - x[ui])/sigma - 1.0/K)/(sqrt(2.0)*K*sigma) + 
                (sqrt(2.0/M_PI)*(x[ui] - m)*exp(-pow(((m - x[ui])/sigma - 1.0/K),2.0)))/(K*sigma*sigma)) *
                (c - ((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*sigma) - (amp*exp(1.0/(2.0*K*K)))/(2.0*K));

        dfdK+=2.0*(((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*K*sigma) + 
                (sqrt(2.0/M_PI)*(x[ui] - m)*exp(-pow(((m - x[ui])/sigma - 1.0/K),2.0))/(K*K*K*sigma)) + 
                (amp*exp(1.0/(2.0*K*K)))/(2.0*K*K) + (amp*exp(1.0/(2.0*K*K)))/(2.0*K*K*K*K))*
                (c - ((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*sigma) - (amp*exp(1.0/(2.0*K*K)))/(2.0*K));
        dfds+=2.0*(((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*sigma*sigma) - 
                (sqrt(2.0/M_PI)*(m - x[ui])*(x[ui] - m)*exp(-pow(((m - x[ui])/sigma - 1.0/K),2.0)))/(K*sigma*sigma*sigma))*
                (c - ((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*sigma) - (amp*exp(1.0/(2.0*K*K)))/(2.0*K));

        dfda+= -(exp(1.0/(2.0*K*K))*(-(amp*exp(1.0/(2.0*K*K)))/(2.0*K) + c - ((x[ui] - m)*gsl_sf_erfc((m - x[ui])/sigma - 1.0/K))/(sqrt(2.0)*K*sigma)))/K;
    }

    df->data[0]=dfdK;
    df->data[1]=dfdm;
    df->data[2]=dfda;
    df->data[3]=dfds;
}


}