#include <iostream>
#include <cstring>
#include <sstream>
#include <numeric>

#include <boost/format.hpp>

// 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 "solinteraction.h"
#include "utils.h"

const float ENERGY[7][7] = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},{0.0,-1.1, -2.1, -2.2,-1.4,-0.9,-0.6},{0.0,-2.1,-2.4,-3.3,-2.1,-2.1,-1.4},{0.0,-2.2,-3.3,-3.4,-2.5,-2.4,-1.5},{0.0,-1.4,-2.1,-2.5,-1.3,-1.3,-0.5},{0.0,-0.9,-2.1,-2.4,-1.3,-1.3,-1.0},{0.0,-0.6,-1.4,-1.5,-0.5,-1.0,-0.3}};

SolInteraction::SolInteraction()
{
}

SolInteraction::SolInteraction(const unsigned int rna1, const unsigned int rna2, const std::string &seq1, const std::string &seq2):
    rna1_(rna1), rna2_(rna2), seq1_(seq1), seq2_(seq2)
{
}

SolInteraction::SolInteraction(const unsigned int rna1, const unsigned int rna2, const std::string &seq1,
                               const std::string &seq2,
                               const std::vector < std::pair < unsigned int, unsigned int > > &listBp,
                               const int ct,
                               const unsigned int nbCt, const unsigned int id,
                               unsigned int energyModel, float energy,
                               const std::vector < float > &probingDataRna1,
                               const std::vector < float > &probingDataRna2,
                               float lowerThresProbing, float upperThresProbing):
    rna1_(rna1), rna2_(rna2), seq1_(seq1), seq2_(seq2), listBp_(listBp), ct_(ct), nbCt_(nbCt), id_(id)
{
    std::sort(listBp_.begin(), listBp_.end());
    motifDetection();

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

    computeProbing(probingDataRna1, probingDataRna2, lowerThresProbing, upperThresProbing);
}

SolInteraction::SolInteraction(const SolInteraction& that)
{
    rna1_ = that.rna1_;
    rna2_ = that.rna2_;
    seq1_ = that.seq1_;
    seq2_ = that.seq2_;
    listBp_ = that.listBp_;
    score_ = that.score_;
    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 (Pseudoknot * P = dynamic_cast < Pseudoknot * > (that.motifs_[i]))
        {
            motifs_[i] = new Pseudoknot (*P);
        }
    }
}


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

unsigned int SolInteraction::get_rna1_() const
{
    return rna1_;
}

unsigned int SolInteraction::get_rna2_() const
{
    return rna2_;
}

std::string SolInteraction::get_seq1_() const
{
    return seq1_;
}

std::string SolInteraction::get_seq2_() const
{
    return seq2_;
}

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

float SolInteraction::get_score_() const
{
    return score_;
}

unsigned int SolInteraction::get_id_() const
{
    return id_;
}

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

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

unsigned int SolInteraction::get_nbCt_() const
{
    return nbCt_;
}

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

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

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

void SolInteraction::set_score_(float score)
{
    score_ = score;
}

void SolInteraction::set_listBp_(const std::vector < std::pair < unsigned int, unsigned int > > &listBp)
{
    listBp_ = listBp;
}

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

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

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

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

std::string SolInteraction::convToDP() const
{
    std::string first = seq1_;
    std::string second = "";
    int i, size;
    for(i = 0, size = int(first.size()); i != size; i++)
    {
        if(find_i_at_first(listBp_, i))
        {
                second  += "[";
        }
        else
        {
            second  += " ";
        }
    }
    std::string third = seq2_;
    std::string last = "";
    for(i = 0, size = int(third.size()); i != size; i++)
    {
        if(find_i_at_second(listBp_, i))
        {
                last  += "]";
        }
        else
        {
            last  += " ";
        }
    }

    return first + "\n" + second + "\n" + third + "\n" + last;
}

std::string SolInteraction::convToDPInvert() const
{

    std::string first = seq1_;
    std::string second = "";
    int i, size;
    for(i = 0, size = int(first.size()); i != size; i++)
    {
        if(find_i_at_first(listBp_, i))
        {
                second  += "[";
        }
        else
        {
            second  += " ";
        }
    }
    std::string third = seq2_;
    std::string last = "";
    for(i = 0, size = int(third.size()); i != size; i++)
    {
        if(find_i_at_second(listBp_, i))
        {
                last  += "]";
        }
        else
        {
            last  += " ";
        }
    }
    return third + "\n" + last  + "\n" + first + "\n" + second;
}

std::string SolInteraction::to_string() const
{
    return convToDP() + "\tScore = " + boost::str(boost::format("%.2f") % score_);
}

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

