#include "coptimizer.h"

#include <cinstance.h>
#include <cinstancereduite.h>
#include <cgraphe.h>
#include <tuning.h>
#include <indices.h>
#include <selecteur_1.h>
#include <selecteur_thermiques.h>
#include <selecteur_scenario.h>
#include <selecteur_scenario_1thermique.h>
#include <filtreAval.h>
#include <filtreavalcplex.h>
#include <filtreavalcplexrestelibre.h>
#include <filtreavalcplexultime.h>
#include <filtreavalcplexultime_1passe.h>
#include <vector>
#include <ctime>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <fstream>
#include <unistd.h>
#include <pthread.h>

#include "DA.h"
#include "cdualcplex.h"
#include "cdualcplexscenariomoyen.h"
#include "cdualcplexunecolonnethermique.h"
#include "cdualcplexcolonnesthermiques.h"
#include "cdualcplexdiagsthermiques.h"

using std::vector;
using namespace std;

COptimizer::COptimizer( bool multithreading )
{
    multithread             = multithreading ;
    pourcentage_de_temps_1  = 0         ;
    pourcentage_de_temps_2  = 0         ;
    duree_cplex             = 5         ;
    heure_debut_passe       = 0         ;
    duree_passe_precedente  = 0         ;
    delta_duree_passe       = 0         ;

    coeff_recharges = 206613465.3984422 ;
    coeff_scenarios = 26000 ;
    if (multithreading) {
        coeff_selecteur = 540 ;
    }
    else {
        coeff_selecteur = 298 ;
    }
    prev_selecteur = (double) I() * T()                 / coeff_selecteur * 1000000.0 / (double) CLOCKS_PER_SEC ;
    prev_recharges = (double) I() * T() * S() * S()     / coeff_recharges * 1000000.0 / (double) CLOCKS_PER_SEC ;
    prev_scenarios = (double) (I() + J()) * T() * S()   / coeff_scenarios * 1000000.0 / (double) CLOCKS_PER_SEC ;

    duree_filtre_aval       = std::min(3 * (double)CInstance::params().getTimeLimit() / 4,
                                        prev_selecteur + prev_recharges + prev_scenarios) ;


    m_passes                        = 0         ;

    nb_colonnes_cout_non_positif    = 0         ;
    nb_rejets_couts_positifs        = 0         ;
    nb_rejets_existants             = 0         ;

    m_dualCalc                      = NULL      ;

    m_cout_solution         = (euros) INFINI ;
    m_prochaine_suppression = 2 * (H() + I()) ;

#ifdef VERBOSE_CPLEX_INTERNAL
    m_cplex_internal_silent = false     ;
#else
    m_cplex_internal_silent = true      ;
#endif
#ifdef VERBOSE_SOLUTION_STATUS
    m_duals_silent = false              ;
#else
    m_duals_silent = true               ;
#endif
#ifdef WRITE_CPLEX_FILE
    m_cplex_writeToFile = true          ;
#else
    m_cplex_writeToFile = false         ;
#endif
}


COptimizer::~COptimizer()
{
    if (m_dualCalc!=NULL)
        delete m_dualCalc ;
}




struct PredAjoutColonne {
    CDualCalc  * m_da ;
    COptimizer * m_optimizer ;
    PredAjoutColonne( COptimizer * opt, CDualCalc * da ) { m_optimizer=opt; m_da=da; }
    void operator()(CGraphe * g);
};

void verbose_rejet_plan_existant( CGraphe * g);
void verbose_rejet_cout_positif( CGraphe * g);




/**
 *  calcul des nouvelles colonnes réparti sur plusieurs threads
 * ( 1 thread par cpu disponible )
 */

// structure partagée par les différents threads
struct CIterThreadGraph {
    vector<CGraphe *>  graphes ;        // liste des graphes à traiter
    vector<CGraphe *>::iterator it ;    // itérateur sur le prochain graphe à traiter
    pthread_mutex_t lock_it ;           // mutex pour bloquer l'itérateur quand un thread "se sert"
    CIterThreadGraph( vector<CGraphe *> & _graphes ) : graphes(_graphes)
    {
        it = graphes.begin();
        pthread_mutex_init(&lock_it, NULL);
    }
};
// fonction "start" d'un thread
void * threaded_creer_colonnes( void * /* (CIterThreadGraph *) */ iter_graphes )
{
    CIterThreadGraph * itg = (CIterThreadGraph *) iter_graphes ; // itérateur sur les graphes à traiter
    while (true) {
        CGraphe * graphe = NULL ;
        // recherche du prochain graphe à traiter
        pthread_mutex_lock( & itg->lock_it );
        if ( itg->it != itg->graphes.end() ) {
            graphe = * itg->it ;
            itg->it ++ ;
        }
        pthread_mutex_unlock( & itg->lock_it );
        // fin du thread s'il n'y a plus de graphes à traiter
        if (graphe==NULL) break ;
        // création de colonne s'il y a un graphe à traiter
        CGraphe::creer_colonnes( graphe );
    }
    return NULL ;
}
// création de colonnes sur plusieurs threads
void COptimizer::creer_colonnes_multithread()
{
    // nombre de cpus
    int  ncpus = sysconf(_SC_NPROCESSORS_CONF);
    // autant de threads que de cpus
    vector<pthread_t> threads(ncpus);

    // lancement des différents threads
    CIterThreadGraph threaded_graphes(graphes) ;
    for (vector<pthread_t>::iterator th=threads.begin(); th!=threads.end(); th++)
    {
        pthread_create( &(*th), NULL, threaded_creer_colonnes, (void *) &threaded_graphes );
    }
    // attente de la fin des threads lancés
    for (vector<pthread_t>::iterator th=threads.begin(); th!=threads.end(); th++)
    {
        pthread_join( *th, NULL );
    }
}




