multiRange.cpp

Go to the documentation of this file.

/*
 * multiRanges.cpp - Atom probe rangefile class, with multiple identities 
 * 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/>.
 */

//Public headers
#include "atomprobe/io/multiRange.h"
#include "atomprobe/helper/version.h"
#include "atomprobe/helper/XMLHelper.h"
#include "atomprobe/isotopes/overlaps.h"

//Local headers
#include "atomprobe/helper/aptAssert.h"
#include "helper/helpFuncs.h"


#include <algorithm>
#include <fstream>
#include <iomanip>
#include <ctime>
#include <memory>
#include <stack>
#include <map>


namespace AtomProbe
{

using std::pair;
using std::map;
using std::string;

const char MULTIRANGE_FORMAT_VERSION[] = "0.0.1";

enum
{
    MULTIRANGE_ERR_NO_PARSER=1,
    MULTIRANGE_ERR_NO_XML_DOC,
    MULTIRANGE_ERR_NO_XML_CONTEXT,
    MULTIRANGE_ERR_NO_NODEPTR,
    MULTIRANGE_NOT_VALID_ROOTNODE,
    MULTIRANGE_ERR_NO_CONTENT,
    MULTIRANGE_ERR_NO_TIMESTAMP,
    MULTIRANGE_ERR_NO_MOLECULES,
    MULTIRANGE_ERR_NO_MOLECULE_ENTRIES,
    MULTIRANGE_ERR_MOL_MISSING_ENTRY,
    MULTIRANGE_ERR_MOL_MISSING_NAME,
    MULTIRANGE_ERR_MOL_MISSING_COMPONENTS,
    MULTIRANGE_ERR_BAD_RANGE,
    MULTIRANGE_ERR_BAD_RANGES,
    MULTIRANGE_ERR_MOL_MISSING_COLOUR,
    MULTIRANGE_ERR_MOL_BAD_COLOUR,
    MULTIRANGE_ERR_NO_RANGES,
    MULTIRANGE_ERR_BAD_GROUP,
    MULTIRANGE_ERR_SELF_INCONSISTENT,
    MULTIRANGE_ERR_RANGE_NOPARENT,
    MULTIRANGE_ERR_ENUM_END,
};


using std::set;
using std::vector;


//Given a set of identifiers with a specific position
// where the value points to another position
// |0 1 2 3 4|
// {0,1,2,1,3}, link all identifiers such that 
// they have the same base value
// eg the above would become {0,0,1,0,0}
// as 0->0, 1->0,2->2, 1->0, 3->1->0
// to give  {0,0,2,0,0}
// and we renumber them so they fill a contigous set of numbers
template<class T>
void linkIdentifiers(vector<T> &link)
{
    //Chain the identifiers to make them unique
    for(auto ui=0; ui<link.size(); ui++)
    {
        if(link[ui] < ui)
            link[ui] = link[link[ui]];
    }

    //Make them contiguous
    map<size_t,size_t> uniqMap;
    for(auto &v : link)
    {
        if(uniqMap.find(v) != uniqMap.end())
            v=uniqMap[v];
        else
        {
            size_t mapSize=uniqMap.size();
            uniqMap[v]=mapSize;
            v=mapSize;
        }
    }
}

bool SIMPLE_SPECIES::operator==(const SIMPLE_SPECIES &oth) const
{ 
    return (atomicNumber == oth.atomicNumber && count == oth.count) ;
}

bool operator<(const SIMPLE_SPECIES &a, const SIMPLE_SPECIES &b)
{
    return std::tie(a.atomicNumber,a.count) < std::tie(b.atomicNumber,b.count);
}

bool pairOverlaps(float aStart, float aEnd, float bStart,float bEnd)
{
    ASSERT(aStart < aEnd);
    ASSERT(bStart < bEnd);

    if(aStart >=bStart && aStart<= bEnd)
        return true;

    if(aEnd >=bStart && aEnd<= bEnd)
        return true;

    return false;
}

MultiRange::MultiRange() :errState(0)
{
}


MultiRange::MultiRange(const RangeFile &r, const AbundanceData &aTable) :errState(0)
{
    set<size_t> badIons;

    for(auto i=0; i <r.getNumIons(); i++)
    {
        //try to decompose the ion name into usable fragements
        vector<pair<string,size_t> > fragments;
        if(!RangeFile::decomposeIonNames(r.getName(i),fragments))
        {
            badIons.insert(i);
            continue;
        }

        //Obtain the symbol index for each fragment
        vector<pair<size_t,unsigned int> > symbolIndices;
        aTable.getSymbolIndices(fragments,symbolIndices);

        bool bad;
        bad=false;

        //work out what this molecule is
        set<SIMPLE_SPECIES> molecule;
        molecule.clear();
        for(const auto &idxPair: symbolIndices)
        {
            if(idxPair.first == -1)
            {
                bad=true;
                break;
            }
            
            SIMPLE_SPECIES species;

            species.atomicNumber= aTable.getAtomicNumber(idxPair.first);
            species.count=idxPair.second;

            molecule.insert(species);

        }

        if(bad)
        {
            badIons.insert(i);
            continue;
        }
        
        moleculeData.push_back(molecule);
        colours.push_back(r.getColour(i));
        ionNames.push_back(r.getName(i));
    }

    //Find each range, check if we were able to convert the ranges to a multirange
    for(unsigned int i=0; i<r.getNumRanges(); i++)
    {
        //Ignore any ranges that belonged to ions we couldnt identify
        unsigned int ionID;
        ionID =  r.getIonID(i);
        if(badIons.find(ionID) != badIons.end())
            continue;
        ranges.push_back(r.getRange(i));
        ionIDs.push_back(r.getIonID(i));
    }

    ASSERT(isSelfConsistent());
}

bool MultiRange::operator==(const MultiRange &other) const
{
    if(ionNames != other.ionNames)
        return false;
    if(moleculeData != other.moleculeData)
        return false;
    if(colours != other.colours)
        return false;

    if(ranges != other.ranges)
        return false;

    if(rangeGrouping != other.rangeGrouping)
        return false;

    if(ionIDs !=other.ionIDs)
        return false;

    return true;
}


bool MultiRange::isSelfConsistent() const
{
    unsigned int ionSize=moleculeData.size();


    if(ionSize !=colours.size())
        return false;

    unsigned int rangeSize=ranges.size();

    if(rangeSize != ionIDs.size())
        return false;


    for(auto &i :moleculeData)
    {
        if(i.empty())
            return false;
    }
    
    for(auto &name : ionNames )
    {
        if(name.empty())
            return false;
    }

    for(auto &rng : ranges)
    {
        if(rng.first >= rng.second)
            return false;
    }


    //Range groups should either be absent, or identically sized to ranges
    if((rangeGrouping.size()) && (rangeGrouping.size() != ranges.size()))
        return false;

    return true;
}


unsigned int MultiRange::addIon(const set<SIMPLE_SPECIES> &ions, 
                const std::string &name, const RGBf &ionCol)
{
    ASSERT(isSelfConsistent());

    //Molecule data must be unique
    for(auto i:moleculeData)
    {
        if( i == ions)
            return -1;
    }

    moleculeData.push_back(ions);   
    colours.push_back(ionCol);
    ionNames.push_back(name);


    ASSERT(isSelfConsistent());

    return moleculeData.size()-1;
}

unsigned int MultiRange::addIon(const SIMPLE_SPECIES &ion, 
                const std::string &name, const RGBf &ionCol)
{
    ASSERT(isSelfConsistent());

    set<SIMPLE_SPECIES> sSet;
    sSet.insert(ion);
    return addIon(sSet,name,ionCol);
}

unsigned int MultiRange::addRange(float start, float end, unsigned int ionID)
{
    ASSERT(ionID < moleculeData.size());
    ASSERT(isSelfConsistent());
    if(end <= start)
        return -1;

    ionIDs.push_back(ionID);
    ranges.push_back(std::make_pair(start,end));

    return ionIDs.size()-1;
}

unsigned int MultiRange::addRange(const pair<float, float > &rng, unsigned int ionID)
{
    return addRange(rng.first,rng.second,ionID);
}

void MultiRange::setRangeGroups(const std::vector<unsigned int> &groups)
{
    ASSERT(groups.empty() || groups.size() == ranges.size());
    rangeGrouping=groups;
}

bool MultiRange::write(const char *fileName, unsigned int format) const
{
    std::ofstream f(fileName);
    if(!f)
        return false;

    ASSERT(isSelfConsistent());
    using std::endl;
    switch(format)
    {
        case MULTIRANGE_FORMAT_XML:
        {

            f << "<multirange version=\"" << MULTIRANGE_FORMAT_VERSION <<"\">" << endl;
            std::time_t t = std::time(nullptr);

            string timeStr;
//Workaround for missing std::put_time (GCC bug 54354). Affects some ubuntu 18.04 installs 
#if __GNU_C__ > 4
            timeStr=std::put_time(std::gmtime(&t), "%c %Z");
#else
            timeStr="";
#endif
            f << tabs(1) << "<created date=\"" <<  timeStr<< 
                "\" writer=\"libatomprobe " << LibVersion::getVersionStr() << "\"/>" << endl;
            f << endl;

            f << tabs(1) << "<molecules>" << endl;

            for(auto ui=0; ui<moleculeData.size(); ui++)
            {
                f << tabs(2) << "<entry>" <<  endl;
                f << tabs(3) << "<name string=\"" << ionNames[ui] << "\"/>" <<  endl;
                f << tabs(3) << "<colour hex=\"" << colours[ui].toHex() << "\"/>" << endl;

                f << tabs(3) << "<components>" << endl;
                for( auto mol : moleculeData[ui])
                {
                    f << tabs(4) << "<isotope atomicnumber=\"" << mol.atomicNumber <<
                            "\" count=\"" << mol.count << "\"/>" << endl;

                }
                f << tabs(3) << "</components>" << endl;
                f << tabs(2) << "</entry>" << endl;

            }

            f << tabs(1) << "</molecules>" << endl;
            f << tabs(1) << "<ranges>" << endl; 
            for(auto ui=0;ui<ranges.size(); ui++)
            {
                f << tabs(2) << "<range start=\"" << ranges[ui].first  << "\" end=\"" << ranges[ui].second << "\" parent=\"" << ionNames[ionIDs[ui]];

                if(rangeGrouping.size())
                    f << " group=\"" << rangeGrouping[ui] << "\" ";
                    
                f << "\"/>" << endl;
            }
            f << tabs(1) << "</ranges>" << endl; 
            f << "</multirange>" <<endl;

            break;  
        }
        default:
            return false;
    }

    return true;
}


bool MultiRange::open(const char *filename, unsigned int format) 
{
    using std::stack; 

    //Check the file opens OK without bothering with an XML
    // parse - this is a bit of a hack
    std::ifstream f(filename);
    if(!f)
        return false;
    clear();

    f.close();




    //parse XML document
    xmlDocPtr doc;
    xmlParserCtxtPtr context;
    context = xmlNewParserCtxt();

        
    if(!context)
    {
        errState=MULTIRANGE_ERR_NO_PARSER;
        return false;
    }
    
    doc = xmlCtxtReadFile(context, filename, NULL, 
                XML_PARSE_NOENT| XML_PARSE_NONET);
    
    if(!doc)
    {
        errState=MULTIRANGE_ERR_NO_XML_DOC;
        return false;
    }

    //Automatic XML free device.
    // - will automatically call free function on exit
    std::unique_ptr<_xmlDoc,void (*)(void *)> cleanup{doc,XMLFreeDoc};

    if(!context->valid)
    {
        xmlFreeParserCtxt(context);
        errState=MULTIRANGE_ERR_NO_XML_CONTEXT;
        return false;
    }
    xmlFreeParserCtxt(context);
    
    //retrieve root node    
    xmlNodePtr nodePtr = xmlDocGetRootElement(doc);
    if(!nodePtr)
    {
        errState=MULTIRANGE_ERR_NO_NODEPTR;
        return false;
    }
    
    //check value of root node
    if(xmlStrcmp(nodePtr->name, (const xmlChar *)"multirange"))
    {
        errState=MULTIRANGE_NOT_VALID_ROOTNODE;
        return false;
    }   
    //FIXME: Do something with the version data.

    if(!nodePtr->xmlChildrenNode)
    {
        errState=MULTIRANGE_ERR_NO_CONTENT;
        return false;
    }
    nodePtr = nodePtr->xmlChildrenNode;

    if(XMLHelpFwdToElem(nodePtr,"created"))
    {
        errState=MULTIRANGE_ERR_NO_TIMESTAMP;
        return false;
    }

    //FIXME: Do someting with the timestamp data

    //Read the molecule section (the species in the rangefile)
    //----
    if(XMLHelpFwdToElem(nodePtr,"molecules"))
    {
        errState=MULTIRANGE_ERR_NO_MOLECULES;
        return false;
    }

    if(!nodePtr->xmlChildrenNode)
    {
        errState=MULTIRANGE_ERR_NO_MOLECULE_ENTRIES;
        return false;
    }

    stack<xmlNodePtr> nodeStack;
    nodeStack.push(nodePtr);
    nodePtr=nodePtr->xmlChildrenNode;

    while(nodePtr)
    {
        if(XMLHelpFwdToElem(nodePtr,"entry"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_ENTRY;
            return false;
        }
        nodeStack.push(nodePtr);
        nodePtr=nodePtr->xmlChildrenNode;
        
        if(XMLHelpFwdToElem(nodePtr,"name"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_NAME;
            return false;
        }

        std::string str;
        if(XMLHelpGetProp(str,nodePtr,"string"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_NAME;
            return false;
        }
        ionNames.push_back(str);
        
        if(XMLHelpFwdToElem(nodePtr,"colour"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_COLOUR;
            return false;
        }

        if(XMLHelpGetProp(str,nodePtr,"hex"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_COLOUR;
            return false;
        }
        unsigned char v[4];
        RGBf col;
        if(!parseColString(str,v[0],v[1],v[2],v[3]))
        {
            errState=MULTIRANGE_ERR_MOL_BAD_COLOUR;
            return false;
        }
        col.red=v[0]/255.0; col.green =v[1]/255.0; col.blue=v[2]/255.0;
        colours.push_back(col);

        if(XMLHelpFwdToElem(nodePtr,"components"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
            return false;
        }
        

        set<SIMPLE_SPECIES> mol;
        nodeStack.push(nodePtr);
        nodePtr = nodePtr->xmlChildrenNode;
        while(nodePtr)
        {
            SIMPLE_SPECIES species;
            if(XMLHelpFwdToElem(nodePtr,"isotope"))
            {
                errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
                return false;
            }
            if(XMLHelpGetProp(species.atomicNumber,nodePtr,"atomicnumber"))
            {
                errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
                return false;
            }
            if(XMLHelpGetProp(species.count,nodePtr,"count"))
            {
                errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
                return false;
            }
            mol.insert(species);

            nodePtr=nodePtr->next;

        
            while(nodePtr && std::string((char*)nodePtr->name) == "text")
                nodePtr=nodePtr->next;
        }

        moleculeData.push_back(mol);
        nodeStack.pop();

        nodePtr=nodeStack.top();
        nodeStack.pop();

        nodePtr=nodePtr->next;
        while(nodePtr && std::string((char*)nodePtr->name) == "text")
            nodePtr=nodePtr->next;
    }

    nodePtr=nodeStack.top();
    nodeStack.pop();
    //----
    
    //Read the ranges section
    //----
    if(XMLHelpFwdToElem(nodePtr,"ranges"))
    {
        errState=MULTIRANGE_ERR_NO_RANGES;
        return false;
    }

    if(!nodePtr->xmlChildrenNode)
    {
        errState=MULTIRANGE_ERR_NO_RANGES;
        return false;
    }

    nodeStack.push(nodePtr);
    nodePtr=nodePtr->xmlChildrenNode;

    while(nodePtr)
    {
        if(XMLHelpFwdToElem(nodePtr,"range"))
        {
            errState=MULTIRANGE_ERR_BAD_RANGES;
            return false;
        }
        std::pair<float,float> p;
        if(XMLHelpGetProp(p.first,nodePtr,"start"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
            return false;
        }
        if(XMLHelpGetProp(p.second,nodePtr,"end"))
        {
            errState=MULTIRANGE_ERR_MOL_MISSING_COMPONENTS;
            return false;
        }

        if(p.first > p.second)
        {
            errState=MULTIRANGE_ERR_BAD_RANGE;
            return false;
        }

        ranges.push_back(p);
        //Look up the parent string from the existing list
        std::string parentStr;
        if(XMLHelpGetProp(parentStr,nodePtr,"parent"))
        {
            errState=MULTIRANGE_ERR_RANGE_NOPARENT;
            return false;
        }

        auto parentIt = std::find(ionNames.begin(),ionNames.end(),parentStr);
        if(parentIt ==ionNames.end())
        {
            errState=MULTIRANGE_ERR_RANGE_NOPARENT;
            return false;
        }

        unsigned int parentId=parentIt-ionNames.begin();
        ionIDs.push_back(parentId);


        //Look up the group identifier, if present
        std::string groupStr;
        if(!XMLHelpGetProp(groupStr,nodePtr,"group"))
        {
            unsigned int rGroup;
            if(stream_cast(rGroup,groupStr))
            {
                errState=MULTIRANGE_ERR_BAD_GROUP;
                return false;
            }
            rangeGrouping.push_back(rGroup);
        }

        //Spin past text nodes
        nodePtr=nodePtr->next;
        while(nodePtr && std::string((char*)nodePtr->name) == "text")
            nodePtr=nodePtr->next;
    }
    if(rangeGrouping.size() && (rangeGrouping.size()  !=ranges.size()))
    {
        errState=MULTIRANGE_ERR_BAD_GROUP;
        return false;
    }
    //----


    //nodePtr=nodeStack.top();
    nodeStack.pop();
    
    if(!isSelfConsistent())
    {
        errState=MULTIRANGE_ERR_SELF_INCONSISTENT;
        return false;
    }

    return true;
}

std::string MultiRange::getErrString() const
{
    const char *errMsg[] = { "Unable to create valid parser",
                 "Unable to create XML document structure",
                 "Unable to create XML document context",
                 "Unable to create node pointer",
                 "Unable to find valid root node",
                 "Unable to find content node",
                 "Unable to find timestamp node",
                 "Unable to find molecules node",
                 "Unable to find molecule entry node",
                 "Missing molecule entry",
                 "Missing molecule name",
                 "Missing molecule components",
                 "Bad range found",
                 "Bad ranges found",
                 "Missing molecule colour",
                 "Bad molecule colour",
                 "Bad group identifier",
                 "Data is not self-consistent",
                 "No parent range found"};

    ASSERT(errState <MULTIRANGE_ERR_ENUM_END);

    return errMsg[errState];
}

unsigned int MultiRange::getNumRanges() const
{
    return ranges.size();
}


set<SIMPLE_SPECIES> MultiRange::getMolecule(unsigned int ionID) const
{
    ASSERT(ionID < moleculeData.size());
    return moleculeData[ionID];
}

std::string MultiRange::getIonName(unsigned int ionID) const
{
    ASSERT(ionID < ionNames.size());
    return ionNames[ionID];
}   

unsigned int MultiRange::getNumRanges(unsigned int ionID) const
{
    ASSERT(ionID < moleculeData.size());
    unsigned int v=std::count(ionIDs.begin(),ionIDs.end(),ionID);

    return v;
}

unsigned int MultiRange::getNumIons() const
{
    return moleculeData.size();
}

unsigned int MultiRange::getIonID(unsigned int range) const
{
    ASSERT(range < ranges.size());

    return ionIDs[range];
}

std::pair<float,float> MultiRange::getRange(unsigned int rng) const
{
    ASSERT(rng<ranges.size());
    return ranges[rng];
}

std::vector<unsigned int> MultiRange::getRanges(float mass) const
{
    vector<unsigned int> tmp;
    for(auto ui=0;ui<ranges.size();ui++)
    {
        if(ranges[ui].first <= mass && ranges[ui].second > mass)
            tmp.push_back(ui);
    }

    return tmp;
}


void MultiRange::setIonID(unsigned int range, unsigned int newIonID)
{
    ASSERT(newIonID < ionIDs.size());
    ionIDs[range] = newIonID;
}


RGBf MultiRange::getColour(unsigned int ionID) const
{
    ASSERT(ionID<moleculeData.size());
    return colours[ionID];
}

void MultiRange::setColour(unsigned int ionID, const RGBf &col)
{
    ASSERT(ionID<moleculeData.size());
    colours[ionID]=col;
}

bool MultiRange::isRanged(float mass) const
{
    for(auto range : ranges)
    {
        if( mass >= range.first  &&
                mass < range.second )
            return true;
    }
    
    return false;
}

bool MultiRange::isRanged(const IonHit &ion) const
{
    return isRanged(ion.getMassToCharge());
}

void MultiRange::range(vector<IonHit> &ions) const
{
    vector<IonHit> rangedVec;

    unsigned int numIons=ions.size();
    rangedVec.reserve(numIons);

    //Range each ion
    for(auto ion :ions)
    {
        if(isRanged(ion))
            rangedVec.push_back(ion);
    }

    ions.swap(rangedVec);
}

void MultiRange::clear()
{
    ionNames.clear();
    moleculeData.clear();
    colours.clear();
    ranges.clear();
    ionIDs.clear();
    rangeGrouping.clear();

    errState=0;

}

void MultiRange::splitOverlapping(std::vector<MultiRange> &splitMR,float tolerance) const
{
    ASSERT(rangeGrouping.size());

    //Obtain a "flattened" (join ranges together in mass blocks if overlapping)
    // representation. 
    vector<FLATTENED_RANGE> flatRange;
    flattenToMassAxis(flatRange);

    //Now group the ranges according to if they sit in a given flattened range rep.
    map<size_t,size_t> groupMapping;
    vector<size_t> link(flatRange.size(),-1);
    
    //Construct graph to identify the edges that link different sections of the range grouping together 
    for(auto ui=0;ui<flatRange.size();ui++)
    {
    
        //Loop over the connected ranges, map back to this flatRange
        for(auto contained  : flatRange[ui].containedRangeIDs)
        {
            size_t gId = rangeGrouping[contained];
            if(groupMapping.find(gId) == groupMapping.end())
                groupMapping[gId] = ui;

            //If we found one, remember the first occurrence
            if(link[ui] == -1)
                link[ui] = groupMapping[gId];
            else
            {
                //always map back to the first occurrence, we will link these later
                link[ui] = std::min((size_t)groupMapping[gId],(size_t)link[ui]);
            }
        }

    }

    ASSERT(std::find(link.begin(),link.end(),-1) == link.end());

    //Remap them such that they refer to a specific group
    //  -we always mapped back to the first occurrence.
    // - this also ensures a contigous grouping
    linkIdentifiers(link);

    //Reserve space for the groups
    splitMR.resize(
        *std::max_element(link.begin(),link.end())+1);

    //examine each group, and add range and other supporting data
    // into each split group
    for(auto ui=0;ui<link.size();ui++)
    {
        size_t dest;
        dest=link[ui];
        for(auto uj=0;uj<ranges.size();uj++)
        {
            //If the range overlaps a given pair, then add it
            if(pairOverlaps(ranges[uj].first,ranges[uj].second,
                flatRange[ui].startMass,flatRange[ui].endMass))
                splitMR[dest].copyDataFromRange(*this,uj);
        }
    }


#ifdef DEBUG
    for(auto &s : splitMR)
    {
        ASSERT(s.isSelfConsistent());
    }
#endif
}

void MultiRange::flattenToMassAxis(std::vector<FLATTENED_RANGE> &ionMapping, float tolerance) const
{
    ionMapping.clear();

    if(ranges.empty())
        return;

    auto rngCopy = ranges;

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

    FLATTENED_RANGE fRng;

    fRng.startMass=rngCopy[0].first;
    fRng.endMass=rngCopy[0].second;

    auto offFirst= std::distance(ranges.begin(),
            std::find(ranges.begin(),ranges.end(),rngCopy[0]));

    fRng.containedRangeIDs.push_back(offFirst);
    fRng.containedIonIDs.push_back(ionIDs[offFirst]);

    for(auto ui=1u;ui<rngCopy.size();ui++)
    {
        if(fRng.endMass > rngCopy[ui].first-tolerance)
        {
            //Overlaps, extend fRng
            fRng.endMass = rngCopy[ui].second+tolerance;


            auto off= std::distance(ranges.begin(),
                    std::find(ranges.begin(),ranges.end(),rngCopy[ui]));
            fRng.containedRangeIDs.push_back(off);
            fRng.containedIonIDs.push_back(ionIDs[off]);
        }
        else
        {
            //does not overlap. Close out group and create new
            ionMapping.push_back(fRng);
            fRng.containedRangeIDs.clear();
            fRng.containedIonIDs.clear();

            fRng.startMass=rngCopy[ui].first;
            fRng.endMass=rngCopy[ui].second;

            //Convert back to ID of original range
            auto off= std::distance(ranges.begin(),
                    std::find(ranges.begin(),ranges.end(),rngCopy[ui]));
            fRng.containedRangeIDs.push_back(off);
            fRng.containedIonIDs.push_back(ionIDs[off]);
        }


    }
    ionMapping.push_back(fRng);


#ifdef DEBUG
    //Check masses are strictly increasing
    for(auto v : ionMapping)
    {
        ASSERT(v.startMass < v.endMass);
    }
#endif
}


void MultiRange::copyDataFromRange(const MultiRange &src,unsigned int srcRngID)
{
    ASSERT(srcRngID < src.ranges.size());

    unsigned int srcIonID;
    srcIonID=src.ionIDs[srcRngID];

    //Check to see if we have the ion name already
    auto it = std::find(ionNames.begin(),ionNames.end(),
                    src.ionNames[srcIonID]);

    //if so, just record ionID data
    if(it !=ionNames.end())
    {
        size_t destID = std::distance(ionNames.begin(),it);
        ionIDs.push_back(destID);

    }
    else
    {
        //copy other other data too 
        ionIDs.push_back(ionNames.size());
        moleculeData.push_back(src.moleculeData[srcIonID]);
        ionNames.push_back(src.ionNames[srcIonID]);
        colours.push_back(src.colours[srcIonID]);
    }
    
    ranges.push_back(src.ranges[srcRngID]);

    ASSERT(isSelfConsistent());
    rangeGrouping.clear();
}

#ifdef DEBUG

set<SIMPLE_SPECIES> createElement(const AbundanceData &abundance,
            const string &symbol, unsigned int count)
{
    SIMPLE_SPECIES ss;
    set<SIMPLE_SPECIES> s;
    auto elemIdx= abundance.symbolIndex(symbol.c_str());
    ss.atomicNumber = abundance.isotope(elemIdx,0).atomicNumber;
    ss.count=count;

    s.insert(ss);
    return s;
}

bool MultiRange::testSplit(const AbundanceData &abundance)
{
    using std::pair;

    //Create an Fe-Ni and Si-N overlap
    MultiRange m;

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

    vector<string> symbols = {"Fe","Ni","Si","N"};
    enum{ FE=0,NI,SI,N,ATOM_END};

    vector<size_t> abundanceIdx,ionID;
    //Add each atom as an ion
    for(auto ui=0;ui<ATOM_END;ui++)
    {
        auto s=createElement(abundance,symbols[ui],1);
        ionID.push_back(m.addIon(s,symbols[ui],ionCol));
    }
    abundance.getSymbolIndices(symbols,abundanceIdx);

    //Construct ranges in both the 1+ and 2+ charge state for all, but N and Si.
    // as these are not molecules, we don't need to group.
    //For N, Si, only 1+ and 2+ respectively
    const float MASS_TOL=0.1;
    vector<unsigned int> rangeGrouping;
    for(auto ui=0;ui<ATOM_END;ui++)
    {
        vector<pair<float,float> > massDist;
        abundance.generateSingleAtomDist(abundanceIdx[ui], 1,massDist);
        for(auto &md : massDist)
        {
            pair<float,float> rng;
            rng = std::make_pair(md.first-MASS_TOL,md.first+MASS_TOL);
            if(ui != N)
            {
                rng.first/=2;
                rng.second/=2;
            }


            //Add the range - we will use the mass distribution
            // to have multiRange recompute the range-> ion mapping for us
            m.addRange(rng,ui);

            //Record the atom's type for the overlap. All ranges for the same ion belong together (as we only have 1 charge state per ion)
            rangeGrouping.push_back(ui);
        }
    }

    m.setRangeGroups(rangeGrouping);

    ASSERT(m.isSelfConsistent());

    //Now attempt to split this. we shoudl end up with an Fe-Ni
    // system and an Si-N system
    vector<MultiRange> splitMulti;
    m.splitOverlapping(splitMulti);

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

    for(auto ui=0;ui<splitMulti.size();ui++)
    {
        ASSERT(splitMulti[ui].getNumIons() ==2);
        ASSERT(splitMulti[ui].getNumRanges() > 2);


        vector<string> names;
        for(auto uj=0;uj<splitMulti[ui].getNumIons();uj++)
        {
            string name = splitMulti[ui].getIonName(uj);
            ASSERT(std::find(names.begin(),names.end(),name) 
                    == names.end())
            names.push_back(name);
        }


        //If we have Fe, we should have Ni
        if(std::find(names.begin(),names.end(),symbols[FE]) !=names.end())
        {
            ASSERT(std::find(names.begin(),names.end(),symbols[NI]) != names.end())
        }
        
        //If we have Si, we should have N
        if(std::find(names.begin(),names.end(),symbols[SI]) !=names.end())
        {
            ASSERT(std::find(names.begin(),names.end(),symbols[N]) != names.end())
        }

    }


    return true;
}

bool MultiRange::test()
{
    MultiRange rng;


    set<SIMPLE_SPECIES> molecule;

    //Make FeH2+
    SIMPLE_SPECIES s;
    s.atomicNumber=26;
    s.count=1;
    molecule.insert(s);
    s.atomicNumber=1;
    s.count=2;
    molecule.insert(s);

    RGBf col;
    col.fromHex("0xFF0000");
    auto feH2Id = rng.addIon(molecule,"FeH",col);
    rng.addRange(28.9,29.3,feH2Id);

    //write test file
    rng.write("test.rngx");


    //Read test file
    {
    MultiRange rngB;

    TEST(rngB.open("test.rngx"), "Multirange load");


    TEST(rngB == rng,"write/load check");
    }

    //Add H2+ to rangefile
    s.atomicNumber=1;
    s.count=2;
    set<SIMPLE_SPECIES> newMol; 
        
    {
        newMol.insert(s);
        
        auto h2Id = rng.addIon(newMol,"H2",col);
        rng.addRange(1.7, 2.5, h2Id);
        
    
        //Set the grouping for the ranges
        vector<unsigned int> rangeGroup= { 0,1 };
        rng.setRangeGroups(rangeGroup);
    
        //Now split the problem, as we have FeH2+ and H2+,
        // we should get two muti-ranges
        vector<MultiRange> mrSplit;
        rng.splitOverlapping(mrSplit);


        TEST(mrSplit.size() == 2,"Split test");
    
        TEST(mrSplit[0].getNumIons() ==1,"split ion count");
        TEST(mrSplit[1].getNumIons() ==1,"split ion count");

        TEST(mrSplit[0].getIonName(0) != mrSplit[1].getIonName(0),"split : unique ion names");
    }



    return true;
}
#endif
}