#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <iterator>
#include <functional>
#include <utility>
#include <algorithm>
#include <string>
#include <cctype>
#include <cstring>

#include <unistd.h>
#include <math.h>

#include "graph.h"
#include "utils.h"

// float compatibility between two structures Peut-etre besoin pour pseudonoeud ?
float compatible(Structure str1, Structure str2) { //, float threshold, int l) {
    float match = 10;

    if( str1.get_rna_() == str2.get_rna_())
    {
       /*int maxConflict = -floor((threshold/100.0 -1.0)*float(l));
       int counter = 0;
       size_t i, size;
       std::vector < std::pair < unsigned int, unsigned int > > l1 = str1.get_listBp_(), l2 = str2.get_listBp_();
       for(i = 0, size = l1.size(); i != size and counter <= maxConflict; i++)
       {
           std::vector < unsigned int > coord = find_bp_with_i(l2, l1[i].first);
           if (!coord.empty()) {
               if ((l2[coord[0]].first == l1[i].first and l2[coord[0]].second != l1[i].second)
                       or (l2[coord[0]].second == l1[i].first and l2[coord[0]].first != l1[i].second))
                   counter++;
           }
           coord = find_bp_with_i(l2, l1[i].second);
           if (!coord.empty()) {
               if ((l2[coord[0]].first == l1[i].second and l2[coord[0]].second != l1[i].first)
                       or (l2[coord[0]].second == l1[i].second and l2[coord[0]].first != l1[i].first))
                   counter++;
           }
       }
       match = 0;
       if (counter <= maxConflict)
           match = float(l-counter)/float(l)*100.0;*/
        match = 0;
    }
    return match;
}

// float compatibility between a structure and an interaction
float compatible(Structure str, SolInteraction inter, float threshold, int l) {
    float match = 10;
    int counter = 0;
    size_t i, size;
    if(uint(str.get_rna_()) == inter.get_rna1_()) {
       int maxConflict = floor(-(threshold/100.0 -1.0)*float(l));
       std::vector < std::pair < unsigned int, unsigned int > > l1 = str.get_listBp_(), l2 = inter.get_listBp_();
       for(i = 0, size = l1.size(); i != size and counter <= maxConflict; i++) {
           bool coord = find_i_at_first(l2, l1[i].first);
           if (coord)
               counter++;
           coord = find_i_at_first(l2, l1[i].second);
           if (coord)
               counter++;
       }
       match = 0;
       if (counter <= maxConflict)
           match = float(l-counter)/float(l)*10.0;
    } else if (str.get_rna_() == inter.get_rna2_()) {
        int maxConflict = floor(-(threshold/100.0 -1.0)*float(l));
        std::vector < std::pair < unsigned int, unsigned int > > l1 = str.get_listBp_(), l2 = inter.get_listBp_();
        for(i = 0, size = l1.size(); i != size and counter <= maxConflict; i++)
        {
            bool coord = find_i_at_second(l2, l1[i].first);
            if (coord)
                counter++;
            coord = find_i_at_second(l2, l1[i].second);
            if (coord)
                counter++;
        }
        match = 0;
        if (counter <= maxConflict)
            match = float(l-counter)/float(l)*10.0;
    }
    return match;
}