/**
 *               Algorithme principal de l'optimisation
 */

void COptimizer::start()
{

  _LBvalue = INFINI ;

#ifdef VERBOSE_TEMPS
    std::cout
            << "\t" << "Magnitude_I " << I() << std::endl
            << "\t" << "Magnitude_J " << J() << std::endl
            << "\t" << "Magnitude_K " << K() << std::endl
            << "\t" << "Magnitude_T " << T() << std::endl
            << "\t" << "Magnitude_H " << H() << std::endl
            << "\t" << "Magnitude_S " << S() << std::endl
            << std::endl ;

    std::cout
            << "Temps pour selecteur + filtre aval :\t" << duree_filtre_aval << std::endl
            << "Temps pour le sélecteur            :\t" << prev_selecteur       << std::endl
            << "Temps pour les recharges           :\t" << prev_recharges       << std::endl
            << "Temps pour les scénarios           :\t" << prev_scenarios       << std::endl
            << std::endl ;
#endif
    // Duals
    switch (CInstance::params().getDualMethod()) {
        case d_dualascent :
            m_dualCalc = new DualAscent( CInstanceReduite::getEDEM(), I() );
//             m_dualCalc = new DualAscent( CInstanceReduite::getEDEM_min(), I() );
            break;
        case d_cplexduals_scenarios :
            m_dualCalc = new CDualCplex(m_cplex_internal_silent, m_duals_silent, m_cplex_writeToFile);
            break;
        case d_cplexduals : // default
            m_dualCalc = new CDualCplexScenarioMoyen(m_cplex_internal_silent, m_duals_silent, m_cplex_writeToFile);
            break;
        case d_cplexduals_1ct :
            m_dualCalc = new CDualCplexUneColonneThermique(m_cplex_internal_silent, m_duals_silent, m_cplex_writeToFile);
            break;
        case d_cplexduals_ct :
            m_dualCalc = new CDualCplexColonnesThermiques(m_cplex_internal_silent, m_duals_silent, m_cplex_writeToFile);
            break;
        case d_cplexduals_dt :
            m_dualCalc = new CDualCplexDiagsThermiques(m_cplex_internal_silent, m_duals_silent, m_cplex_writeToFile);
            break;
    }
    u = u_initial();

    /**
     *              CREATION DES GRAPHES
     */

    debut_optimization = time(NULL);
    /* pour chaque centrale, créer un graphe */
    temps1 = clock() ;
    for (powerplants2_t i=(powerplants2_t) 0; i<I(); i = (powerplants2_t) (i+1))
    {
        graphes.push_back( new CGraphe(i) );
    }
    temps2 = clock() ;
    verbose_temps_graphes();


    /**
     *              CREATION DES COLONNES
     */

    /* tant qu'il reste du temps */
    while (prepare_strategie())
    {
        verbose_debut_passe();

        if (m_passes++ > 0 and strategie_dual_ascent)
        {
            verbose_avant_dual_ascent();
            vector<euros_par_kWh> new_u = m_dualCalc->calculeDuals(); // voir DA.cpp
            double temp = m_dualCalc->getObjectiveValue() ;
            std::cout << "m_dualCalc->getObjectiveValue() : " << temp << std::endl ;
            if (_LBvalue >= temp+0.1) // pour le log
            {
                _LBvalue = temp;
                _nbIterToLB = m_passes;
            }
            m_cr = m_dualCalc->getReducedCosts();
            verbose_apres_dual_ascent(new_u);
            u = new_u;
            /* suppression des colonnes à coût réduit positif */
            /*vector<euros> cr = m_dualCalc->getReducedCosts();
            for (int p=CPlan::getNbPlans()-1; p>=0; p--)
            {
//                 euros & _cout_reduit = CPlan::getPlan((size_t)p)->m_cout_reduit ;
                euros _cout_reduit = CPlan::getPlan((size_t)p)->m_cout_reduit ;
                cr[p] = CPlan::getPlan((size_t)p)->recalcule_cout_reduit();
                if ( (euros)  -1 < _cout_reduit  and cr[p] <= (euros) -1 and _cout_reduit.v != cr[p].v)
                    cout << "le coût réduit de la colonne du plan " << CPlan::getPlan((size_t)p)->m_id
                            << " est devenu négatif : (" << cr[p] << " < " << _cout_reduit << ")" << endl;

                if (_cout_reduit <= (euros) 0 and (euros)1 <= cr[p] and _cout_reduit.v != cr[p].v)
                    cout << "coût réduit devenu positif sur la colonne du plan " << CPlan::getPlan((size_t)p)->m_id <<
                            " (" << cr[p] << " > " << _cout_reduit << ")" << endl;

                _cout_reduit = cr[p] ;
            }
            */
            //suppression_colonnes(graphes, maintenant);
        }

        /* passage des variables duales aux graphes */
        verbose_avant_viterbi();
        CGraphe::prepare_viterbi(&u, /*true*/strategie_dual_ascent, strategie_utilise_contraintes_de_cliques) ;

        /* générer les nouvelles colonnes pour chaque centrale */
        verbose_avant_creer_colonnes();
        if (multithread)
            creer_colonnes_multithread() ;
        else
            std::for_each( graphes.begin(), graphes.end(), CGraphe::creer_colonnes );
        verbose_apres_creer_colonnes();

        /* ajout des colonnes à la liste des colonnes */
        size_t nbPlansAvantAjouts = CPlan::getNbPlans() ;
        PredAjoutColonne ajoute(this, m_dualCalc);
        std::for_each( graphes.begin(), graphes.end(), ajoute );
        for (size_t p=nbPlansAvantAjouts; p<CPlan::getNbPlans(); p++)
            CConstraint::ajoute_plan_contraintes_13_a_21(CPlan::getPlan(p));

        /* production éventuelle d'un fichier solution */
        vector<const CPlan *> solution;
        produit_solution(graphes, maintenant, solution, true, m_passes);

        /* gestion du temps */
        calcule_temps_fin_passe();
    }
    _nbTotalIter = m_passes ;

    /**
     * calcul et sortie d'une solution
     */
    _nbPlans = CPlan::getNbPlans();
    vector<const CPlan *> solution;
//     produit_solution(graphes, maintenant, solution, true);
    _FinalValue = produit_meilleure_solution(graphes, maintenant, solution, true).v;

    // libération de mémoire
    std::for_each(graphes.begin(), graphes.end(), CGraphe::detruit);
    CPlan::libere_plans();
}


