#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 <sys/stat.h>
#include <string.h>
#include <boost/algorithm/string.hpp>
#include <math.h>

#include "Extract_data/fa.h"
#include "rna.h"
#include "structure.h"
#include "solinteraction.h"
#include "complexe.h"
#include "predictor.h"
#include "IP/ipcomp.h"
#include "Extract_data/extract_user_ct.h"
#include "Extract_data/extract_str_inter.h"
#include "Extract_data/extract_probing_data.h"
#include "IP/findkbestparetoset.h"
#include "utils.h"
#include "graph.h"

//#define _DEBUG

//Probleme de la clique et decomposition de graph pour les pseudonoeuds
//Heuristique multicritere
// TODO Verifier la seed de predictor
// TODO Enlever variables calculs temps

// Remove clashes from base pair list and return several lists of base pairs with no clashes
std::vector < std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > > removeClashes(
        std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > listBp, size_t i) {

    std::vector < std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > > listBps, listBps2;

    if (i != listBp.size() -1) {
        // Check if there are clashes
        std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > lbp, listBp2;
        std::vector < int > clashes, clashes2;
        // Look in the listBp[i+1:end] if there are clashes with the base pair i
        std::copy(listBp.begin() + i + 1, listBp.end(),  std::back_inserter(lbp));
        clashes = find_bp_with_i(lbp, listBp[i].first);
        clashes2 = find_bp_with_i(lbp, listBp[i].second);
        clashes.insert(clashes.end(), clashes2.begin(), clashes2.end());
        // If there are clashes
        if (!clashes.empty()) {
            bool consistent = true;
            // Create for each clash a new listBp without the other clashes
            // clashes.push_back(-1);
            for (int j = int(clashes.size())-1, size = -1; j != size; j--) {
                // The bp is consistent if it is part of an helix
                consistent = false;
                for (ulong k = 0, size2 = listBp.size(); k != size2 or !consistent; k++) {
                    if ((listBp[k].first.first == listBp[ulong(clashes[ulong(j)])].first.first
                         and listBp[k].second.first == listBp[ulong(clashes[ulong(j)])].second.first)
                            and (
                         (listBp[k].first.second > listBp[ulong(clashes[ulong(j)])].first.second
                         and listBp[k].second.second < listBp[ulong(clashes[ulong(j)])].second.second)
                            or (listBp[k].first.second < listBp[ulong(clashes[ulong(j)])].first.second
                                and listBp[k].second.second > listBp[ulong(clashes[ulong(j)])].second.second)
                                )) {
                        consistent = true;
                    }
                }
                // If the remaining base pair is consistent, remove clashes in the following
                if (consistent) {
                    listBp2 = listBp;
                    for(int k = int(clashes.size())-1, size2 = -1; k != size2; k--)
                        if (k != j) {
                            listBp2.erase(listBp2.begin() + clashes[ulong(k)] + i+1);
                        }
                    // Remove clashes in the following
                    listBps2 = removeClashes(listBp2, i+1);
                    listBps.insert(listBps.end(), listBps2.begin(), listBps2.end());
                }
            }
        } else {
            listBps = removeClashes(listBp, i+1);
        }
    } else {
        listBps.push_back(listBp);
    }

    return listBps;
}

