xref: /dports/biology/lamarc/lamarc-2.1.8/src/tools/rangex.cpp

// $Id: rangex.cpp,v 1.23 2011/12/29 23:44:19 ewalkup Exp $

/*
  Copyright 2002  Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein

  This software is distributed free of charge for non-commercial use
  and is copyrighted.  Of course, we do not guarantee that the software
  works, and are not responsible for any damage you may cause or have.
*/

#include <cassert>

#include "local_build.h"

#include "rangex.h"
#include "stringx.h"
#include "registry.h"

using std::string;
using std::make_pair;

//------------------------------------------------------------------------------------

string ToString(rangepair rpair)
{
    //Note!  This prints stuff out in 'user units' instead of 'internal units'.
    // In other words, if the range is <0,34>, this is an 'open-ended' interval
    // and the actual active sites are 0:33
    //
    //Furthermore, if the user is expecting to have no zeroes in the output, we
    // need to convert the 0 to -1, and display "-1:33".
    rpair.second--;
    rpair = ToNoZeroesIfNeeded(rpair);
    string retval = ToString(rpair.first);
    if (rpair.second > rpair.first)
    {
        retval += ":" + ToString(rpair.second);
    }
    return retval;
}

//------------------------------------------------------------------------------------

string ToString(rangeset rset)
{
    rangeset::iterator rpair = rset.begin();
    if (rpair == rset.end())
    {
        return "(none)";
    }
    string retval = ToString(*rpair);
    rpair++;
    for ( ; rpair != rset.end(); rpair++)
    {
        retval += ", " + ToString(*rpair);
    }
    return retval;
}

//------------------------------------------------------------------------------------
// Returns a rangeset as a set containing a single pair, that defined by its two arguments.

rangeset MakeRangeset(long int low, long int high)
{
    rangeset retset;
    retset.insert(make_pair(low, high));

    return retset;
}

//------------------------------------------------------------------------------------
// Returns a rangeset as the set union of the first arg (as a rangepair) plus the second arg (as a rangeset).

rangeset AddPairToRange(const rangepair& addpart, const rangeset& rset)
{
    rangeset retset = rset;

    rangepair low  = make_pair(addpart.first, addpart.first);
    rangepair high = make_pair(addpart.second, addpart.second);

    long int newlow = low.first;
    long int newhigh = high.second;

    rangesetiter early = retset.lower_bound(low);
    if (early != retset.begin())
    {
        early--;
        if (early->second >= low.first)
        {
            //'low' falls within early's interval
            newlow = early->first;
        }
    }

    rangesetiter late = retset.upper_bound(high);
    //We need to increment this late.first == high.first
    if (late != retset.end())
    {
        if (late->first == high.first)
        {
            late++;
        }
    }
    if (late != retset.begin())
    {
        late--;
        if (late->second > high.first)
        {
            //'high' falls within late's interval
            newhigh = late->second;
        }
    }

    early = retset.lower_bound(make_pair(newlow, newlow+1));
    late  = retset.upper_bound(make_pair(newhigh-1, newhigh));

    retset.erase(early, late);
    retset.insert(make_pair(newlow, newhigh));

    return retset;
}

//------------------------------------------------------------------------------------
// Returns a rangeset as the set difference of the second arg (as a rangeset) minus the first arg (as a rangepair).

rangeset RemovePairFromRange(const rangepair& removepart, const rangeset& rset)
{
    rangeset retset;
    for (rangeset::iterator range=rset.begin(); range != rset.end(); range++)
    {
        if ((range->first >= removepart.first) && (range->second <= removepart.second))
        {
            //Don't add it.
        }
        else
        {
            if ((range->first < removepart.first) && (range->second > removepart.second))
            {
                //Add two outside halves.
                retset.insert(make_pair(range->first, removepart.first));
                retset.insert(make_pair(removepart.second, range->second));
            }
            else if ((range->second > removepart.first) && (range->second <= removepart.second))
            {
                //Add only the first half.
                retset.insert(make_pair(range->first, removepart.first));
            }
            else if ((range->first >= removepart.first) && (range->first <= removepart.second))
            {
                //Add only the second half.
                retset.insert(make_pair(removepart.second, range->second));
            }
            else
            {
                //Add the whole thing.
                retset.insert(*range);
            }
        }
    }
    return retset;
}

//------------------------------------------------------------------------------------
// Returns a rangeset as the set difference of the second arg (as a rangeset) minus the first arg (as a rangeset).

rangeset RemoveRangeFromRange(const rangeset& removerange, const rangeset& rset)
{
    rangeset retset = rset;

    for (rangeset::iterator removepair = removerange.begin(); removepair != removerange.end(); removepair++)
    {
        retset = RemovePairFromRange(*removepair, retset);
    }

    return retset;
}

//------------------------------------------------------------------------------------
// LS NOTE:
// These functions originally written for range.cpp, but ended up being needed elsewhere.
//
// Note that 'Union' (used to be 'OR' but I got confused too often) actually does the exact same thing that
// 'AddRangeToRange' used to, but much faster, so I took out AddRangeToRange entirely.  Perhaps there's a
// similar way to speed up RemoveRangeFromRange?  We'll see if it shows up in the profiler.

rangeset Union(const rangeset& set1, const rangeset& set2)
{
    if (set1.empty()) return set2;
    if (set2.empty()) return set1;

    rangeset mergeset;
    merge(set1.begin(),set1.end(),set2.begin(),set2.end(),
          inserter(mergeset,mergeset.begin()));

    rangeset newset;
    rangeset::iterator rit;

    for(rit = mergeset.begin(); rit != mergeset.end(); ++rit)
    {
        ORAppend(*rit,newset);
    }

    return newset;
} // OR