/**
 *              suppression de colonnes à coût réduit positif
 */

void COptimizer::suppression_colonnes(vector<CGraphe *> & graphes, time_t & maintenant)
{
    if (CPlan::getNbPlans() < m_prochaine_suppression)
        return ;

    m_prochaine_suppression = CPlan::getNbPlans() * 5 ; //+ H() + I() ;

    cout << "Suppression de colonnes. Avant : " << CPlan::getNbPlans() << endl;
    // calcul d'une solution
    vector<const CPlan *> plans;
    euros ancien_cout_solution = m_cout_solution;
    produit_solution(graphes, maintenant, plans, false) ;

    cout << "Suppression de colonnes. Coût solution = " << m_cout_solution << " (ancien cout = "<< ancien_cout_solution << ")" << endl;
    if (ancien_cout_solution <= m_cout_solution  or m_cout_solution<=(euros)0)
    {
        return ;
    }


    // calcul et enregistrement des nouveaux coûts réduits
    vector<euros> cr = m_dualCalc->getReducedCosts();
    for (size_t p=0; p < CPlan::getNbPlans() ; p++)
        CPlan::getPlan(p)->m_cout_reduit = cr[p] ;

    // suppression de colonnes à coût réduit positif ne faisant pas partie de la solution
    for (size_t p=CPlan::getNbPlans()-1; (signed) p >= 0 /*and CPlan::getNbPlans()>=(unsigned)H()+I()*/; p--)
    {
        if (m_plans_solution.find(CPlan::getPlan(p))!=m_plans_solution.end()) continue ;
//         if (CPlan::getPlan(p)->recalcule_cout_reduit()>(euros)0)
        if (CPlan::getPlan(p)->m_cout_reduit>(euros)0)
        {
            m_dualCalc->retraitColonne(CPlan::getPlan(p));
            CPlan::libere_plan(p);
        }
    }
    cout << "Suppression de colonnes. Après : " << CPlan::getNbPlans() << endl;
    size_t nb_cr_neg = 0 ;
    for (size_t p=0; p < CPlan::getNbPlans() ; p++)
        nb_cr_neg += CPlan::getPlan(p)->m_cout_reduit <= (euros)0 ? 1 : 0 ;
    cout << "Suppression de colonnes. Nb de colonnes à c.r.<=0 : " << nb_cr_neg << endl;


    m_prochaine_suppression = CPlan::getNbPlans() * 5 ; //+ H() + I() ;
}