SolInteraction& SolInteraction::operator=(const SolInteraction& that)
{
    rna1_ = that.rna1_;
    rna2_ = that.rna2_;
    seq1_ = that.seq1_;
    seq2_ = that.seq2_;
    listBp_ = that.listBp_;
    score_ = that.score_;
    ct_ = that.ct_;
    ic_ = that.ic_;
    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 (Pseudoknot * P = dynamic_cast < Pseudoknot * > (that.motifs_[i]))
        {
            local_motifs_[i] = new Pseudoknot (*P);
        }
    }
    for(size_t i = 0, size = motifs_.size(); i != size; i++)
        delete motifs_[i];
    motifs_ = local_motifs_;

    return *this;
}

bool SolInteraction::operator<(const SolInteraction& s)
{
    return score_ < s.score_;
}

bool SolInteraction::operator<(const SolInteraction& s) const
{
    return score_ < s.score_;
}

bool SolInteraction::operator>(const SolInteraction& s)
{
    return score_ > s.score_;
}

bool SolInteraction::operator>(const SolInteraction& s) const
{
    return score_ > s.score_;
}

bool SolInteraction::operator<=(const SolInteraction& s)
{
    return score_ <= s.score_;
}

bool SolInteraction::operator<=(const SolInteraction& s) const
{
    return score_ <= s.score_;
}

bool SolInteraction::operator>=(const SolInteraction& s)
{
    return score_ >= s.score_;
}

bool SolInteraction::operator>=(const SolInteraction& s) const
{
    return score_ >= s.score_;
}

bool SolInteraction::operator==(const SolInteraction& s)
{
    /*bool res = true;
    if( listBp_.size() == s.listBp_.size())
    {
        if(s.get_rna1_() == get_rna1_() and s.get_rna2_() == get_rna2_())
        {
            for(unsigned int i = 0; i < listBp_.size(); i++)
            {
                if(listBp_.get_coord(i) != s.listBp_.get_coord(i))
                {
                    res = false;
                    break;
                }
            }
        }
        else if (s.get_rna1_() == get_rna2_() and s.get_rna2_() == get_rna1_())
        {
            for(unsigned int i = 0; i < listBp_.size(); i++)
            {
                if(listBp_.get_coordF(i) != s.listBp_.get_coordS(i) or listBp_.get_coordS(i) != s.listBp_.get_coordF(i))
                {
                    res = false;
                    break;
                }
            }
        }
        else
        {
            res = false;
        }
    }
    else
    {
        res = false;
    }*/

    return s.rna1_ == rna1_ and s.rna2_ == rna2_ and s.id_ == id_;
}


bool SolInteraction::operator==(const SolInteraction& s) const
{
    /*bool res = true;
    if( listBp_.size() == s.listBp_.size())
    {
        if(s.get_rna1_() == get_rna1_() and s.get_rna2_() == get_rna2_())
        {
            for(unsigned int i = 0; i < listBp_.size(); i++)
            {
                if(listBp_.get_coord(i) != s.listBp_.get_coord(i))
                {
                    res = false;
                    break;
                }
            }
        }
        else if (s.get_rna1_() == get_rna2_() and s.get_rna2_() == get_rna1_())
        {
            for(unsigned int i = 0; i < listBp_.size(); i++)
            {
                if(listBp_.get_coordF(i) != s.listBp_.get_coordS(i) or listBp_.get_coordS(i) != s.listBp_.get_coordF(i))
                {
                    res = false;
                    break;
                }
            }
        }
        else
        {
            res = false;
        }
    }
    else
    {
        res = false;
    }*/

    return s.rna1_ == rna1_ and s.rna2_ == rna2_ and s.id_ == id_;
}

std::vector < std::pair < unsigned int, unsigned int > > SolInteraction::convToBP(std::string structure)
{
    std::vector<std::pair<int,int>> opening;
    std::vector < std::pair < unsigned int, unsigned int > > listBp;
    int index = 0;
    int rna1 = 0;
    int rna2 = 1;
    int pass = false;
    for(size_t i = 0, size = structure.size(); i != size; i++)
    {
        if (structure[i] == '(')
        {
            if (!pass)
            {
                opening.push_back(std::make_pair(rna1, index));
            }
            else
            {
                opening.push_back(std::make_pair(rna2, index));
            }
        }
        else if (structure[i] == ')')
        {
            if(pass and std::get<0>(opening[opening.size() -1]) != rna2)
            {
                listBp.push_back(std::make_pair(std::get<1>(opening[opening.size()-1]), index));
            }
            opening.pop_back();
        }
        else if (structure[i] == '&')
        {
            pass = true;
            index = -1;
        }
        index++;
    }

    return listBp;
}

