#include <iostream>
#include <vector>
#include <utility>
#include <string>
#include <boost/format.hpp>
#include <boost/algorithm/algorithm.hpp>
#include <cmath>
#include <sstream>
#include <numeric>

// Import NUPACK energy computation
extern "C" {
#include "thermo/core.h"
}

// Import ViennaRNA energy computation
extern "C"{
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/data_structures.h"
#include "ViennaRNA/model.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/utils.h"
#include "ViennaRNA/read_epars.h"
#include "ViennaRNA/constraints.h"
#include "ViennaRNA/eval.h"
#include "ViennaRNA/cofold.h"
#include "ViennaRNA/file_formats.h"
#include "ViennaRNA/subopt.h"
#include "RNAeval_cmdl.h"
}

#include "structure.h"
#include "Motifs/hairpinloop.h"
#include "Motifs/interloop.h"
#include "Motifs/multiloop.h"
#include "Motifs/helix.h"
#include "Motifs/pseudoknot.h"
#include "utils.h"

Structure::Structure()
{
}

Structure::Structure (const int rna, const std::string &seq,
                      const std::vector < std::pair < unsigned int, unsigned int > > &listBp,
                      const int ct, const int nbCt, const int id,
                      unsigned int energyModel, float energy,
                      const std::vector < float > &probingData,
                      float lowerThresProbing,
                      float upperThresProbing):
    rna_(rna), seq_(seq), listBp_(listBp), ct_(ct), nbCt_(nbCt), id_(id)
{
    motifDetection();

    if (energyModel == 0)
        obj1_ = energy;
    else if (energyModel == 1)
        computeEnergyVienna();
    else if (energyModel == 2)
        computeEnergyNUPACK();

    computeProbing(probingData, lowerThresProbing, upperThresProbing);
}

Structure::Structure(const Structure& that)
{
    rna_ = that.rna_;
    seq_ = that.seq_;
    obj1_ = that.obj1_;
    listBp_ = that.listBp_;
    ct_ = that.ct_;
    ic_ = that.ic_;
    nbCt_ = that.nbCt_;
    id_ = that.id_;
    probing_ = that.probing_;
    motifs_ = std::vector < Motif * > (that.motifs_.size());
    for(size_t i = 0, size = that.motifs_.size(); i != size; i++)
    {
        if (Helix * S = dynamic_cast < Helix * > (that.motifs_[i]))
        {
            motifs_[i] = new Helix (*S);
        }
        else if (Interloop * I = dynamic_cast < Interloop * > (that.motifs_[i]))
        {
            motifs_[i] = new Interloop (*I);
        }
        else if (Hairpinloop * H = dynamic_cast < Hairpinloop * > (that.motifs_[i]))
        {
            motifs_[i] = new Hairpinloop (*H);
        }
        else if (Pseudoknot * P = dynamic_cast < Pseudoknot * > (that.motifs_[i]))
        {
            motifs_[i] = new Pseudoknot (*P);
        }
        else if (Multiloop * M = dynamic_cast < Multiloop * > (that.motifs_[i]))
        {
            motifs_[i] = new Multiloop (*M);
        }
    }
}

Structure::~Structure()
{
    for(size_t i = 0, size = motifs_.size(); i != size; i++)
        delete motifs_[i];
}

int Structure::get_rna_() const
{
    return rna_;
}

std::string Structure::get_seq_() const
{
    return seq_;
}

float Structure::get_obj1_() const {
    return obj1_;
}

std::vector < std::pair < unsigned int, unsigned int > > Structure::get_listBp_() const
{
    return listBp_;
}

int Structure::get_ct_() const
{
    return ct_;
}

float Structure::get_ic_() const
{
    return ic_;
}

int Structure::get_nbCt_() const
{
    return nbCt_;
}

int Structure::get_id_() const
{
    return id_;
}

float Structure::get_probing_() const
{
    return probing_;
}

std::vector < Motif * > Structure::get_motifs_() const
{
    return motifs_;
}

void Structure::set_obj1_(float obj1){
    obj1_ = obj1;
}

void Structure::set_id_(int id)
{
    id_ = id;
}

void Structure::set_nbCt_(int nbCt)
{
    nbCt_ = nbCt;
}

void Structure::set_ct_(int ct)
{
    ct_ = ct;
}

void Structure::set_ic_(float ic)
{
    ic_ = ic;
}

void Structure::set_probing_(float probing)
{
    probing_ = probing;
}

std::string Structure::convToDP() const
{
    return seq_ + "\n" + convToDPonly();
}

std::string Structure::convToDP(char c1, std::vector<char> forbid, bool first) const
{
    return seq_ + "\n" + convToDPonly(c1, forbid, first);
}

// From RNAstructure
std::string Structure::convToDPonly() const
{
    /*std::string res(uint(seq_.size()), '.');

    char  c1 = 'A';
    for (int i = 0; i < int(listBp_.size()); i++)
    {
        if(i > 0 and (listBp_[i].first - listBp_[i-1].first > 1 or
                listBp_[i].second - listBp_[i-1].second > 1) )
            c1++;
        res[listBp_[i].first] = c1;
        res[listBp_[i].second] = c1;

    }
    return res;*/

    std::string res = "";

    std::vector < uint > ranks(uint(seq_.size())+1);
    const char*const brackets = "()<>{}[]AaBbCcDd";
    const uint MAX_LEVEL = strlen(brackets)/2; // Highest pseudoknot level allowed.

    // Initialize basepr
    std::vector < uint > basepr = std::vector < uint > (uint(seq_.size()+1), 0);
    for (size_t i = 0, size = listBp_.size(); i != size; i++) {
        basepr[listBp_[i].first+1] = listBp_[i].second+1;
        basepr[listBp_[i].second+1] = listBp_[i].first+1;
    }

    // Write out the STRUCTURE LINE
    getPseudoknotRanks(ranks, basepr);

    for (uint i=1;i<=uint(seq_.size());i++) {
            // Normally pairs are encoded with "()" but psedoknots (and higher-order knots) are encoded with alternate brackets.
            int level = std::min(ranks[i], MAX_LEVEL) - 1;  // Get the knot level. 0 = no pseudoknot, 1 = first-order pseudoknots, 2 = second-order knots (i.e. which crossed other pseudoknots)
            if (basepr[i]>i) { res += brackets[2*level]; } // The opening bracket, e.g. '('
            else if (basepr[i]==0) { res += "."; }
            else { res += brackets[2*level+1]; } // The closing bracket, e.g. ')'
    }
    return res;
}