std::vector < Complexe > convertCliquesToComplexes(
        const std::vector < std::tuple < std::vector < uint >, std::vector < float > > > &cliques0,
        const std::vector < unsigned int > &vertices0,
        const std::vector < Rna > &rnaList_,
        const std::vector < Structure > &strs_,
        const std::vector < SolInteraction > &inters_,
        int lastStr) {

    std::vector < unsigned int > rnaId;
    std::vector < std::string > rnaName;
    std::vector < std::string > seqList;

    std::vector < Complexe > comps;
    std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > listBP;
    std::vector <
            std::vector < std::pair < std::pair < uint, uint > , std::pair < uint, uint > > > > listBPs;
    std::vector < std::pair < uint, uint > > lbp;
    std::vector < uint > clique0;
    size_t i, size, j, size2, k, size3;

    for(i = 0, size = uint(rnaList_.size()); i != size; i++) {
        rnaId.push_back(uint(i));
        rnaName.push_back(rnaList_[i].get_name_());
        seqList.push_back(rnaList_[i].get_seq_());
    }

    for(i = 0, size = uint(cliques0.size()); i != size; i++) {
        listBP.clear();
        lbp.clear();

        clique0 = std::get<0>(cliques0[i]);

        if(!clique0.empty()){
            for(j=0, size2 = uint(clique0.size()); j != size2; j++) {
                if (int(clique0[j]) <= lastStr) {
                    lbp = strs_[vertices0[clique0[j]]].get_listBp_();
                    for (k = 0, size3 = lbp.size(); k != size3; k++)
                        listBP.push_back(
                                    std::make_pair(
                                        std::make_pair(strs_[vertices0[clique0[j]]].get_rna_(),
                                                       lbp[k].first),
                                    std::make_pair(strs_[vertices0[clique0[j]]].get_rna_(),
                                                   lbp[k].second) ));
                } else {
                    if(lastStr != -1){
                        lbp = inters_[vertices0[clique0[j]]-lastStr-1].get_listBp_();
                    }else {
                        lbp = inters_[vertices0[clique0[j]]].get_listBp_();
                    }
                    for (k = 0, size3 = lbp.size(); k != size3; k++) {
                        if(lastStr != -1){
                        listBP.push_back(
                                    std::make_pair(
                                        std::make_pair(inters_[vertices0[clique0[j]]-uint(lastStr)-1].get_rna1_(),
                                                       lbp[k].first),
                                    std::make_pair(inters_[vertices0[clique0[j]]-uint(lastStr)-1].get_rna2_(),
                                                   lbp[k].second) ));
                        } else {
                            listBP.push_back(
                                        std::make_pair(
                                            std::make_pair(inters_[vertices0[clique0[j]]].get_rna1_(),
                                                           lbp[k].first),
                                        std::make_pair(inters_[vertices0[clique0[j]]].get_rna2_(),
                                                       lbp[k].second) ));
                        }
                    }

                }
            }
            std::sort(listBP.begin(), listBP.end());
            // Remove redundance in listBP
            for(j = 0, size2 = listBP.size()-1; j != size2; j++) {
                if (listBP[j] == listBP[j+1]) {
                    listBP.erase(listBP.begin() + j);
                    size2--;
                    j--;
                }
            }

            listBPs = removeClashes(listBP, 0);

            for (j = 0, size2 = listBPs.size(); j != size2; j++)
                comps.push_back(
                            Complexe (rnaId, rnaName, seqList, listBPs[j], std::get<1>(cliques0[i])));
        }
    }

    return comps;
}

// Delete pointers
void del(std::vector < Interloop * > ctIntLoop_,
         std::vector < Hairpinloop * > ctHairLoop_,
         std::vector < Helix * > ctHelix_,
         std::vector < Multiloop * > ctMultiLoop_,
         std::vector < Pseudoknot * > ctPseudo_)
{
    size_t j, size2;
    for (j = 0, size2 = ctIntLoop_.size(); j != size2; j++)
        delete ctIntLoop_[j];
    for (j = 0, size2 = ctHairLoop_.size(); j != size2; j++)
        delete ctHairLoop_[j];
    for (j = 0, size2 = ctHelix_.size(); j != size2; j++)
        delete ctHelix_[j];
    for (j = 0, size2 = ctMultiLoop_.size(); j != size2; j++)
        delete ctMultiLoop_[j];
    for (j = 0, size2 = ctPseudo_.size(); j != size2; j++)
        delete ctPseudo_[j];
}