/**
 *
 *              condition pour continuer l'optimisation
 *
 */
bool COptimizer::continue_passes()
{
    if (CInstance::params().getLimitType()==CParams::TIME_LIMIT)
    {
#ifdef VERBOSE_PASSES
    std::cout
        << "planification de passes : duree_passe_precedente=" << duree_passe_precedente
        << " + delta_duree_passe="<< delta_duree_passe << " + duree_cplex="<< duree_cplex
        << " + duree_filtre_aval="<< duree_filtre_aval
        << " = " << duree_passe_precedente + delta_duree_passe + duree_cplex + duree_filtre_aval << std::endl ;
    std::cout << "planification de passes : duree_restante=" << CInstance::params().getDeadline() - time(NULL) << std::endl;
#endif
        return (time(NULL) + duree_passe_precedente + delta_duree_passe + duree_cplex + duree_filtre_aval
                < CInstance::params().getDeadline());
    }
    else
        return (m_passes < (signed) CInstance::params().getPassLimit());
}


/**
 *
 *          ajout d'une colonne, sous conditions
 *
 */

void PredAjoutColonne::operator()(CGraphe * g)
{
    // comptage du nombre de colonnes à coup réduit positif
    std::cout << " plan centrale " << g->m_i << ",\tc.r.=" << g->m_plan->m_cout_reduit << std::endl;
    // test sur le coût du plan à ajouter
    if (m_optimizer->m_passes!=1
        and m_optimizer->strategie_rejete_couts_reduits_positifs
        and m_optimizer->strategie_dual_ascent
        and g->m_plan->m_cout_reduit > (euros)10
        )
    {
        // rejet
        m_optimizer->nb_rejets_couts_positifs ++ ;
        verbose_rejet_cout_positif(g);
        return ;
    }
    // test sur l'existence du même plan dans la base
    if (CPlan::retiens_plan(g->m_plan, m_optimizer->strategie_accepte_plan_existant))
    {
        if (g->m_plan->m_cout_reduit < (euros) 10)
            m_optimizer->nb_colonnes_cout_non_positif ++ ;
        m_da->ajoutcolonne(g->m_plan);
        g->accepte_colonne();
    }
    else
    {
        // rejet
        m_optimizer->nb_rejets_existants ++ ;
        verbose_rejet_plan_existant( g );
    }
}


/**
 *
 *          calcul des u_i initiaux
 *
 */

vector<euros_par_kWh> & COptimizer::u_initial() const
{
    static vector<euros_par_kWh> s_u( H()+I(), (euros_par_kWh) 0 );
#ifndef U_INITIAL_ZERO
    euros_par_kWh cout_max(0);
    for (powerplants2_t i=(powerplants2_t) 0; i<I(); ++i)
    {
        for (campaign_t k(0); k<type2(i).K; ++k) {
            cout_max = std::max( cout_max, C(i,k) );
        }
    }
    for (week_t h(0); h<H(); ++h)
        s_u[h] = cout_max + (euros_par_kWh)1 ;
#endif
    return s_u ;
}



/**
 *
 *          Production d'une solution
 *
 */