// From RNAstructure
void Structure::getPseudoknotRanks(std::vector < uint > &ranks, std::vector < uint > basepr) const {

    // `list` first contains ALL basepairs.  We then call findPseudoknots to fill it with just the pseudoknots.
    //  The pseudoknots found in each round become the "normal pairs" for each subsequent round of finding pseudoknots.
    std::vector < int > list(basepr.size());
    std::copy(basepr.begin(), basepr.end(), list.begin());

    for(unsigned int i=0;i<ranks.size();i++)
            ranks[i]=list[i]==0 ? 0 : 1;
    // Results now has a 0 for all positions without basepairs and a 1 at each position with a pair (regardless of whether it is crossing or not)

    while(hasPseudoknots(list)) {
            findPseudoknots(list, &list);  // Determine which are the cossing (pseudoknot) bonds.
            // We've already accounted for the optimal set of non-crossing bonds.
            // so now just increment results[i] for each bond in the set of crossing bonds.
            for(unsigned int i=0;i<ranks.size();i++)
                    if (list[i]!=0) ranks[i]++;
    }
}

// From RNAstructure
bool Structure::hasPseudoknots(const std::vector < int > &pairs) const {
    const int length = pairs.size();
    IntervalStack stack(std::min(8,length/4));
    const int FIRST_NUC = 1; // In RNAstructure the convention is that index 1 is the first nucleotide (instead of 0)
    stack.push(FIRST_NUC, length-1);
    while (stack.pop()) {
        int k; // Represents the 3' end of the basepair for which the 5' end is at stack.i
        // increase stack.i until it is either == stack.j or stack.i points to a basepair.
        while(stack.i <= stack.j && 0==(k = pairs[stack.i])) // find out if there is a bond at position i (and make sure that i is the 5' end).
            stack.i++;

        // If no basepair was found between i and j, pop that interval. we are done looking inside it.
        if (stack.i > stack.j)
            continue; // pop next item from stack

        // whenever a bond is encountered, two new intervals are pushed -- the interval INSIDE the bond and the interval AFTER it.
        // The 5' end is always encountered first (even for crossing bonds) because the lower interval is processed before the upper interval.
        // So stack.i will ALWAYS be the 5' position and bases[stack.i] will ALWAYS be the 3' position.
        if (k < stack.i) std::cerr << "Programming logic error. 5' end encountered in ::hasPseudoknots" << std::endl;

        if (k > stack.j) return true; // If the 3' end is outside of the interval, this represents a crossing bond. So return true.
        if (k + 1 <= stack.j)  // If there are any nucleotides between the 3'-end and j, push that interval
        stack.push(k+1, stack.j);
        if (stack.i + 1 <= k-1) // If there are any nucleotides between the 5' and 3' end, push that interval
        stack.push(stack.i+1, k-1);
    }
    return false;
}

// From RNAstructure
void Structure::findPseudoknots(const std::vector < int > &pairs,
                                std::vector < int > *knotted,
                                std::vector < int > *optimal) const {

    const unsigned int length = pairs.size();
    const unsigned int FIRST_NUC = 1; // In RNAstructure the convention is that index 1 is the first nucleotide (instead of 0)
    if (length==0) return;

    if (optimal==NULL && knotted==NULL) return; // At least one result vector should be passed in or there is no point in doing the calculation.

    // Create an upper-triangular matrix for both outer and trace.
    // Since the first base is at position 1, the 0th column and row are never used, so we can create smaller
    //   arrays and do pointer math to make outer[1][1] actually point to the 0th element (i.e. first memory location).
    // Additionally note that for trace, we can make the array even smaller because i>j for all trace[i][j] so
    //   we can delete the center diagonal. So trace[1][2] points to the 0th element. (i.e. rows are shifted by 1, columns by 2)
    // In `outer`: rows go from i=1 to N; cols from i to N  (where N=length-1 is the index of the last base)
    short** outer=(new short*[length-FIRST_NUC])-FIRST_NUC;  for(int i=FIRST_NUC; i<length; ++i) outer[i]=(new short[length-i])-i;
    // In `trace`: // rows go from i=1 to N-1; cols from i+1 to N
    // The size of trace is length-2, but it is but it is only shifted by 1. This is because it is one row shorter than outer.
    bool**  trace=(new bool*[length-FIRST_NUC-1])-FIRST_NUC; for(int i=FIRST_NUC; i<length-1; ++i) trace[i]=(new bool[length-i-1])-i-1;

    // Sanity test for triangular array and pointer math: (these will NOT necessarily result in access violations or segfaults, even if the pointer math is wrong)
    // DEBUG: for(int i=FIRST_NUC; i<length; i++){for(int j=i; j<length; j++){outer[i][j]=i*j; cout<<outer[i][j]<<" ";}cout<<endl;}
    // DEBUG: for(int i=FIRST_NUC; i<length-1; i++){for(int j=i+1; j<length; j++){outer[i][j]=i*j; cout<<outer[i][j]<<" ";}cout<<endl;}

    for (int i=FIRST_NUC; i<length; i++) outer[i][i]=0; // initialize diagonal of `outer` to 0 (this ensures proper start conditions)

    // The center diagonal of `outer` has been set to 0.
    // Each iteration of `n` fills the next diagonal to the (upper) right.
    // The first iteration represents all intervals [i, i+1]. The second is [i, i+2] etc. so the size of the interval
    // increases with each iteration. (Note: "interval" means a segment of the nucleobase sequence. The "size" is simply `j - i`)
    // Iteration 1:  [1,2], [2,3], [3,4]...[N-2,N-1] [N-1,N]  (size: 1)       (where N = length-1 is the last valid index in `bases`.)
    // Iteration 2:  [1,3], [2,4], [3,5]...[N-2, N]           (size: 2)
    // Ieration N-1: [1,N-1],[2,N]                            (size: N-1)
    // Iteration N:  [1,N]                                    (size: N)
    // The value at outer[i, j] represents the maximum number of non-crossing bonds **fully contained** within the sequence-interval [i, j].
    // (Fully contained means  bond.i >=i and bond.j <= j for every bond)
    // trace[i, j] is true if the bond at bases[i] should be included (i.e. it maximizes the number of bonds in the interval [i,j])

    for (int n = 1; n < length; n++) {
        for (int i = FIRST_NUC; i < length - n; i++) {
            int j = i + n;
            // First assume that either no bond starts here. So the total non-crossing bonds within the interval [i, j] is unchanged compared
            // to the (smaller) interval [i+1, j] which would have been set in the previous iteration.
            outer[i][j] = outer[i + 1][j];
            trace[i][j] = false; // Default value assume there is no bond here.
            int k = pairs[i];   // k == 0 if there is no basepair.
            // if k > i, then i is the 5' end, and k is the 3' end of the pair.
            // if k < i then k is the 5' end and i is the 3' end.
            if (k != 0 && k > i && k <= j) { // i.e. If this is the 3' end of a basepair and the 5' end (k) is inside the current interval [i, j]
                int tmp = 1; // Add 1 for THIS bond

                if (i + 1 <= k - 1)
                    tmp += outer[i + 1][k - 1]; // Add all the bonds **fully contained** by this one (i.e. which have i' > i and k' < k)

                if (k + 1 <= j)
                    tmp += outer[k + 1][j];  // Add in all the bonds AFTER this one but still **fully contained** by the interval [k+1, j]

                if (tmp >= outer[i][j]) {     // If inclusion of this bond maximizes the total number of non-crossing bonds in the interval [i, j]
                    outer[i][j] = (short) tmp;
                    trace[i][j] = true;
                }
            }
        }
    }
    // DEBUG:
    //cout<<endl<<"pairs: ";  for(int i=FIRST_NUC; i<length; i++) cout<<pairs[i]<<"\t";  cout<<endl;
    //for(int i=FIRST_NUC; i<length; i++) {
    //      for(int j=FIRST_NUC; j<=i; j++) cout<<"\t";
    //      for(int j=i+1; j<length; j++)
    //              cout<<outer[i][j]<<(trace[i][j]?"**\t":"\t");
    //      cout<<endl;     }
    // We no longer need the "outer" data, so we reuse it here:
    // Note that the user may pass in the same vector for pairs AND either optimal or knotted,
    // so we have to make a copy of the pairing data in case the pairs array is modified when setting values in `optimal` or `knotted`.
    short *results = outer[FIRST_NUC]; for(int i = FIRST_NUC; i < length; i++) results[i]=pairs[i]; // Fill in results with the "pairs" information.
    // Backtrace
    IntervalStack stack(std::min(uint(8), length / 4));
    stack.push(1, length - 1);
    while (stack.pop()) {
        //DEBUG: cout << "stack: " << stack.i << "," << stack.j << endl;
        while(stack.i < stack.j && !trace[stack.i][stack.j]) // Trace[i][j] is true if the pair with its 5' end at i is in the optimal set.
            stack.i++;
        if (stack.i >= stack.j)
            continue; // Pop next item from stack
        //DEBUG: cout << "found: " << stack.i << "," << stack.j << ": " << trace[stack.i][stack.j] << endl;

        // here k != -1 and i < j
        int k = pairs[stack.i];
        results[stack.i] = -k; // make the 5' end negative to indicate it's in the optimal set.
        results[k] = -pairs[k]; // make the 3' end negative
        if (stack.i + 1 < k-1)
            stack.push(stack.i+1, k-1);
        if (k + 1 < stack.j)
            stack.push(k+1, stack.j);
    }
    //DEBUG: cout<<endl<<"Results: ";for(int i=FIRST_NUC; i<length; i++) cout<<results[i]<<"   ";cout<<endl;

    if (optimal != NULL) {
        if (optimal->size() < length) optimal->resize(length);
        for (int i = FIRST_NUC; i < length; i++)
            (*optimal)[i]=results[i]<0 ? -results[i] : 0; // Fill in pair information if this pair IS in the optimal results list.
    }
    if (knotted != NULL) {
        if (knotted->size() < length) knotted->resize(length);
        for (int i = FIRST_NUC; i < length; i++)
            (*knotted)[i]=results[i]>0 ? results[i] : 0; // Fill in pair information if this pair is NOT in the optimal results list.
    }

    // Note: This should perform the inverse of the pointer math done at allocation.
    for(int i=FIRST_NUC;i<length;i++)   delete[] (outer[i]+i);   delete[] (outer+FIRST_NUC);
    for(int i=FIRST_NUC;i<length-1;i++) delete[] (trace[i]+i+1); delete[] (trace+FIRST_NUC);
    // FYI - do not delete `results` because it is just a pointer to outer[0].
}

