massconvert.cpp

Go to the documentation of this file.

/* 
 * abundance.cpp : module to load, manipulate and compute natural abundance information
 * Copyright (C) 2016  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 <algorithm>
#include <iostream>

#include <cmath>
#include <limits>


#include "atomprobe/isotopes/abundance.h"
#include "atomprobe/helper/stringFuncs.h"

#include "atomprobe/helper/aptAssert.h"

#include "atomprobe/isotopes/massconvert.h"


namespace AtomProbe
{
using std::map;
using std::vector;
using std::string;


void convertMolToMass(const AbundanceData &atomTable,
    const map<unsigned int, double> &compositionMols,
    map<unsigned int, double> &compositionMass)
{
    map<unsigned int, double> massValues;

    for(map<unsigned int, double>::const_iterator it=compositionMols.begin(); it!=compositionMols.end();++it)
        massValues[it->first] = atomTable.getNominalMass(it->first);
    

    //compute the total summed molar values
    double totalMols=0;
    for(map<unsigned int, double>::const_iterator it=compositionMols.begin();
        it!=compositionMols.end(); ++it)
        totalMols+=it->second;


    double totalMass=0;
    for(map<unsigned int, double>::const_iterator it=compositionMols.begin();
        it!=compositionMols.end(); ++it)
        totalMass+=massValues.at(it->first)*it->second/totalMols;

    for(map<unsigned int, double>::const_iterator it=compositionMols.begin();
        it!=compositionMols.end(); ++it)
    {
        ASSERT(compositionMass.find(it->first) == compositionMass.end());
        compositionMass[it->first] = it->second/totalMols*massValues.at(it->first)/totalMass;
    }

    
}

//FIXME : Error handling? What if user supplies an atom that does not exist?
void convertMassToMol(const AbundanceData &atomTable,
    const map<unsigned int, double> &compositionMass,
    map<unsigned int, double> &compositionMols)
{
    map<unsigned int, double> massValues;

    for(map<unsigned int,double>::const_iterator it=compositionMass.begin(); it!=compositionMass.end();++it)
        massValues[it->first] = atomTable.getNominalMass(it->first);

    double totalMol=0;
    for(map<unsigned int,double>::const_iterator it=compositionMass.begin(); it!=compositionMass.end();++it)
    {
        double thisMol;
        thisMol= (it->second/massValues.at(it->first));
        compositionMols[it->first] = thisMol;
        totalMol+=thisMol;
    }

    for(map<unsigned int,double>::const_iterator it=compositionMass.begin(); it!=compositionMass.end();++it)
    {
        compositionMols[it->first]/=totalMol;
    }
    
}


unsigned int parseCompositionData(const std::vector<std::string> &atomList, 
    const AbundanceData &massTable, std::map<unsigned int, double> &atomData)
{
    //Dummy composition value to use as a placeholder when computing composition data
    const double BALANCE_TOKEN = -128.0f;
    
    double totalComp=0.0f;
    //Number of times we have seen the "bal" marker before
    unsigned int balanceNumber=0;
    
    for(unsigned int ui=0;ui<atomList.size();ui++)
    {
        std::string s;
        s=stripWhite(atomList[ui]);
        
        vector<string> strs;
        splitStrsRef(s.c_str(),",\t ",strs);

        if(strs.empty())
            break;

        stripZeroEntries(strs);

        if(strs.size() !=2)
            return COMPPARSE_BAD_INPUT_FORMAT;
        
        //Convert bal or numeric string to lowercase
        strs[1]=lowercase(strs[1]);
    
        //Ensure bal value is used the right number of times    
        double value;
        if(stream_cast(value,strs[1]) && strs[1] != "bal")
            return COMPPARSE_BAD_COMP_VALUE;
        else if( strs[1] == "bal" && balanceNumber)
            return COMPPARSE_BAD_BAL_USAGE; 

        //parse atomic data 
        size_t atomicNumber;
        atomicNumber= massTable.symbolIndex(strs[0].c_str())+1;
        if(!atomicNumber)
            return COMPPARSE_BAD_SYMBOL;
        
        if(atomData.find(atomicNumber)  != atomData.end())
            return COMPPARSE_DUPLICATE_SYMBOL;

        //Record atomic number->composition pairing
        if(strs[1] != "bal")
        {
            atomData[atomicNumber-1] = value; 
            totalComp+=value;
        }
        else
        {
            balanceNumber=atomicNumber-1;
            atomData[atomicNumber-1] = BALANCE_TOKEN;
        }
        
    }
    if(balanceNumber)
    {
        if(totalComp > 1.0f )
        {
            return COMPPARSE_BAD_TOTAL_COMP;
        }
        
        atomData[balanceNumber] = 1.0-totalComp;
        
    }

    return 0;
}

#ifdef DEBUG
#include "atomprobe/helper/aptAssert.h"
#include "helper/helpFuncs.h"


bool testMassToMolConversion()
{
    AbundanceData natData;
    if(natData.open("naturalAbundance.xml"))
        return false;

    map<unsigned int, double> molComp,massComp;

    unsigned int idxFe,idxC;
    idxFe=natData.symbolIndex("Fe");
    massComp[idxFe]=0.96;
    idxC=natData.symbolIndex("C");
    massComp[idxC]=0.04;

    convertMassToMol(natData,massComp, molComp);
    
    TEST(massComp.size() ==2, "Mass Composition data size");
    float total=0;
    for(auto it=molComp.begin();it!=molComp.end(); it++)    
    {
        TEST(massComp.find(it->first) != molComp.end(),"Symbol present in both mass and molar comp");
        TEST(it->second <=1 && it->second >=0, "Composition should be in 0->1");    
        total+=it->second;
    }

    TEST(TOL_EQ(total,1.0f),"total composition sums to 1");

    TEST(TOL_EQ(molComp[idxFe],0.8377079),"Fe composition");
    TEST(TOL_EQ(molComp[idxC],0.162292),"C composition");

    return true;
}

bool testMolToMassConversion()
{
    AbundanceData natData;
    if(natData.open("naturalAbundance.xml"))
        return false;

    map<unsigned int, double> molComp,massComp;

    unsigned int idxFe,idxC;
    idxFe=natData.symbolIndex("Fe");
    molComp[idxFe]=0.96;
    idxC=natData.symbolIndex("C");
    molComp[idxC]=0.04;

    convertMolToMass(natData,molComp, massComp);
    
    TEST(massComp.size() ==2, "Mass Composition data size");
    float total=0;
    for(auto it=massComp.begin();it!=massComp.end(); it++)  
    {
        TEST(molComp.find(it->first) != molComp.end(),"Symbol present in both mass and molar comp");    
        total+=it->second;
    }

    TEST(TOL_EQ(total,1.0f),"total composition sums to 1");

    TEST(TOL_EQ(massComp[idxFe],0.99115),"Fe composition");
    TEST(TOL_EQ(massComp[idxC],0.008849557),"C composition");

    return true;
}

bool testMassMolConversion()
{
    TEST(testMolToMassConversion(),"Mol->mass conversion");
    TEST(testMassToMolConversion(),"Mass->Mol conversion");

    return true;
}
#endif
}