euros COptimizer::produit_solution( vector<CGraphe *> & graphes, time_t & maintenant, vector<const CPlan *> & plans, bool sur_disque, int passes )
{
#ifndef PRODUIT_SOLUTION_A_CHAQUE_PASSE
    if (passes!=-1) return (euros) -1;
#endif

    euros cout_solution(-1);

    /**
     *              APPEL A CPLEX
     */
    temps1 = clock() ;
    Selecteur * selecteur = NULL ;
    switch (CInstance::params().getSelecteur()) {
        case sel_1ct :
            selecteur = new Selecteur_1;
            break;
        case sel_ct :
            selecteur = new Selecteur_Thermiques;
            break;
        case sel_1cts :
            selecteur = new Selecteur_Scenario_1Thermique;
            break;
        case sel_cts :
            selecteur = new Selecteur_Scenario;
            break;
    }
    selecteur->setOption( multithread ? CPlexModel::Opportunistic : CPlexModel::None ) ;
//     selecteur->setOption( CPlexModel::Time, prev_selecteur ) ;

    vector<size_t> x ;

#ifndef NO_CPLEX
    try {
        x = selecteur->resolution();
        temps2 = clock() ;
        _MIPvalue = selecteur->m_cplex.getObjectiveValue() ; // pour le log
        _nbContraintesClique = selecteur->getNbContraintesDeClique() ; // pour le log

#ifdef VERBOSE_TEMPS
        std::cout << "time(selecteur)" << "\t"
                << time(NULL) - CInstance::params().getTimeStart()
                << "\t" << time(NULL) - maintenant
                << "\tnb.plans\t" << CPlan::getNbPlans()
                << "\tt.steps\t"      << T()
                << "\tcliques\t"    << selecteur->m_cliques.size()
                << std::endl;
        std::cout << "clock(selecteur)" << "\t"
                << (temps2 - temps1) / (double) CLOCKS_PER_SEC
                << std::endl ;
        maintenant = time(NULL);
#endif

        delete selecteur;
        selecteur = NULL ;
    }
    catch (CPlexFailed e) {
        if (selecteur)
            delete selecteur ;
        if (passes!=-1) return (euros)-1;
        else
            throw (string("Erreur dans l'appel à CPlex : ") + e);
    }
#else
    x = vector<size_t>(I());
    for (size_t p=0; p<CPlan::getNbPlans(); p++)
        x[ CPlan::getPlan(p)->m_i ] = p ;
#endif

    /**
     *          libération de mémoire pour le filtre aval
     */
//     for_each(graphes.begin(), graphes.end(), Delete());


    /**
     *          appel au filtre aval
     */
    {
        temps1 = clock() ;
        std::ostringstream oss;
        if (passes!=-1) oss << "_passe_" << printIntSurQuatreChiffres(passes) ;
        FiltreAval * filtreAval = NULL;
//         try {
        switch (CInstance::params().getFiltreAval()) {
            case fa_manuel :
                filtreAval = new FiltreAval();
                break;
            case fa_cplex :
                filtreAval = new FiltreAvalCplex();
                break;
            case fa_cplex_reste_libre :
                filtreAval = new FiltreAvalCplexResteLibre();
                break;
            case fa_cplex_ultime :
                filtreAval = new FiltreAvalCplexUltime();
                break;
            case fa_cplex_ultime_1passe :
                filtreAval = new FiltreAvalCplexUltime1Passe();
                break;
        }
        cout_solution = filtreAval->calcule_solution(x, multithread, prev_recharges);
        if (CInstance::params().getFiltreAval()==fa_cplex_ultime)
            _MinFiltreAvalPhase1Value = ((FiltreAvalCplexUltime *)filtreAval)->minPhase1Value() ;       // pour le log

        for (size_t p=0; p<x.size(); p++)
            plans.push_back( CPlan::getPlan(x[p]) );

        if (cout_solution < m_cout_solution)
        {
	  // enregistrement des colonnes de la solution
            m_cout_solution = cout_solution ;
            m_plans_solution.clear();
            for (vector<const CPlan *>::iterator ptp=plans.begin(); ptp!=plans.end(); ptp++)
            {
                m_plans_solution.insert(*ptp);
            }
        } else if (sur_disque) {
            x = vector<size_t>(I());
            for (set<const CPlan *>::iterator itp=m_plans_solution.begin(); itp!=m_plans_solution.end(); itp++)
                x[ (*itp)->m_i ] = CPlan::findIndex(*itp) ;
            cout_solution = filtreAval->calcule_solution(x);
        }
        if (sur_disque) {
            ofstream ofs((CInstance::params().getSolutionName()+oss.str()).c_str());
            filtreAval->printOutput(ofs, CInstance::params().getTimeStart());
            ofs.close();
        }

#ifdef EXPORTE_PLANS_PAR_SCENARIO
        if (passes==-1) // dernière passe
            filtreAval->exportePlans();
#endif
/*        } catch (...) {
        if (filtreAval!=NULL) delete filtreAval ;
        filtreAval = NULL ;
        throw;
    }*/
        if (filtreAval!=NULL) delete filtreAval;
        filtreAval = NULL ;

        temps2 = clock() ;
#ifdef VERBOSE_TEMPS
        std::cout << "time(f.aval) " << "\t"
                << time(NULL) - CInstance::params().getTimeStart()
                << "\t" << time(NULL) - maintenant
                << std::endl;
        maintenant = time(NULL);
        std::cout << "clock(f.aval)" << "\t"
                << (temps2 - temps1) / (double) CLOCKS_PER_SEC
                << std::endl ;
#endif
    }

    return cout_solution;
}



/**
 *
 *          Production d'une solution avec plusieurs essais sur la partie 1 du filtre aval ultime
 *
 */