std::string Structure::convToDPonly(char c1, std::vector<char> forbid, bool first) const
{
    std::string res(int(seq_.size()), '.');

        bool F = false;

        for (int i = 0; i < int(listBp_.size()); i++)
        {
            if(i > 0 and (listBp_[i].first - listBp_[i-1].first > 1 or
                       listBp_[i].second - listBp_[i-1].second > 1) )
            {
                c1++;
                while(std::find(forbid.begin(), forbid.end(), c1) != forbid.end())
                {
                    if (c1 == 127 and first)
                    {
                        c1 = 33;
                        first = false;
                    }
                    else if (c1 == 127 and not first)
                    {
                        F = true;
                        break;
                    }
                    else
                    {
                        c1++;
                    }
                }
            }
            res[listBp_[i].first] = c1;
            res[listBp_[i].second] = c1;
            forbid.push_back(c1);
            if(F)
                break;
        }
        if(F)
            res = "Not dot parenthesis display available.";

        return res;
}

std::string Structure::to_string() const
{
    return convToDP() + "\t" + boost::str(boost::format("%.2f") % obj1_);
}

std::string Structure::to_Json() const
{
    return "{\"seq\":\"" + convToDP() + "\", \"obj1\":\" " + boost::str(boost::format("%.2f") % obj1_) + "\"},";
}

/*int Structure::find_bp(std::pair< unsigned int, unsigned int > p) const
{
    std::vector< std::pair< unsigned int, unsigned int > >::const_iterator it = std::find(listBp_.begin(), listBp_.end(), p);
    int r = -1;
    if(it != listBp_.end())
        r = std::distance(listBp_.begin(), it);
    return r;
}*/

Structure& Structure::operator=(const Structure& that)
{
    rna_ = that.rna_;
    seq_ = that.seq_;
    obj1_ = that.obj1_;
    listBp_ = that.listBp_;
    ct_ = that.ct_;
    nbCt_ = that.nbCt_;
    id_ = that.id_;
    probing_ = that.probing_;
    std::vector < Motif * > local_motifs_ = std::vector < Motif * > (that.motifs_.size());
    for(size_t i = 0, size = that.motifs_.size(); i != size; i++)
    {
        if (Helix * S = dynamic_cast < Helix * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Helix (*S);
        }
        else if (Interloop * I = dynamic_cast < Interloop * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Interloop (*I);
        }
        else if (Hairpinloop * H = dynamic_cast < Hairpinloop * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Hairpinloop (*H);
        }
        else if (Pseudoknot * P = dynamic_cast < Pseudoknot * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Pseudoknot (*P);
        }
        else if (Multiloop * M = dynamic_cast < Multiloop * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Multiloop (*M);
        }
    }
    for(size_t i = 0, size = motifs_.size(); i != size; i++)
        delete motifs_[i];
    motifs_ = local_motifs_;

    return *this;
}

bool Structure::operator<(const Structure& s)
{
    return obj1_ < s.obj1_;
}

bool Structure::operator<(const Structure& s) const
{
    return obj1_ < s.obj1_;
}

bool Structure::operator>(const Structure& s)
{
    return obj1_ > s.obj1_;
}

bool Structure::operator>(const Structure& s) const
{
    return obj1_ > s.obj1_;
}

bool Structure::operator<=(const Structure& s)
{
    return obj1_ <= s.obj1_;
}

bool Structure::operator<=(const Structure& s) const
{
    return obj1_ <= s.obj1_;
}