void outputCliques (
        std::string CLIQUESfile,
        const std::vector < std::tuple < std::vector < uint >, std::vector < float > > > &cliques,
        const std::vector < unsigned int > &vertices,
        const std::vector < float > &verticesWE,
        const std::vector < float > &verticesWP,
        const std::vector < float > &verticesWCt,
        const uint &nbEdges,
        const double &density,
        const double &total_time5
        ) {

    std::ofstream cliquesFile;
    std::vector < std::string > vertices_str;
    std::vector < std::string > cliques_str;

    for(size_t a=0, size = cliques.size(); a != size; a++) {
        vertices_str.clear();
        for(size_t b = 0, size2 = std::get<0>(cliques[a]).size(); b != size2; b++){
            vertices_str.push_back( "(" + std::to_string(std::get<0>(cliques[a])[b] + 1)
                            + ", {" + std::to_string(-verticesWE[std::get<0>(cliques[a])[b]])
                            + ", " +  std::to_string(verticesWP[std::get<0>(cliques[a])[b]])
                            + ", " + std::to_string(verticesWCt[std::get<0>(cliques[a])[b]])
                            + "})" );
        }
        cliques_str.push_back( "{ " + join(vertices_str, " | ") + "}\tweight: {"
                      + std::to_string(-std::get<1>(cliques[a])[1]) + ", "
                      + std::to_string(std::get<1>(cliques[a])[2]) + ", "
                      + std::to_string(std::get<1>(cliques[a])[0]) + "}" );
    }

    cliquesFile.open(CLIQUESfile);
    cliquesFile << "|V|=" << vertices.size() << "\t|E|=" << nbEdges << "\td=" << density << std::endl;
    for(size_t a=0, size = cliques_str.size(); a != size; a++) {
        cliquesFile << cliques_str[a] << std::endl;
    }
    cliquesFile << "found: " << cliques.size() << " cliques\ttook: " << total_time5 << "s" << std::endl;
    cliquesFile.close();
}


// Man RCPred
void usage()
{
    std::cout << "C-RCPred (Constrained RNA Complex Prediction)" << std::endl
              << " -f       Input fasta file" << std::endl
              << " -s       Input secondary structure files," << std::endl
              << " -i       Input interaction secondary structure files," << std::endl
              << " -c       Constraint files where constraints are in the same order than the secondary structure and interaction files" << std::endl
              << " -p       Probing data file for each RNA, in the same order than the fasta files (if no data for an RNA, an empty file is needed)" << std::endl
              << " -o       Format of output \"d\" (dot-parenthesis, default), \"j\" (JSON) or \"b\" (base pair list format)" << std::endl
              << " -e       Energy model (0: No energy model; use energies computed by upstream tools, 1: ViennaRNA package model; default, 2: NUPACK model)" << std::endl
              << " -t       Threshold for compatibility RNA (0->100, default 100)" << std::endl
              << " -l       Lower probing threshold (0->100, if not specified, probing values between 0% and the upper threshold once normalized will not be taken into account)" << std::endl
              << " -u       Upper probing threshold  (0->100, if not specified, probing values between the lower threshold and 100% once normalized will not be taken into account)" << std::endl
              << " -g       Output graph file" << std::endl
              << " -k       Output cliques predicted" << std::endl
              << " -q       Output cliques predicted without duplicates" << std::endl
              << " -m       Maximum number of structures to output" << std::endl
              << " -h       Shows this message" << std::endl;
}

