errorEstimate.cpp

Go to the documentation of this file.

#include "atomprobe/atomprobe.h"

using namespace std;
using namespace AtomProbe;


//This uses poisson ratio error to determine error bars on a line profile from IVAS
int main(int argc, char *argv[])
{
    if(argc !=2)
    {
        cerr << "Usage: " << argv[0] << " CSV_FILE" << endl;
        return 1;
    }


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

    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;
    }

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


    vector<vector<float > > counts;

    counts.resize(intensityData.size());
    for(auto &v : counts)
    {
        v.resize(intensityData[0].size());
    }

    unsigned int distanceCol=-1;
    unsigned int countCol=-1;
    vector<unsigned int> compCols;


    for(auto ui=0;ui<headings.size();ui++)
    {
        if(headings[ui].find("Distance") != std::string::npos)
        {
            distanceCol=ui;
            continue;
        }
        
        if(headings[ui].find("Atom Count") != std::string::npos)
        {
            countCol=ui;
            continue;
        }

        //Check for a column we care about
        if(headings[ui].find("%") != std::string::npos &&
                headings[ui].find("Sigma") == std::string::npos)
        {
            cerr << "Using :" << headings[ui] << endl;
            compCols.push_back(ui);
        }
    }

    if(distanceCol == -1 || countCol == -1)
    {
        cerr << "Unable to locate distance or atom count columns. Aborting" << endl;
        return 1;
    }

    
    //Perform selection from large table to make subtable

    vector<string> perColumnHeading;
    vector<vector<float> > perColumnCounts;
    for(auto v : compCols)
    {
        perColumnHeading.push_back(headings[v]);
        perColumnCounts.push_back(intensityData[v]);
    }


    vector<float> distances;
    string distHeader = headings[distanceCol];
    for(auto v : intensityData[distanceCol])
        distances.push_back(v);

    const float ALPHA=0.99;
    const unsigned int NTRIALS=1000;

    //Now extract the counts across the row, and feed into error routine
    vector<vector<pair<float,float> > > errBounds; 

    vector<vector<float> > perRowCounts=perColumnCounts;
    transposeVector(perRowCounts);

    for(auto vec: perRowCounts)
    {
        //Denominator
        float dCounts=0;
        for(auto v : vec)
        {
            if(!isnan(v))
                dCounts+=v;
        }

        
        vector<pair<float,float> > errs;
        if(!dCounts)
        {
            errs.resize(vec.size(),make_pair(-1,-1));
            errBounds.push_back(errs);
            errs.clear();
            continue;
        }

        float lBound,uBound;
        for(auto ui=0;ui<vec.size();ui++)
        {
            lBound=-1;uBound=-1;
            if(!isnan(vec[ui]))
            {
                if(!numericalEstimatePoissRatioConf(vec[ui],dCounts,
                        ALPHA,NTRIALS,randGen.getRng(),lBound,uBound))
                {
                    //could not determine confidence
                    lBound=0;uBound=1;
                }

                //Cap lBound,uBound
                if(uBound > 1)
                    uBound=1;

            }
            errs.push_back(make_pair(lBound,uBound));
        }
        cerr <<" Error bound size:" << errs.size() << endl; 
        errBounds.push_back(errs);
        errs.clear();
    }

    //Dump error estimates
    cout <<distHeader << "\t" ;  
    for(auto s : perColumnHeading)
        cout << s << "(low)\t" << s << "(hi)" << "\t";
    cout << endl;

    for(auto ui=0u;ui<errBounds.size();ui++)
    {
        cout << distances[ui] << "\t";
        for(auto e : errBounds[ui])
        {
            cout << e.first << "\t" << e.second << "\t";
        }
        cout << endl;
    }

}