#include "filtreavalcplexultime.h"


#include <cinstance.h>
#include <cinstancereduite.h>
#include <cplan.h>
#include <cgraphe.h>
#include <fstream>
using std::ofstream;

// #define WRITE_CPLEX_FILE

FiltreAvalCplexUltime::FiltreAvalCplexUltime()
    : FiltreAval(),
                 m_cplex( "filtreAval", CPlexModel::Minimisation, CPlexModel::Linear,
#ifndef VERBOSE_CPLEX_INTERNAL
                          true,
#else
                          false,
#endif
#ifndef VERBOSE_SOLUTION_STATUS
                          true,
#else
                          false,
#endif
#ifdef WRITE_CPLEX_FILE
                          true,
                          true
#else
                          false,
                          false
#endif
                        )
{
// std::cout << "**** FiltreAvalCplexUltime::FiltreAvalCplexUltime()" << std::endl ;

    m_type2_prod_min.resize(I());
    m_type2_prod_max.resize(I());
}



// structure partagée par les différents threads
struct CIterThreadScenario {
    scenario_t      m_s ;           // dernier scenario pris en charge
    FiltreAvalCplexUltime * m_fa ;  // filtre aval contenant la solution à compléter
    euros           m_cout ;        // coût cumulé des solutions
    pthread_mutex_t lock_it ;       // mutex pour bloquer m_s et m_cout
    CIterThreadScenario( FiltreAvalCplexUltime * fa )
    {
        m_fa = fa ;
        pthread_mutex_init(&lock_it, NULL);
    }
};
// fonction "start" d'un thread
void * threaded_solution_scenario( void * iter_scenarios )
{
    CIterThreadScenario * its = (CIterThreadScenario *) iter_scenarios ; // itérateur sur les scenarios à traiter
    CPlexMutable cplex_model( "filtreAval", CPlexModel::Minimisation, CPlexModel::Linear,
#ifndef VERBOSE_CPLEX_INTERNAL
                              true,
#else
                              false,
#endif
#ifndef VERBOSE_SOLUTION_STATUS
                              true,
#else
                              false,
#endif
#ifdef WRITE_CPLEX_FILE
                              true,
                              true
#else
                                      false,
                                      false
#endif
                            );
    /* initialisation du modèle cplex avec les parties indépendantes du scénario */
    its->m_fa->setproblemdata( cplex_model ) ;

    /* parcours des scénarii et calcul des solutions */
    while (true) {
        scenario_t s = S() ;
        // recherche du prochain scenario à traiter
        pthread_mutex_lock( & its->lock_it );
        if ( its->m_s != S() ) {
            s = its->m_s ;
            ++ its->m_s ;
        }
        pthread_mutex_unlock( & its->lock_it );
        // fin du thread s'il n'y a plus de graphes à traiter
        if (s==S()) break ;

        // calcul de la solution du scenario à traiter
        euros _cout ;
        try {
            _cout = its->m_fa->calcule_solution( s, cplex_model );
        } catch (CPlexFailed e) {
            _cout = (euros) INFINI ;
        }
        // ajout du coût à la solution globale
        pthread_mutex_lock( & its->lock_it );
        its->m_cout += _cout ;
        pthread_mutex_unlock( & its->lock_it );
    }
    return NULL ;
}
// création de colonnes sur plusieurs threads
void FiltreAvalCplexUltime::calcule_solutions_scenarios_multithread()
{
//     std::cout << "**** FiltreAvalCplexUltime::calcule_solutions_scenarios_multithread()" << std::endl ;

    // 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
    CIterThreadScenario threaded_scenario(this) ;
    for (vector<pthread_t>::iterator th=threads.begin(); th!=threads.end(); th++)
    {
        pthread_create( &(*th), NULL, threaded_solution_scenario, (void *) &threaded_scenario );
    }
    // attente de la fin des threads lancés
    for (vector<pthread_t>::iterator th=threads.begin(); th!=threads.end(); th++)
    {
        pthread_join( *th, NULL );
    }
    if ((euros) INFINI <= threaded_scenario.m_cout)
        throw string("Au moins un scénario n'a pas trouvé de solution.");
}

euros FiltreAvalCplexUltime::calcule_solution( scenario_t s, CPlexMutable & cplex )
{
//     std::cout << "**** FiltreAvalCplexUltime::calcule_solution( scenario_t s, CPlexMutable & cplex )" << std::endl ;

    vector<double> deCplex;

    setproblemdata(s, cplex);
    cplex.resolution();


    reponse[s].costPourScenario = (euros) 0;

        // coût des recharges
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
        for (campaign_t k; k<K(); ++k)
            reponse[s].costPourScenario += C(i,k) * recharges[i][k] ;

        // remplissage des productions
        // thermiques
    for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j)
    {
        deCplex = cplex.getResult( string("production thermique centrale ")+j );
        for (timestep_t t; t<T(); ++t) {
            reponse[s].prod_thermique[j][t] = (kWh) deCplex[t];//std::max( EMIN(j,t,s), std::min( EMAX(j,t,s), (kWh) deCplex[t]));
            reponse[s].prod_thermique_totale[t] += reponse[s].prod_thermique[j][t];
            reponse[s].costPourScenario += C(j,t,s) * reponse[s].prod_thermique[j][t];
        }
    }
        // nucléaires
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        deCplex = cplex.getResult( string("production nucléaire centrale ")+i );
        for (timestep_t t; t<T(); ++t) {
            reponse[s].prod_nucleaire[i][t] = std::max( (kWh)0, std::min( m_type2_prod_max[i][t], (kWh) deCplex[t]));
            reponse[s].prod_nucleaire_totale[t] += reponse[s].prod_nucleaire[i][t];
        }
        reponse[s].stock[i] = m_plans[i]->remplis_stock(reponse[s].prod_nucleaire[i], recharges[i]);
    }
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        deCplex = cplex.getResult( string("fuel fin  , centrale ")+i );
        reponse[s].costPourScenario -= C(i,K_Tp1) * (kWh) deCplex[m_plans[i]->m_arcs.size()-1];
    }
    return reponse[s].costPourScenario ;
}

