rangeCheck.cpp

Go to the documentation of this file.

/*
 * rangeCheck.cpp : Rangefile sanity checking routines 
 * Copyright (C) 2017  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 "atomprobe/algorithm/rangeCheck.h"
#include "atomprobe/helper/aptAssert.h"

#include <string>
#include <utility>

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

namespace AtomProbe {

//This holdes information about decomposed range information 
struct RANGE_MOLECULE
{
    bool isOK;

    string name;
    //First is element symbol offset (eg H->0).
    // second is multiplicity (eg H2 has value 2)
    vector<pair<size_t,size_t> > components;

};

//Check to see if a value is contained within a given pair
bool pairContains(const pair<float,float> &p, float m)
{
    ASSERT(p.first < p.second);
    return (p.first <= m && p.second >=m);
}

//returns a range "molecule" - a wrapper for decomposition of a range
// if an ion cannot be broken down into components (eg unknown ion), then
// an empty molecule is returned, with isOK=false.
RANGE_MOLECULE getRangeMolecule(const AbundanceData &massTable,
        const RangeFile &rng, unsigned int rangeId)
{
    RANGE_MOLECULE mol;
    auto ionId=rng.getIonID(rangeId);
    mol.name=rng.getName(ionId);
    mol.components.clear();
    mol.isOK=true;

    //  
    size_t symbolIdx;
    symbolIdx=massTable.symbolIndex(mol.name.c_str());

    if(symbolIdx ==(size_t)-1)
    {
        vector<pair<string,size_t> > namedComponent;
        //try decomposing. If it fails, return a bad molecule
        if(!RangeFile::decomposeIonNames(mol.name,namedComponent))
        {
            mol.isOK=false;
            return mol;
        }

        //Check decomposed elements are identifiable
        for(size_t uj=0;uj<namedComponent.size();uj++)
        {
            symbolIdx=massTable.symbolIndex(namedComponent[uj].first.c_str());
            if( symbolIdx == (size_t)-1)
            {
                mol.isOK=false;
                return mol;
            }
        
        
            //Store component
            mol.components.push_back(make_pair(symbolIdx,namedComponent[uj].second));
        }
    }
    else
    {
        mol.components.push_back(make_pair(symbolIdx,1));
    }

    return mol;
}

void checkMassRangingCorrectness(const RangeFile &rng,AbundanceData &massTable,
        float massTolerance, unsigned int maxChargeState, 
        unsigned int maxComponents,std::vector<bool> &badRanges)
{
    badRanges.resize(rng.getNumRanges(),false);
        
    for(size_t ui=0;ui<rng.getNumRanges();ui++)
    {
        //Check we can decompose a range
        RANGE_MOLECULE rngMol;
        rngMol=getRangeMolecule(massTable,rng,ui);
        if(rngMol.components.empty() || !rngMol.isOK)
            continue;
        
        //Work out the number of atoms in the molecule
        size_t totalComponents;
        totalComponents=0;
        for(size_t uk=0;uk<rngMol.components.size();uk++)
            totalComponents+=rngMol.components[uk].second;
    
        //Skip check if there are too many components, that will cause
        // problems when checking   
        if(totalComponents > maxComponents)
            continue;


        pair<float,float> thisRange;
        thisRange=rng.getRange(ui);
        
        //widen range by mass 
        thisRange.first-=massTolerance;
        thisRange.second+=massTolerance;

        //Convert rngMol to more "friendly" form for isotope distribution
        // generation
        vector<size_t> elements;
        vector<size_t> composition;
        for(size_t uj=0;uj<rngMol.components.size();uj++)
        {
            for(size_t uk=0;uk<rngMol.components[uj].second; uk++)
            {
                elements.push_back(rngMol.components[uj].first);
                composition.push_back(1);
            }
        }
        
        //Generate isotope distribution
        vector<pair<float,float> > isotopeDist;
        massTable.generateIsotopeDist(elements,composition,isotopeDist);

        //The isotopeDist ranged must match a theoretical range
        // Check each range for several differing chargestates
        bool foundMass;
        foundMass=false;
        for(size_t uj=1;uj<=maxChargeState;uj++)
        {
            for(size_t um=0;um<isotopeDist.size();um++)
            {
                if(pairContains(thisRange,isotopeDist[um].first/(float)uj))
                {
                    foundMass=true;
                    break;
                }
            }
            if(foundMass)
                break;
        }

        //If we were not able to find an isotope that fits within
        // the specified range, then this is an incorrect range 
        if(!foundMass)
            badRanges[ui]=true;
    }

}


#ifdef DEBUG
bool testRangeChecking()
{
    //Perform direct range-test
    // should be no incorrectranges
    // then add a bad range
    AbundanceData natData;
    if(natData.open("naturalAbundance.xml"))
        return false;

    RangeFile rng;
    RGBf c;
    c.red=c.green=c.blue=1.0f;

    auto ionIdTi=rng.addIon("Ti","Ti",c);
    auto ionIdH = rng.addIon("H","H",c);

    //Add Ranges
    auto rngH=rng.addRange(0.95,1.05,ionIdH);
    auto rngTi=rng.addRange(23.9,24.1,ionIdTi);


    vector<bool> badRanges;
    checkMassRangingCorrectness(rng,natData,
            0.05,3,3,badRanges);

    TEST(badRanges.size() == 2,"Check bad range count");

    TEST((badRanges[0] == false) && (badRanges[1] == false)
        ,"Ranges should be OK");

    //Move Ti to a bad mass
    rng.moveBothRanges(rngTi, 17,18);
    checkMassRangingCorrectness(rng,natData,
            0.05,3,3,badRanges);
    TEST((badRanges[rngH] == false) && (badRanges[rngTi] == true)
        ,"Ranges should be wrong");

    return true;
}
#endif
}