abundance.cpp

Go to the documentation of this file.

/* 
 * abundance.cpp : module to load, manipulate and compute natural abundance information
 * Copyright (C) 2015  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 <cstdlib>
#include <stack>
#include <map>
#include <cmath>
#include <limits>
#include <algorithm>
#include <numeric>

#include "atomprobe/isotopes/abundance.h"
#include "atomprobe/helper/XMLHelper.h"
#include "atomprobe/helper/stringFuncs.h"
#include "atomprobe/helper/misc.h"
#include "atomprobe/helper/maths/misc.h"


#include "helper/helpFuncs.h"

#include "atomprobe/helper/aptAssert.h"

namespace AtomProbe
{


using std::vector;
using std::pair;
using std::make_pair;
using std::map;
using std::string;

const char *ABUNDANCE_ERROR[] = { "Unable to read abundance data (opening file)",
                 "Unable to create XML reader.",
                 "Bad property found in XML file",
                 "XML document did not match expected layout (DTD validation)",
                 "Unable to find required node during parse",
                 "Root node missing, expect <atomic-mass-table>!",
                 "Found incorrect root node. Expected <atomic-mass-table>"
};


const char *AbundanceData::getErrorText(size_t errorCode)
{
    ASSERT(errorCode < ABUNDANCE_ERROR_ENUM_END);
    return ABUNDANCE_ERROR[errorCode];
}

size_t AbundanceData::numIsotopes() const
{
    size_t v=0;
    for(size_t ui=0;ui<isotopeData.size();ui++)
        v+=isotopeData[ui].size();
    
    return v;
}

size_t AbundanceData::numElements() const
{
    return isotopeData.size();
}

size_t AbundanceData::open(const char *file, bool strict)
{

    xmlDocPtr doc;
    xmlParserCtxtPtr context;

    context =xmlNewParserCtxt();

    if(!context)
        return ABUNDANCE_ERR_NO_CONTEXT;

    //Open the XML file
    doc = xmlCtxtReadFile(context, file, NULL, XML_PARSE_DTDVALID| XML_PARSE_NOENT | XML_PARSE_NONET);

    if(!doc)
    {
        xmlFreeParserCtxt(context);
        return ABUNDANCE_ERR_BAD_DOC;
    }
    else
    {
        //Check for context validity
        if(!context->valid)
        {
            if(strict)
            {
                xmlFreeDoc(doc);
                xmlFreeParserCtxt(context);
                return ABUNDANCE_ERROR_FAILED_VALIDATION;
            }
        }
    }

    try
    {
    
    //retrieve root node    
    std::stack<xmlNodePtr> nodeStack;
    xmlNodePtr nodePtr = xmlDocGetRootElement(doc);

    if(!nodePtr)
        throw ABUNDANCE_ERROR_MISSING_ROOT_NODE;
    
    //This *should* be an abundance file
    if(xmlStrcmp(nodePtr->name, (const xmlChar *)"atomic-mass-table"))
        throw ABUNDANCE_ERROR_WRONG_ROOT_NODE;

    nodeStack.push(nodePtr);

    nodePtr=nodePtr->xmlChildrenNode;
    while(!XMLHelpFwdToElem(nodePtr,"entry"))
    {
        ISOTOPE_ENTRY curIsoData;
        
        if(XMLHelpGetProp(curIsoData.symbol,nodePtr,"symbol"))
            throw ABUNDANCE_ERROR_BAD_VALUE;
        
        if(XMLHelpGetProp(curIsoData.atomicNumber,nodePtr,"atomic-number"))
            throw ABUNDANCE_ERROR_BAD_VALUE;


        nodeStack.push(nodePtr);
        nodePtr=nodePtr->xmlChildrenNode;
        
        //Move to natural-abundance child
        if(XMLHelpFwdToElem(nodePtr,"natural-abundance"))
            throw ABUNDANCE_ERROR_MISSING_NODE;
        
        nodePtr=nodePtr->xmlChildrenNode;

        vector<ISOTOPE_ENTRY> curIsotopes;  
        //TODO: value checking  
        while(!XMLHelpFwdToElem(nodePtr,"isotope"))
        {
            //Spin to mass node
            if(XMLHelpGetProp(curIsoData.massNumber,nodePtr,"mass-number"))
                throw ABUNDANCE_ERROR_MISSING_NODE;
            
            nodeStack.push(nodePtr);
            nodePtr=nodePtr->xmlChildrenNode;

            //Spin to mass node
            if(XMLHelpFwdToElem(nodePtr,"mass"))
                throw ABUNDANCE_ERROR_MISSING_NODE;

            if(XMLHelpGetProp(curIsoData.mass,nodePtr,"value"))
                throw ABUNDANCE_ERROR_BAD_VALUE;
            
            if(XMLHelpGetProp(curIsoData.massError,nodePtr,"error"))
                throw ABUNDANCE_ERROR_BAD_VALUE;
            else
                curIsoData.massError=0;


            //spin to abundance node
            if(XMLHelpFwdToElem(nodePtr,"abundance"))
                throw ABUNDANCE_ERROR_MISSING_NODE;
            
            if(XMLHelpGetProp(curIsoData.abundance,nodePtr,"value"))
                throw ABUNDANCE_ERROR_BAD_VALUE;
            
            if(XMLHelpGetProp(curIsoData.abundanceError,nodePtr,"error"))
                throw ABUNDANCE_ERROR_BAD_VALUE;
            else
                curIsoData.abundanceError=0;

            curIsotopes.push_back(curIsoData);

            nodePtr=nodeStack.top();
            nodePtr=nodePtr->next;
            nodeStack.pop();
        }

        isotopeData.push_back(curIsotopes);
        curIsotopes.clear();    

        nodePtr=nodeStack.top();
        nodeStack.pop();
        nodePtr=nodePtr->next;
        
    }
    }
    catch( int &excep)
    {
        xmlFreeDoc(doc);
        xmlFreeParserCtxt(context);
        return excep;
    }
        
    xmlFreeDoc(doc);
    xmlFreeParserCtxt(context);
    return 0;
}


size_t AbundanceData::symbolIndex(const char *symbol, bool caseSensitive) const
{
    ASSERT(isotopeData.size());
    if(!caseSensitive)
    {
        for(size_t ui=0;ui<isotopeData.size();ui++)
        {
            if(isotopeData[ui].size() && lowercase(isotopeData[ui][0].symbol) == lowercase(symbol) )
                return ui;
        }
    }
    else
    {
        for(size_t ui=0;ui<isotopeData.size();ui++)
        {
            if(isotopeData[ui].size() && isotopeData[ui][0].symbol == symbol)
                return ui;
        }
    }

    return (size_t)-1;
}


void AbundanceData::isotopeIndex(size_t elemIn, float mass ,size_t &isoIndex) const
{
    ASSERT(isotopeData.size());
    for(size_t uj=0;uj<isotopeData[elemIn].size();uj++)
    {
        if(isotopeData[elemIn][uj].mass == mass)
        {
            isoIndex=uj;
            return;
        }
    }

    isoIndex=(size_t)-1;
    return;
}

void AbundanceData::isotopeIndex(size_t elem, float mass, float tolerance, size_t &isoIndex) const
{
    ASSERT(isotopeData.size());
    for(size_t uj=0;uj<isotopeData[elem].size();uj++)
    {
        if(fabs(isotopeData[elem][uj].mass - mass) < tolerance)
        {
            isoIndex=uj;
            return;
        }
    }

    isoIndex=(size_t)-1;
    return;
}


const ISOTOPE_ENTRY &AbundanceData::isotope(size_t elemIdx,size_t isotopeIdx) const
{
    ASSERT(isotopeData.size());
    return isotopeData[elemIdx][isotopeIdx];
}

void AbundanceData::generateIsotopeDist(const vector<size_t> &elementIdx,
                    const vector<size_t> &frequency,
                vector<pair<float,float> > &massDist) const
{
    ASSERT(isotopeData.size());
    ASSERT(frequency.size() == elementIdx.size());
    //Search out the given isotopes, and compute the peaks
    // that would be seen for this isotope combination

    //map of vectors for the available isotopes
    map<unsigned int, vector< pair<float,float> > > isotopeMassDist;

    //For each isotope, retrieve  its (mass,probability dist) values
    // placing into a map
    //--
    for(size_t ui=0;ui<elementIdx.size();ui++)
    {
        size_t curIso;
        curIso = elementIdx[ui];

        //If we have seen isotope before, move on
        if(isotopeMassDist.find(curIso) != isotopeMassDist.end())
            continue;

        vector<pair<float,float> > thisIsotope;
        for(size_t uj=0;uj<isotopeData[curIso].size();uj++)
        {
            //Compute the probability of this isotope's presence in
            // the sampled concentration
            thisIsotope.push_back(make_pair(isotopeData[curIso][uj].mass,
                isotopeData[curIso][uj].abundance));
        }

        isotopeMassDist[curIso] = thisIsotope;
        thisIsotope.clear();

    }
    //--


    //Given the isotopes we have, permute the mass spectra
    vector<pair<float,float> > peakProbs;
    for(size_t ui=0;ui<elementIdx.size();ui++) //Loop over element type
    {

        for(size_t repeat=0;repeat<frequency[ui];repeat++) //Loop over the numberof times element occurs
        {
            vector< pair<float,float> >::const_iterator  isoBegin,isoEnd;
            isoBegin = isotopeMassDist[elementIdx[ui]].begin();
            isoEnd = isotopeMassDist[elementIdx[ui]].end();

            vector<pair<float,float> > newProbs;
            //If this is the very first item in our list,
            // simply push on the value, rather than modifying the
            // distribution
            if(peakProbs.empty())
            {
                //The masses will be added to, and the probabilities multiplied 
                for(vector<pair<float,float> >::const_iterator it=isoBegin;
                    it!=isoEnd;++it)
                    peakProbs.push_back(*it);
            }
            else
            {
                const unsigned int peakProbSize = peakProbs.size();
                newProbs.resize(peakProbSize*isotopeMassDist[elementIdx[ui]].size());
                unsigned int offset=0;

                for(size_t uj=0;uj<peakProbSize;uj++) //Loop over mass distribution
                {
            
                    //Convolve it with other peaks
                    //The masses will be added to, and the probabilities multiplied 
                    for(vector<pair<float,float> >::const_iterator it=isoBegin;
                        it!=isoEnd;++it) 
                    {
                        pair<float,float> newMass;

                        newMass.first=peakProbs[uj].first+it->first;
                        newMass.second=peakProbs[uj].second*it->second;
                        newProbs[offset]=newMass;
                        offset++;
                    }

                }
                
                peakProbs.swap(newProbs);
                newProbs.clear();
            }
        }
    }


    float tolerance=sqrt(std::numeric_limits<float>::epsilon());

    vector<bool> killPeaks(peakProbs.size(),false);
    const unsigned int peakProbSize=peakProbs.size();
#pragma omp parallel for
    //Find the non-unique peaks and sum them
    for(size_t ui=0;ui<peakProbSize;ui++)
    {
        for(size_t uj=ui+1;uj<peakProbSize;uj++)
        {
            if(fabs(peakProbs[ui].first-peakProbs[uj].first) <tolerance &&
                peakProbs[uj].second > 0.0f)
            {

                peakProbs[ui].second+=peakProbs[uj].second;
                peakProbs[uj].second=0.0f;
                killPeaks[uj]=true;
            }
        }
    }
    
    vectorMultiErase(peakProbs,killPeaks);
    massDist.swap(peakProbs);
}

void AbundanceData::generateGroupedIsotopeDist(const vector<size_t> &elementIdx,
                        const vector<size_t> &frequency, vector<pair<float,float> > &massDist, 
                        float massTolerance) const
{
    //First, generate the isotope distribtion
    generateIsotopeDist(elementIdx,frequency, massDist);

    //We can't group zero or 1 results
    if(massDist.size() < 2)
        return;

    //Now sort the isotope distribution by mass
    ComparePairFirst cmp;
    std::sort(massDist.begin(),massDist.end(),cmp);

    //Now we have to accumulate overlapping entries.

    vector<vector<pair<float,float> > > overlappingEntries;

    vector<pair<float,float> > curGroup;
    curGroup.push_back(massDist[0]);

    float threshPos=massDist[0].first + massTolerance;
    for(unsigned int ui=1;ui<massDist.size(); ui++)
    {
        if(massDist[ui].first > threshPos)
        {
            //Terminate the current group
            overlappingEntries.push_back(curGroup);
            curGroup.clear();


        }
        curGroup.push_back(massDist[ui]);
        threshPos=massDist[ui].first + massTolerance;
    }

    overlappingEntries.push_back(curGroup);

    massDist.clear();
    //OK, so now we have each overlap set as a group, average each
    // group's mass to provide a good mass estimate, but sum the abundance intensities
    for(unsigned int ui=0;ui<overlappingEntries.size(); ui++)
    {
        vector<float> x,y;

        for(unsigned int uj=0;uj<overlappingEntries[ui].size();uj++)
        {
            x.push_back(overlappingEntries[ui][uj].first);
            y.push_back(overlappingEntries[ui][uj].second);
        }

        float wMean;
        wMean=weightedMean(x,y);

        float summedAbundance;
        summedAbundance=std::accumulate(y.begin(),y.end(),0.0f);

        massDist.push_back(make_pair(wMean,summedAbundance));


    }




}

float AbundanceData::abundanceBetweenLimits( const std::vector<size_t> &elemIdx,
                        const std::vector<size_t> &frequency,
                        float massStart, float massEnd) const
{
    vector<pair<float,float> > massDist;
    generateIsotopeDist(elemIdx,frequency,massDist);
    float abundance=0;
    for(const auto &isotope : massDist)
    {
        if(isotope.first >=massStart &&  isotope.first < massEnd)
            abundance+=isotope.second;


    }

    ASSERT(abundance >=0 && abundance <=1);
    return abundance;
}

size_t AbundanceData::getMajorIsotopeFromElemIdx(size_t elementIdx) const
{
    ASSERT(elementIdx < isotopeData.size());
    const std::vector<ISOTOPE_ENTRY> &curIsotopes= isotopeData[elementIdx];

    ASSERT(isotopeData.size());
    
    size_t maxIdx=0;
    float maxV=curIsotopes[0].abundance;
    for(size_t ui=1; ui<curIsotopes.size();ui++)
    {
        if(curIsotopes[ui].abundance > maxV )
        {
            maxIdx=ui;
            maxV=curIsotopes[ui].abundance;
        }
    }

    return maxIdx;
}

unsigned int AbundanceData::getAtomicNumber(size_t elementIdx) const
{
    ASSERT(elementIdx < isotopeData.size());
    const std::vector<ISOTOPE_ENTRY> &curIsotopes= isotopeData[elementIdx];

#ifdef DEBUG
    for(size_t ui=1;ui<curIsotopes.size();ui++)
    {
        ASSERT(curIsotopes[ui].atomicNumber == curIsotopes[0].atomicNumber);
    }
#endif
    return curIsotopes[0].atomicNumber;
}

void AbundanceData::generateSingleAtomDist(size_t atomIdx, unsigned int repeatCount, 
                        std::vector<std::pair<float,float> > &massDist) const
{
    massDist.clear();
        
    //Obtain the single Isotope distribution
    vector<pair<float,float> > singleDist;
    const std::vector<ISOTOPE_ENTRY> &curIso=isotopeData[atomIdx];
    for(size_t ui=0;ui<curIso.size();ui++)
    {
        singleDist.push_back(make_pair(curIso[ui].mass,curIso[ui].abundance));
    }

    //Duplicate the single dist for the first atom
    massDist=singleDist;

    //For the number of elements, branch the distribution
    // out from this using the single mass distribution to continually
    // create new members 
    for(size_t ui=1;ui<repeatCount;ui++)
    {
        vector<pair<float,float> >  newDist;
        //Add the new element to the multiDist
        for(size_t uj=0;uj<singleDist.size();uj++)
        {
            for(size_t uk=0;uk<massDist.size();uk++)
            {
                pair<float,float> newMass;
                newMass.first = massDist[uk].first + singleDist[uj].first;
                newMass.second = massDist[uk].second*singleDist[uj].second;
                newDist.push_back(newMass);
            }
        }
        massDist.swap(newDist);
        newDist.clear();
    }


}



void AbundanceData::getSymbolIndices(const vector<string> &symbols,vector<size_t> &indices) const
{
    indices.resize(symbols.size());
    for(size_t ui=0;ui<symbols.size(); ui++)
        indices[ui]=symbolIndex(symbols[ui].c_str());

}

void AbundanceData::getSymbolIndices(const vector<pair<string,size_t> > &symbols,
            vector<pair<size_t,unsigned int >> &indices) const
{
    //Extract first
    vector<string> sV;
    vector<size_t> iV;
    sV.resize(symbols.size()); iV.resize(symbols.size());
    for(unsigned int ui=0;ui<symbols.size(); ui++)
    {
        sV[ui]= symbols[ui].first;
        iV[ui]= symbols[ui].second;
    }

    vector<size_t> indicesV;
    getSymbolIndices(sV,indicesV);

    //Repackage into pair
    indices.resize(sV.size());
    for(unsigned int ui=0; ui<indices.size();ui++)
    {
        indices[ui].first = indicesV[ui];  //Index
        indices[ui].second = iV[ui]; //Payload
    }

}

std::string AbundanceData::elementNames(size_t start,size_t end, char separator) const
{
    ASSERT(start < end &&  isotopeData.size());
    std::string retStr,sep;
    sep=separator;
    for(size_t ui=start; ui<end; ui++)  
    {
        size_t symbolIdx;
        symbolIdx=symbolIdxFromAtomicNumber(ui);
        if(symbolIdx !=(size_t)-1)
            retStr+=elementName(symbolIdx) + sep;
    }

    return retStr;
}

size_t AbundanceData::symbolIdxFromAtomicNumber(size_t atomicNumber) const
{
    for(size_t ui=0;ui<isotopeData.size(); ui++)
    {
        if(isotopeData[ui].size() && 
            isotopeData[ui][0].atomicNumber == atomicNumber)
            return ui;
    }

    return (size_t)-1;
}

const std::vector<ISOTOPE_ENTRY> &AbundanceData::isotopes(size_t offset) const 
{
    ASSERT(offset < isotopeData.size()); return isotopeData[offset];
}

float AbundanceData::getNominalMass(size_t elementIdx) const
{

    const vector<ISOTOPE_ENTRY> &isoDat=isotopeData[elementIdx];

    float total=0;
    for(size_t ui=0;ui<isoDat.size();ui++)
    {
        total+=isoDat[ui].abundance*isoDat[ui].mass;
    }

    //return mean
    return total;
    
}

size_t AbundanceData::getNearestCharge(const vector<pair<size_t,size_t> > &element,float targetMass,size_t maxCharge) const
{
    //find the nearest isotope
    vector<pair<float,float> > isotopeDist;

    vector<size_t> elemIdx, elemFreq;
    elemIdx.resize(element.size());
    elemFreq.resize(element.size());

    for(size_t ui=0;ui<element.size();ui++)
    {
        elemIdx[ui] = element[ui].first;
        elemFreq[ui] = element[ui].second;
    }

    generateIsotopeDist(elemIdx,elemFreq,isotopeDist);

    float minDelta=std::numeric_limits<float>::max();
    //this is the position of the isotope in the dist
    size_t bestCharge=0;
    for(size_t curCharge=1;curCharge<maxCharge; curCharge++)
    {
        for(size_t ui=0;ui<isotopeDist.size();ui++)
        {
            float curDelta;
            curDelta=fabs(isotopeDist[ui].first/curCharge-targetMass);
            if( curDelta< minDelta)
            {
                minDelta = curDelta;
                bestCharge=curCharge;
            }
        }
    
    }

    return bestCharge;
}

#ifdef DEBUG
void AbundanceData::checkErrors() const
{
    //Ensure all isotopes sum to 1-ish
    // Rounding errors limit our correctness here.
    for(size_t ui=0;ui<isotopeData.size();ui++)
    {
        if(!isotopeData[ui].size())
            continue;

        float sum;
        sum=0.0f;
        for(size_t uj=0; uj<isotopeData[ui].size();uj++)
        {
            sum+=isotopeData[ui][uj].abundance;
        }

        ASSERT(fabs(sum -1.0f) < 0.000001);

    }


    //Ensure Ti has 5 isotopes (original data file was missing)
    ASSERT(isotopeData[symbolIndex("Ti")].size() == 5);
}

bool AbundanceData::runUnitTests()
{
    AbundanceData massTable;
    TEST(massTable.open("../data/naturalAbundance.xml") == 0,"load table");
    //FIXME: Getting the isotope dis

    size_t ironIndex=massTable.symbolIndex("Fe");
    TEST(ironIndex != (size_t)-1,"symbol lookup");

    //Generate the mass peak dist for iron
    vector<size_t> elements;
    vector<size_t> concentrations;
    elements.push_back(ironIndex);
    concentrations.push_back(1);

    std::vector<std::pair<float,float> > massDist;
    massTable.generateIsotopeDist(elements,concentrations,massDist);

    TEST(massDist.size() == 4, "Iron has 4 isotopes");


    elements.clear();
    massDist.clear();

    //Check the isotope dist (grouped) for Mo2
    //==========
    size_t molyIdx = massTable.symbolIndex("Mo");
    concentrations.clear();
    //Mo2
    elements.push_back(molyIdx);
    concentrations.push_back(2);

    
    massTable.generateGroupedIsotopeDist(elements,concentrations,massDist,0.005f);

    // AJL Reports the following
    const float MO2_ISOTOPES[15][2] ={
        {183.8136,0.022},
        {185.8119,0.0275},
        {186.8126,0.0473},
        {187.811,0.0581},
        {188.8119,0.0578},
        {189.8111,0.1278},
        {190.8108,0.0708},
        {191.8118,0.1315},
        {192.811,0.1087},
        {193.8115,0.1074},
        {194.8124,0.0768},
        {195.8117,0.0904},
        {196.8135,0.0184},
        {197.8129,0.0465},
        {199.8149,0.0093},
    };

    ComparePairFirst cmp;
    std::sort(massDist.begin(),massDist.end(),cmp);

    const unsigned int NUM_ISOTOPES_MO2=15;
    TEST(massDist.size() == NUM_ISOTOPES_MO2,"Mo2 dist test");

    for(unsigned int ui=0;ui<NUM_ISOTOPES_MO2; ui++)
    {
        float m,i;
        m=MO2_ISOTOPES[ui][0];
        i=MO2_ISOTOPES[ui][1];

        float dm,di;
        dm = fabs(m-massDist[ui].first);
        di = fabs(i-massDist[ui].second);
        TEST(dm < 0.01,"Mo2 Isotope position");
        TEST(di < 0.01,"Mo2 Isotope intensity");

    }
    //==========


    return true;    

}
#endif
}