void FiltreAvalCplexUltime::setPlans( const vector<size_t> & x )
{
//     std::cout << "**** FiltreAvalCplexUltime::setPlans( const vector<size_t> & x )" << std::endl ;

    //Récupération des colonnes sélectionnées
    for( size_t i =0 ; i<x.size() ; i++) {
        m_plans[CPlan::getPlan(x[i])->m_i] = CPlan::getPlan(x[i]);
    }

    // vérification : on doit avoir un plan par centrale nucléaire
    for (size_t i=0; i<m_plans.size(); i++) {
        if (m_plans[i]==NULL)
            throw string("Erreur : une centrale n'a pas de plan de production");
    }

    // remplissage des min et max par pas de temps
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
        m_plans[i]->remplis_min_max(m_type2_prod_min[i],m_type2_prod_max[i]);
}





euros FiltreAvalCplexUltime::calcule_solution_phase_1(
        std::vector<size_t> x,
        std::vector<std::vector<kWh> > & _recharges,
        double temps_max,
        bool multithread )
{
//     std::cout << "**** FiltreAvalCplexUltime::calcule_solution_phase_1" << std::endl ;


    euros _cout = (euros) INFINI ;

    setPlans( x ) ;

    // création du modèle et appel à cplex pour le calcul des recharges
    // base de temps : la semaine

#ifdef VERBOSE_TEMPS
    std::cout << "FiltreAvalCplexUltime::calcule_solution : définition du problème des recharges" << std::endl ;
#endif

    setproblemdata_recharges();

//     std::cout << "FiltreAvalCplexUltime::calcule_solution : écriture sur disque pour DEBUG" << std::endl ;
//     m_cplex.write_pb();
//     m_cplex.write_names();
#ifdef VERBOSE_TEMPS
    std::cout << "FiltreAvalCplexUltime::calcule_solution : résolution du problème des recharges" << std::endl ;
#endif

    if (temps_max > 0)
        m_cplex.setOption( CPlexModel::Time, temps_max );

    _cout = (euros) m_cplex.resolution();

#ifdef VERBOSE_TEMPS
    std::cout << "FiltreAvalCplexUltime::calcule_solution : problème des recharges résolu" << std::endl ;
    if (m_cplex.getInfeasible()) std::cout << "\t\t\tINFAISABLE" <<std::endl ;
#endif

    vector<double> deCplex;

    // récupération des recharges
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        deCplex = m_cplex.getResult( string("recharges ")+i );
        for (campaign_t k; k<K(); ++k) {
            _recharges[i][k] = (kWh) deCplex[k] ;
        }
    }

    return (euros) _cout ;
}




euros FiltreAvalCplexUltime::calcule_solution_phase_2(bool multithread)
{
//     std::cout << "**** FiltreAvalCplexUltime::calcule_solution_phase_2" << std::endl ;

    valeurObjectif = (euros) INFINI ;

    if ( multithread ) {
        this->calcule_solutions_scenarios_multithread() ;
    }
    else {
        // création du modèle, partie indépendante des scénarios
        // base de temps : unité du problème initial
        setproblemdata(m_cplex);

        // Boucle sur les scénarios et récupération de la solution
        for (scenario_t s = (scenario_t) 0 ; s < S() ; ++s)
        {
#ifdef VERBOSE_TEMPS
            clock_t temps3 = clock() ;
#endif
            calcule_solution( s, m_cplex );
#ifdef VERBOSE_TEMPS
            clock_t temps4 = clock() ;
            std::cout << "clock(scenario_" << s << ")\t"
                    << (temps4 - temps3) / (double) CLOCKS_PER_SEC
                    << std::endl ;
#endif
        }
    }

    // calcul du coût moyen de la solution
    valeurObjectif = (euros)0 ;
    for (scenario_t s = (scenario_t) 0 ; s < S() ; ++s)
        valeurObjectif += reponse[s].costPourScenario ;
    valeurObjectif /= S();

    return valeurObjectif ;
}

/**
 *
 * Calcul d'une solution. La première phase du filtre est lancée plusieurs fois
 * s'il y a le temps et s'il y a plusieurs jeux de colonnes à essayer
 *
 */