bool Structure::operator>=(const Structure& s)
{
    return obj1_ >= s.obj1_;
}

bool Structure::operator>=(const Structure& s) const
{
    return obj1_ >= s.obj1_;
}

bool Structure::operator==(const Structure& s)
{
    /*bool res = true;
    if( ! ( rna_ == s.rna_ and obj1_ == s.obj1_ and convToDP().compare(s.convToDP()) == 0 ))
    {
        res = false;
    }*/
    return rna_ == s.rna_ and s.id_ == id_;
}

bool Structure::operator==(const Structure& s) const
{
    /*bool res = true;
    if( ! ( rna_ == s.rna_ and obj1_ == s.obj1_ and convToDP().compare(s.convToDP()) == 0 ))
    {
        res = false;
    }*/
    return rna_ == s.rna_ and s.id_ == id_;
}

void Structure::makeListBp(std::string structureDP, std::vector < std::pair < unsigned int, unsigned int > > &listBp)
{
    std::vector < unsigned int > openingPar;
    std::vector < unsigned int > openingBra;
    for (int i = 0; i < int(structureDP.size()); i++)
    {
        if(structureDP[i] == '(' )
        {
            openingPar.push_back(i);
        }
        else if(structureDP[i] == ')')
        {
            listBp.push_back(std::make_pair(openingPar.back(), i));
            openingPar.pop_back();
        }
        else if(structureDP[i] == '[' )
        {
            openingBra.push_back(i);
        }
        else if(structureDP[i] == ']')
        {
            listBp.push_back(std::make_pair(openingBra.back(), i));
            openingBra.pop_back();
        }
    }
    std::sort(listBp.begin(), listBp.end());
}

bool Structure::checkStructure(std::string structure)
{
    bool res = true;
    size_t oPar = std::count(structure.begin(), structure.end(), '(');
    size_t cPar = std::count(structure.begin(), structure.end(), ')');
    size_t oBra = std::count(structure.begin(), structure.end(), '[');
    size_t pBra = std::count(structure.begin(), structure.end(), ']');
    if(oPar != cPar or oBra != pBra) {
        res = false;
        throw (std::string("Invalid number of base pairs."));
        return res;
    }
    return res;
}

// Method to print the tree
void Structure::print(helix *tmp)
{
    /*std::cout << "root (" << tmp->pos_start.first << "," << tmp->pos_start.second
              << ") (" << tmp->pos_end.first << "," << tmp->pos_end.second << ") pseudo=" << tmp->pseudo << " crossingHelixs= " << tmp->crossingHelixs.size() << std::endl;*/
    if (!tmp->crossingHelixs.empty())
    {
        /*std::cout << "crossingHelixs ";
        for(size_t j = 0, size2 = tmp->crossingHelixs.size(); j != size2; j++)
            std::cout << " (" << tmp->crossingHelixs[j]->pos_start.first << "," <<  tmp->crossingHelixs[j]->pos_start.second << ") (" << tmp->crossingHelixs[j]->pos_end.first << "," <<  tmp->crossingHelixs[j]->pos_end.second;
        std::cout << std::endl;*/
    }
    if (!tmp->child.empty())
    {
        for(size_t j = 0, size2 = tmp->child.size(); j != size2; j++)
        {
            /*std::cout << "child de root [" << tmp->pos_start.first << "," << tmp->pos_start.second
                      << "][" << tmp->pos_end.first << "," << tmp->pos_end.second << "] (" << tmp->child[j]->pos_start.first << "," << tmp->child[j]->pos_start.second
                      << ") (" << tmp->child[j]->pos_end.first << "," << tmp->child[j]->pos_end.second << ") pseudo=" << tmp->child[j]->pseudo << " crossingHelixs= " << tmp->child[j]->crossingHelixs.size() << std::endl;*/
            print(tmp->child[j]);
        }
    }
}

// Method to check the children for the multiloop motif
bool Structure::checkChildrenMulti(helix *child, helix *parent)
{
    bool crossingHelix = true;
    if(child->pseudo)
    {
        crossingHelix = false;
        size_t i = 0, size2 = child->crossingHelixs.size();
        int counter = 0;
        while(!crossingHelix and i != size2 and counter < 2)
        {
            if(!child->crossingHelixs[i]->pseudo and child->crossingHelixs[i]->parent == parent)
            {
                counter++;
                if (counter == 2)
                    crossingHelix = true;
            }
            i++;
        }
    }
    if(crossingHelix and !child->child.empty())
    {
        bool checkChild = true;
        size_t i = 0, size2 = child->child.size();
        while(checkChild and i != size2)
        {
            checkChild = checkChildrenMulti(child->child[i], parent);
            i++;
        }
        crossingHelix = crossingHelix and checkChild;
    }

    return crossingHelix;
}

// Method to check the children for the Internal loop motif
bool Structure::checkChildrenIntern(helix *child, helix *parent)
{
    bool internal = true;
    if (!child->child.empty())
        for(size_t i = 0, size = child->child.size(); i != size and internal; i++)
            if(child->child[i]->pseudo)
                if(std::find(child->child[i]->crossingHelixs.begin(), child->child[i]->crossingHelixs.end(), child) == child->child[i]->crossingHelixs.end()
                        or std::find(child->child[i]->crossingHelixs.begin(), child->child[i]->crossingHelixs.end(), parent) == child->child[i]->crossingHelixs.end())
                    internal = false;

    return internal;
}

