/***************************************************************************
 *   Copyright (C) 2006 by BUI Quang Minh, Steffen Klaere, Arndt von Haeseler   *
 *   minh.bui@univie.ac.at   *
 *                                                                         *
 *   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 2 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, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "splitgraph.h"
#include <iostream>
#include <fstream>
#include <cctype>
#include <algorithm>
#include "tree/node.h"
#include "ncl/ncl.h"
#include "nclextra/myreader.h"
#include "tree/mtree.h"
#include "tree/mtreeset.h"


bool compareSplit(Split* sp1, Split* sp2) {
    if (sp1->countTaxa() != sp2->countTaxa())
        return sp1->countTaxa() < sp2->countTaxa();
    else
        return sp1->firstTaxon() < sp2->firstTaxon();
}

//#define MY_DEBUG
/********************************************************
    Defining SplitGraph methods
********************************************************/
SplitGraph::SplitGraph()
        : vector<Split*>()
{
    pda = NULL;
    taxa = NULL;
    splits = NULL;
    sets = NULL;
    trees = NULL;
    mtrees = NULL;
    areas_boundary = NULL;
}

SplitGraph::SplitGraph(Params &params) : vector<Split*>() {
    init(params);
}

void SplitGraph::createBlocks() {
    taxa = new NxsTaxaBlock();
    splits = new MSplitsBlock(this);
    pda = new MPdaBlock(this);
    sets = new MSetsBlock();
    trees = new TreesBlock(taxa);
    //mtrees = NULL;
}


void SplitGraph::init(Params &params)
{
    mtrees = NULL;
    if (params.intype == IN_NEWICK) {
        // read the input file, can contain more than 1 tree
        mtrees = new MTreeSet(params.user_file, params.is_rooted, params.tree_burnin, params.tree_max_count);
        //mtree = new MTree(params.user_file, params.is_rooted);

        if (params.is_rooted) {
            params.sub_size++;
            params.min_size++;
        }
        if (mtrees->isRooted() && params.root != NULL)
            outError(ERR_CONFLICT_ROOT);
        //SplitIntMap hash_ss;
        mtrees->convertSplits(*this, params.split_threshold, params.split_weight_summary, params.split_weight_threshold);

        if (verbose_mode >= VB_DEBUG)
            saveFileStarDot(cout);
    } else {
        createBlocks();
//        if (params.is_rooted)
//            outError(ERR_ROOT_NET);
    
         cout << "Reading input file " << params.user_file << "..." << endl;

        MyReader nexus(params.user_file);
    
        nexus.Add(taxa);
        nexus.Add(splits);
        nexus.Add(pda);
        nexus.Add(sets);
        nexus.Add(trees);

        MyToken token(nexus.inf);
        nexus.Execute(token);
        if (trees->GetNumTrees() > 0) { 
            if (getNSplits() > 0) 
                outError("Ambiguous input file, pls only specify either SPLITS block or TREES block");
            convertFromTreesBlock(params.tree_burnin, params.tree_max_count, params.split_threshold, 
                params.split_weight_summary, params.split_weight_threshold, params.tree_weight_file);
        }

    }
    
    if (verbose_mode >= VB_DEBUG)
        taxa->Report(cout);
    //splits->Report(cout);
    //reportConflict(cout);
    if (params.pdtaxa_file != NULL) {
        if (sets->getNSets() > 0)
            outError("Taxa sets were already specified in the input file");
        cout << "Reading taxa sets in file " << params.pdtaxa_file << "..." << endl;
    
        bool nexus_formated = (detectInputFile(params.pdtaxa_file) == IN_NEXUS);
        if (nexus_formated) {
            MyReader nexus(params.pdtaxa_file);
            nexus.Add(sets);
            MyToken token(nexus.inf);
            nexus.Execute(token);
        } else {
            readTaxaSets(params.pdtaxa_file, sets);
        }
        if (sets->getNSets() == 0)
            outError("No taxa sets found");
    }

    areas_boundary = NULL;
    if (params.areas_boundary_file) {
        if (sets->getNSets() == 0) outError("No taxon sets defined yet");
        areas_boundary = new double [sets->getNSets() * sets->getNSets()];
        cout << "Reading sets relation file " << params.areas_boundary_file << "..." << endl;
        readAreasBoundary(params.areas_boundary_file, sets, areas_boundary);
    }

    if (verbose_mode >= VB_DEBUG && sets->getNSets() > 0)
        sets->Report(cout);

    if (sets->getNSets() > 0 && taxa->GetNumTaxonLabels() == 0) {
        AddTaxaFromSets();
    }
    if (taxa->GetNumTaxonLabels() == 0)
        outError("No taxa found");
    if (getNSplits() == 0) {
        //outError(ERR_NO_SPLITS);
        createStarTree();
    }
    cout << getNTaxa()-params.is_rooted <<
        " taxa and " << getNSplits()-params.is_rooted << " splits." << endl;

}