euros FiltreAvalCplexUltime::calcule_solution(std::vector<std::vector<size_t> > x, double max_time_passes_1, bool multithread)
{

//     std::cout << "**** FiltreAvalCplexUltime::calcule_solution(std::vector<std::vector<size_t> > x, double max_time_passes_1, bool multithread)" << std::endl ;



    time_t debut = time(NULL) ;
    time_t maintenant = debut ;
    time_t temps_recharges = maintenant ;
    time_t max_duree_recharges  = 0 ;
    time_t duree_ecoulee        = 0 ;

    size_t sol_id = 0 , meilleure_sol_id = 0 ;
    bool _infeasible = true ;
    euros _cout_min = (euros) INFINI ;
    euros _cout ;
    vector< vector<kWh> > _recharges( I(), vector<kWh>(K()) ) ;
    while ( (_infeasible or (duree_ecoulee + max_duree_recharges < max_time_passes_1)) and sol_id<x.size() )
    {
        _cout = calcule_solution_phase_1( x[sol_id], _recharges, (_infeasible ? -1 : max_time_passes_1 - duree_ecoulee), multithread );

        bool this_infeasible = m_cplex.getInfeasible() ;
        /* libération mémoire */
        m_cplex.reinit();

        if ( (_cout < _cout_min and (this_infeasible == _infeasible)) or (_infeasible and not this_infeasible) ) {
            _cout_min = _cout ;
            recharges = _recharges ;
            meilleure_sol_id = sol_id ;
            _infeasible = _infeasible and this_infeasible ;
        }

        maintenant = time(NULL) ;

#ifdef VERBOSE_TEMPS
        std::cout << "time(solution_" << sol_id << "/" << x.size() - 1 << ")\t"
                << maintenant - temps_recharges
                << std::endl ;
        std::cout << "cout(solution_" << sol_id << "/" << x.size() - 1 << ")\t"
                << _cout
                << std::endl ;
#endif

        max_duree_recharges = std::max( max_duree_recharges, maintenant-temps_recharges );
        temps_recharges = maintenant ;
        duree_ecoulee = maintenant - debut ;
#ifdef VERBOSE_TEMPS
        std::cout << "duree_ecoulee " << duree_ecoulee << std::endl ;
        std::cout << "duree_max_recharge " << max_duree_recharges << std::endl ;
        std::cout << "max_time_passes_1 : " << max_time_passes_1 << std::endl ;
#endif
        sol_id ++ ;
    }

    _minPh1v = (double) _cout_min.v ; // pour le log

    /* calcul des solutions par scénario */
#ifdef VERBOSE_TEMPS
    clock_t temps1 = clock() ;
#endif

//     std::cout << "    solution " << meilleure_sol_id << " retenue pour phase 2" << std::endl;

    setPlans( x[meilleure_sol_id] );
    calcule_solution_phase_2(multithread) ;

#ifdef VERBOSE_TEMPS
    clock_t temps2 = clock() ;

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


    return valeurObjectif;
}

/**
 *
 * calcule d'une solution en utilisant un seul jeu de colonnes.
 * La phase 1 du filtre est appelée une seule fois.
 *
 */

euros FiltreAvalCplexUltime::calcule_solution(std::vector<size_t> x, double prev_recharges, bool multithread)
{
//     std::cout << "**** FiltreAvalCplexUltime::calcule_solution(std::vector<size_t> x, double prev_recharges, bool multithread)" << std::endl ;

    clock_t temps1 = clock() ;
    time_t maintenant = time(NULL) ;
    time_t temps_recharges = maintenant, temps_scenario = maintenant ;

    euros _cout1 = calcule_solution_phase_1( x, recharges, prev_recharges, multithread ) ;
    /* libération mémoire */
    m_cplex.reinit();

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

    /* calcul des solutions par scénario */
    temps1 = clock() ;

    calcule_solution_phase_2(multithread) ;

    temps2 = clock() ;

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


    return valeurObjectif;
}


FiltreAvalCplexUltime::~FiltreAvalCplexUltime()
{
//     std::cout << "**** FiltreAvalCplexUltime::~FiltreAvalCplexUltime()" << std::endl ;

}