// Method to generate the arborescence of all the motifs with the Object corresponding to the different motifs
std::vector < Motif * > Structure::readTree(helix *tmp)
{
    std::vector < Motif * >  res;

    if(!tmp->child.empty())
    {
        if(tmp->pos_start.first < tmp->pos_end.first) {
            try {
                res.push_back(new Helix(rna_, rna_, tmp->pos_start.first, tmp->pos_end.first,
                                        tmp->pos_start.second, tmp->pos_end.second,
                                        std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                        seq_, seq_));
            } catch (std::string e) {
                std::cout << e << std::endl;
            }
        }


        if(!tmp->pseudo)
        {
            bool crossingHelix = true;
            std::vector < std::pair < int, int > > closingBp;

            for(unsigned int i(0); i < tmp->child.size(); i++)
            {
                bool bulge_i = (std::abs(std::get<0>(tmp->pos_end) - std::get<0>(tmp->child[i]->pos_start)) > 1) ;
                bool bulge_j =(std::abs(std::get<1>(tmp->pos_end) - std::get<1>(tmp->child[i]->pos_start)) > 1);

                // Internal loops
                if (((bulge_i and !bulge_j) or (!bulge_i and bulge_j)) and tmp->child.size() == 1 and !tmp->child[i]->pseudo)
                    if (checkChildrenIntern(tmp->child[i], tmp)) {
                        try {
                            res.push_back(new Interloop(rna_, seq_, tmp->pos_end.first, tmp->child[i]->pos_start.first,
                                                        tmp->pos_end.second, tmp->child[i]->pos_start.second,
                                                        std::abs(tmp->pos_end.first - tmp->child[i]->pos_start.first),
                                                        std::abs(tmp->pos_end.second - tmp->child[i]->pos_start.second)));
                        } catch (std::string e) {
                            std::cout << e << std::endl;
                        }
                    }
                if(bulge_i && bulge_j && tmp->child.size()==1 && !tmp->child[i]->pseudo)
                    if (checkChildrenIntern(tmp->child[i], tmp)) {
                        try {
                            res.push_back(new Interloop(rna_, seq_, tmp->pos_end.first, tmp->child[i]->pos_start.first,
                                                        tmp->pos_end.second, tmp->child[i]->pos_start.second,
                                                        std::abs(tmp->pos_end.first - tmp->child[i]->pos_start.first),
                                                        std::abs(tmp->pos_end.second - tmp->child[i]->pos_start.second)));
                        } catch (std::string e) {
                            std::cout << e << std::endl;
                        }
                    }

                // Multiloop
                if(tmp->child.size() > 1 && tmp->crossingHelixs.size() == 0 and crossingHelix)
                {
                        crossingHelix = checkChildrenMulti(tmp->child[i], tmp);
                        closingBp.push_back(std::make_pair(tmp->child[i]->pos_start.first, tmp->child[i]->pos_start.second));
                }

                std::vector < Motif * > tmp_res = readTree(tmp->child[i]);
                res.insert(res.end(),tmp_res.begin(),tmp_res.end());
            }

            if(tmp->child.size() > 1 and tmp->crossingHelixs.size() == 0 and crossingHelix) {
                try {
                    res.push_back(new Multiloop(rna_, seq_, closingBp));
                } catch (std::string e) {
                    std::cout << e << std::endl;
                }
            }
        }
        else
        {
            for(size_t i = 0, size = tmp->crossingHelixs.size(); i != size; i++) {
                try {
                    res.push_back(new Pseudoknot(rna_, rna_,
                                         Helix(rna_, rna_, tmp->pos_start.first, tmp->pos_end.first,
                                               tmp->pos_start.second, tmp->pos_end.second,
                                               std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                               seq_, seq_),
                                         Helix(rna_, rna_,
                                               tmp->crossingHelixs[i]->pos_start.first, tmp->crossingHelixs[i]->pos_end.first,
                                               tmp->crossingHelixs[i]->pos_start.second, tmp->crossingHelixs[i]->pos_end.second,
                                               std::abs(tmp->crossingHelixs[i]->pos_end.first - tmp->crossingHelixs[i]->pos_start.first),
                                               seq_, seq_)));
                } catch (std::string e) {
                    std::cout << e << std::endl;
                }
            }

            for(unsigned int i(0); i < tmp->child.size(); i++)
            {
                if(tmp->pos_end.first < tmp->child[i]->pos_start.first
                        and tmp->child[i]->pos_start.first < tmp->child[i]->pos_start.second
                        and tmp->child[i]->pos_start.second < tmp->pos_start.second)
                {
                    bool bulge_i = (std::abs(std::get<0>(tmp->pos_end) - std::get<0>(tmp->child[i]->pos_start)) > 1) ;
                    bool bulge_j =(std::abs(std::get<1>(tmp->pos_end) - std::get<1>(tmp->child[i]->pos_start)) > 1);

                    if(((bulge_i and !bulge_j) or (!bulge_i and bulge_j)) and tmp->child.size()==1 and tmp->child[i]->pseudo)
                        if (checkChildrenIntern(tmp->child[i], tmp)) {
                            try {
                                res.push_back(new Interloop(rna_, seq_, tmp->pos_end.first, tmp->child[i]->pos_start.first,
                                                        tmp->pos_end.second, tmp->child[i]->pos_start.second,
                                                        std::abs(tmp->pos_end.first - tmp->child[i]->pos_start.first),
                                                        std::abs(tmp->pos_end.second - tmp->child[i]->pos_start.second)));
                            } catch (std::string e) {
                                std::cout << e << std::endl;
                            }
                        }

                    if(bulge_i and bulge_j and tmp->child.size()==1 and tmp->child[i]->pseudo)
                        if (checkChildrenIntern(tmp->child[i], tmp)) {
                            try {
                                res.push_back(new Interloop(rna_, seq_, tmp->pos_end.first, tmp->child[i]->pos_start.first,
                                                            tmp->pos_end.second, tmp->child[i]->pos_start.second,
                                                            std::abs(tmp->pos_end.first - tmp->child[i]->pos_start.first),
                                                            std::abs(tmp->pos_end.second - tmp->child[i]->pos_start.second)));
                            } catch (std::string e) {
                                std::cout << e << std::endl;
                            }
                        }
                }
                std::vector < Motif * > tmp_res = readTree(tmp->child[i]);
                res.insert(res.end(),tmp_res.begin(),tmp_res.end());

            }
        }

        return res;
    }
    else if (tmp->crossingHelixs.empty())
    {
        if(tmp->pos_start.first < tmp->pos_end.first) {
            try {
                res.push_back(new Helix(rna_, rna_, tmp->pos_start.first, tmp->pos_end.first,
                                    tmp->pos_start.second, tmp->pos_end.second,
                                    std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                    seq_, seq_));
            } catch (std::string e) {
                std::cout << e << std::endl;
            }
        }

        if(tmp->pseudo)
            for(size_t i = 0, size = tmp->crossingHelixs.size(); i != size; i++) {
                try {
                    res.push_back(new Pseudoknot(rna_, rna_,
                                         Helix(rna_, rna_,
                                               tmp->pos_start.first, tmp->pos_end.first,
                                               tmp->pos_end.second, tmp->pos_start.second,
                                               std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                               seq_, seq_),
                                         Helix(rna_, rna_,
                                               tmp->crossingHelixs[i]->pos_start.first, tmp->crossingHelixs[i]->pos_end.first,
                                               tmp->crossingHelixs[i]->pos_end.second, tmp->crossingHelixs[i]->pos_start.second,
                                               std::abs(tmp->crossingHelixs[i]->pos_end.first - tmp->crossingHelixs[i]->pos_start.first),
                                               seq_, seq_)));
                } catch (std::string e) {
                    std::cout << e << std::endl;
                }
            }
        else
            res.push_back(new Hairpinloop(rna_, seq_, tmp->pos_end.first, tmp->pos_end.second, std::abs(tmp->pos_end.second - tmp->pos_end.first)));

        return res;
    }
    else
    {
        if(tmp->pos_start.first < tmp->pos_end.first) {
            try {
                res.push_back(new Helix(rna_, rna_,
                                    tmp->pos_start.first, tmp->pos_end.first,
                                    tmp->pos_start.second, tmp->pos_end.second,
                                    std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                    seq_, seq_));
            } catch (std::string e) {
                std::cout << e << std::endl;
            }
        }
        if(tmp->pseudo)
            for(size_t i = 0, size = tmp->crossingHelixs.size(); i != size; i++) {
                try {
                    res.push_back(new Pseudoknot(rna_, rna_,
                                         Helix(rna_, rna_,
                                               tmp->pos_start.first, tmp->pos_end.first,
                                               tmp->pos_start.second, tmp->pos_end.second,
                                               std::abs(tmp->pos_end.first - tmp->pos_start.first),
                                               seq_, seq_),
                                         Helix(rna_, rna_,
                                               tmp->crossingHelixs[i]->pos_start.first, tmp->crossingHelixs[i]->pos_end.first,
                                               tmp->crossingHelixs[i]->pos_start.second, tmp->crossingHelixs[i]->pos_end.second,
                                               std::abs(tmp->crossingHelixs[i]->pos_end.first - tmp->crossingHelixs[i]->pos_start.first),
                                               seq_, seq_)));
                } catch (std::string e) {
                    std::cout << e << std::endl;
                }
            }
        return res;
    }
}