void SplitGraph::saveCheckpoint() {
    if (empty()) return;
    int ntax = getNTaxa();
//    checkpoint->startStruct("S");
    CKP_SAVE(ntax);
    int nsplits = size();
    CKP_SAVE(nsplits);
    checkpoint->startList(size());
    for (iterator it = begin(); it != end(); it++) {
        checkpoint->addListElement();
        stringstream ss;
        ss << (*it)->getWeight();
        for (int i = 0; i < ntax; i++)
            if ((*it)->containTaxon(i))
                ss << " " << i;
        checkpoint->put("", ss.str());
    }
    checkpoint->endList();
//    checkpoint->endStruct();
    CheckpointFactory::saveCheckpoint();
}

void SplitGraph::restoreCheckpoint() {
    int ntax, nsplits;
    CheckpointFactory::restoreCheckpoint();
//    checkpoint->startStruct("S");

    if (!CKP_RESTORE(ntax)) return;
    CKP_RESTORE(nsplits);
    checkpoint->startList(nsplits);
    for (int split = 0; split < nsplits; split++) {
        checkpoint->addListElement();
        string str;
        bool found = checkpoint->getString("", str);
        ASSERT(found);
        stringstream ss(str);
        double weight;
        ss >> weight;
        Split *sp = new Split(ntax, weight);
        for (int i = 0; i < ntax; i++) {
            int tax;
            if (ss >> tax) {
                sp->addTaxon(tax);
            } else
                break;
        }
        push_back(sp);
    }
    checkpoint->endList();
//    checkpoint->endStruct();
}

int SplitGraph::getNTrivialSplits() {
    int count = 0;
    for (iterator it = begin(); it != end(); it++)
        if ((*it)->trivial() >= 0)
            count++;
    return count;
}


void SplitGraph::createStarTree() {
    cout << "No splits found, creating a star tree with branch length of 1..." << endl;
    int ntaxa = taxa->GetNumTaxonLabels();
    for (int i = 0; i < ntaxa; i++) {
        Split *sp = new Split(ntaxa, 1.0);
        sp->addTaxon(i);
        push_back(sp);
    }
    cout << "NOTE: subsequent PD will correspond to species richness." << endl;
}


void SplitGraph::AddTaxaFromSets() {
    cout << "Taking taxa from SETS block..." << endl;
    for (int i = 0; i < sets->getNSets(); i++)
        for(vector<string>::iterator it = sets->getSet(i)->taxlist.begin(); 
            it != sets->getSet(i)->taxlist.end(); it++) 
            if (!taxa->IsAlreadyDefined(NxsString(it->c_str()))) {
                taxa->AddTaxonLabel(NxsString(it->c_str()));
            }    
}

void SplitGraph::freeMem() {
    for (reverse_iterator it = rbegin(); it != rend(); it++) {
        //(*it)->report(cout);
        delete *it;
    }
    clear();
    if (areas_boundary) delete areas_boundary;
    if (trees) delete trees;
    if (sets) delete sets;
    if (pda) delete pda;
    if (splits) delete splits;
    if (taxa) delete taxa;
    if (mtrees) delete mtrees;
}

SplitGraph::~SplitGraph()
{
    freeMem();
}