void FiltreAvalCplexUltime::setproblemdata_recharges()
{
//     std::cout << "**** FiltreAvalCplexUltime::setproblemdata_recharges()" << std::endl ;

    // matrices de production min et max par centrale et par semaine
    vector<vector<kWh> > _type2_prod_min_h;
    vector<vector<kWh> > _type2_prod_max_h;
    vector<kWh> production_nucleaire_maximale;
    vector<vector<kWh> > demande_si_type2_max;
    vector<vector<kWh> > modulation_a_garantir;
    vector<int> nb_moduls ;

    _type2_prod_min_h.resize( I(), vector<kWh>(H()) );
    _type2_prod_max_h.resize( I(), vector<kWh>(H()) );

    // calcul préalable des bornes de production par pas de temps puis par semaine, pour chaque centrale
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i) {
        // déjà fait dans setPlans
//         m_plans[i]->remplis_min_max(m_type2_prod_min[i],m_type2_prod_max[i]);
        for (week_t h=(week_t) 0; h<H(); ++h) {
            kWh prodmin = (kWh) 0, prodmax = (kWh) 0;
            for (timestep_t t=week2timestep(h); t<week2timestep(h+(wlen_t)1); ++t) {
                prodmin += m_type2_prod_min[i][t];
                prodmax += m_type2_prod_max[i][t];
            }
            _type2_prod_min_h[i][h] = prodmin ;
            _type2_prod_max_h[i][h] = prodmax ;
        }
    }

    // vecteur de la production nucléaire maximale
    production_nucleaire_maximale.resize( T() );
    for ( timestep_t t; t<T(); ++t) {
        for (powerplants2_t i; i<I(); ++i) {
            production_nucleaire_maximale[t] += m_type2_prod_max[i][t] ;
        }
    }

    // matrice de la demande "diminuée" par la production nucléaire max et
    // par la production thermique min,
    // pour chaque scenario
    demande_si_type2_max.resize( S(), vector<kWh>(T()) );
    for (scenario_t s; s<S(); ++s) {
        for (timestep_t t; t<T(); ++t) {
            demande_si_type2_max[s][t] = EDEM(t,s) - production_nucleaire_maximale[t] ;
            for (powerplants1_t j; j<J(); ++j) {
                demande_si_type2_max[s][t] -= EMIN(j,t,s);
            }
        }
    }

    // matrice de la modulation totale à garantir, pour chaque scenario et par
    // semaine. Il faut de la modulation (sur les nucléaires) lorsqu'on a un pas
    // de temps durant lequel la production nucléaire max
    // dépasse (demande - prod_thermique_min)
    modulation_a_garantir.resize( S(), vector<kWh>(T()) );
    nb_moduls.resize( S(), 0 ) ;

    int nb_moduls_total = 0;
    for (scenario_t s; s<S(); ++s) {
        for (timestep_t t; t<T(); ++t) {
            if (demande_si_type2_max[s][t]<(kWh)0) {
                modulation_a_garantir[s][t] = demande_si_type2_max[s][t] ;
                nb_moduls[s] ++ ;
            }
        }
#ifdef VERBOSE_AJOUT_PAS_DE_TEMPS
        std::cout << "nb modulations timestep, scenario " << s << " = " << nb_moduls[s] << std::endl;
#endif
        nb_moduls_total += nb_moduls[s] ;
    }
#ifdef VERBOSE_AJOUT_PAS_DE_TEMPS
    std::cout << "nb modulations timestep, TOTAL = " << nb_moduls_total << std::endl;
#endif
    /** définition des zones horizontales */

    // zones horizontales = variables

    // 1 - production par semaine pour chaque scénario

//     std::cout << "   1 - production par semaine pour chaque scénario" << std::endl;

    for (scenario_t s=(scenario_t)0; s<S(); ++s)
    {
        for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j) {
            m_cplex.addZoneHorizontale( string("ptc")+j+"s"+s, H()+nb_moduls[s]);
            m_cplex.setZoneHorizontale( string("ptc")+j+"s"+s );
            for ( week_t h=(week_t) 0; h<H(); ++h) {
                kWh _emin ;
                kWh _emax ;
                // séparation de la production sur les timesteps à modulation imposée
                for (timestep_t t=week2timestep(h); t<week2timestep(h+wlen_t(1)); ++t) {
                    if (modulation_a_garantir[s][t]>=(kWh)0) {
                        _emin += EMIN(j,t,s);
                        _emax += EMAX(j,t,s);
                    }
                }
                m_cplex.addToObj(   h,
                                    _emin.v,
                                        _emax.v,
                                            CPlexModel::Continuous,
                                            C(j,s,h).v / S() );
            }
            // timesteps à modulation imposée
            int tt = H() ;
            for (timestep_t t; t<T(); ++t) {
                if (modulation_a_garantir[s][t]<(kWh)0) {
                    m_cplex.addToObj(   tt,
                                        EMIN(j,t,s).v,
                                        EMAX(j,t,s).v,
                                        CPlexModel::Continuous,
                                        C(j,t,s).v / S() );
                    tt++;
                }
            }
        }
        for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i) {
            m_cplex.addZoneHorizontale( string("pnc")+i+"s"+s, H()+nb_moduls[s]);
            m_cplex.setZoneHorizontale( string("pnc")+i+"s"+s);
            for ( week_t h=(week_t) 0; h<H(); ++h) {
                kWh _emin ;
                kWh _emax ;
                // séparation de la production sur les timesteps à modulation imposée
                for (timestep_t t=week2timestep(h); t<week2timestep(h+wlen_t(1)); ++t) {
                    if (modulation_a_garantir[s][t]>=(kWh)0) {
                        _emin += m_type2_prod_min[i][t];
                        _emax += m_type2_prod_max[i][t];
                    }
                }
                m_cplex.addToObj(   h,
                                    _emin.v,
                                    _emax.v,
                                    CPlexModel::Continuous,
                                    0 );
            }
            // timesteps à modulation imposée
            int tt = H() ;
            for (timestep_t t; t<T(); ++t) {
                if (modulation_a_garantir[s][t]<(kWh)0) {
                    m_cplex.addToObj(   tt,
                                        m_type2_prod_min[i][t].v,
                                        m_type2_prod_max[i][t].v,
                                        CPlexModel::Continuous,
                                        0 );
                    tt++;
                }
            }
        }
    }

    // 2 - recharges de chaque campagne

