deconvolution.cpp

Go to the documentation of this file.

/* deconvolution.cpp :  Mass overlap deconvolution algorithm
 * Copyright (C) 2017  Daniel 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/atomprobe.h"
#include "atomprobe/helper/maths/misc.h"
#include "atomprobe/helper/aptAssert.h"
#include "helper/helpFuncs.h"

#include <utility>
#include <vector>
#include <map>

using std::pair;
using std::map;
using std::vector;


#ifdef DEBUG
//DEBUG
using std::endl;
using std::cout;
#include <string>
#endif

namespace AtomProbe
{

float nullBackground(float rangeStart,float rangeEnd)
{
    return 0;
}


bool leastSquaresOverlapSolve(const AtomProbe::MultiRange &rangeData, const AtomProbe::AbundanceData &abundance, 
                const std::vector<IonHit> &hits, float (*backgroundEstimator)(float,float),
                std::vector<std::pair<unsigned int,float> > &decomposedHits, float rangeTolerance)
{

    if(!rangeData.getNumIons() || 
            !rangeData.getNumRanges() || hits.empty() || rangeTolerance < 0)
        return false;


    vector<MultiRange> decomposedRanges;

    //FIXME: Currently mass tolerance not taken into account
    // when splitting overlapping ranges
    rangeData.splitOverlapping(decomposedRanges,rangeTolerance);

    if(decomposedRanges.empty())
        return false;

    vector<vector<float> > solutions;
    vector<vector<unsigned int> > solvedIonIds;
    vector<float> residuals;
    //Loop through each disconnected overlap problem
    for(const auto &mRng : decomposedRanges)
    {
        //"Flatten" the problem onto the mass axis
        std::vector<FLATTENED_RANGE> flatRange;
        mRng.flattenToMassAxis(flatRange);

        //A unique set of ionIDs in this problem
        vector<unsigned int> containedIonIds;
        containedIonIds.clear();
        {
        std::set<unsigned int> idSet;
        for(const auto &e : flatRange)
        {
            //Add all ionIDs from this range set
            for(const auto id : e.containedIonIDs)
                idSet.insert(id);
        }

        //Move unique identifiers into vector
        containedIonIds.assign(idSet.begin(),idSet.end());

        }

        

        gsl_matrix *m;
        m=gsl_matrix_alloc(flatRange.size(),containedIonIds.size());
        gsl_matrix_set_zero(m);

        //Construct the abundance matrix from the flattened range data
        // for this problem
        for(auto ui=0u; ui< containedIonIds.size();ui++)
        {
            std::set<SIMPLE_SPECIES> ss;
            ss = mRng.getMolecule(containedIonIds[ui]);
            vector<size_t> elemIdx,frequency;

            elemIdx.clear(); frequency.clear();
            for(const auto &mol : ss)
            {
                frequency.push_back(mol.count);
                elemIdx.push_back(abundance.symbolIdxFromAtomicNumber(
                            mol.atomicNumber));

                //Check we were able to identify the atomic number
                if(elemIdx.back() == -1)
                    return false;
            }


            for(auto uj=0u; uj< flatRange.size();uj++)
            {
                float thisAbundance;
                thisAbundance=abundance.abundanceBetweenLimits(elemIdx,frequency,
                        flatRange[uj].startMass,flatRange[uj].endMass);
                gsl_matrix_set(m,uj,ui,thisAbundance);
            }
        }

        //TODO: We could speed this up, as we actually are performing the SVD twice.
        // once is to estimate rank, and the other is to perform the least-squares.
        // this is unnecessary
        unsigned int rank;
        rank = estimateRank(m);

        //FIXME: Warn if exactly constrained (invalid residual)
        if(rank < mRng.getNumIons())
            return false;

        gsl_vector *b;
        b=gsl_vector_alloc(flatRange.size());

        vector<unsigned int> counts;
        counts.resize(flatRange.size());
        for(auto &h : hits)
        {
            for(auto uj=0;uj<flatRange.size();uj++)
            {
                if(flatRange[uj].startMass < h.getMassToCharge()
                    && flatRange[uj].endMass> h.getMassToCharge())
                {
                    counts[uj]++;
                    break;
                }
            }

        }

        if(backgroundEstimator)
        {
            vector<float> background;
            background.resize(counts.size(),0);
            for(auto ui=0;ui<flatRange.size();ui++)
            {
                background[ui]=(*backgroundEstimator)(
                    flatRange[ui].startMass,flatRange[ui].endMass);
            }

            //Preform background subtraction
            for(auto ui=0;ui<counts.size();ui++)
                gsl_vector_set(b,ui,counts[ui]-background[ui]);
        }
        else
        {
            for(auto ui=0;ui<counts.size();ui++)
                gsl_vector_set(b,ui,counts[ui]);
        }


        //Now we have the matrix and vector,
        // perform the actual solve step
        gsl_vector *x;
        x=gsl_vector_alloc(mRng.getNumIons());

        solveLeastSquares(m,b,x);   

        //Copy out data from GSL structures to vectors
        vector<float> soln;
        soln.resize(x->size);
        for(auto ui=0;ui<x->size;ui++)
            soln[ui]=x->data[ui];

        solutions.push_back(soln);
        //FIXME: These are the ionIDs of the child multirange. We need to
        // track them back to their parents. Use the mapping from MultiRange::splitOverlapping to solve this
        solvedIonIds.push_back(containedIonIds);

        //FIXME: residual not tracked

        gsl_matrix_free(m);
        gsl_vector_free(b);
        gsl_vector_free(x);
    }

    //Squish each sub-problem solution back into a global solution
    map<unsigned int,float> decomposeMap;
    ASSERT(solutions.size() == solvedIonIds.size());
    for(auto ui=0;ui<solvedIonIds.size();ui++)
    {
        ASSERT(solutions[ui].size() == solvedIonIds[ui].size());
        for(auto uj=0;uj<solvedIonIds[ui].size();uj++)
        {
            auto it =decomposeMap.find(solvedIonIds[ui][uj]);

            if( it == decomposeMap.end())
                decomposeMap[solvedIonIds[ui][uj]] = solutions[ui][uj];
            else
                it->second+=solutions[ui][uj];
        }
    }

    for(auto &v : decomposeMap)
    {
        decomposedHits.push_back(std::make_pair(v.first,v.second));
    }

    return true;
}

#ifdef DEBUG

void  runFe_Ni_Sim(MultiRange &mRng,const AbundanceData &natTable,
        const unsigned int numIron,
        const unsigned int numNickel,
        vector<ISOTOPE_ENTRY> entriesFe, vector<ISOTOPE_ENTRY> &entriesNi,
        vector<float> ironCDF, vector<float> &nickelCDF,
        vector<float> &feCounts, vector<float> &niCounts)
{
    gsl_rng *rng = AtomProbe::randGen.getRng();

    vector<IonHit> ions;

    //Add iron atoms, randomly from spectra
    unsigned int actualFeCount=0,actualNiCount=0;
    for(auto ui=0;ui<numIron;ui++)
    {
        float r;
        r=gsl_rng_uniform(rng);

        unsigned int atom=-1;
        for(auto uj=0;uj < ironCDF.size();uj++)
        {
            if(r <=ironCDF[uj])
            {
                atom=uj;
                break;
            }
        }

        if(atom == (unsigned int) -1)
            continue;
        ions.push_back(IonHit(0,0,0,entriesFe[atom].mass));
        actualFeCount++;
    }

    //Add Nickel atoms, randomly from spectra
    for(auto ui=0;ui<numNickel;ui++)
    {
        float r;
        r=gsl_rng_uniform(rng);

        unsigned int atom=-1;
        for(auto uj=0;uj < nickelCDF.size();uj++)
        {
            if(r <=nickelCDF[uj])
            {
                atom=uj;
                break;
            }
        }

        if(atom == (unsigned int) -1)
            continue;
        ions.push_back(IonHit(0,0,0,entriesNi[atom].mass));
        actualNiCount++;
    }

    //So we now have overlap, pass to solver
    vector<pair<unsigned int,float> > decomposedResults;
    if(!leastSquaresOverlapSolve(mRng, natTable,
                ions,nullptr,decomposedResults,0.05))
    {
        ASSERT(false);
    }

    ASSERT(decomposedResults.size() ==2);

    for(auto &v : decomposedResults)
    {
        std::string name;
        name=mRng.getIonName(v.first);
        if( name == "Fe")
            feCounts.push_back(v.second);
        else if(name == "Ni")
            niCounts.push_back(v.second);
        else
        {
            ASSERT(false);
        }
    }
}

bool testDeconvolve()
{

    //Construct an Fe-Ni overlap
    MultiRange mRng;

    RGBf col;
    col.red=col.green=col.blue=1.0f;


    SIMPLE_SPECIES ss;
    ss.atomicNumber=26;
    ss.count=1;
    auto ionIdFe=mRng.addIon(ss,"Fe",col);

    
    ss.atomicNumber=28;
    ss.count=1;
    auto ionIdNi=mRng.addIon(ss,"Ni",col);


    AbundanceData natTable;
    auto errCode=natTable.open("../data/naturalAbundance.xml");

    if(errCode)
    {
        WARN(false,"Unable to read abundance table. Skipping test");
        return true;
    }

    //synthesise some ions
    const unsigned int NUM_IRON =100;
    const unsigned int NUM_NICKEL=200;
    vector<float> ironCDF,nickelCDF;

    unsigned int ironIdx,nickelIdx;
    ironIdx=natTable.symbolIndex("Fe");
    nickelIdx=natTable.symbolIndex("Ni");

    ASSERT(ironIdx != (unsigned int) -1);
    ASSERT(nickelIdx != (unsigned int) -1);

    vector<ISOTOPE_ENTRY> entriesFe = natTable.isotopes(ironIdx);
    vector<ISOTOPE_ENTRY> entriesNi = natTable.isotopes(nickelIdx);

    vector<unsigned int> rangeGrouping;

    //Fe+    
    const float MASSTOL =0.1;
    float accum=0;
    for(auto v: entriesFe)
    {
        accum+=v.abundance;
        ironCDF.push_back(accum);
        mRng.addRange(v.mass-MASSTOL,v.mass+MASSTOL,ionIdFe);

        rangeGrouping.push_back(0);
    }

    //Ni+
    accum=0;
    for(auto v: entriesNi)
    {
        accum+=v.abundance;
        nickelCDF.push_back(accum);
        mRng.addRange(v.mass-MASSTOL,v.mass+MASSTOL,ionIdNi);
        rangeGrouping.push_back(1);
    }

    ASSERT( (ironCDF.back() <= 1.0f) );
    ASSERT( (nickelCDF.back() <= 1.0f) );


    mRng.setRangeGroups(rangeGrouping);

    //Fe and Ni count vectors.

    vector<float> counts[2];

    //Accumulate a series of different trials, and statistically check results
    const unsigned int NTRIALS=100;
    for(auto ui=0;ui<NTRIALS;ui++)
    {
        runFe_Ni_Sim(mRng,natTable,NUM_IRON,NUM_NICKEL,
            entriesFe,entriesNi,ironCDF, nickelCDF, counts[0],counts[1]);
    }

    float meanFe,stdFe;
    float meanNi,stdNi;
    meanAndStdev(counts[0],meanFe,stdFe);
    meanAndStdev(counts[1],meanNi,stdNi);

    float errFe=fabs(meanFe-NUM_IRON)/(float)NUM_IRON ;
    float errNi=fabs(meanNi-NUM_NICKEL)/(float)NUM_NICKEL;

    //Check for percentage error. At most 20% error, averaged over all trials
    // if NUM_IRON is low, this will be bad per-trial
    // this is set conservatively, so that a poor user doesn't trip it randomly
    ASSERT(errFe< 0.20);
    ASSERT(errNi< 0.20);

    return true;
}

#endif

}