curveFit.cpp

Go to the documentation of this file.

/* curveFit.cpp : Sample program to fit a simple mass spectrum 
 * Copyright (C) 2015  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */



#include <iostream>

#include "atomprobe/atomprobe.h"

using namespace std;
using namespace AtomProbe;

int main(int argc, char *argv[]) 
{
    if(argc !=  2)
    {
        cerr << "USAGE: " << argv[0] << " spectrum" << endl;

        cerr << " an example file, spectrum-example.txt, is distributed with libatomprobe" << endl;
        return 1;
    }

    //Load the file 
    std::string filename = argv[1];
    cerr << "Processing file :" << filename << endl;

    vector<vector<float> > data;
    vector<string> headers;
    if(loadTextData(filename.c_str(),data,headers))
    {
        cerr << "failed to load data" << endl;
        return 1;
    }

    if(data.empty())
    {
        cerr << "0 rows found. Aborting" << endl;
        return 2;
    }

    transposeVector(data);

    if(data[0].size() !=2)
    {
        cerr << "Expected nx2 data size, got nx" << data[0].size() << endl;
        return 3;
    }

    cerr << "Loaded data : " << data.size() << " rows" << endl;

    //Flip to column format (?)
    transposeVector(data);
    
    

    //Square-root all  x-values, which moves the data from mass to TOF space
    for(auto &v : data[0])
        v=sqrt(v);

    vector<double> xData,yData;
    xData.resize(data[0].size());
    yData.resize(data[0].size());
    for(auto ui=0u;ui<data[0].size();ui++)
    {
        xData[ui] = data[0][ui];
        yData[ui] = data[1][ui];
    }
    //Fit the voigt function
    double sigma,gamma,mu,amp;
    if(!fitVoigt(xData,yData, sigma,gamma,mu,amp))
    {
        cerr << "Curve fit failed. " << endl;
    }


    cerr << "VOIGT " << endl;
    cerr << "=========" << endl;
    cerr << "sigma:\t" << sigma <<endl;
    cerr << "gamma:\t" << gamma <<endl;
    cerr << "mu:\t" << mu << " ( Sqrt : " << mu*mu<< ")" << endl;
    cerr << "amplitude:\t" << amp <<endl;
    cerr << "=========" << endl;
    


    vector<double> fittedVoigtProfile;
    voigtProfile(xData,sigma,gamma,mu,amp,fittedVoigtProfile);
    
    //Now feed this into provide a log-Gaussian like function
    double K;

    cerr << "Attempting semi-random fits of exponential-normal distribution." << endl;
    cerr << "This can be unstable if too far from initial conditions" << endl;
    cerr << " This is effectively a budget simulated anneal, and can fail. Try running a few times." << endl;

    vector<double> bestV = {K,sigma,mu,amp};

    const unsigned int N_ATTEMPTS=200;
    unsigned int iter=0;
    float bestLSQ=std::numeric_limits<float>::max();
    auto r = randGen.getRng();

    //FIXME: proper simulated anneal. This is a shitty random re-search method
    auto handler=gsl_set_error_handler_off();
    bool firstIt=true;
    do
    {
        //do our best to fit the curve. If first iteration, guess values heuristically
        fitExpNorm(xData,yData,K,mu,sigma,amp,firstIt);
        firstIt=false;

        //Compute residual^2
        double lsqThis;
        vector<double> fittedExpNormProfile;
        expNorm(xData,K,mu,sigma,amp,fittedExpNormProfile);
        lsqThis=lsq(yData,fittedExpNormProfile);
        if(lsqThis < bestLSQ)
        {
            bestV[0]=K;
            bestV[1]=sigma;
            bestV[2]=mu;
            bestV[3]=amp;
            bestLSQ=lsqThis;
            cerr << "Best: (K,Sigma,Mu,Amp) : (" << K << "," << sigma << "," << mu << "," << amp << ")" << endl;
        }

        const float SCALE=0.35;
        K = bestV[0] + gsl_ran_gaussian(r,bestV[0]*SCALE);
        sigma = bestV[1] + gsl_ran_gaussian(r,bestV[1]*SCALE);
        mu= bestV[2] + gsl_ran_gaussian(r,bestV[2]*SCALE);
        amp = bestV[3] + gsl_ran_gaussian(r,bestV[3]*SCALE);

        //Force positivity for amplitude & sigma
        amp=fabs(amp);
        sigma=fabs(sigma);

        
        iter++;
    }while (iter < N_ATTEMPTS);
    gsl_set_error_handler(handler);


    K = bestV[0] ;
    sigma = bestV[1] ;
    mu= bestV[2] ;
    amp = bestV[3] ;


    cerr << "ExpNorm" << endl;
    cerr << "=========" << endl;
    cerr << "K:\t" << K <<endl;
    cerr << "mu:\t" << mu << " ( Sqrt : " << mu*mu<< ")" << endl;
    cerr << "sigma:\t" << sigma <<endl;
    cerr << "amplitude:\t" << amp <<endl;
    cerr << "---------" << endl;
    cerr << "LSQ:" << bestLSQ << endl;
    cerr << "=========" << endl;

    vector<double> fittedExpNormProfile;
    expNorm(xData,K,mu,sigma,amp,fittedExpNormProfile);


    ofstream fOut("fittedProfile.txt");
    if(!fOut)
    {
        cerr << "Failed to open fittedProfile.txt for output" << endl;
        return 6;
    }

    fOut << "x\torig\tvoigt\tllg" << endl;
    for(auto ui=0u; ui<xData.size();ui++)
    {
        fOut << (xData[ui]) << "\t" << yData[ui] << "\t" << fittedVoigtProfile[ui] << "\t" << fittedExpNormProfile[ui]<<endl;
    }


    cerr << "Wrote curves to : fittedProfile.txt" << endl;

}