//     std::cout << "   2 - recharges de chaque campagne" << std::endl;

    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        m_cplex.addZoneHorizontale( string("recharges ")+i, K());
        m_cplex.setZoneHorizontale( string("recharges ")+i );
        // 0 par défaut
        for (campaign_t k; k<K(); ++k)
        {
            m_cplex.addToObj( k, 0, 0, CPlexModel::Continuous, 0 );
        }
        // entre RMIN et RMAX si la recharge est programmée
        for (size_t a=0; a<m_plans[i]->m_arcs.size(); a++)
        {
            const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
            campaign_t k = arc.m_depart->m_campaign ;
            if (arc.m_depart->m_type == TYPE_OUTAGE) {
                m_cplex.addToObj(   k, RMIN(i,k).v, RMAX(i,k).v, CPlexModel::Continuous, C(i,k).v );
            }
        }
    }

    // 3 - fuel au début et à la fin de chaque arc, et après recharge

//     std::cout << "   3 - fuel au début et à la fin de chaque arc, et après recharge" << std::endl;

    for (scenario_t s=(scenario_t)0; s<S(); ++s)
    {
//         std::cout << "       s = " << s << std::endl ;
        for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
        {
//             std::cout << "       i = " << i << std::endl ;
//             std::cout << "       m_plans[i]->m_arcs.size() = " << m_plans[i]->m_arcs.size() << std::endl ;

            m_cplex.addZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s, m_plans[i]->m_arcs.size());
            m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
            for (size_t a=0;a<m_plans[i]->m_arcs.size();a++) m_cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);

            m_cplex.addZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s, m_plans[i]->m_arcs.size());
            m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
            for (size_t a=0;a<m_plans[i]->m_arcs.size();a++) m_cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);

            m_cplex.addZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s, m_plans[i]->m_arcs.size());
            m_cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s);
            for (size_t a=0;a<m_plans[i]->m_arcs.size()-1;a++) m_cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);

//             std::cout << "       C(i,K_Tp1).v = " << C(i,K_Tp1).v << std::endl ;

            m_cplex.addToObj(   m_plans[i]->m_arcs.size()-1, 0, INFINI, CPlexModel::Continuous, - C(i,K_Tp1).v / S() );
        }
    }




    /** contraintes */

    /// A - production = demande

//     std::cout << "   A - production = demande" << std::endl;

    for (scenario_t s=(scenario_t)0; s<S(); ++s)
    {
        m_cplex.addZoneVerticale( string("dem_s_")+s, H()+nb_moduls[s] ) ;
        m_cplex.setZoneVerticale( string("dem_s_")+s);
        // partie droite : demande
        for ( week_t h=(week_t) 0; h<H(); ++h) {
            kWh dem = EDEM(s,h) ;
            // séparation de la production sur les timesteps à modulation imposée
            for (timestep_t t=week2timestep(h); t<week2timestep(h+wlen_t(1)); ++t) {
                if (modulation_a_garantir[s][t]<(kWh)0) {
                    dem -= EDEM(t,s);
                }
            }
            m_cplex.addToRHS(h, CPlexModel::Equals, dem.v);
        }
        // timesteps à modulation imposée
        int tt = H() ;
        for (timestep_t t; t<T(); ++t) {
            if (modulation_a_garantir[s][t]<(kWh)0) {
                m_cplex.addToRHS(   tt,CPlexModel::Equals, EDEM(t,s).v );
                tt++;
            }
        }
        // partie gauche : production thermique et nucléaire
        for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j) {
            m_cplex.setZoneHorizontale( string("ptc")+j+"s"+s);
            for ( week_t h=(week_t) 0; h<H()+nb_moduls[s]; ++h)
                m_cplex.addToZ(h,h,1);
        }
        for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i) {
            m_cplex.setZoneHorizontale( string("pnc")+i+"s"+s);
            for ( week_t h=(week_t) 0; h<H()+nb_moduls[s]; ++h)
                m_cplex.addToZ(h,h,1);
        }
    }

    /// B - productions et stocks pour chaque arc