//------------------------------------------------------------------------------------
// Helper function for 'Union'

void ORAppend(rangepair newrange, rangeset& oldranges)
{
    if (oldranges.empty())
    {
        oldranges.insert(newrange);
        return;
    }

    rangeset::iterator last = --oldranges.end();

    assert(newrange.first >= last->first); // Expect sequential, sorted input.

    if (newrange.second <= last->second) return; // New is contained in old.

    if (newrange.first > last->second)  // New is after old.
    {
        oldranges.insert(oldranges.end(),newrange);
        return;
    }

    newrange.first = last->first;       // New starts within old and extends past it.
    oldranges.erase(last);
    oldranges.insert(oldranges.end(),newrange);

    return;
} // ORAppend

//------------------------------------------------------------------------------------
// Used to be 'AND' until I got confused too often.

rangeset Intersection(const rangeset& mother, const rangeset& father)
{
    rangeset::const_iterator m = mother.begin();
    rangeset::const_iterator f = father.begin();

    rangeset result;
    rangepair newpair;

    if (mother.empty() || father.empty()) return result;

    while (true)
    {
        newpair.first = std::max((*m).first, (*f).first);
        newpair.second = std::min((*m).second, (*f).second);

        if (newpair.first < newpair.second)
            result.insert(result.end(),newpair);

        if ((*m).second < (*f).second) ++m;
        else ++f;

        if (m == mother.end() || f == father.end()) return result;
    }
} // Intersection

//------------------------------------------------------------------------------------
// Counts the sites (actually, any object; could be links) in a rangeset.

long int CountSites(rangeset rset)
{
    long int retval(0);

    for (rangeset::iterator rpair = rset.begin(); rpair != rset.end(); rpair++)
    {
        retval += rpair->second - rpair->first;
    }

    return retval;
}

//------------------------------------------------------------------------------------
// Measures the distance, in links, between the two ends of the sites in its argument.

long CountLinks(const rangeset& sites)
{
    if (sites.empty()) return 0L;

    long start = sites.begin()->first;
    long end = sites.rbegin()->second;
    assert(start < end);

    return end - start - 1;
} // CountLinks

//------------------------------------------------------------------------------------

rangeset TargetLinksFrom(const rangeset& targetsites)
{
    rangeset targetlinks;

    // If the targetsites are empty we really do want to return an empty rangeset.
    if (!targetsites.empty())
    {
        long firstlink(targetsites.begin()->first);
        long lastlink(targetsites.rbegin()->second-1);

        if (firstlink != lastlink)
        {
            rangepair newlinks(firstlink,lastlink);
            targetlinks = AddPairToRange(newlinks,targetlinks);
        }
    }

    return targetlinks;
} // TargetLinksFrom

//------------------------------------------------------------------------------------

rangeset AddToRangeset(const rangeset& rset, long offset)
{
    rangeset::const_iterator range;
    rangeset newRangeSet;

    for(range = rset.begin(); range != rset.end(); ++range)
    {
        long newFirst = range->first + offset;
        long newSecond = range->second + offset;

        newRangeSet.insert(make_pair(newFirst,newSecond));
    }
    return newRangeSet;
}

rangeset SubtractFromRangeset(const rangeset& rset, long offset)
{
    long invOffset = 0 - offset;
    return AddToRangeset(rset,invOffset);
}

//------------------------------------------------------------------------------------

bool IsInRangeset(const rangeset& targetset, long target)
{
    if (targetset.empty())
        return false;

    rangeset::const_iterator range;

    for(range = targetset.begin(); range != targetset.end(); ++range)
    {
        if (range->second <= target) continue;
        if (range->first <= target) return true;
        break;
    }

    return false;
} // IsInRangeset

//------------------------------------------------------------------------------------
// Used in Range::AreLowSitesOnInactiveBranch_Rg to diagnose whether sites of interest
// occur on both sides of the cutpoint.

bool Surrounds(const rangeset& targetsites, long targetlink)
{
    if (targetsites.empty()) return false;

    return (targetlink >= targetsites.begin()->first && targetlink < targetsites.rbegin()->second - 1);
} // Surrounds

//------------------------------------------------------------------------------------
//These functions are to be used when converting site labels for display to
// the user, and from *menu* input from the user.  It is assumed that in the
// XML, the 'user input' has already been converted to the 'sequential' version.

long int ToSequentialIfNeeded(long int input)
{
    //Note:  input might equal zero if it was the upper end of a rangepair,
    // since rangepairs are open-ended by default.  In other words, if the user
    // types in "-20:-1" in the menu, the first thing that happens is that these
    // numbers are converted to the pair <-20, 0>.  To convert these values to
    // the 'sequential' numbering system, this needs to be converted to the
    // pair <-19, 1>.

    if (registry.GetConvertOutputToEliminateZeroes() && (input <= 0))
    {
        input++;
    }

    return input;
}

//------------------------------------------------------------------------------------

long int ToNoZeroesIfNeeded(long int input)
{
    if (registry.GetConvertOutputToEliminateZeroes() && (input <= 0))
    {
        input--;
    }

    return input;
}

//------------------------------------------------------------------------------------

rangepair ToSequentialIfNeeded(rangepair input)
{
    if (registry.GetConvertOutputToEliminateZeroes())
    {
        input.first = ToSequentialIfNeeded(input.first);
        input.second = ToSequentialIfNeeded(input.second);
    }

    return input;
}

//------------------------------------------------------------------------------------

rangepair ToNoZeroesIfNeeded(rangepair input)
{
    if (registry.GetConvertOutputToEliminateZeroes())
    {
        input.first = ToNoZeroesIfNeeded(input.first);
        input.second = ToNoZeroesIfNeeded(input.second);
    }

    return input;
}

//____________________________________________________________________________________