specsim.cpp

Go to the documentation of this file.

#include "atomprobe/atomprobe.h"

#include <tuple>
#include <numeric>

using namespace std;
using namespace AtomProbe;

const float MASS_TOL=0.05;

typedef vector<pair<string,size_t> > FRAGMENT;


//Obtain the mass distributions for the given species fragments
// this does not take into account the  charge state nor the relative intensities of the charges
void getMassDistributions(const vector<FRAGMENT> &fragmentVec, const vector<string> &parentSpecies,
    const AbundanceData &natTable, vector<vector<pair<float,float> > > &massDistVec);

void normaliseVec(vector<float> &v)
{
    float normTotal=0;
    normTotal=std::accumulate(v.begin(),v.end(),0.0f);
    for(auto &frac : v)
        frac/=normTotal;
}

int main(int argc, char *argv[])
{

    //Obtain the natural abundance file
    //---
    //Use default unless overridden
    string sAbundance = "naturalAbundance.xml";
    if(argc >1)
        sAbundance =argv[1];

    AbundanceData natTable;
    
    if(natTable.open(sAbundance.c_str()))
    {
        cerr << "Failed to open abundance table : " << sAbundance << endl;
        return 1;
    }
    //---

    //Obtain the species from the user
    //---
    std::vector<string> species;

    std::string specStr;
    cerr << "Enter species names, space separated" << endl;
    getline(cin,specStr);

    //Split the string into a vector, around spaces
    splitStrsRef(specStr.c_str()," ",species);

    if(species.empty())
    {
        cerr << "no species specified, aborting" << endl;
        return 1;
    }
    
    
    std::string compStr;
    cerr << "Enter species atomic compositions, space separated. This will be normalised" << endl;
    getline(cin,compStr);

    vector<float> compositions;
    vector<string> compStrVec;
    //Split the string into a vector, around spaces
    splitStrsRef(compStr.c_str()," ",compStrVec);
    stripZeroEntries(compStrVec);

    for(auto &i : compStrVec)
    {
        float f;
        if(stream_cast(f,i))
        {
            cerr << "Unable to understand composition : \"" << i << "\", aborting";
            return 1;
        }

        if(f < 0)
        {
            cerr << "Composition cannot be negative" << endl;
            return 1;
        }
        compositions.push_back(f);
    }

    if(compositions.size() != species.size())
    {
        cerr << "Number of compositions did not match number of species, aborting" << endl;
        return 1;
    }

    for(auto frac : compositions)
    {
        if(frac < 0)
        {
            cerr << "Negative compositions not allowed. Aborting" << endl;
            return 1;
        }
    }

    //Normalise compositions
    normaliseVec(compositions); 
    //---

    
    //request the charge state
    //--
    vector<vector<int> >  chargeStates;
    cout << "Enter charge states, as numbers, e.g. 1 3"<<  endl;
    for(auto &sS: species)
    {
        string chargeStr;
        vector<string> chargesStrVec;

        cout << sS << "? ";
        getline(cin,chargeStr);


        splitStrsRef(chargeStr.c_str()," ",chargesStrVec);
        stripZeroEntries(chargesStrVec);

        vector<int> charges;
        charges.clear();

        for(auto &s : chargesStrVec)
        {
            int i;
            if(stream_cast(i,s))
            {
                cerr << "Unable to understand charge value, aborting" << endl;
                return 1;
            }
            if(i < 1)
            {
                cerr << "Charge states must be positive" << endl;
                return 1;
            }

            charges.push_back(i);
        }


        chargeStates.push_back(charges);

    }
    //--

    //Check we have some charge states
    //--
    unsigned int nCharges=0;
    for(auto &v : chargeStates)
        nCharges+=v.size();

    if(!nCharges)
    {
        cerr << "No charges specified, aborting" << endl;
        return 1;
    }

    //--


    //Request the charge state ratios from the user
    //----
    vector<vector<float> > chargeIntensities;
    bool haveMultiCharge=false;
    for( auto &v : chargeStates)
    {
        if(v.size()  >1)
        {
            haveMultiCharge=true;
            break;
        }

    }


    if(haveMultiCharge)
        cout << "Enter charge state intensities (0->1), will be normalised " << endl;
    for(auto i=0;i<chargeStates.size();i++)
    {
        vector<float> cR;
        cR.clear();

        if(chargeStates[i].size() > 1)
        {
            for(auto c : chargeStates[i])
            {
                float f;
                cout << species[i] << "(" <<c << ") ?";
                cin >> f;

                if(f < 0)
                {
                    cerr << "Charge state intensity must be greater than 0. Aborting" << endl;
                    return 1;
                }

                cR.push_back(f);
            }
                
            normaliseVec(cR);   
        }
        else
            cR.push_back(1);

        chargeIntensities.push_back(cR);
    }

    //----

    vector<FRAGMENT> fragmentVec;
    
    //Break each species, such as Fe2H2 into components {Fe,2},{H,2}    
    for(auto &s : species)
    {
        vector<pair<string,size_t> > frags;
        if(!RangeFile::decomposeIonNames(s,frags))
        {
            cerr << "WARN : Failed to decompose fragment : " << s << ", skipping" << endl;
        }

        fragmentVec.push_back(frags);
    }

    
    //Generate the mass distributions for each species
    vector<vector<pair<float,float> > > massDistVec;
    getMassDistributions(fragmentVec,species,natTable,massDistVec);

    //ASSERT(compositions.size() == massDistVec.size());


    cerr << "Output is (mass, relative intensity, charge) " << endl;
    for(auto ui=0;ui<massDistVec.size();ui++)
    {
        cerr << "--- " <<  species[ui] << "---" << endl;
        vector<tuple<float,float,int> > speciesPeaks;
        for(auto j=0;j<chargeStates[ui].size();j++)
        {
            int charge = chargeStates[ui][j];

            for(auto &m : massDistVec[ui])
            {
                float chargeIntensity = chargeIntensities[ui][j];
                speciesPeaks.push_back(
                std::make_tuple(m.first/charge,m.second*compositions[ui]*chargeIntensity, charge)); 
            }
        }

        //merge the peaks that are within tolerance
        //FIXME: This algorithm doesn't work if multiple peaks have 
        // transitivity (ie a links to b links to c)
        // its also a massive hack.
        //=========
        //Find any species that links to this one
        vector<pair<int,int> > linked;
        for(auto i=0;i<speciesPeaks.size();i++)
        {
        
            //Skip values already linked    
            bool preLinked;
            preLinked=false;
            for(auto p : linked)
            {
                if(p.second == i)
                {
                    preLinked=true;
                    break;
                }
            }

            if(preLinked)
                continue;
            
            for(auto j=i+1;j<speciesPeaks.size();j++)
            {
                float mass1,mass2;
                mass1=std::get<0>(speciesPeaks[j]); 
                mass2=std::get<0>(speciesPeaks[i]);
                int charge1,charge2;
                charge1=std::get<2>(speciesPeaks[j]);
                charge2=std::get<2>(speciesPeaks[i]);
                

                if(fabs( mass1/charge1- mass2/charge2) < MASS_TOL)
                    linked.push_back(make_pair(i,j));
            }


        }

        //Merge any peaks that are linked
        vector<std::tuple<float,float,vector<int> > > mergedPeaks;
        for(auto i : linked)
        {
            tuple<float,float,int> t1,t2;
            t1 = speciesPeaks[i.first];
            t2 = speciesPeaks[i.second];

            std::vector<int> charges;
            charges.push_back(std::get<2>(t1));
            charges.push_back(std::get<2>(t2));


            mergedPeaks.push_back(
                std::make_tuple((std::get<0>(t1) + std::get<0>(t2))*0.5,
                    std::get<1>(t1)+std::get<1>(t2),charges));

        }

        vector<bool> unlinked;
        unlinked.resize(speciesPeaks.size(),true);
        for(auto i : linked)
            unlinked[i.first]=unlinked[i.second]=false;

        for(auto i=0;i<speciesPeaks.size();i++)
        {
            if(unlinked[i])
            {
                auto t = speciesPeaks[i];
                vector<int> v;
                v.push_back(std::get<2>(t));
                mergedPeaks.push_back(
                    make_tuple(std::get<0>(t),std::get<1>(t),v));
            }
        }
        //=========

        for(auto &t : mergedPeaks)
        {
            string sCharges;
            sCharges.clear();

            for(auto charge : std::get<2>(t))
            {
                string tmp;
                stream_cast(tmp,charge);
                sCharges+=tmp + " ";
            }

             cerr << std::get<0>(t) << "," << std::get<1>(t) << 
                "," << sCharges << endl;
        }

        cerr << "-------------------------------" << endl << endl;
    }


    
}