//     std::cout << "   B - productions et stocks pour chaque arc" << std::endl;

    // fuel avant recharge (a+1) = fuel après recharge(a) + somme productions(a)
    for (scenario_t s=(scenario_t)0; s<S(); ++s)
    {
        for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
        {
            // fuel début des arcs - fuel fin des arcs précédents = 0
            size_t nb_arcs = m_plans[i]->m_arcs.size() ;
            m_cplex.addZoneVerticale( string("fd=ff ")+s+" "+i, nb_arcs-1 );
            m_cplex.setZoneVerticale( string("fd=ff ")+s+" "+i);
            m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
            for (size_t a=1; a<nb_arcs; a++)    m_cplex.addToZ(a-1, a, 1);
            m_cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s);
            for (size_t a=0; a<nb_arcs-1; a++)  m_cplex.addToZ(a,   a,-1);
            for (size_t a=0; a<nb_arcs-1; a++)  m_cplex.addToRHS(a, CPlexModel::Equals, 0);

            // fuel après refuel - B*fuel avant refuel  - recharge = C
            m_cplex.addZoneVerticale( string("far_is_f(f)")+s+" "+i+" ", nb_arcs );
            m_cplex.setZoneVerticale( string("far_is_f(f)")+s+" "+i+" ");
            for (size_t a=0; a<nb_arcs; a++)
            {
                const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
                campaign_t k = arc.m_depart->m_campaign;
                switch ( arc.m_depart->m_type ) {

                    case TYPE_OUTAGE :
                        m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
                        m_cplex.addToZ(a, a, ((coeff_t)1-(coeff_t)1/Q(i,k)).v);
                        m_cplex.setZoneHorizontale( string("recharges ")+i);
                        m_cplex.addToZ(a, k, 1);
                        m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
                        m_cplex.addToZ(a, a, -1);
                        m_cplex.addToRHS(a, CPlexModel::Equals, (BO(i,(campaign_t)(k-1)) * ((Q(i,k)-(coeff_t)1)/Q(i,k)) - BO(i,k)).v);
                        break;

                    default :
                        m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
                        m_cplex.addToZ(a, a, 1);
                        m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
                        m_cplex.addToZ(a, a,-1);
                        m_cplex.addToRHS(a, CPlexModel::Equals, 0);
                        break;
                }
            }

            for (size_t a=0; a<nb_arcs; a++)
            {
                const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
                campaign_t k = arc.m_depart->m_campaign;
                if ( arc.m_depart->m_type == TYPE_OUTAGE ) {
                    // fuel après refuel < SMAX
                    m_cplex.addZoneVerticale( string("far<smax")+s+" "+i+" "+a, 1 );
                    m_cplex.setZoneVerticale( string("far<smax")+s+" "+i+" "+a);

                    m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a, 1);
                    m_cplex.addToRHS(0, CPlexModel::LowerThan, SMAX(i,k).v);

                    // fuel avant refuel < AMAX
                    m_cplex.addZoneVerticale( string("far<amax")+s+" "+i+" "+a, 1 );
                    m_cplex.setZoneVerticale( string("far<amax")+s+" "+i+" "+a);

                    m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a, 1);
                    m_cplex.addToRHS(0, CPlexModel::LowerThan, AMAX(i,k).v);
                }
                if ( arc.m_depart->m_type != TYPE_BO and arc.m_fin->m_date!=0 ) {
                    // fuel fin > BO
                    m_cplex.addZoneVerticale( string("ff>BO")+s+" "+i+" "+a, 1 );
                    m_cplex.setZoneVerticale( string("ff>BO")+s+" "+i+" "+a);

                    m_cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a, 1);
                    m_cplex.addToRHS(0, CPlexModel::GreaterThan, sup(BO(i,k)).v);
                }
                if ( arc.m_depart->m_type == TYPE_BO or arc.m_depart->m_date == (week_t)0 ) {
                    // fuel debut fixé
                    m_cplex.addZoneVerticale( string("fuel début fixe, scenario ")+s+" centrale "+i+" arc "+a, 1 );
                    m_cplex.setZoneVerticale( string("fuel début fixe, scenario ")+s+" centrale "+i+" arc "+a);

                    m_cplex.setZoneHorizontale( string("fuel début, centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a, 1);
                    m_cplex.addToRHS(0, CPlexModel::Equals, arc.m_depart->m_fuel.v);
                }
            }
            // fuel après - fuel fin - production = 0
            m_cplex.addZoneVerticale( string("fa_-_ff_-_p_=_0_")+s+"_"+i, nb_arcs );
            for (size_t a=0; a<nb_arcs; a++)
            {
                const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
                campaign_t k = arc.m_depart->m_campaign;

                // fuel après - fuel fin - production = 0
                m_cplex.setZoneVerticale( string("fa_-_ff_-_p_=_0_")+s+"_"+i );
                m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
                m_cplex.addToZ(a, a, 1);
                m_cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s);
                m_cplex.addToZ(a, a,-1);
                m_cplex.setZoneHorizontale( string("pnc")+i+"s"+s);
                kWh prod_max = (kWh) 0;
                for ( week_t h=arc.m_depart->m_date_campaign; h<arc.m_fin->m_date; ++h) {
                    m_cplex.addToZ(a,h,-1);
                    prod_max += EMAX(i,h);
                }
                // timesteps à modulation imposée
                int tt = H() ;
                for (timestep_t t; t<week2timestep(arc.m_fin->m_date); ++t) {
                    if (modulation_a_garantir[s][t]<(kWh)0) {
                        if (t>=week2timestep(arc.m_depart->m_date_campaign))
                            m_cplex.addToZ(a,tt,-1);
                        tt++;
                    }
                }
                m_cplex.addToRHS(a, CPlexModel::Equals, 0);

                // modulation
                if (arc.m_depart->m_type != TYPE_BO) {
                    m_cplex.addZoneVerticale( string("fa_m_ff_greater_pmax_m_mmax")+s+" "+i+" "+a, 1 );
                    m_cplex.setZoneVerticale( string("fa_m_ff_greater_pmax_m_mmax")+s+" "+i+" "+a );
                    m_cplex.setZoneHorizontale( string("fuel après, centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a, 1);
                    m_cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i+" scenario "+s);
                    m_cplex.addToZ(0, a,-1);
                    m_cplex.addToRHS(0, CPlexModel::GreaterThan, (prod_max-MMAX(i,k)).v);
                }
            }
        }
    }
//     std::cout << "--- FiltreAvalCplexUltime::setproblemdata_recharges terminé" << std::endl;
}



