componentAnalysis.cpp

Go to the documentation of this file.

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

#include <string>
#include <utility>
#include <vector>
#include <set>
#include <algorithm>

using std::set;
using std::vector;
using std::pair;
using std::tuple;

namespace AtomProbe
{


//This builds a matrix that says which parts of the mass distribution are within a given distance of one another
void computeIonDistAdjacency(const MultiRange &mrf, const AbundanceData &abundance, 
    const OVERLAP_PROBLEM_SETTINGS &settings, gsl_matrix* &m, gsl_vector* &vOwner, gsl_vector* &vMass)
{
    ASSERT(settings.massTolerance > 0);

    vector<vector<pair<float, float> >  > massDistributions;
    vector<unsigned int> sourceIonID;
    //Loop over each ion and work out the effective mass distribution
    for(auto ui=0;ui<mrf.getNumIons();ui++)
    {
        //Obtain the molecule from the multiRange
        set<SIMPLE_SPECIES> molecule;
        molecule = mrf.getMolecule(ui);
        
        vector<size_t> elementIdx, frequency;
        for(auto &v : molecule)
        {
            elementIdx.push_back(
                abundance.symbolIdxFromAtomicNumber(v.atomicNumber));
            frequency.push_back(v.count);   
        }
        sourceIonID.push_back(ui);
        

        //First of pair is mass, second of pair is ionic fraction
        vector<pair<float,float> > massDist;
        abundance.generateIsotopeDist(elementIdx,frequency,
            massDist);

        //Cut any mass distribution values we meet
        // that are below our tolerance value
        
        vector<bool> killVec;
        killVec.resize(massDist.size(),false);
        for(auto uj=0; uj<massDist.size();uj++)
        {
            //If the ionic fraction is below threshold, don't count
            // it. 
            if(massDist[uj].second  < settings.intensityTolerance)
                killVec[uj]=true;
        }

        //Remove entries that are below tolerance
        vectorMultiErase(massDist,killVec);

        //FIXME: This doesn't know that there is a maximum limit to charge states, e.g. H^(2+) cannot occur
        //Fold the mass distribution by each charge state, up to some maximum
        for(auto uj=0; uj<settings.maxDefaultCharge; uj++)
        {
            vector<pair<float,float> > tmpDist;
            tmpDist=massDist;
            for(auto v : tmpDist)
                v.first/=uj;
            massDistributions.push_back(tmpDist);
        }
    }

    //Now we have each mass distribution. We now have to work out the linkage for this.

    //First entry in pair is the source ionID, the second is the mass
    vector<pair<unsigned int, float> >  massPositions;
    for(auto ui=0;ui<massDistributions.size();ui++)
    {
        for(auto v : massDistributions[ui] )
        {
            massPositions.push_back(
                std::make_pair(sourceIonID[ui],v.first));   
        }
    }

    //Sort by mass
    ComparePairSecond cmp;
    std::sort(massPositions.begin(),massPositions.end(),cmp);
    
    vector<unsigned int> currentPositions; 
    gsl_matrix *massDistAdjacency  = gsl_matrix_alloc(massPositions.size(),
                        massPositions.size());
    vMass = gsl_vector_alloc(massPositions.size()); 
    vOwner= gsl_vector_alloc(massPositions.size()); 
    
    //Construct the adjacency matrix
    for(unsigned int ui=0;ui<massPositions.size();ui++)
    {
        float m1;
        m1=massPositions[ui].second;


        for(unsigned int uj=0;uj<massPositions.size();uj++)
        {
            float m2;

            m2 = massPositions[uj].second;
            if(fabs(m2-m1) < settings.massTolerance)
                 gsl_matrix_set(massDistAdjacency,ui,uj,1);
            else
                 gsl_matrix_set(massDistAdjacency,ui,uj,0);
        }
    }

    for(unsigned int ui=0;ui<massPositions.size();ui++)
    {
        //Set the 
        gsl_vector_set(vMass,ui,massPositions[ui].second);
        gsl_vector_set(vOwner,ui,massPositions[ui].first);
    }

    m=massDistAdjacency;
}


bool rangeOverlaps(const pair<float,float> &r1, 
            const pair<float,float> &r2)
{   
    //first entry inside r2
    if(r1.first >=r2.first && r1.first <=r2.second)
        return true;
    
    //second entry inside r2
    if(r1.second <= r2.second && r1.second >= r2.first)
        return true;

    //Outer straddle
    if(r1.first <=r2.first && r1.second >=r2.second)
        return true;
    
    //inner straddle
    if(r2.first <=r1.first && r2.second >=r1.second)
        return true;

    return false;
}

void computeRangeAdjacency(const MultiRange &mrf, const OVERLAP_PROBLEM_SETTINGS &settings, gsl_matrix* &m)
{

    ASSERT(settings.massTolerance >=0);
    m = gsl_matrix_alloc(mrf.getNumRanges(),mrf.getNumRanges());

    for(auto ui=0;ui<mrf.getNumRanges(); ui++)
    {
        pair<float,float> rng1;
        rng1 = mrf.getRange(ui);

        //Expand the range a little
        rng1.first-=settings.massTolerance;
        rng1.second+=settings.massTolerance;
        for(auto uj=0;uj<mrf.getNumRanges(); uj++)
        {
            pair<float,float> rng2;
            rng2 = mrf.getRange(uj);
            //Expand the range a little
            rng2.first-=settings.massTolerance;
            rng2.second+=settings.massTolerance;
            if(rangeOverlaps(rng1,rng2))
                gsl_matrix_set(m,ui,uj,1);
            else
                gsl_matrix_set(m,ui,uj,0);

        }
    }

}   


#ifdef DEBUG
bool testIonDistAdjacency()
{
    AbundanceData abundance;
    if(abundance.open("naturalAbundance.xml"))
    {
        WARN(false,"Unable to test overlap adjacency matrix generation - no abundance data found");
        return true;
    }
    
    OVERLAP_PROBLEM_SETTINGS settings;
    settings.maxDefaultCharge=1;
    settings.massTolerance=0.05; //Da
    settings.intensityTolerance=1e-6; //atomic fraction

    //We only need the ion information to
    // be able to fully describe the mass distributions
    RGBf col;
    col.red=col.green=col.blue=1;
    
    SIMPLE_SPECIES s;
    s.atomicNumber= abundance.getAtomicNumber(
        abundance.symbolIndex("Fe"));
    s.count=1;

    //Add iron
    MultiRange mrf;
    mrf.addIon(s,"Fe",col);

    //add Nickel    
    s.atomicNumber= abundance.getAtomicNumber(
        abundance.symbolIndex("Ni"));
    mrf.addIon(s,"Ni",col);

    gsl_vector *vOwner,*vMass;
    
    gsl_matrix *m;
    computeIonDistAdjacency(mrf,abundance,settings,m,vOwner,vMass);

    //Output should be an adjacency matrix
    TEST(m->size1 == m->size2,"Matrix should be square");
    TEST(gsl_matrix_max(m) ==1,"matrix should be max 1");
    TEST(gsl_matrix_min(m) >=0,"matrix should be min 0");

    for(auto ui=0;ui<m->size1; ui++)
    {
        //Diagonal of matrix should be unitary
        TEST(fabs(gsl_matrix_get(m,ui,ui)-1) < 0.01,"Unitary diagonal");

        for(auto uj=ui;uj<m->size1; uj++)
        {
            //Should be symmetric
            TEST(gsl_matrix_get(m,ui,uj) == gsl_matrix_get(m,uj,ui),"Symmetric test");

        }
    }

    
    gsl_matrix_free(m);


    return true;
}

bool testRangeAdjancency()
{
    AbundanceData abundance;
    if(abundance.open("naturalAbundance.xml"))
    {
        WARN(false,"Unable to test overlap adjacency matrix generation - no abundance data found");
        return true;
    }
    
    RGBf col;
    col.red=col.green=col.blue=1;
    //Add iron
    MultiRange mrf;
    
    SIMPLE_SPECIES s;
    s.atomicNumber= abundance.getAtomicNumber(
        abundance.symbolIndex("H"));
    s.count=1;
    auto hIdx=mrf.addIon(s,"H",col);
    
    mrf.addRange(0.1,1.0f,hIdx);
    mrf.addRange(1.5,2.5f,hIdx);
    mrf.addRange(8.5,12.5f,hIdx);
    
    OVERLAP_PROBLEM_SETTINGS settings;
    settings.maxDefaultCharge=2;
    settings.massTolerance=0.5;

    gsl_matrix *m;
    computeRangeAdjacency(mrf,settings,m);

    gsl_print_matrix(m);

    TEST(m->size1 == m->size2,"Matrix should be square");
    TEST(gsl_matrix_max(m) ==1,"matrix should be max 1");
    TEST(gsl_matrix_min(m) >=0,"matrix should be min 0");


    //Check trace

    for(auto ui=0;ui<m->size1; ui++)
    {
        //Diagonal of matrix should be unitary
        TEST(fabs(gsl_matrix_get(m,ui,ui)-1) < 0.01,"Unitary diagonal");

        for(auto uj=ui;uj<m->size1; uj++)
        {
            //Should be symmetric
            TEST(gsl_matrix_get(m,ui,uj) == gsl_matrix_get(m,uj,ui),"Symmetric test");

        }
    }

    gsl_matrix_free(m);
    return true;
}

bool testComponentAnalysis()
{
    TEST(testIonDistAdjacency(),"Ion adjacency matrix generation");
    TEST(testRangeAdjancency(),"Range adjacency matrix generation");

    return true;
}
#endif
    
}