void getMassDistributions(const vector<FRAGMENT> &fragmentVec, const vector<string> &parentSpecies,
    const AbundanceData &natTable, vector<vector<pair<float,float> > > &massDistVec)
{
    //ASSERT(fragmentVec.size() == parentSpecies.size());

    for(auto ui=0;ui<fragmentVec.size();ui++)
    {
        //Look up this fragment set (e.g. {{Fe,2},{H,2}}) in the 
        // abundance table
        vector<pair<size_t,unsigned int> > outputVec;
        natTable.getSymbolIndices(fragmentVec[ui],outputVec);

        //Check it was OK, and record as needed
        bool fragmentLookupOK;
        fragmentLookupOK=true;
        vector<size_t> indexVec,freqVec;
        for( auto &p : outputVec)
        {
            //Spit out error if lookup failed
            if(p.first == (size_t)-1)
            {
                cerr << "Unable to lookup mass for component:" << ui<< endl;
                fragmentLookupOK=false;
                break;
            }

            //Record the symbol index for ths component
            indexVec.push_back(p.first);
            //Record the count
            freqVec.push_back(p.second);
        }

        //Lookup failed, provide more error data
        if(!fragmentLookupOK)
        {
            cerr << "Component lookup for Fragment :" << parentSpecies[ui] << "failed" << endl;
            continue;
        }



        //Obtain the distribution for this fragment set
        vector<pair<float,float> > massDist;

        natTable.generateGroupedIsotopeDist(indexVec,freqVec,massDist,MASS_TOL);
        massDistVec.emplace_back(massDist);
    }
}