std::string SolInteraction::convToDPRnaEval()
{
    std::string second = "";
    for(unsigned int i = 0; i < uint(seq1_.size()); i++) {
        if(find_i_at_first(listBp_, i)) {
                second  += "(";
        } else {
            second  += ".";
        }
    }
    std::string last = "";
    for(int i = 0; i < uint(seq2_.size()); i++) {
        if(find_i_at_second(listBp_, i)){
                last  += ")";
        } else {
            last  += ".";
        }
    }

    return second + "&" + last;
}

bool SolInteraction::checkStructure(std::string structure)
{
    bool res = true;
    long oPar = std::count(structure.begin(), structure.end(), '(');
    long cPar = std::count(structure.begin(), structure.end(), ')');
    long oBra = std::count(structure.begin(), structure.end(), '[');
    long pBra = std::count(structure.begin(), structure.end(), ']');

    long link = std::count(structure.begin(), structure.end(), '&');
    if(oPar != cPar or oBra != pBra or int(link) != 1) {
        res = false;
        throw(std::string("Invalid number of base pairs."));
    }

    return res;
}

// Simple model to compute energy
void SolInteraction::computeEnergy()
{
    score_ = 0.0;
    int nti1, nti2, ntj1, ntj2, typei, typej, j;
    for(uint i = 0, size = uint(listBp_.size()); i != size; i++)
    {
        nti1 = listBp_[i].first;
        nti2 = listBp_[i].second;
        ntj1 = nti1 +1;
        ntj2 = nti2 -1;
        j = find_bp(listBp_, std::make_pair(ntj1, ntj2) );
        if(j != -1 )
        {
            std::string bpi, bpj;
            bpi = seq1_[nti1];
            bpj = seq1_[ntj1];
            bpi += seq2_[nti2];
            bpj += seq2_[ntj2];

            typei = type(bpi);
            typej = type(bpj);
            score_ += ENERGY[typei][typej];
        }
    }
}

int SolInteraction::type(std::string bp)
{
    int type = 0;
    if(bp.compare("AU") == 0)
    {
        type = 1;
    }
    else if(bp.compare("CG") == 0 )
    {
        type = 2;
    }
    else if(bp.compare("GC") == 0 )
    {
        type = 3;
    }
    else if(bp.compare("GU") == 0 )
    {
        type = 4;
    }
    else if(bp.compare("UG") == 0 )
    {
        type = 5;
    }
    else if(bp.compare("UA") == 0 )
    {
        type = 6;
    }

    return type;
}

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

    if(tmp->pos_start.first < tmp->pos_end.first) {
        try {
            res.push_back(new Helix(rna1_, rna2_,
                                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),
                                seq1_, seq2_));
        } catch (std::string e) {
            std::cout << e << std::endl;
        }
    }
    for(size_t i = 0, size = tmp->crossingHelixs.size(); i != size; i++) {
        try {
            res.push_back(new Pseudoknot(rna1_, rna2_,
                                 Helix(rna1_, rna2_,
                                       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),
                                       seq1_, seq2_),
                                 Helix(rna1_, rna2_,
                                       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),
                                       seq1_, seq2_)));
        } catch (std::string e) {
            std::cout << e << std::endl;
        }
    }

    return res;
}

void SolInteraction::motifDetection()
{
    // Tree generation of the helixI positions (only with the this->parenthesis vector)
    std::vector < helixI * > roots;
    helixI *root = new helixI();

    std::vector < std::pair < unsigned int, unsigned int > > listBp = listBp_;

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

    size_t i, size;

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

            if(diff)
            {
                root->pos_end=std::make_pair(std::get<0>(listBp[val]), std::get<1>(listBp[val]));
                for (i = 0, size = roots.size(); i != size; i++)
                    if(listBp[val].second > roots[i]->pos_start.second )
                        root->crossingHelixs.push_back(roots[i]);

                helixI *tmp= new helixI ();
                tmp->pos_start=std::make_pair(std::get<0>(listBp[val+1]), std::get<1>(listBp[val+1]));
                roots.push_back(tmp);
                root=tmp;
            }
        }
    }

    if(!listBp.empty()){
        root->pos_end=std::make_pair(std::get<0>(listBp.back()), std::get<1>(listBp.back()));
        for (i = 0, size = roots.size(); i != size; i++)
            if(listBp.back().second > roots[i]->pos_start.second)
                root->crossingHelixs.push_back(roots[i]);
    }
    else
    {
        root->pos_end=std::make_pair(0,0);
    }

    // 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(auto val: roots)
        delete val;
}