euros COptimizer::produit_meilleure_solution(
                                             vector<CGraphe *> & graphes,
                                             time_t & maintenant,
                                             vector<const CPlan *> & plans,
                                             bool sur_disque )
{
    euros cout_solution(-1);

    /**
     *              APPEL A CPLEX
     */
    temps1 = clock() ;
    Selecteur * selecteur = NULL ;
    switch (CInstance::params().getSelecteur()) {
        case sel_1ct :
            selecteur = new Selecteur_1;
            break;
        case sel_ct :
            selecteur = new Selecteur_Thermiques;
            break;
        case sel_1cts :
            selecteur = new Selecteur_Scenario_1Thermique;
            break;
        case sel_cts :
            selecteur = new Selecteur_Scenario;
            break;
    }
    selecteur->setOption( multithread ? CPlexModel::Opportunistic : CPlexModel::None ) ;
    selecteur->setOption( CPlexModel::Time, prev_selecteur ) ;

    vector< vector<size_t> > x ;

    try {
        x = selecteur->populate();
        temps2 = clock() ;

	_MIPvalue = selecteur->m_cplex.getObjectiveValue() ; // pour le log
        _nbContraintesClique = selecteur->getNbContraintesDeClique() ; // pour le log

#ifdef VERBOSE_TEMPS
        std::cout << "time(selecteur)" << "\t"
                << time(NULL) - CInstance::params().getTimeStart()
                << "\t" << time(NULL) - maintenant
                << "\tnb.plans\t" << CPlan::getNbPlans()
                << "\tt.steps\t"      << T()
                << "\tcliques\t"    << selecteur->m_cliques.size()
                << std::endl;
        std::cout << "clock(selecteur)" << "\t"
                << (temps2 - temps1) / (double) CLOCKS_PER_SEC
                << std::endl ;
        maintenant = time(NULL);
#endif

        delete selecteur;
        selecteur = NULL ;
    }
    catch (CPlexFailed e) {
        if (selecteur)
            delete selecteur ;
        throw (string("Erreur dans l'appel à CPlex : ") + e);
    }

    /**
     *          libération de mémoire pour le filtre aval
     */
//     for_each(graphes.begin(), graphes.end(), Delete());


    /**
     *          appel au filtre aval
     */
    try {
            temps1 = clock() ;
            FiltreAvalCplexUltime * filtreAval = new FiltreAvalCplexUltime();

            // calcul de la durée autorisée pour la première phase du filtre aval
            time_t duree_phase_1 = 0 ;
            if (CInstance::params().getLimitType()==CParams::TIME_LIMIT)
                duree_phase_1 = CInstance::params().getDeadline() - time(NULL) - prev_scenarios * 2 ;
            else
                duree_phase_1 = 3600;

            duree_phase_1 = std::max( duree_phase_1, (time_t) 1 ) ;

            cout_solution = filtreAval->calcule_solution(x, duree_phase_1, multithread);
            _MinFiltreAvalPhase1Value = filtreAval->minPhase1Value() ;       // pour le log
            _FinalValue = (double) cout_solution.v ;   // pour le log

            if (sur_disque) {
                ofstream ofs(CInstance::params().getSolutionName().c_str());
                filtreAval->printOutput(ofs, CInstance::params().getTimeStart());
                ofs.close();
            }

            if (filtreAval!=NULL) delete filtreAval;
            filtreAval = NULL ;

            temps2 = clock() ;
    #ifdef VERBOSE_TEMPS
            std::cout << "time(f.aval) " << "\t"
                    << time(NULL) - CInstance::params().getTimeStart()
                    << "\t" << time(NULL) - maintenant
                    << std::endl;
            maintenant = time(NULL);
            std::cout << "clock(f.aval)" << "\t"
                    << (temps2 - temps1) / (double) CLOCKS_PER_SEC
                    << std::endl ;
    #endif


        return cout_solution;
    } catch (CPlexFailed e) {
        throw string( "Erreur dans le filtre aval, appel à CPlex : ") + e ;
    }
}



/**
 *
 *          CALCULS DES TEMPS et VERBOSITE
 *
 */

void COptimizer::calcule_temps_fin_passe()
{
    maintenant = time(NULL);
    if (m_passes>1)
        delta_duree_passe = std::max( maintenant - heure_debut_passe - duree_passe_precedente, delta_duree_passe );
    duree_passe_precedente = maintenant - heure_debut_passe ;
}



void verbose_rejet_plan_existant( CGraphe * g)
{
#ifdef VERBOSE_PLANS_REJETES
    cout << "rejet du plan " << g->m_plan->m_id
            << ", centrale " << g->m_i
            << " - plan existant" << endl;
#endif
}


void verbose_rejet_cout_positif( CGraphe * g)
{
#ifdef VERBOSE_PLANS_REJETES
    cout << "rejet du plan " << g->m_plan->m_id
            << ", centrale " << g->m_i
            << " - cout_reduit(" << g->m_plan->m_cout_reduit
            << ") > 0"
            << endl;
#endif
}


void COptimizer::verbose_temps_graphes()
{
    fin_creation_graphe = time(NULL);
#ifdef VERBOSE_TEMPS
     std::cout << "time(graphes) " << "\t"
             << fin_creation_graphe - CInstance::params().getTimeStart()
             << "\t" << fin_creation_graphe - debut_optimization
             << std::endl;
     maintenant = time(NULL);
     std::cout << "clock(graphes)" << "\t"
             << (temps2 - temps1) / (double) CLOCKS_PER_SEC
             << std::endl ;
#endif
}