void SplitGraph::convertFromTreesBlock(int burnin, int max_count, double split_threshold, 
    int split_weight_summary, double weight_threshold, const char *tree_weight_file) {
    cout << trees->GetNumTrees() << " tree(s) loaded" << endl;
    if (burnin >= trees->GetNumTrees())
        outError("Burnin value is too large");
    if (burnin > 0)
    cout << burnin << " beginning tree(s) discarded" << endl;
    mtrees = new MTreeSet();
    
    for (int i = burnin; i < trees->GetNumTrees() && (i < burnin+max_count); i++) {
        stringstream strs(trees->GetTranslatedTreeDescription(i), ios::in | ios::out | ios::app);
        strs << ";";
        MTree *tree = mtrees->newTree();
        bool myrooted = trees->IsRootedTree(i);
        tree->readTree(strs, myrooted);
        mtrees->push_back(tree);
        mtrees->tree_weights.push_back(1);
    }
    mtrees->checkConsistency();
    //SplitIntMap hash_ss;
    
    if (tree_weight_file) 
        readIntVector(tree_weight_file, burnin, max_count, mtrees->tree_weights);
/*    else if (!weights)
        tree_weights.resize(size(), 1);*/

    if (mtrees->size() != mtrees->tree_weights.size())
        outError("Tree file and tree weight file have different number of entries");    
    mtrees->convertSplits(*this, split_threshold, split_weight_summary, weight_threshold);
}



void SplitGraph::report(ostream &out)
{

    out << endl;
    out << "Split network contains ";

    if (size() == 0)
    {
        out << "no split" << endl;
    }
    else if (size() == 1)
        out << "one split" << endl;
    else
        out << size() << " splits" << endl;

    if (size() == 0)
        return;

    sort(begin(), end(), compareSplit);
    int k = 0;
    for (iterator it = begin(); it != end(); it++,k++)
    {
        out << '\t' << (k+1) << '\t';
        (*it)->report(out);
    }
}

void SplitGraph::reportConflict(ostream &out)
{
    int k = 0;
    out << "Compatible splits: " << endl;
    for (iterator i = begin(); i != end(); i++, k++)
    {
        out << (k+1) << '\t';
        int k2 = 1;
        for (iterator j = begin(); j != end(); j++, k2++)
            if ( j != i && (*i)->compatible(*(*j)))
            {
                out << k2 << " ";
            }
        out << endl;
    }
}

/**
    calculate sum of weights of preserved splits in the taxa_set
    @param taxa_set a set of taxa
*/
double SplitGraph::calcWeight(Split &taxa_set)
{
    double sum = 0.0;
    for (iterator it = begin(); it != end(); it++)
        if ((*it)->preserved(taxa_set))
            sum += (*it)->getWeight();
    return sum;
}

int SplitGraph::countSplits(Split &taxa_set)
{
    int cnt = 0;
    for (iterator it = begin(); it != end(); it++)
        if ((*it)->preserved(taxa_set))
            cnt++;
    return cnt;
}

int SplitGraph::countInternalSplits(Split &taxa_set)
{
    int cnt = 0;
    for (iterator it = begin(); it != end(); it++)
        if ((*it)->trivial() < 0 && (*it)->preserved(taxa_set))
            cnt++;
    return cnt;
}

/**
    calculate sum of weights of all splits
*/
double SplitGraph::calcWeight() {
    double sum = 0.0;
    for (iterator it = begin(); it != end(); it++)
        sum += (*it)->weight;
    return sum;
}

double SplitGraph::calcTrivialWeight() {
    double sum = 0.0;
    for (iterator it = begin(); it != end(); it++)
        if ((*it)->trivial() >= 0)
            sum += (*it)->weight;
    return sum;
}

double SplitGraph::maxWeight() {
    double m = -1e6;
    for (iterator it = begin(); it != end(); it++)
        if (m < (*it)->weight) m = (*it)->weight;
    return m;
}

void SplitGraph::generateTaxaSet(char *filename, int size, int overlap, int times) {
    ofstream out(filename);
    if (!out.is_open())
        outError(ERR_WRITE_OUTPUT, filename);
    ASSERT(overlap <= size);
    int total = 2*size - overlap;
    int ntaxa = getNTaxa();
    for (int cnt = 0; cnt < times; cnt++) {
        // generate random taxon index 
        IntVector ranvec;
        BoolVector occur(ntaxa, false);
        int i;
        for (i = 0; i < total; i++) {
            int rnum;
            do { rnum = random_int(ntaxa); } while (occur[rnum]);
            ranvec.push_back(rnum);
            occur[rnum] = true;
        }
        // now write the first set
        out << size << endl;
        for (i = 0; i < size; i++) 
            out << taxa->GetTaxonLabel(ranvec[i]) << endl;
        out << endl;
        // now write the second set
        out << size << endl;
        for (i = size-overlap; i < total; i++) 
            out << taxa->GetTaxonLabel(ranvec[i]) << endl;
        out << endl;
    }
    out.close();
}