bool SolInteraction::checkListBp(std::string seq1, std::string seq2, 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(seq1[listBp[i].first], seq2[listBp[i].second])) {
            res = false;
            throw (std::string("Error in the motif or constraint, base pair " + std::to_string(seq1[listBp[i].first])
                + "-" + std::to_string(seq2[listBp[i].second]) + "isn't allow."));
        }
    }

    return res;
}

std::string SolInteraction::convToNUPACK() const
{
    int i,j, 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(seq1_.size()) + int(seq2_.size()+1);
    char *thefold;
    int thepairs[int(seq1_.size())+ int(seq2_.size())+1];

    // Initialize the pairs
    for (int i = 0, size = int(seq1_.size())+ int(seq2_.size())+1; i != size; i++) {
        thepairs[i] = -1;
    }
    for (size_t i = 0, size = listBp_.size(); i != size; i++) {
        thepairs[listBp_[i].first] = listBp_[i].second + int(seq1_.size()) + 1;
        thepairs[listBp_[i].second + int(seq1_.size()) + 1] = listBp_[i].first;
    }

    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 (i = 0; i < seqlength; i++) {
    pairlist[i] = (int *) malloc(2 * sizeof(int));
    }

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

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

    //offSet = 0;
    lastL = -1;
    lastR = seqlength;
    for( i = 0; i < seqlength; i++) {
        if( thepairs[i] != -1 && thepairs[i] > i) {
            if( i > lastR || thepairs[i] > lastR) {
                for( 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];
    }
    thefold[int(seq1_.size())] = '+';

    free(parensString); parensString = NULL;

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

    return thefold;
}

void SolInteraction::computeEnergyNUPACK() {

    size_t i, size;

    // Sequence initialisation
    std::string s = seq1_ + '+' + seq2_;
    char * seq = &(s[0]);

    // Structure initialisation
    std::string structNupack = convToNUPACK();
    char parens[structNupack.size()]; // Structure in dot-paren notation of NUPACK
    for(i = 0, size = structNupack.size(); i != size; i++)
        parens[i] = structNupack[i];

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

    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;
    }*/

    score_ = float(ene);
}

void SolInteraction::computeEnergyVienna() {

    char *structure, *tmp, *rec_sequence, *rec_rest;
    std::string str = convToDPRnaEval();
    rec_rest = &(str[0]);
    char **rec_rest2 = &(rec_rest);
    std::string seq = seq1_ + '&' + seq2_;
    rec_sequence = &(seq[0]);
    vrna_md_t md;

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

    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);

    // Computations
    vc->cutpoint = int(seq1_.size());
    score_ = vrna_eval_structure_v(vc, structure, 0, NULL); // Erreur à cette ligne

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

// Compute the F1-score between the probing data and the interaction and make it the probing weight
void SolInteraction::computeProbing(const std::vector < float > &probingDataRna1, const std::vector < float > &probingDataRna2, float lowerThresProbing, float upperThresProbing)
{
    probing_ = 0.0;

    if (std::accumulate(probingDataRna1.begin(), probingDataRna1.end(), 0) != -1*int(probingDataRna1.size()) and
        std::accumulate(probingDataRna2.begin(), probingDataRna2.end(), 0) != -1*int(probingDataRna2.size())) {

        int TP = 0, TN = 0, FP = 0, FN = 0;
        float precision = 0, recall = 0;
        bool bp;

        for(size_t i = 0, size = seq1_.size(); i != size; i++) {
            bp = find_i_at_first(listBp_, i);
            if (probingDataRna1[i] <= (lowerThresProbing/100) and bp == true)
                TP++;
            if (probingDataRna1[i] >= (upperThresProbing/100) and bp == true)
                FP++;
            if (probingDataRna1[i] >= (upperThresProbing/100) and bp == false)
                TN++;
            if (probingDataRna1[i] <= (lowerThresProbing/100) and bp == false)
                FN++;
        }

        for(size_t i = 0, size = seq2_.size(); i != size; i++) {
            bp = find_i_at_second(listBp_, i);
            if (probingDataRna2[i] <= (lowerThresProbing/100) and bp == true)
                TP++;
            if (probingDataRna2[i] >= (upperThresProbing/100) and bp == true)
                FP++;
            if (probingDataRna2[i] >= (upperThresProbing/100) and bp == false)
                TN++;
            if (probingDataRna2[i] <= (lowerThresProbing/100) and bp == false)
                FN++;
        }

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

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

    }
}