// float compatibility between two interactions
float compatible(SolInteraction inter1, SolInteraction inter2, float threshold) {
    float match = 10;
    int len1, len2;

    if(inter1.get_rna1_() == inter2.get_rna1_() and inter1.get_rna2_() == inter2.get_rna2_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size())*2;
        len2 = int(l2.size())*2;
       int maxConflict = floor(-double(threshold/100.0 -1.0)*(double(len1) + double(len2)));
       int counter = 0;
       size_t i, size;
       for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
           int coord = find_i_at_first_int(l2, l1[i].first);
           if (coord != -1)
               //if (l2[ulong(coord)].second != l1[i].second)
                   counter++;
           coord = find_i_at_second_int(l2, l1[i].second);
           if (coord != -1)
               //if (l2[ulong(coord)].first != l1[i].first)
                   counter++;
       }
       match = 0;
       if (counter <= maxConflict)
           match = float(len1 + len2 - counter)/float(len1 + len2)*10.0;
    } else if(inter1.get_rna1_() == inter2.get_rna2_() and inter1.get_rna2_() == inter2.get_rna1_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size())*2;
        len2 = int(l2.size())*2;
        int maxConflict = floor(-(threshold/100.0 -1.0)*float(len1 + len2));
        int counter = 0;
        size_t i, size;
        for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
            int coord = find_i_at_second(l2, l1[i].first);
            if (coord != -1)
                //if (l2[ulong(coord)].second == l1[i].first and l2[ulong(coord)].first != l1[i].second)
                    counter++;
            coord = find_i_at_first(l2, l1[i].second);
            if (coord != -1)
                //if (l2[ulong(coord)].first == l1[i].second and l2[ulong(coord)].second != l1[i].first)
                    counter++;
        }
        match = 0;
        if (counter <= maxConflict)
            match = float(len1 + len2 - counter)/float(len1 + len2)*10.0;
     } else if(inter1.get_rna1_() == inter2.get_rna1_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size());
        len2 = int(l2.size());
        int maxConflict = floor(-(threshold/100.0 -1.0)*float(len1 + len2));
        int counter = 0;
        size_t i, size;
        for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
            bool coord = find_i_at_first(l2, l1[i].first);
            if (coord)
                counter++;
        }
        match = 0;
        if (counter <= maxConflict)
            match = float(len1 + len2 -counter)/float(len1 + len2)*10.0;
     } else if (inter1.get_rna2_() == inter2.get_rna2_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size());
        len2 = int(l2.size());
         int maxConflict = floor(-(threshold/100.0 -1.0)*float(len1 + len2));
         int counter = 0;
         size_t i, size;
         for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
             bool coord = find_i_at_second(l2, l1[i].second);
             if (coord)
                 counter++;
         }
         match = 0;
         if (counter <= maxConflict)
             match = float(len1 + len2-counter)/float(len1 + len2)*10.0;
     } else if(inter1.get_rna1_() == inter2.get_rna2_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size());
        len2 = int(l2.size());
        int maxConflict = floor(-(threshold/100.0 -1.0)*float(len1 + len2));
        int counter = 0;
        size_t i, size;
        for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
            bool coord = find_i_at_second(l2, l1[i].first);
            if (coord)
                counter++;
        }
        match = 0;
        if (counter <= maxConflict)
            match = float(len1 + len2 -counter)/float(len1 + len2)*10.0;
     } else if (inter1.get_rna2_() == inter2.get_rna1_()) {
        std::vector < std::pair < unsigned int, unsigned int > > l1 = inter1.get_listBp_(), l2 = inter2.get_listBp_();
        len1 = int(l1.size());
        len2 = int(l2.size());
         int maxConflict = floor(-(threshold/100.0 -1.0)*float(len1 + len2));
         int counter = 0;
         size_t i, size;
         for(i = 0, size = l1.size(); i != size and counter <= maxConflict+1; i++) {
             bool coord = find_i_at_first(l2, l1[i].second);
             if (coord)
                 counter++;
         }
         match = 0;
         if (counter <= maxConflict)
             match = float(len1 + len2-counter)/float(len1 + len2)*10.0;
     }
    return match;
}

// Two intervals are linked ?
bool areLinked (Helix h1, Helix h2, const std::vector < Rna > &rnaList) {
    bool match = false;

    std::vector< std::pair < unsigned int, unsigned int > > bps1 = h1.getListBp_(), bps2 = h2.getListBp_();

    if(!bps1.empty() and !bps2.empty()) {
        int i1, j1, i1bis, j1bis;
        if(h1.getRna_() == h2.getRna_() and h1.getRna2_() == h2.getRna2_()) {
            j1 = int(h1.getJ1_()) + int(rnaList[ulong(h1.getRna_())].get_n_());
            j1bis = int(h2.getJ1_()) + int(rnaList[ulong(h2.getRna_())].get_n_());
            i1 = h1.getI1_();
            i1bis = h2.getI1_();

            if((i1 < i1bis and i1bis < j1 and j1 < j1bis)
                   or (i1bis < i1 and i1 < j1bis and j1bis < j1))
                match = true;
        }
        else if(h1.getRna_() == h2.getRna2_() and h1.getRna2_() == h2.getRna_()) {
            j1 = int(h1.getJ1_()) + int(rnaList[ulong(h1.getRna_())].get_n_());
            j1bis = int(h2.getI1_()) + int(rnaList[ulong(h1.getRna_())].get_n_());
            i1 = h1.getI1_();
            i1bis = h2.getJ1_();

            if((i1 < i1bis and i1bis < j1 and j1 < j1bis)
                    or (i1bis < i1 and i1 < j1bis and j1bis < j1))
                match = true;
        }

       // std::cout << "h1 " << h1.getId_() << " h2 " << h2.getId_() << " " << match <<
                    // " .... " << i1 << " " << i1bis << " " << j1 << " " << j1bis << std::endl;
    }
    return match;
}