void COptimizer::verbose_debut_passe()
{
    maintenant = time(NULL);
    heure_debut_passe = maintenant ;

    /* calcul des variables duales */
#ifdef VERBOSE_PASSES
    std::cout << std::endl << "passe n°" << m_passes + 1
        << (strategie_utilise_contraintes_de_cliques ? " - cliques" : "")
        << (strategie_dual_ascent ? " - duales" : "")
        << (strategie_accepte_plan_existant ? " - accepte plans existants" : "")
        << (strategie_rejete_couts_reduits_positifs ? " - rejette coûts réduits positifs" : " - accepte coûts réduits positifs")
        << "#############" << std::endl ;
#endif
#ifdef _RELEASE_
    pourcentage_de_temps_2 = ((heure_debut_passe - CInstance::params().getTimeStart()) * 100)
                            / CInstance::params().getTimeLimit() ;
    if (pourcentage_de_temps_2 != pourcentage_de_temps_1)
        std::cout << "elapsed time : " << pourcentage_de_temps_2 << "%" << std::endl;
    pourcentage_de_temps_1 = pourcentage_de_temps_2;
#endif
}



void COptimizer::verbose_avant_dual_ascent()
{
#ifdef VERBOSE_PASSES
        std::cout << "  calcul des variables duales" << std::endl ;
#endif
        temps1 = clock();
}



void COptimizer::verbose_apres_dual_ascent(vector<euros_par_kWh> new_u)
{
    temps2 = clock();
#ifdef VERBOSE_TEMPS
        std::cout << "time(DA) " << "\t"
                << time(NULL) - CInstance::params().getTimeStart()
                << "\t" << time(NULL) - maintenant
                << std::endl;
        maintenant = time(NULL);
        std::cout << "clock(DA)" << "\t"
                << (temps2 - temps1) / (double) CLOCKS_PER_SEC
                << std::endl ;
#endif
#ifdef VERBOSE_DUALS
        cout << "duals : " ;
        for (size_t t = 0; t<u.size(); t++)
            cout << new_u[t] << "\t" ;
        cout << endl;
        cout << "diffs : " ;
        for (size_t t = 0; t<u.size(); t++)
            cout << new_u[t]-u[t] << "\t" ;
        cout << endl;
#endif
}



void COptimizer::verbose_avant_viterbi()
{
#ifdef VERBOSE_PASSES
        std::cout << "  préparation viterbi" << std::endl ;
#endif
}



void COptimizer::verbose_avant_creer_colonnes()
{
#ifdef VERBOSE_PASSES
        std::cout << "  génération des colonnes" << std::endl ;
#endif
        temps1 = clock() ;
}



struct PredAffichePlanDeProduction {
    void operator()(CGraphe * g) {
        std::cout << *(g->m_plan) << std::endl;
    }
} ;
void COptimizer::verbose_apres_creer_colonnes()
{
    temps2 = clock() ;
#ifdef VERBOSE_TEMPS
    std::cout << "time(colonnes) " << "\t"
            << time(NULL) - CInstance::params().getTimeStart()
            << "\t" << time(NULL) - maintenant
            << std::endl;
    std::cout << "clock(colonnes)" << "\t"
            << (temps2 - temps1) / (double) CLOCKS_PER_SEC
            << std::endl ;
#endif

    maintenant = time(NULL);

#ifdef VERBOSE_PLANS
    PredAffichePlanDeProduction affichePlan;
    std::for_each( graphes.begin(), graphes.end(), affichePlan );
#endif

#ifdef VERBOSE_PASSES
    std::cout << "  ajout des nouvelles colonnes" << std::endl ;
#endif
}



/* écriture du fichier daniel.txt : résumé des infos à mettre dans le tableau */
void COptimizer::log_daniel_txt()
{
    std::ofstream daniel("daniel.txt") ;
    daniel
            << "data\tRelaxedMPLB\tIterToLB\tTotalIter\tnbPlans\tnbContraintes\tMPsol\tF.A.Phase1\tFinalSol\tMaxTime\tElapsedTime"
            << std::endl
            << CInstance::params().getInstanceName() << "\t"
            << _LBvalue << "\t"
            << _nbIterToLB << "\t"
            << _nbTotalIter << "\t"
            << _nbPlans << "\t"
            << _nbContraintesClique << "\t"
            << _MIPvalue << "\t"
            << _MinFiltreAvalPhase1Value << "\t"
            << _FinalValue << "\t"
            << CInstance::params().getTimeLimit() << "\t"
            << time(NULL) - CInstance::params().getTimeStart()
            << std::endl ;
    daniel.close();
}

/**
 *                  STRATEGIES
 */