void SplitGraph::calcDistance(char *filename) {
    ofstream out(filename);
    if (!out.is_open())
        outError(ERR_WRITE_OUTPUT, filename);
    mmatrix(double) dist;
    int i, j;    
    calcDistance(dist);

    int ntaxa = getNTaxa();

    // now write the distances in phylip .dist format
    out << ntaxa << endl;
    
    for (i = 0; i < ntaxa; i++) {
        out << taxa->GetTaxonLabel(i) << "   ";
        for (j = 0; j < ntaxa; j++) {
            out << dist[i][j] << "  ";
        }
        out << endl;
    }
    out.close();
}

void SplitGraph::calcDistance(mmatrix(double) &dist) {
    int ntaxa = getNTaxa();
    iterator it;
    vector<int> vi, vj;
    vector<int>::iterator i, j;

    dist.resize(ntaxa);
    for (mmatrix(double)::iterator di = dist.begin(); di != dist.end(); di++)
        (*di).resize(ntaxa, 0);

    for (it = begin(); it != end(); it++) {
        (*it)->getTaxaList(vi, vj);
        for (i = vi.begin(); i != vi.end(); i++)
            for (j = vj.begin(); j < vj.end(); j++) {
                dist[*i][*j] += (*it)->weight;
                dist[*j][*i] += (*it)->weight;
            }
    }

}


void SplitGraph::calcDistance(mmatrix(double) &dist, vector<int> &taxa_order) {
    int ntaxa = getNTaxa();
    int i, j;

    mmatrix(double) my_dist;
    calcDistance(my_dist);
    dist.resize(ntaxa);
    for (i = 0; i < ntaxa; i++) {
        dist[i].resize(ntaxa);
        for (j = 0; j < ntaxa; j++)
            dist[i][j] = my_dist[taxa_order[i]][taxa_order[j]];
    }
}

bool SplitGraph::checkCircular(mmatrix(double) &dist) {
    return true;
    int ntaxa = getNTaxa();
    Split taxa_set(ntaxa, 0.0);
    for (int i = 0; i < ntaxa-2; i++)
        for (int j = i+1; j < ntaxa-1; j++)
            for (int k = j+1; k < ntaxa; k++) {
                taxa_set.addTaxon(i);
                taxa_set.addTaxon(j);
                taxa_set.addTaxon(k);
                taxa_set.weight = calcWeight(taxa_set);
                if (fabs(2 * taxa_set.weight - (dist[i][j] + dist[i][k] + dist[j][k])) > 0.0001) {
                    cout << "Taxa " << i << " " << j << " " << k;
                    cout << " do not satisfy circular equation!" << endl;
                    cout << "Weight = " << taxa_set.weight << endl;
                    cout << "Sum dist/2 = " << (dist[i][j] + dist[i][k] + dist[j][k]) / 2.0 << endl;
                    cout << "dist = " << dist[i][j] << " " << dist[i][k] << " "
                         << dist[j][k] << endl;
                    return false;
                }
                taxa_set.removeTaxon(i);
                taxa_set.removeTaxon(j);
                taxa_set.removeTaxon(k);
            }
    return true;
}