void modelGraph(const std::vector < Rna > &rnaList_,
                const std::vector < Structure > &strs_, const std::vector < SolInteraction > &inters_,
                const std::vector < std::vector < Helix > > &helices_, std::vector < unsigned int > &vertices0,
                std::vector < float > &verticesWCt, std::vector < float > &verticesWP,
                std::vector < float > &verticesWE, std::vector < std::vector < unsigned int > > &Nv0,
                std::vector < std::vector < unsigned int > > &NvB0, std::vector < std::vector < unsigned int > > &G2vertices0,
                std::vector < std::vector < std::vector < unsigned int > > > &G2Nv0, std::vector < std::vector < unsigned int > > &G2G1correspondance0,
                //std::vector < unsigned int > &ctVertices, std::vector < unsigned int > &rmVertices,
                float threshold, std::string SEQfile, uint &nbEdges, int &lastStr) {

    unsigned int i, size, j, size2;

    // Recover vertices GRAPH 1
    // Secondary structures
    for(i=0, size = uint(strs_.size()); i != size; i++) {
            vertices0.push_back(uint(vertices0.size()));
            verticesWE.push_back(strs_[i].get_obj1_());
            verticesWCt.push_back(strs_[i].get_ic_());
            verticesWP.push_back(strs_[i].get_probing_());
            // Constraints vertices
            /*if(strs_[i].get_ct_() == 100) {
                ctVertices.push_back(uint(vertices0.size())-1);
                rmVertices.push_back(uint(vertices0.size())-1);
            }*/
    }
    lastStr = -1;
    if (vertices0.size() != 0)
        lastStr = int(vertices0.size()-1);
    // Interaction
    for(i=0, size = uint(inters_.size()); i != size; i++) {
        vertices0.push_back(uint(vertices0.size()));
        verticesWE.push_back(inters_[i].get_score_());
        verticesWCt.push_back(inters_[i].get_ic_());
        verticesWP.push_back(inters_[i].get_probing_());

        // Constraints vertices
        /*if(inters_[i].get_ct_() == 100) {
            ctVertices.push_back(uint(vertices0.size())-1);
            rmVertices.push_back(uint(vertices0.size())-1);
        }*/
    }
    // Recover vertices GRAPH 2
    for(i = 0, size = helices_.size(); i != size; i++) {
        G2vertices0.push_back(std::vector < uint > ());
        G2Nv0.push_back(std::vector < std::vector < uint > > ());
        G2G1correspondance0.push_back(std::vector < uint > ());
        for (j = 0, size2 = helices_[i].size(); j != size2; j++) {
            G2vertices0[i].push_back(uint(j));
            if(lastStr != -1)
                G2G1correspondance0[i].push_back(helices_[i][j].getIdSolInter_()+lastStr+1);
            else
                G2G1correspondance0[i].push_back(helices_[i][j].getIdSolInter_());
        }
    }
#ifdef _DEBUG
    std::cout << "avant recover edge G1" << std::endl;
#endif
    // G1 Recover edges, find if match between interaction/interaction and interaction/structure
    float wCp;
    for (i = 0, size = uint(vertices0.size()); i != size; i++) {
        Nv0.push_back(std::vector< unsigned int > ());
        //NvWCp.push_back(std::vector< float > ());
        NvB0.push_back(std::vector< unsigned int > ());

        for (j = 0, size2 = uint(vertices0.size()); j != size2; j++) {
            if(i != j) {
                if (int(i) <= lastStr and int(j) <= lastStr) {
                    wCp = compatible(strs_[vertices0[i]], strs_[vertices0[j]]);
                } else if (int(i) <= lastStr and int(j) > lastStr) {
                    wCp = compatible(strs_[vertices0[i]], inters_[vertices0[j]-uint(lastStr)-1],
                            threshold, int(rnaList_[ulong(strs_[vertices0[i]].get_rna_())].get_n_()));
                } else if (int(i) > lastStr and int(j) <= lastStr) {
                    wCp = compatible(strs_[vertices0[j]], inters_[vertices0[i]-uint(lastStr)-1],
                                threshold, int(rnaList_[ulong(strs_[vertices0[j]].get_rna_())].get_n_()));
                } else {
                    if (lastStr != -1)
                        wCp = compatible(inters_[vertices0[i]-uint(lastStr)-1], inters_[vertices0[j]-uint(lastStr)-1],
                                threshold);
                    else
                        wCp = compatible(inters_[vertices0[i]], inters_[vertices0[j]],
                                threshold);
                }

                //std::cout << i << " " << j << " " << wCp << std::endl;
                if (wCp != 0.0) {
                    //NvWCp[i].push_back(wCp);
                    Nv0[i].push_back(uint(j));
                } else {
                    //NvWCp[i].push_back(0);
                    NvB0[i].push_back(uint(j));
                }
            }
        }
    }
#ifdef _DEBUG
    std::cout << "avant recover edge G2" << std::endl;
#endif
    // G2 Recover edges, find if match between intervals
    bool match;
    for (i = 0, size = uint(G2vertices0.size()); i != size; i++) {
        for (j = 0, size2 = uint(G2vertices0[i].size()); j != size2; j++) {
            G2Nv0[i].push_back(std::vector< unsigned int > ());

#ifdef _DEBUG
            std::cout << "combi " << i << ", " << j << " : ";
#endif
            // Pourquoi c'est le graph g2_2 (1-2) qui est rempli et pas g2_1 (0-2)
            for (size_t k = 0, size3 = uint(G2vertices0[i].size()); k != size3; k++) {
                if(j != k and std::find(Nv0[G2G1correspondance0[i][j]].begin(),
                                        Nv0[G2G1correspondance0[i][j]].end(), G2G1correspondance0[i][k])
                        != Nv0[G2G1correspondance0[i][j]].end()) {
                        // The two SolInteraction are compatible, then check if the helices are pseudoknotted
                        match = areLinked(helices_[i][j], helices_[i][k],
                                rnaList_);
                    if (match) {
#ifdef _DEBUG
                        std::cout << k << " ";
#endif
                        G2Nv0[i][j].push_back(uint(k));
                    }
                }
            }
#ifdef _DEBUG
            std::cout << std::endl;
#endif
        }
    }

#ifdef _DEBUG
    std::cout << "SANS CT nbV=" << vertices0.size() << std::endl;
#endif
}

