zechBackground.cpp

Go to the documentation of this file.

/* zechBackground.cpp : Sample program to perform a max-likelihood estimation on mass spectra
 * Copyright (C) 2019  D 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;

#ifndef ASSERT
    #include <cassert>
    #define ASSERT(f) assert(f)
#endif

template<class T>
bool dumpHistogramToFile(std::vector<std::vector<T> > &hist, const char *filename)
{
    ofstream f(filename);

    if(!f)
        return false;

    for(auto ui=0u;ui<hist.size();ui++)
    {
        for(auto uj=0u;uj<hist[ui].size();uj++)
            f << hist[ui][uj] << "\t";

        f << endl;
    }

    return true;
}


void zechCorrect(vector<float> &background, vector<float> &observation, float alpha,
            vector<float> &corrected)

{
    //Peform subtraction
    ASSERT(background.size() == observation.size());

    corrected=observation;
    const float THRESHOLD=550; // Counts
    for(auto i=0u;i<background.size();i++)
    {
        if(observation[i] > THRESHOLD || background[i] > THRESHOLD)
            corrected[i]=observation[i]-background[i];
        else if(!observation[i])
        {
            //pass : estimator doesn't work with obs=0
        }
        else
        {
            

            float estimate;
            if(AtomProbe::zechConfidenceLimits(background[i],observation[i],alpha,estimate))
                corrected[i]=estimate;
            else
                corrected[i]=observation[i]-background[i];

        }
    }
}

int main(int argc, char *argv[]) 
{
    //Check for correct program usage
    if(argc < 3)
    {
        //print error message, and exit
        cerr << "USAGE : " << argv[0] << " SPECTRUM_FILE [alpha1 alpha2 ....] OUTFILE" << endl;

        cerr << "In the spectrum file, the columns should be ordered as mass, count, background" << endl;
        return 1;
    }
            

    //Load CSV Data
    vector<string> headings;
    vector<vector<float > > intensityData;

    //M/c, Uncorrected, Corrected, Background Estimate
    auto errCode=loadTextData(argv[1],intensityData,headings,",\t");

    if(errCode)
    {
        cerr << "Error loading text file :" << errCode <<endl;
        return 1;
    }

    if(!intensityData.size())
    {
        cerr << "input data appears to be empty" << endl;
        return 1;
    
    }
    vector<float> alpha;
    alpha.resize(argc-3,-1);

    if(!alpha.size())
    {
        cerr << "Alpha not specified, using 0.5 (most probable corretion) " << endl;
        alpha.push_back(0.5);
    }
    else
    {
        unsigned int OFFSET=2;
        for(auto ui=0u;ui<alpha.size();ui++)
        {
            float tmp;
            if(stream_cast(tmp,argv[ui+OFFSET]))
            {
                cerr << "Unable to understand specified alpha :" << argv[ui+OFFSET] << endl;
                return 1;
            }

            alpha[ui]=tmp;
        }
    }

    


    transposeVector(intensityData);
    cerr << "Loaded :" << intensityData.size() << "Rows" << endl;


    vector<float> mass,observed,background;
    mass.resize(intensityData.size());
    observed.resize(intensityData.size());
    background.resize(intensityData.size());

    for(auto ui=0;ui<intensityData.size();ui++)
    {
        ASSERT(intensityData[ui].size() > 3);

        //mass[ui] = intensityData[ui][0];
        observed[ui] = intensityData[ui][1];
        background[ui] = intensityData[ui][3];

    }

    vector<vector<float> > corrected;
    corrected.resize(alpha.size());
    for(unsigned int ui=0;ui<alpha.size();ui++)
    {
        ASSERT(alpha[ui] > 0 && alpha[ui] < 1.0f);
        zechCorrect(background,observed,alpha[ui],corrected[ui]);
    }

    transposeVector(corrected);


    const char *outFile=argv[argc-1];
    if(!dumpHistogramToFile(corrected,outFile))
    {
        cerr << "Unable to write output to file" << endl;
        return 1;
    }

    cerr << "Complete. output written to  :" << outFile << endl;
}