void SplitGraph::generateCircular(Params &params) {
    int i, j;
    int ntaxa = params.sub_size;
    int num_splits = (params.num_splits > 0) ? params.num_splits : 3 * ntaxa;
    if (num_splits < ntaxa) 
        outError(ERR_FEW_SPLITS); 

    taxa = new NxsTaxaBlock();
    splits = new MSplitsBlock(this);

    double threshold = (ntaxa > 3) ? (double)(num_splits - ntaxa) / (ntaxa*(ntaxa-3)/2) : 0.0;

    // insert all trivial splits
    for (i = 0; i < ntaxa; i++) {
        double weight = randomLen(params);
        Split *sp = new Split(ntaxa, weight);
        sp->addTaxon(i);
        push_back(sp);
        ostringstream str;
        str << "T" << (i+1);
        taxa->AddTaxonLabel(NxsString(str.str().c_str()));
        splits->cycle.push_back(i);
    }

    // randomly insert internal splits
    for (i = 0; i < ntaxa-2 && getNSplits() < num_splits; i++)
        for (j = i+1; j < ntaxa && j < ntaxa-3+i; j++) {
            double choice = random_double();
            if (choice > threshold) continue;
            double weight = randomLen(params);
            Split *sp = new Split(ntaxa, weight);
            for (int k = i; k <= j; k++)
                sp->addTaxon(k);
            push_back(sp);
            if (getNSplits() >= num_splits) break;
        }

    ofstream out(params.user_file);
    if (!out.is_open()) {
        outError(ERR_WRITE_OUTPUT, params.user_file);
    }

    saveFileNexus(out);
    out.close();
} 

void SplitGraph::saveFileNexus(ostream &out, bool omit_trivial) {
    int ntaxa = getNTaxa();
    int i;
    out << "#nexus" << endl << endl;
    out << "BEGIN Taxa;" << endl;
    out << "DIMENSIONS ntax=" << ntaxa << ";" << endl;
    out << "TAXLABELS" << endl;
    for (i = 0; i < ntaxa; i++)
        out << "[" << i+1 << "] '" << taxa->GetTaxonLabel(i) << "'" << endl;
    out << ";" << endl << "END; [Taxa]" << endl << endl;
    out << "BEGIN Splits;" << endl;
    out << "DIMENSIONS ntax=" << ntaxa << " nsplits=" << ((omit_trivial) ? getNSplits() - getNTrivialSplits() : getNSplits()) << ";" << endl;
    out << "FORMAT labels=no weights=yes confidences=no intervals=no;" << endl;
    if (isCircular()) {
        out << "CYCLE";
        for (i = 0; i < ntaxa; i++) 
            out << " " << splits->cycle[i] + 1;
        out << ";" << endl;
    }
    out << "MATRIX" << endl;
    int near_zeros = 0;
    int zeros = 0;
    for (iterator it = begin(); it != end(); it++) {
        if (omit_trivial && (*it)->trivial() >= 0) continue;
        if ((*it)->weight == 0.0) zeros ++;
        if ((*it)->weight <= 1e-6) near_zeros ++;
        out << "\t" << (*it)->weight << "\t";
        for (i = 0; i < ntaxa; i++)
            if ((*it)->containTaxon(i))
                out << " " << i+1;
        out << "," << endl;
    }
    out << ";" << endl << "END; [Splits]" << endl << endl;
    if (near_zeros) {
        //outWarning("Some nearly-zero split weights observed");
        //cout << zeros << " zero-weights and " << near_zeros << " near zero weights!" << endl;
    }
}

void SplitGraph::saveFileStarDot(ostream &out, bool omit_trivial) {
    int ntaxa = getNTaxa();
    int i;
    for (iterator it = begin(); it != end(); it++) {
        if (omit_trivial && (*it)->trivial() >= 0) continue;
        bool swap_code = !(*it)->containTaxon(0);
        if (swap_code) {
            for (i = 0; i < ntaxa; i++)
                out << (((*it)->containTaxon(i)) ? '.' : '*');
        } else {
            for (i = 0; i < ntaxa; i++)
                out << (((*it)->containTaxon(i)) ? '*' : '.');
        }
        out << "\t" << (*it)->weight << endl;
    }
}

void SplitGraph::saveFile(const char* out_file, InputType file_format, bool omit_trivial) {
    try {
        ofstream out;
        out.exceptions(ios::failbit | ios::badbit);
        out.open(out_file);
        if (file_format == IN_NEXUS) 
            saveFileNexus(out, omit_trivial);
        else
            saveFileStarDot(out, omit_trivial);
        out.close();
    } catch (ios::failure) {
        outError(ERR_WRITE_OUTPUT, out_file);
    }
}