/*void addHardCt(std::vector < unsigned int > &vertices,
               std::vector < std::vector < unsigned int > > &Nv,
               std::vector < std::vector < unsigned int > > &NvB,
               const std::vector < std::vector < unsigned int > > &NvB0,
               std::vector < std::vector < unsigned int > > &G2vertices,
               std::vector < std::vector < std::vector < unsigned int > > > &G2Nv,
               std::vector < std::vector < unsigned int > > &G2G1correspondance,
               std::vector < unsigned int > &ctVertices, std::vector < unsigned int > &rmVertices) {

    uint i, size, j, size2;
    std::cout << "Il faut recuperer les sommets (pour pas les implémenter) qui ne sont compatibles avec aucune des ct qui sont demandées ? " << ctVertices.size() << std::endl;
    std::sort(ctVertices.begin(), ctVertices.end());

    std::vector < unsigned int > dest, intersec;

    // Vertices to remove: is computed over G1
    intersec = NvB0[ctVertices[0]];
    for(size_t i = 1, size = ctVertices.size(); i != size; i++) {
        std::cout << "ctVertices[i] : " << ctVertices[i] << std::endl;
        dest.clear();
        std::set_intersection(intersec.begin(), intersec.end(), NvB0[ctVertices[i]].begin(), NvB0[ctVertices[i]].end(), std::back_inserter(dest));
        intersec = dest;
    }

    rmVertices = ctVertices;
    rmVertices.insert(rmVertices.end(), intersec.begin(), intersec.end());
    std::sort(rmVertices.begin(), rmVertices.end());


    // Remove the constrained and to be removed vertices
    for(i = 0, size = uint(vertices.size()); i != size; i++)
        if (std::find(ctVertices.begin(), ctVertices.end(), vertices[i]) != ctVertices.end()) {
            vertices.erase(vertices.begin() + i);
            i--;
            size--;
        }

#ifdef _DEBUG
    std::cout << "vertices" << std::endl;
#endif
    for(i = 0, size = uint(Nv.size()); i != size; i++)
        for(j = 0, size2 = uint(Nv[i].size()); j != size2; j++)
            if (std::find(ctVertices.begin(), ctVertices.end(), Nv[i][j]) != ctVertices.end()) {
                Nv[i].erase(Nv[i].begin() + j);
                j--;
                size2--;
            }
#ifdef _DEBUG
    std::cout << "Nv" << std::endl;
#endif
    for(i = 0, size = uint(NvB.size()); i != size; i++)
        for(j = 0, size2 = uint(NvB[i].size()); j != size2; j++)
            if (std::find(ctVertices.begin(), ctVertices.end(), NvB[i][j]) != ctVertices.end()) {
                NvB[i].erase(NvB[i].begin() + j);
                j--;
                size2--;
            }
#ifdef _DEBUG
    std::cout << "NvB" << std::endl;
#endif

    // G2
    for(i = 0, size = uint(G2vertices.size()); i != size; i++) {
        for(j = 0, size2 = uint(G2vertices[i].size()); j != size2; j++)
            if (std::find(ctVertices.begin(), ctVertices.end(), G2vertices[i][j]) != ctVertices.end()) {
                G2vertices[i].erase(G2vertices[i].begin() + j);
                G2G1correspondance[i].erase(G2G1correspondance[i].begin() + j);
                j--;
                size2--;
            }
#ifdef _DEBUG
        std::cout << "G2 vertices" << std::endl;
#endif
        for(j = 0, size2 = uint(G2Nv[i].size()); j != size2; j++)
            for(size_t k = 0, size3 = uint(G2Nv[i][j].size()); k != size3; k++)
                if (std::find(ctVertices.begin(), ctVertices.end(), G2Nv[i][j][k]) != ctVertices.end()) {
                    G2Nv[i][j].erase(G2Nv[i][j].begin() + k);
                    k--;
                    size3--;
                }
#ifdef _DEBUG
        std::cout << "G2G1correspondance" << std::endl;
#endif
    }
}*/