void FiltreAvalCplexUltime::setproblemdata(scenario_t s, CPlexMutable & cplex)
{
//     std::cout << "**** FiltreAvalCplexUltime::setproblemdata(scenario_t s, CPlexMutable & cplex)" << std::endl ;

    for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j) {
        cplex.setZoneHorizontale( string("production thermique centrale ")+j );
        for ( timestep_t t=(timestep_t) 0; t<T(); ++t)
            cplex.addToObj(   t,
                                EMIN(j,t,s).v,
                                     EMAX(j,t,s).v,
                                          CPlexModel::Continuous,
                                          C(j,t,s).v );
    }

    cplex.setZoneVerticale( string("production=demande"));
    // partie droite : demande
    for ( timestep_t t=(timestep_t) 0; t<T(); ++t) {
        cplex.addToRHS(t, CPlexModel::Equals, EDEM(t,s).v);
    }

}


void FiltreAvalCplexUltime::setproblemdata(CPlexMutable & cplex)
{
//     std::cout << "**** FiltreAvalCplexUltime::setproblemdata(CPlexMutable & cplex)" << std::endl ;

    // calcul préalable des bornes de production par pas de temps, pour chaque centrale
    // déjà fait dans setproblemdata_recharges
//     for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
//         m_plans[i]->remplis_min_max(m_type2_prod_min[i],m_type2_prod_max[i]);

    /** définition des zones horizontales */

    // zones horizontales = variables

    // 1 - production par pas de temps pour chaque scénario


    for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j) {
        cplex.addZoneHorizontale( string("production thermique centrale ")+j, T());
    }
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i) {
        cplex.addZoneHorizontale( string("production nucléaire centrale ")+i, T());
        cplex.setZoneHorizontale( string("production nucléaire centrale ")+i);
        for ( timestep_t t=(timestep_t) 0; t<T(); ++t)
            cplex.addToObj( t,
                                m_type2_prod_min[i][t].v,
                                m_type2_prod_max[i][t].v,
                                CPlexModel::Continuous,
                                0 );
    }
    // 2 - recharges de chaque campagne
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        cplex.addZoneHorizontale( string("recharges ")+i, K());
        cplex.setZoneHorizontale( string("recharges ")+i );
        // 0 par défaut
        for (campaign_t k; k<K(); ++k)
        {
            cplex.addToObj( k, recharges[i][k].v, recharges[i][k].v, CPlexModel::Continuous, C(i,k).v );
        }
    }
    // 3 - fuel au début et à la fin de chaque arc, et après recharge
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        cplex.addZoneHorizontale( string("fuel début, centrale ")+i, m_plans[i]->m_arcs.size());
        cplex.setZoneHorizontale( string("fuel début, centrale ")+i);
        for (size_t a=0;a<m_plans[i]->m_arcs.size();a++) cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);

        cplex.addZoneHorizontale( string("fuel après, centrale ")+i, m_plans[i]->m_arcs.size());
        cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
        for (size_t a=0;a<m_plans[i]->m_arcs.size();a++) cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);

        cplex.addZoneHorizontale( string("fuel fin  , centrale ")+i, m_plans[i]->m_arcs.size());
        cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i);
        for (size_t a=0;a<m_plans[i]->m_arcs.size()-1;a++) cplex.addToObj(a,0,INFINI,CPlexModel::Continuous,0);
        cplex.addToObj(   m_plans[i]->m_arcs.size()-1, 0, INFINI, CPlexModel::Continuous, - C(i,K_Tp1).v );
    }




    /** contraintes */

    /// A - production = demande
    cplex.addZoneVerticale( string("production=demande"), T() ) ;
    cplex.setZoneVerticale( string("production=demande"));

    // partie gauche : production thermique et nucléaire
    for (powerplants1_t j=(powerplants1_t)0; j<J(); ++j) {
        cplex.setZoneHorizontale( string("production thermique centrale ")+j);
        for ( timestep_t t=(timestep_t) 0; t<T(); ++t)
            cplex.addToZ(t,t,1);
    }
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i) {
        cplex.setZoneHorizontale( string("production nucléaire centrale ")+i);
        for ( timestep_t t=(timestep_t) 0; t<T(); ++t)
            cplex.addToZ(t,t,1);
    }

    /// B - productions et stocks pour chaque arc

    // fuel avant recharge (a+1) = fuel après recharge(a) + somme productions(a)
    for (powerplants2_t i=(powerplants2_t)0; i<I(); ++i)
    {
        // fuel début des arcs - fuel fin des arcs précédents = 0
        size_t nb_arcs = m_plans[i]->m_arcs.size() ;
        cplex.addZoneVerticale( string("fd=ff ")+i, nb_arcs-1 );
        cplex.setZoneVerticale( string("fd=ff ")+i);
        cplex.setZoneHorizontale( string("fuel début, centrale ")+i);
        for (size_t a=1; a<nb_arcs; a++)    cplex.addToZ(a-1, a, 1);
        cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i);
        for (size_t a=0; a<nb_arcs-1; a++)  cplex.addToZ(a,   a,-1);
        for (size_t a=0; a<nb_arcs-1; a++)  cplex.addToRHS(a, CPlexModel::Equals, 0);

        // fuel après refuel - B*fuel avant refuel  - recharge = C
        cplex.addZoneVerticale( string("far_is_f(f)")+i+" ", nb_arcs );
        cplex.setZoneVerticale( string("far_is_f(f)")+i+" ");
        for (size_t a=0; a<nb_arcs; a++)
        {
            const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
            campaign_t k = arc.m_depart->m_campaign;
            switch ( arc.m_depart->m_type ) {

                case TYPE_OUTAGE :
                    cplex.setZoneHorizontale( string("fuel début, centrale ")+i );
                    cplex.addToZ(a, a, ((coeff_t)1-(coeff_t)1/Q(i,k)).v);
                    cplex.setZoneHorizontale( string("recharges ")+i);
                    cplex.addToZ(a, k, 1);
                    cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
                    cplex.addToZ(a, a, -1);
                    cplex.addToRHS(a, CPlexModel::Equals, (BO(i,(campaign_t)(k-1)) * ((Q(i,k)-(coeff_t)1)/Q(i,k)) - BO(i,k)).v);
                    break;

                default :
                    cplex.setZoneHorizontale( string("fuel début, centrale ")+i);
                    cplex.addToZ(a, a, 1);
                    cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
                    cplex.addToZ(a, a,-1);
                    cplex.addToRHS(a, CPlexModel::Equals, 0);
                    break;
            }
        }

        for (size_t a=0; a<nb_arcs; a++)
        {
            const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
            campaign_t k = arc.m_depart->m_campaign;
            if ( arc.m_depart->m_type == TYPE_OUTAGE ) {
                // fuel après refuel < SMAX
                cplex.addZoneVerticale( string("far<smax")+" "+i+" "+a, 1 );
                cplex.setZoneVerticale( string("far<smax")+" "+i+" "+a);

                cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
                cplex.addToZ(0, a, 1);
                cplex.addToRHS(0, CPlexModel::LowerThan, SMAX(i,k).v);

                // fuel avant refuel < AMAX
                cplex.addZoneVerticale( string("far<amax")+" "+i+" "+a, 1 );
                cplex.setZoneVerticale( string("far<amax")+" "+i+" "+a);

                cplex.setZoneHorizontale( string("fuel début, centrale ")+i);
                cplex.addToZ(0, a, 1);
                cplex.addToRHS(0, CPlexModel::LowerThan, AMAX(i,k).v);
            }
            if ( arc.m_depart->m_type != TYPE_BO and arc.m_fin->m_date!=0 ) {
                // fuel fin > BO
                cplex.addZoneVerticale( string("ff>BO")+" "+i+" "+a, 1 );
                cplex.setZoneVerticale( string("ff>BO")+" "+i+" "+a);

                cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i);
                cplex.addToZ(0, a, 1);
                cplex.addToRHS(0, CPlexModel::GreaterThan, sup(BO(i,k)).v);
            }
            if ( arc.m_depart->m_type == TYPE_BO or arc.m_depart->m_date == (week_t)0 ) {
                // fuel debut fixé
                cplex.addZoneVerticale( string("fuel début fixe,")+" centrale "+i+" arc "+a, 1 );
                cplex.setZoneVerticale( string("fuel début fixe,")+" centrale "+i+" arc "+a);

                cplex.setZoneHorizontale( string("fuel début, centrale ")+i);
                cplex.addToZ(0, a, 1);
                cplex.addToRHS(0, CPlexModel::Equals, arc.m_depart->m_fuel.v);
            }
        }
        // fuel après - fuel fin - production = 0
        cplex.addZoneVerticale( string("fa_-_ff_-_p_=_0_")+"_"+i, nb_arcs );
        for (size_t a=0; a<nb_arcs; a++)
        {
            const CPlan::CArc & arc = m_plans[i]->m_arcs[a];
            campaign_t k = arc.m_depart->m_campaign;

            // fuel après - fuel fin - production = 0
            cplex.setZoneVerticale( string("fa_-_ff_-_p_=_0_")+"_"+i );
            cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
            cplex.addToZ(a, a, 1);
            cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i);
            cplex.addToZ(a, a,-1);
            cplex.setZoneHorizontale( string("production nucléaire centrale ")+i);
            kWh prod_max = (kWh) 0;
            for ( timestep_t t=week2timestep(arc.m_depart->m_date_campaign); t<week2timestep(arc.m_fin->m_date); ++t) {
                cplex.addToZ(a,t,-1);
                prod_max += EMAX(i,t);
            }
            cplex.addToRHS(a, CPlexModel::Equals, 0);

            // modulation
            if (arc.m_depart->m_type != TYPE_BO) {
                cplex.addZoneVerticale( string("fa_m_ff_greater_pmax_m_mmax")+" "+i+" "+a, 1 );
                cplex.setZoneVerticale( string("fa_m_ff_greater_pmax_m_mmax")+" "+i+" "+a );
                cplex.setZoneHorizontale( string("fuel après, centrale ")+i);
                cplex.addToZ(0, a, 1);
                cplex.setZoneHorizontale( string("fuel fin  , centrale ")+i);
                cplex.addToZ(0, a,-1);
                cplex.addToRHS(0, CPlexModel::GreaterThan, (prod_max-MMAX(i,k)).v);
            }
        }
    }
}