bool COptimizer::prepare_strategie()
{
    static enum {cliques, dual_ascent} etat = cliques ;
    static bool premiere_fois = true ;
    static bool option_rejete_couts_reduits_positifs = false ;
    static bool option_accepte_plan_existant = true ;
    static char no_changes = 0 ;
    static int ancien_nb_colonnes_cout_non_positif = -1 ;

    if (premiere_fois) {
        option_rejete_couts_reduits_positifs = CInstance::params().getRejetCoutsReduitsPositifs() ;
        option_accepte_plan_existant = not CInstance::params().getRejetPlansExistants() ;
        premiere_fois = false  ;
        strategie_stop = false ;
    }
    switch (CInstance::params().getStrategie())
    {
        case toujoursDualAscent_sansCliques :

            strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
            strategie_accepte_plan_existant = option_accepte_plan_existant ;
            strategie_dual_ascent = true ;
            strategie_utilise_contraintes_de_cliques = false ;
            if (nb_rejets_couts_positifs+nb_rejets_existants==(unsigned)I())
                strategie_stop = true ;
            break;

        case jamaisDualAscent_avecCliques :

            strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
            strategie_accepte_plan_existant = option_accepte_plan_existant ;
            strategie_dual_ascent = false ;
            strategie_utilise_contraintes_de_cliques = true ;
            if (nb_rejets_couts_positifs+nb_rejets_existants==(unsigned)I())
                strategie_stop = true ;
            break;


        case Cliques_puis_DualAscent_AvecCliques_quand_plus_de_nouveaux_plans :

            switch (etat)
            {
                case cliques :
                    if (nb_rejets_couts_positifs+nb_rejets_existants<(unsigned)I())
                    {
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = false ;
                        strategie_dual_ascent = false ;
                        strategie_utilise_contraintes_de_cliques = true ;
                    }
                    else
                    {
                        etat = dual_ascent ;
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = option_accepte_plan_existant ;
                        strategie_dual_ascent = true ;
                        strategie_utilise_contraintes_de_cliques = true ;
                    }
                    break;

                case dual_ascent :
                    if (nb_rejets_couts_positifs+nb_rejets_existants==(unsigned)I())
                        strategie_stop = true ;
                    break;
            }
            break;


        case Cliques_puis_DualAscent_SansCliques_quand_plus_de_nouveaux_plans :

            switch (etat)
            {
                case cliques :
                    if (nb_rejets_couts_positifs+nb_rejets_existants<(unsigned)I())
                    {
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = false ;
                        strategie_dual_ascent = false ;
                        strategie_utilise_contraintes_de_cliques = true ;
                    }
                    else
                    {
                        etat = dual_ascent ;
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = option_accepte_plan_existant ;
                        strategie_dual_ascent = true ;
                        strategie_utilise_contraintes_de_cliques = false ;
                    }
                    break;

                case dual_ascent :
                    if (nb_rejets_couts_positifs+nb_rejets_existants==(unsigned)I())
                        strategie_stop = true ;
                    break;
            }
            break;


        case Cliques_et_DualAscent_alternes :

            switch (etat)
            {
                case cliques :
                    if (nb_colonnes_cout_non_positif > ancien_nb_colonnes_cout_non_positif)
                    {
                        etat    = dual_ascent ;
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = option_accepte_plan_existant ;
                        strategie_dual_ascent = true ;
                        strategie_utilise_contraintes_de_cliques = false ;
                        ancien_nb_colonnes_cout_non_positif = nb_colonnes_cout_non_positif ;
                    }
                    else // on reste dans l'état cliques sans recalculer les variables duales
                    {
                        etat    = cliques ;
                        strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                        strategie_accepte_plan_existant = false ;
                        strategie_dual_ascent = false ;
                        strategie_utilise_contraintes_de_cliques = true ;
                    }
                    break;

                case dual_ascent :
                    etat    = cliques ;
                    strategie_rejete_couts_reduits_positifs = option_rejete_couts_reduits_positifs ;
                    strategie_accepte_plan_existant = false ;
                    strategie_dual_ascent = (nb_colonnes_cout_non_positif > ancien_nb_colonnes_cout_non_positif) ;
                    strategie_utilise_contraintes_de_cliques = true ;
                    break;
            }
	    if (nb_rejets_couts_positifs+nb_rejets_existants==(unsigned)I())
	      no_changes ++ ;
	    else
	      no_changes = 0 ;
	    if (no_changes == 2)
	      strategie_stop = true ;
	    else
	      strategie_stop = false ;
            break;
    }
    nb_rejets_couts_positifs = 0 ;
    nb_rejets_existants = 0 ;
#ifdef VERBOSE_PASSES
    std::cout << "planification de passes : strategie_stop=" << (strategie_stop?"true":"false") << std::endl;
#endif
    return (not strategie_stop and continue_passes());
}