// Main
int main(int argc, char *argv[]) {
    std::string SEQfile = "", STRfile = "", INTfile = "", PROBINGfile = "", CTfile = "", printing = "d", GRAPHfile = "", CLIQUESfile = "", CLIQUESfile2 = "";
    unsigned int T = 1000, I = 2000, energyModel = 1;
    float threshold = 100, lowerThresProbing = 20, upperThresProbing = 80;
    float alphas = 0.7, alphar = 0.92, Po = 0.75;
    uint maxSol = 50;


    if(argc > 1) {
        // Parse options
        char ch;
        while ((ch=char(getopt(argc, argv, "hf:s:i:c:o:t:e:p:l:u:g:k:q:m:?")))!=-1) {
            switch (ch) {
            case 'f':
                SEQfile = optarg;
                break;
            case 's':
                STRfile = optarg;
                break;
            case 'i':
                INTfile = optarg;
                break;
            case 'c':
                CTfile = optarg;
                break;
            case 'o':
                printing = optarg;
                break;
            case 'g':
                GRAPHfile = optarg;
                break;
            case 'k':
                CLIQUESfile = optarg;
                break;
            case 'q':
                CLIQUESfile2 = optarg;
                break;
            case 'm':
                maxSol = uint(atoi(optarg));
                break;
            case 'e':
                energyModel = uint(atoi(optarg));
                if (energyModel != 0 and energyModel != 1 and energyModel != 2 and energyModel != 3) {
                    std::cout << "-e: value must be 0, 1, 2 or 3." << std::endl;
                    usage();
                    return EXIT_SUCCESS;
                }
                break;
            case 't':
                threshold = atof(optarg);
                if (threshold < 0 or threshold > 100) {
                    std::cout << "-t: value must be between 0 and 100." << std::endl;
                    usage();
                    return EXIT_SUCCESS;
                }
                break;
            case 'p':
                PROBINGfile = optarg;
                break;
            case 'l':
                if (atof(optarg) < 0 or atof(optarg) > 100) {
                    std::cout << "-l: value must be between 0 and 100." << std::endl;
                    usage();
                    return EXIT_SUCCESS;
                }
                else {
                  lowerThresProbing = atof(optarg) ;
                }
                break;
            case 'u':
                if (atof(optarg) < 0 or atof(optarg) > 100) {
                    std::cout << "-u: value must be between 0 and 100." << std::endl;
                    usage();
                    return EXIT_SUCCESS;
                }
                else {
                  upperThresProbing = atof(optarg);
                }
                if (upperThresProbing <= lowerThresProbing) {
                    std::cout << "-u: value must be greater than -l value." << std::endl;
                    usage();
                    return EXIT_SUCCESS;
                }
                break;
            case 'h': case '?': default:
                std::cout << "help" << std::endl;
                usage();
                return EXIT_SUCCESS;
            }
        }
    } else {
        usage();
        return EXIT_SUCCESS;
    }



    if(SEQfile == "") {
        std::cout << "Fasta file is missing." << std::endl;
        return EXIT_FAILURE;
    } else {
        std::vector < Rna > rnaList_; // Vector of RNA sequences
        std::vector < Structure > strs_; // Vector of structures
        std::vector < SolInteraction > inters_; // Vector of interactions
        std::vector < std::vector < Helix > > helices_; // Vector of helices for each pair of RNA (G2 graph)

        std::vector< std::tuple < std::pair < uint, uint >, std::pair < uint, uint >, uint > > ctAppTwo_; // Vector of Weak constraints of pairs of bases that must be paired together
        std::vector< std::tuple < uint, uint, uint > > ctApp_;   // Vector of Weak constraints of bases that must be paired of a complex
        std::vector< std::tuple < uint, uint, uint > > ctSingle_;   // Vector of Weak constraints of bases that must not be paired of a complex
        std::vector < Interloop * > ctIntLoop_; // Vector of Weak constraints of internal loops
        std::vector < Hairpinloop * > ctHairLoop_; // Vector of Weak constraints of hairpin loops
        std::vector < Helix * > ctHelix_;   // Vector of Weak constraints of helixs
        std::vector < Multiloop * > ctMultiLoop_;   // Vector of Weak constraints of multiloops
        std::vector < Pseudoknot * > ctPseudo_;   // Vector of Weak constraints of pseudoknots
        unsigned long i, j, size, size2;

        // Extracting rnas from the fasta file
        Fasta fa;
        try {
            fa = Fasta();
            fa.load(SEQfile.c_str());
            for(j = 0, size2 = uint(fa.seq().size()); j != size2; j++) {

                std::string seq = Rna::format(fa.seq()[j]);
                if (Rna::check_seq(seq)) {
                    rnaList_.push_back(Rna(fa.name()[j], seq));
                } else {
                    throw std::string("The sequence is not supported by the software (allowed nucleosides: A, U, G, C, I, Q, D, T, S, X, B, O, N, R, Y).");
                }
            }

            // Set id of each RNA
            for(i = 0, size = uint(rnaList_.size()); i != size; i++)
                rnaList_[i].set_id_(uint(i));

            // Check number of structure files
            if(STRfile == "") {
                std::cout << "Secondary structure file is missing.";
                usage();
                return EXIT_FAILURE;
            } else {
                // Check number of interaction files
                if(INTfile == "") {
                    std::cout << "Interaction file is missing." << std::endl;
                    usage();
                    return EXIT_FAILURE;
                } else {

                    // Extracting constraints
                    uint nbHardCT =0;
                    if (CTfile != "") {
                        extract_ct(CTfile, rnaList_, ctAppTwo_, ctApp_,
                                        ctSingle_, ctIntLoop_, ctHairLoop_,
                                        ctHelix_, ctMultiLoop_, ctPseudo_, nbHardCT);
                    }

                    if(!PROBINGfile.empty()) {
                        extract_probing(PROBINGfile, rnaList_);
                        if (lowerThresProbing == -10 and upperThresProbing == 110) {
                            // If there is probing data but the user did not specify any probing threshold,
                            // a default threshold is applied
                            lowerThresProbing = 5;
                            upperThresProbing = 95;
                        }
                    }

                    // Extracting secondary structures
                    extract_str(STRfile, rnaList_, strs_, ctAppTwo_, ctApp_,
                                ctSingle_,  ctIntLoop_, ctHairLoop_, ctHelix_,
                                ctMultiLoop_,  ctPseudo_, energyModel, nbHardCT,
                                lowerThresProbing, upperThresProbing);

                    // Extracting interactions
                    extract_inter(INTfile, rnaList_, inters_, helices_,
                                  ctAppTwo_, ctHelix_, ctPseudo_, energyModel, nbHardCT,
                                  lowerThresProbing, upperThresProbing);


                    // MODELIZE GRAPHS
                    if(int(rnaList_.size()) > 1) {
                        // G1 GRAPH
                        std::vector < unsigned int > vertices0;
                        std::vector < float > verticesWCt; // User constraint weight
                        std::vector < float > verticesWP; // Probing weight
                        std::vector < float > verticesWE; // Energy weight
                        // Vector of linked vertices for any vertex
                        std::vector < std::vector < unsigned int > > Nv0;
                        // Vector of no linked vertices for any vertex
                        std::vector < std::vector < unsigned int > > NvB0;

                        // G2 GRAPHS
                        // Vertex vectors of each G2 graphs
                        // Vertices are sorted by the beginning of each interval
                        std::vector < std::vector < unsigned int > > G2vertices0;
                        // Nv vectors of each G2 graphs
                        std::vector < std::vector < std::vector < unsigned int > > > G2Nv0;
                        // G1 vertex correspondance vector of each G2 graph
                        std::vector < std::vector < unsigned int > > G2G1correspondance0;

                        // Constraints vertices for the two graphs
                        std::vector < unsigned int > ctVertices, rmVertices;

                        uint nbEdges;
                        int lastStr;
                        modelGraph(rnaList_, strs_, inters_, helices_, vertices0,
                                        verticesWCt, verticesWP, verticesWE, Nv0,
                                        NvB0, G2vertices0, G2Nv0, G2G1correspondance0,
                                        threshold, SEQfile, nbEdges, lastStr);

                        // G1 GRAPH
                        std::vector < unsigned int > vertices = vertices0;
                        std::vector < std::vector < unsigned int > > Nv = Nv0;
                        std::vector < std::vector < unsigned int > > NvB = NvB0;
                        // G2 GRAPHS
                        std::vector < std::vector < unsigned int > > G2vertices = G2vertices0;
                        std::vector < std::vector < std::vector < unsigned int > > > G2Nv = G2Nv0;
                        std::vector < std::vector < unsigned int > > G2G1correspondance = G2G1correspondance0;

                        if (GRAPHfile != "") {
                            // OUTPUT GRAPH IN FILE
                            std::ofstream graphFile;
                            graphFile.open(GRAPHfile);

                            // Nb vertices
                            graphFile << "i\t" << vertices.size() << "\tw\t3" << std::endl;

                            // Structures vertices
                            for(size_t c = 0, size = strs_.size(); c != size; c++) {
                                graphFile << std::to_string(-verticesWE[c]) << "\t";
                                graphFile << verticesWP[c] << "\t";
                                graphFile << verticesWCt[c];
                                graphFile << std::endl;
                            }

                            // Interactions vertices
                            for(size_t c = 0, size = inters_.size(); c != size; c++) {
                                graphFile << std::to_string(-verticesWE[c+ (lastStr + 1)]) << "\t";
                                graphFile << verticesWP[c+ (lastStr + 1)] << "\t";
                                graphFile << verticesWCt[c+ (lastStr + 1)];
                                graphFile << std::endl;
                            }

                            // Edges
                            for(size_t c = 0, size = Nv.size(); c != size; c++) {
                                for(size_t d = 0, size2 = Nv[c].size(); d != size2; d++) {
                                    if(Nv[c][d] > c) {
                                        graphFile << "e\t" << c + 1 << "\t" << Nv[c][d] + 1 << std::endl; // + 1 to make vertices start with 1
                                    }
                                }
                            }

                            graphFile.close();
                        }


                        if(int(vertices.size()) /*>= int(rnaList_.size())*/ > 0) {

                            time_t s5, e5;
                            double total_time5;
                            s5 = clock();

                            // HEURISTIC
                            Predictor predictor = Predictor(vertices, verticesWCt, verticesWP, verticesWE,
                                                            Nv, NvB,
                                                            G2vertices, G2Nv, G2G1correspondance);

                            double density = 2*double(nbEdges) /( double(vertices.size()) * (double(vertices.size()) - 1.0)); // False because nbEdges always equal to 0, nbEdges not updated in modelGraph
                            I = uint(40488.27*density); // I not used in predictor.predict(), not necessary and False because density is False

                            std::vector < std::tuple < std::vector < uint >, std::vector < float > > > cliques;

                            // Default pseudoknot Level = 3
                            bool paretoRespectHardCT = false;
                            uint loop = 0;
                            while(!paretoRespectHardCT and loop<2) {
                                predictor = Predictor(vertices, verticesWCt, verticesWP, verticesWE,
                                                            Nv, NvB,
                                                            G2vertices, G2Nv, G2G1correspondance);
                                predictor.predict(T, alphas, alphar, Po, 3, I, uint(time(NULL)), nbHardCT);
                                cliques = predictor.get_best_cliques_();
                                for(uint i = 0, size = cliques.size(); i != size; i++) {
                                    if (std::get<1>(cliques[i])[0] >= 100*nbHardCT)
                                        paretoRespectHardCT=true;
                                }
                                loop++;
                            }

                            e5 = clock();
                            total_time5 = double(e5-s5)/double(CLOCKS_PER_SEC);

                            cliques = predictor.get_best_cliques_();

                            nbEdges = 0;
                            for(size_t a = 0, size = Nv.size(); a != size; a++) {
                                for(size_t b = 0, size2 = Nv[a].size(); b != size2; b++) {
                                    if(Nv[a][b] > a) {
                                        nbEdges = nbEdges + 1;
                                    }
                                }
                            }
                            density = 2*double(nbEdges) / (double(vertices.size())*(double(vertices.size())-1.0));


                            if(!cliques.empty()) {

                                if(CLIQUESfile != "") {
                                    outputCliques(CLIQUESfile, cliques, vertices, verticesWE, verticesWP, verticesWCt, nbEdges, density, total_time5);
                                }

                                // CONVERT CLIQUES TO COMPLEXES
                                std::vector < std::tuple < std::vector < uint >, std::vector < float > > > cliques0 = cliques;
                                std::vector < unsigned int > clique, clique0, intersec;
                                std::vector < float > weights;

                                std::vector < Complexe > comps = convertCliquesToComplexes(cliques0, vertices0,
                                        rnaList_, strs_, inters_, lastStr);


                                if (comps.size() != 0) {

                                    // REMOVE DUPLICATE COMPLEXES
                                    Complexe c;
                                    std::vector < Complexe > comps2; // comps without duplicates
                                    std::vector < std::tuple < std::vector < uint >, std::vector < float > > > cliques2; // cliques corresponding to comps2
                                    comps2.clear();
                                    bool duplicate;

                                    for(i = 0; i != comps.size(); i++) { // load listsBP_ for each complex
                                        comps[i].loadListsBP_();
                                    }
                                    for(i = 0; i != comps.size(); i++) {
                                        duplicate = false;
                                        for(j = 0; j != comps2.size() and duplicate == false; j++) {
                                            if (c.isEquivalent(comps[i], comps2[j])) { // if the 2 complexes are the same
                                                duplicate = true;
                                            }
                                        }
                                        if (duplicate == false) {
                                            comps2.push_back(comps[i]);
                                            cliques2.push_back(cliques0[i]);
                                        }
                                    }

                                    if(CLIQUESfile2 != "") {
                                        outputCliques(CLIQUESfile2, cliques2, vertices, verticesWE, verticesWP, verticesWCt, nbEdges, density, total_time5);
                                    }

                                        if(printing == "j") {
                                            for(i = 0; i != comps2.size() and i < maxSol; i++) {
                                                std::cout << comps2[i].to_Json() << std::endl;
                                            }
                                        } else if(printing == "b") {
                                            for(i = 0; i != comps2.size() and i < maxSol; i++) {
                                                std::cout << comps2[i].to_forna() << std::endl;
                                            }
                                        } else {
                                            for(i = 0; i != comps2.size() and i < maxSol; i++) {
                                                std::cout << comps2[i].to_string() << std::endl;
                                            }
                                        }

                                } else {
                                    std::cout << "No solution available." << std::endl;
                                }

                                del(ctIntLoop_, ctHairLoop_, ctHelix_, ctMultiLoop_, ctPseudo_);
                                return EXIT_SUCCESS;
                            } else {
                                del(ctIntLoop_, ctHairLoop_, ctHelix_, ctMultiLoop_, ctPseudo_);
                                std::cout << "No solution available for this seed with the given structures and interactions. Please try again." << std::endl;
                                return EXIT_FAILURE;
                            }

                        } else {
                            del(ctIntLoop_, ctHairLoop_, ctHelix_, ctMultiLoop_, ctPseudo_);
                            std::cout << "No solution available with the given structures and interactions." << std::endl;
                            return EXIT_FAILURE;
                        }
                    } // End of modelize graphs
                } // End check of number of interaction files
            } // End check of number of structure files

            del(ctIntLoop_, ctHairLoop_, ctHelix_, ctMultiLoop_, ctPseudo_);

        } catch (const std::string& str) {
            std::cout << str << std::endl;
            del(ctIntLoop_, ctHairLoop_, ctHelix_, ctMultiLoop_, ctPseudo_);
            return EXIT_FAILURE;
        }
    } // End of fasta file detected

    return EXIT_SUCCESS;
}