void Structure::freeHelixSS(helix *val)
{
    for(auto child: val->child){
        freeHelixSS(child);
    }
    val->child.clear();
    delete val;
}

void Structure::motifDetection()
{
    std::vector< std::pair< unsigned int, unsigned int > > parenthesisClose = listBp_;
    std::vector<helix*> roots;

    std::sort(parenthesisClose.begin(), parenthesisClose.end());

    // Tree generation of the helix positions (only with the this->parenthesis vector)
    helix *root = new helix();
    helix *root2;

    if(!parenthesisClose.empty()){
        root->pos_start=std::make_pair(std::get<0>(parenthesisClose[0]), std::get<1>(parenthesisClose[0]));
    }else{
        root->pos_start=std::make_pair(0,0);
    }
    root->pseudo = false;
    roots.push_back(root);

    size_t i, size;
    std::vector < helix * > openedPseudo;
    std::vector < helix * > crossingHelixsTmp;

    if(!parenthesisClose.empty()){
        for(unsigned int val(0); val<(parenthesisClose.size()-1); val++){
            bool diff = (std::get<0>(parenthesisClose[val+1]) - std::get<0>(parenthesisClose[val]) > 1
                         or std::get<1>(parenthesisClose[val+1]) - std::get<1>(parenthesisClose[val]) > 1);

            if(diff)
            {
                if (!root->pseudo and std::get<1>(parenthesisClose[val]) > std::get<1>(parenthesisClose[val+1])) // Normal parent and child
                {
                    // same child = same branch so same root
                    root->pos_end=std::make_pair(std::get<0>(parenthesisClose[val]), std::get<1>(parenthesisClose[val]));
                    for (i = 0, size = openedPseudo.size(); i != size and root != NULL; i++)
                    {
                        if(((parenthesisClose[val].first < openedPseudo[i]->pos_start.second
                             and openedPseudo[i]->pos_start.second < parenthesisClose[val].second)
                                or (parenthesisClose[val].first < openedPseudo[i]->pos_start.first
                                    and openedPseudo[i]->pos_start.first < parenthesisClose[val].second))
                                and std::find(openedPseudo[i]->crossingHelixs.begin(), openedPseudo[i]->crossingHelixs.end(), root)
                                == openedPseudo[i]->crossingHelixs.end())
                        {
                            openedPseudo[i]->crossingHelixs.push_back(root);
                            root->crossingHelixs.push_back(openedPseudo[i]);
                        }
                        else if (openedPseudo[i]->pos_start.second < parenthesisClose[val].first)
                        {
                            openedPseudo.erase(openedPseudo.begin() + i);
                            size--;
                            i--;
                        }
                    }

                    helix *tmp= new helix ();
                    tmp->pos_start=std::make_pair(std::get<0>(parenthesisClose[val+1]), std::get<1>(parenthesisClose[val+1]));
                    tmp->parent=root;
                    tmp->pseudo = false;
                    root->child.push_back(tmp);
                    root=tmp;
                }
                else if ( (!root->pseudo and std::get<1>(parenthesisClose[val]) < std::get<1>(parenthesisClose[val+1]) and std::get<0>(parenthesisClose[val+1]) < std::get<1>(parenthesisClose[val])) // Normal parent and pseudoknot child
                          or (root->pseudo and std::get<1>(parenthesisClose[val]) > std::get<1>(parenthesisClose[val+1]) and std::get<0>(parenthesisClose[val+1]) < root->parent->pos_start.second and std::get<1>(parenthesisClose[val+1]) > root->parent->pos_start.second)) // Pseudoknot parent and child
                {
                    // Same child = same branch so same root
                    root->pos_end=std::make_pair(std::get<0>(parenthesisClose[val]), std::get<1>(parenthesisClose[val]));
                    for (i = 0, size = openedPseudo.size(); i != size and root != NULL; i++)
                    {
                        if(((parenthesisClose[val].first < openedPseudo[i]->pos_start.second
                             and openedPseudo[i]->pos_start.second < parenthesisClose[val].second)
                                or (parenthesisClose[val].first < openedPseudo[i]->pos_start.first
                                    and openedPseudo[i]->pos_start.first < parenthesisClose[val].second))
                            and std::find(openedPseudo[i]->crossingHelixs.begin(), openedPseudo[i]->crossingHelixs.end(), root)
                                == openedPseudo[i]->crossingHelixs.end())
                        {
                            openedPseudo[i]->crossingHelixs.push_back(root);
                            root->crossingHelixs.push_back(openedPseudo[i]);
                        }
                        else if (openedPseudo[i]->pos_start.second < parenthesisClose[val].first)
                        {
                            openedPseudo.erase(openedPseudo.begin() + i);
                            size--;
                            i--;
                        }
                    }

                    helix *tmp= new helix ();
                    tmp->pos_start=std::make_pair(std::get<0>(parenthesisClose[val+1]), std::get<1>(parenthesisClose[val+1]));
                    tmp->parent=root;
                    tmp->pseudo = true;
                    root2 = root;
                    while (root2 != NULL)
                    {
                        if(root2->pos_start.first < tmp->pos_start.first
                                and tmp->pos_start.first < root2->pos_start.second
                                and root2->pos_start.second < tmp->pos_start.second)
                        {
                            tmp->crossingHelixs.push_back(root2);
                            root2->crossingHelixs.push_back(tmp);
                        }
                        root2=root2->parent;
                    }

                    openedPseudo.push_back(tmp);
                    root->child.push_back(tmp);
                    root=tmp;
                }
                else
                {
                    // New child
                    root->pos_end=std::make_pair(std::get<0>(parenthesisClose[val]), std::get<1>(parenthesisClose[val]));
                    for (i = 0, size = openedPseudo.size(); i != size and root != NULL; i++)
                    {
                        if(((parenthesisClose[val].first < openedPseudo[i]->pos_start.second
                             and openedPseudo[i]->pos_start.second < parenthesisClose[val].second)
                                or (parenthesisClose[val].first < openedPseudo[i]->pos_start.first
                                    and openedPseudo[i]->pos_start.first < parenthesisClose[val].second
                                    and root != openedPseudo[i]->parent))
                            and std::find(openedPseudo[i]->crossingHelixs.begin(), openedPseudo[i]->crossingHelixs.end(), root)
                                == openedPseudo[i]->crossingHelixs.end())
                        {
                            openedPseudo[i]->crossingHelixs.push_back(root);
                            root->crossingHelixs.push_back(openedPseudo[i]);
                        }
                        else if (openedPseudo[i]->pos_start.second < parenthesisClose[val].first)
                        {
                            openedPseudo.erase(openedPseudo.begin() + i);
                            size--;
                            i--;
                        }
                    }

                    helix *tmp = new helix();
                    tmp->pos_start=std::make_pair(std::get<0>(parenthesisClose[val+1]), std::get<1>(parenthesisClose[val+1]));

                    crossingHelixsTmp.clear();
                    while ((root != NULL) and
                            !(!root->pseudo and std::get<1>(root->pos_end) > std::get<1>(parenthesisClose[val+1])
                              and std::get<1>(root->pos_end) > std::get<0>(parenthesisClose[val+1])) // Normal child
                            and !(!root->pseudo and std::get<1>(root->pos_end) < std::get<1>(parenthesisClose[val+1])
                                  and std::get<0>(parenthesisClose[val+1]) < root->pos_end.second
                                  and root->pos_end.first < parenthesisClose[val+1].first)) // Pseudoknot child
                    {
                        if(root->pos_start.first < parenthesisClose[val+1].first and parenthesisClose[val+1].first < root->pos_start.second and root->pos_start.second < parenthesisClose[val+1].second)
                        {
                            crossingHelixsTmp.push_back(root);
                            root->crossingHelixs.push_back(tmp);
                        }
                        root=root->parent;
                    }


                    tmp->crossingHelixs = crossingHelixsTmp;
                    if (root!=NULL)
                    {
                        // New child = new branch with the same previous root
                        tmp->parent=root;
                        if (std::get<1>(root->pos_end) < std::get<1>(parenthesisClose[val+1]) and std::get<0>(parenthesisClose[val+1]) < root->pos_end.second and root->pos_end.first < parenthesisClose[val+1].first)
                        {
                            tmp->pseudo = true;
                            openedPseudo.push_back(tmp);
                        }
                        root2 = root;
                        while (root2 != NULL)
                        {
                            if(root2->pos_start.first < tmp->pos_start.first
                                    and tmp->pos_start.first < root2->pos_start.second
                                    and root2->pos_start.second < tmp->pos_start.second)
                            {
                                tmp->crossingHelixs.push_back(root2);
                                root2->crossingHelixs.push_back(tmp);
                            }
                            root2=root2->parent;
                        }
                        root->child.push_back(tmp);
                        root=tmp;

                    }
                    else
                    {
                        // New child = new branch with new root
                        // Add this new root in roots vector
                        roots.push_back(tmp);
                        root=tmp;
                        root->parent=NULL;
                        root->pseudo = false;
                    }
                }
            }
        }
    }

    if(!parenthesisClose.empty()){
        root->pos_end=std::make_pair(std::get<0>(parenthesisClose.back()), std::get<1>(parenthesisClose.back()));
        for (i = 0, size = openedPseudo.size(); i != size and root != NULL; i++)
        {
            if(((parenthesisClose.back().first < openedPseudo[i]->pos_start.second
                 and openedPseudo[i]->pos_start.second < parenthesisClose.back().second)
                    or (parenthesisClose.back().first < openedPseudo[i]->pos_start.first
                        and openedPseudo[i]->pos_start.first < parenthesisClose.back().second))
                and std::find(openedPseudo[i]->crossingHelixs.begin(), openedPseudo[i]->crossingHelixs.end(), root)
                    == openedPseudo[i]->crossingHelixs.end())
            {
                openedPseudo[i]->crossingHelixs.push_back(root);
                root->crossingHelixs.push_back(openedPseudo[i]);
            }
            else if (openedPseudo[i]->pos_start.second < parenthesisClose.back().first)
            {
                openedPseudo.erase(openedPseudo.begin() + i);
                size--;
                i--;
            }
        }
    }
    else
    {
        root->pos_end=std::make_pair(0,0);
    }

    openedPseudo.clear();

    //print the tree
    /*for(size_t i = 0, size = roots.size(); i != size; i++)
    {
        print(roots[i]);
    }*/

    // To execute the motifs_detection method and to stock the ensemble motif result in motif_storage
    for(unsigned int i(0); i < roots.size(); i++)
    {
        std::vector < Motif * > tmp3 = readTree(roots[i]);
        motifs_.insert(motifs_.end(), tmp3.begin(), tmp3.end());
    }

    /*for(size_t i = 0, size = motif_storage.size(); i != size; i++)
        for(size_t j = 0, size2 = motif_storage[i].size(); j != size2; j++)
            std::cout << motif_storage[i][j] << " ";
    std::cout << std::endl;*/

    for(auto val: roots)
        freeHelixSS(val);

}