// Find complementary graph G2 of a sub-graph from G1
/*void g2ComplementaryG1 (const std::vector < unsigned int >& V1, // G1 sub-graph vertices
                        const std::vector < std::pair < unsigned int, unsigned int > >& E2, // E2 edges
                        std::vector < unsigned int >& V2Comp, // complementary sub-graph vertices
                        std::vector < std::pair < unsigned int, unsigned int > >& E2Comp) // complementary sub-graph edges
{
    V2Comp = V1; // Vertices are the same
    for(unsigned int i = 0, size = uint(E2.size()); i != size; i++)
    {
        if(std::find(V2Comp.begin(), V2Comp.end(), E2[i].first) != V2Comp.end()
                or std::find(V2Comp.begin(), V2Comp.end(), E2[i].second) != V2Comp.end())
        {
            E2Comp.push_back(E2[i]);
        }
    }
}*/

// bool compatibility between two structures
/*bool compatibleTF(Structure str1, Structure str2) {
    bool match = true;

    if( str1.get_rna_() == str2.get_rna_())
    {
       size_t i, size;
       std::vector < std::pair < unsigned int, unsigned int > > l1 = str1.get_listBp_(), l2 = str2.get_listBp_();
       for(i = 0, size = l1.size(); i != size and match; i++)
       {
           std::vector < unsigned int > coord = find_bp_with_i(l2, l1[i].first);
           if (!coord.empty())
           {
               match = false;
           }
           else
           {
               coord = find_bp_with_i(l2, l1[i].second);
               if (!coord.empty())
                   match = false;
           }
       }
    }
    return match;
}*/


// Match between helix and a structure
/*bool compatible(Structure str, Helix helix) {
    bool res = true;
    size_t i = 0, size = str.get_motifs_().size();
    while (res and i != size) {
        if (Helix * h = dynamic_cast < Helix * > (str.get_motifs_()[i])) {
            res = Helix::match(helix, *h);
        }
        i++;
    }
    return res;
}*/

// Two structures form a pseudoknot ?
/*bool arePseudo(Structure str1, Structure str2) {
    bool res = true;
    size_t i = 0, size1 = str1.get_motifs_().size(), j = 0, size2 = str2.get_motifs_().size();
    while (res and i != size1) {
        j = 0;
        while (res and j != size2) {
            if (Helix * h1 = dynamic_cast < Helix * > (str1.get_motifs_()[i]))
                if (Helix * h2 = dynamic_cast < Helix * > (str2.get_motifs_()[j]))
                    if(Helix::match(*h1, *h2))
                        res = Helix::arePseudo(*h1, *h2);
            j++;
        }
        i++;
    }
    return res;
}*/

// Helix and a structure form a pseudoknot ?
/*bool arePseudo(Structure str, Helix helix) {
    bool res = true;
    size_t i = 0, size = str.get_motifs_().size();
    while (res and i != size) {
        if (Helix * h = dynamic_cast < Helix * > (str.get_motifs_()[i])) {
            res = Helix::arePseudo(helix, *h);
        }
        i++;
    }
    return res;
}*/