void SplitGraph::scaleWeight(double norm, bool make_int, int precision) {
    for (iterator itg = begin(); itg != end(); itg ++ )
        if (make_int)
            (*itg)->setWeight( round((*itg)->getWeight()*norm) );
        else if (precision < 0)
            (*itg)->setWeight( (*itg)->getWeight()*norm);
        else 
            (*itg)->setWeight( round((*itg)->getWeight()*norm*pow((double)10.0,precision))/pow((double)10.0,precision));
}
// TODO Implement a more efficient function using Hash Table
bool SplitGraph::containSplit(Split &sp) {
    Split invert_sp(sp);
    invert_sp.invert();
    for (iterator it = begin(); it != end(); it++)
        if ((*(*it)) == sp || (*(*it)) == invert_sp)
            return true;
    return false;
}

double SplitGraph::computeBoundary(Split &area) {
    if (!areas_boundary) return 0.0;
    int nareas = sets->getNSets();
    double boundary = 0.0;
    for (int i = 0; i < nareas; i++) 
    if (area.containTaxon(i)) {
        boundary += areas_boundary[i*nareas+i];
        for (int j = i+1; j < nareas; j++) 
            if (area.containTaxon(j))
                boundary -= 2.0 * areas_boundary[i*nareas+j];
    }
    return boundary;
}

bool SplitGraph::compatible(Split *sp) {
    for (iterator it = begin(); it != end(); it++)
        if (!(*it)->compatible(*sp))
            return false;
    return true;
}

void SplitGraph::findMaxCompatibleSplits(SplitGraph &maxsg) {

    // maximum number of compatible splits = 2n-3!
    int max_splits = getNTaxa() * 2 - 3;

    // myset will be sorted by weight in descending order
    SplitSet myset;
    myset.insert(myset.end(), begin(), end());
    sort(myset.begin(), myset.end(), splitweightcmp);

    // now build the spset
    if (!maxsg.taxa)
        maxsg.taxa = new NxsTaxaBlock();
    if (!maxsg.splits)
        maxsg.splits = new MSplitsBlock(&maxsg);
    if (!maxsg.pda)
        maxsg.pda = new MPdaBlock(&maxsg);

    for (int i = 0; i < getNTaxa(); i++)
        maxsg.taxa->AddTaxonLabel(taxa->GetTaxonLabel(i));
    
    // make the cycle
    maxsg.splits->cycle = splits->cycle;
    // make the splits

    for (SplitSet::iterator it = myset.begin(); it != myset.end(); it++) 
        if (maxsg.compatible(*it)){
            maxsg.push_back(new Split(*(*it)));
            //(*it)->report(cout);
            if (maxsg.size() >= max_splits)
                break;
        }
    myset.clear();
}

bool SplitGraph::isWeaklyCompatible() {
    if (getNSplits() < 3) return true;
    for (iterator it1 = begin(); it1+2 != end(); it1++)
        for (iterator it2 = it1+1; it2+1 != end(); it2++)
            for (iterator it3 = it2+1; it3 != end(); it3++) {
                Split sp1(*(*it1));
                Split sp2(*(*it2));
                Split sp3(*(*it3));
                Split sp(sp1);
                sp *= sp2;
                sp *= sp3;
                if (sp.isEmpty()) continue;
                sp1.invert();
                sp2.invert();
                sp = sp1;
                sp *= sp2;
                sp *= sp3;
                if (sp.isEmpty()) continue;
                sp2.invert();
                sp3.invert();
                sp = sp1;
                sp *= sp2;
                sp *= sp3;
                if (sp.isEmpty()) continue;
                sp1.invert();
                sp2.invert();
                sp = sp1;
                sp *= sp2;
                sp *= sp3;
                if (sp.isEmpty()) continue;
                return false;
            }
    return true;
}


void SplitGraph::getTaxaName(vector<string> &taxname) {
    taxname.clear();
    for (int i = 0; i < getNTaxa(); i++)
        taxname.push_back(taxa->GetTaxonLabel(i));
}

int SplitGraph::findLeafName(string &name) {
    for (int i = 0; i < getNTaxa(); i++)
        if (taxa->GetTaxonLabel(i) == name)
            return i;
    return -1;
}

int SplitGraph::removeTrivialSplits() {
    int removed = 0;
    for (iterator itg = begin(); itg != end(); )  {
        if ((*itg)->trivial() >= 0) {
            removed++;
            delete (*itg);
            (*itg) = back();
            pop_back(); 
        } else itg++;
    }
    return removed;
}