std::string Structure::convToNUPACK() const {

    std::string res;
    int pos;
    char pairSymbols[] = { '(', ')', '{','}', '[', ']', '<', '>' };
    int type = 0;
    int nTypes = 4;
    int **pairlist; // Each row is i,j pair
    int npairs; // number of pairs in structure
    int seqlength = int(seq_.size());
    char *thefold;
    int *thepairs;
    thepairs = (int *) malloc(seqlength * sizeof(int *));

    //initialize the pairs
    int i = 0;
    for(int j = 0; j < seqlength; j++) {
        thepairs[j] = -1;
    }
    for(int j = 0; j < seqlength; j++) {
        if(i < int(listBp_.size())) {
            if(listBp_[i].first == uint(j))
                thepairs[j] = int(listBp_[i].second);
        }
        i++;
    }

    char *parensString;
    parensString = ( char*) malloc( (seqlength+1)*sizeof( char) );
    int lastL, lastR;

    // Allocate memory for pairlist (this is more than we need, but be safe)
    pairlist = (int **) malloc(seqlength * sizeof(int *));
    for (int i = 0; i < seqlength; i++) {
      pairlist[i] = (int *) malloc(2 * sizeof(int));
    }

    // Create pairlist from thepairs
    npairs = 0;
    for(int j = 0; j < seqlength; j++) {
      if(thepairs[j] > j) {
        pairlist[npairs][0] = j;
        pairlist[npairs++][1] = thepairs[j];
      }
    }

    // Creat dot-paren structure
    for( i = 0; i < seqlength; i++) {
      parensString[i] = '.';
    }

    //offSet = 0;
    lastL = -1;
    lastR = seqlength;
    for(int i = 0; i < seqlength; i++) {
      if( thepairs[i] != -1 && thepairs[i] > i) {
        if( i > lastR || thepairs[i] > lastR) {
          for(int j = 0; j < i; j++) {
            if( thepairs[j] > i && thepairs[j] < thepairs[i]) {
              type = (type + 1) % nTypes;
              break;
            }
          }
        }

        parensString[i] = pairSymbols[ 2*type];
        parensString[ thepairs[i]] = pairSymbols[2*type + 1];
        lastL = i;
        lastR = thepairs[i];
      }
    }

    thefold = (char*)malloc(seqlength+1 * sizeof(char));
    for( i = 0; i < seqlength; i++) {
      thefold[i] = '.';
    }
    thefold[seqlength] = '\0';

    pos = 0;
    for( i = 0; i < seqlength; i++) {
      thefold[pos++] = parensString[i];
    }

    free( parensString); parensString = NULL;

    for (i = 0; i < seqlength; i++) {
      free(pairlist[i]);
    }
    free(pairlist);
    free(thepairs);

    res = thefold;
    free(thefold); thefold = NULL;

    return res;
}

void Structure::computeEnergyNUPACK() {

    size_t i, j, size, size2;

    // Sequence initialisation
    std::string s = seq_;
    char * seq = &(s[0]);

    // Structure initialisation
    std::string str = convToNUPACK();
    char* parens = &(str[0]); // Structure in dot-paren notation of NUPACK

    int seqNum[MAXSEQLENGTH+1]; // Sequence into ASCII code

    //int thepairs[MAXSEQLENGTH];
    DBL_TYPE ene; // Energy of the structure
    int vs = 1; // Number of permutation, 1 for structure
    int tmpLength; // Sequence length into char

    tmpLength = strlen(seq);
    convertSeq(seq, seqNum, tmpLength); // Convert seq into ASCII code

    // Parameters initialisation
    // Initialize global parameters
    DNARNACOUNT = RNA;
    DANGLETYPE = 1;
    TEMP_K = 37.0 + ZERO_C_IN_KELVIN;

    DO_PSEUDOKNOTS = 1; // Enable pseudoknot
    SODIUM_CONC = 1.0;
    MAGNESIUM_CONC = 0.0;
    USE_LONG_HELIX_FOR_SALT_CORRECTION = 0;

    ene = naEnergyPairsOrParensFullWithSym( NULL, parens, seqNum, DNARNACOUNT, DANGLETYPE,
                      TEMP_K - ZERO_C_IN_KELVIN, vs,
                      SODIUM_CONC, MAGNESIUM_CONC,
                      USE_LONG_HELIX_FOR_SALT_CORRECTION);

    // Check to see if the results is close to NAD_INFINITY and report
    // error if it is
    /*if (ABS_FUNC(1.0 - ene/NAD_INFINITY) < INF_CUTOFF) {
      std::cout << "\n\n*** Error: invalid base pair(s) or disconnected complex. Check your inputs. ***\n\n" << std::endl;
    }*/

    obj1_ = float(ene);
}

void Structure::computeEnergyVienna()
{
    // Vienna RNA
    char *structure, *tmp, *rec_sequence, *rec_rest;
    std::string str = convToDPonly().c_str();
    vrna_md_t                 md;

    // Apply default model details
    vrna_md_set_default(&md);

    rec_sequence = new char[seq_.size() + 1];
    std::copy(seq_.begin(), seq_.end(), rec_sequence);
    rec_sequence[seq_.size()] = '\0';

    rec_rest = new char[str.size() + 1];
    std::copy(str.begin(), str.end(), rec_rest);
    rec_rest[str.size()] = '\0';

    char ** rec_rest2 = &rec_rest;

    vrna_fold_compound_t *vc = vrna_fold_compound(rec_sequence, &md, VRNA_OPTION_MFE | VRNA_OPTION_EVAL_ONLY);

    tmp = vrna_extract_record_rest_structure((const char **)rec_rest2, 0, 0);    // I had check until read_epars included


    int cp = -1;
    structure = vrna_cut_point_remove(tmp, &cp);

    free(tmp);

    // Computations
    obj1_ = vrna_eval_structure_v(vc, structure, 0, NULL);

    // Free
    free(rec_sequence);
    free(structure);
    free(rec_rest);
    rec_sequence = structure = NULL;
    rec_rest  = NULL;
    rec_rest2 = NULL;
    vrna_fold_compound_free(vc);
}

// Compute the F1-score between the probing data and the secondary structure
void Structure::computeProbing(const std::vector < float > &probingData, float lowerThresProbing, float upperThresProbing)
{
    probing_ = 0.0;

    if (std::accumulate(probingData.begin(), probingData.end(), 0) != -1*int(probingData.size())) {
        int TP = 0, TN = 0, FP = 0, FN = 0;
        float precision = 0, recall = 0;

        std::vector < uint > bp;

        for(size_t i = 0, size = seq_.size(); i != size; i++) {
            bp = find_bp_with_i(listBp_, i);
            if (probingData[i] <= (lowerThresProbing/100.0) and !bp.empty())
                TP++;
            if (probingData[i] >= (upperThresProbing/100.0) and !bp.empty())
                FP++;
            if (probingData[i] >= (upperThresProbing/100.0) and bp.empty())
                TN++;
            if (probingData[i] <= (lowerThresProbing/100.0) and bp.empty())
                FN++;
        }

        // MCC computation
        /*if (float((TP*TN - FP*FN)) - 0 > PRECISION
                and std::sqrt(float((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN))) - 0 > PRECISION)
            probing_ = float((TP*TN - FP*FN))/std::sqrt(float((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)));*/

        // Precision and recall computation
        if (TP+FP != 0) {
          precision = float(TP)/(TP+FP);
        }
        if (TP+FN != 0) {
          recall = float(TP)/(TP+FN);
        }

        // F1-score computation
        if (recall+precision != 0) {
          probing_ = 2*((recall*precision)/(recall+precision));
        }
    }
}

bool Structure::checkListBp(std::string seq, std::vector < std::pair < unsigned int, unsigned int > > &listBp)
{
    bool res = true;
    for(size_t i = 0, size = listBp.size(); i != size; i++)
    {
        if(!app(seq[listBp[i].first], seq[listBp[i].second]))
        {
            std::stringstream ss;
            ss << seq[listBp[i].first] << "-" << seq[listBp[i].second];
            std::string bp;
            ss >> bp;
            res = false;
            throw (std::string("Error the base pair " + bp + "(" + std::to_string(listBp[i].first) + "-" + std::to_string(listBp[i].second) + ") isn't allowed."));
        }
    }

    return res;
}