#include "cgraphe.h"

#include "../instance/cinstance.h"
#include "../instance/cinstancereduite.h"
#include "viterbi.h"
#include "cplan.h"

#include <algorithm>
#include <functional>

using std::pair ;




/**
 * perte de fuel lors d'un refuel
 * @param i index de la centrale
 * @param kmoins1 campagne en cours avant l'outage
 * @param k campagne de l'outage
 * @param x niveau de fuel avant l'outage
 * @return
 */
inline double perte_refuel( powerplants2_t i, campaign_t k, double x )
{
    //return 0 ;
    return ((Q(i,k)-1)/Q(i,k)) * BO(i,(campaign_t)(k-1)) - BO(i,k) + x/Q(i,k) ;
}

/**
 * niveaux possibles après refuel
 * @param i centrale
 * @param k campagne de l'outage
 * @param scourant niveau avant le refuel, auquel on a déjà ôté la perte
 * @param smin niveau minimum après refuel
 * @param smax niveau maximum après refuel
 */
inline void refuel(powerplants2_t i, campaign_t k, double scourant, double & smin, double & smax)
{
    smin = std::min(scourant + RMIN(i, k), SMAX(i,k));
    smax = std::min(scourant + RMAX(i, k), SMAX(i,k));
}

static double __defaut__ ;
/**
 * calcul de la production pendant une semaine
 * @param i         index de la centrale
 * @param k         campagne de production
 * @param h         semaine
 * @param type      type du noeud de départ (TYPE_PREMIER, TYPE_OUTAGE ou TYPE_BO)
 * @param smin      stock minimal au début de la semaine
 * @param smax      stock maximal au début de la semaine
 * @param modul     modulation pouvant encore être utilisée sur la campagne
 * @param emin      production minimale sur la campagne (début de semaine : entrée, fin de semaine : sortie)
 * @param emax      production maximale sur la campagne (début de semaine : entrée, fin de semaine : sortie)
 * @param smin_d    stock minimal au début de la campagne
 * @param smax_d    stock maximal au début de la campagne
 * @param smin_n    stock minimal à la fin de la semaine                                            (sortie)
 * @param smax_n    stock maximal à la fin de la semaine                                            (sortie)
 * @param modul_n   modulation pouvant encore être utilisée sur la campagne à la fin de la semaine  (sortie)
 * @return 
 */
inline bool consommation_semaine( powerplants2_t i, campaign_t k, week_t h, noeud_t type,
                                  double smin, double smax,
                                  double modul,
                                  double & emin,    double &    emax,
                                  double   smin_d,  double      smax_d,
                                  double & smin_n,  double &    smax_n,
                                  double & modul_n = __defaut__ )
{
    if (type == TYPE_BO)
    {
        emax += std::min(smin, (1+epsilon()) * PB(i, k, smin, h, (wlen_t) 1));
        emin += std::min(smax, (1-epsilon()) * PB(i, k, smax, h, (wlen_t) 1));
        smax_n = std::max((double) 0, smax_d - emin );
        smin_n = std::max((double) 0, smin_d - emax );
    } else {
        // test sur le seuil BO
        if ( smax < BO(i, k) )
            return true ;
        emax += EMAX(i,h) ;
        emin += EMAX(i, h) - std::min( EMAX(i,h), modul) ;
        modul_n = modul - std::min( EMAX(i,h), modul );
        smax_n = smax_d - emin ;
        smin_n = smin_d - emax ;
    }
    return false ;
}

inline void consommation_semaine_BO( powerplants2_t i, campaign_t k, week_t h,
                                  double & smin, double & smax)
{
    smin = std::max((double) 0, smin - (1+epsilon()) * PB(i, k, smin, h, (wlen_t) 1) );
    smax = std::max((double) 0, smax - (1-epsilon()) * PB(i, k, smax, h, (wlen_t) 1) );
}

/* attributs de classe */
std::vector<double>  CGraphe::m_somme_ui_EMAX ;
const std::vector<double> * CGraphe::m_u ;
std::vector<Ui> CGraphe::U ;

/* constructeur */
CGraphe::CGraphe(powerplants2_t i)
{
    m_i = i ;
    /* création du vecteur m_somme_EMAX */
    prepareCreation() ;

    /* création des noeuds */
    creerNoeuds();

    /* création des arcs */
    creerArcs();

    /* initialisation du plan de production */
    m_plan = NULL ;
}

CGraphe::~CGraphe()
{
    if (m_plan) delete m_plan ;
    for_each(m_noeuds.begin(), m_noeuds.end(), Delete());
}

/**
 * Prépare la création du graphe
 */
void CGraphe::prepareCreation()
{
    m_somme_EMAX.resize(H()+1,0);
    double somme = 0 ;
    for (week_t h=(week_t)0; h<H(); ++h) {
        somme+=EMAX(m_i,h);
        m_somme_EMAX[h+1] = somme;
    }
}

/**
 * Création des noeuds du graphe
 */

void CGraphe::creerNoeuds()
{
    //
    // Création de la liste des semaines d'outage
    //

    m_campaigns.resize(H());
    for (campaign_t k=(campaign_t)0; k<K(); ++k)
    {
        for (week_t h=TO(m_i,k); h<=TA(m_i,k); ++h)
        {
            m_campaigns[h].push_back(k) ;
        }
    }

    //
    // CREATION du premier noeud
    //

    creer_noeud( (week_t)0,
                  (XI(m_i)-EMAX(m_i,(week_t)0)>=BO(m_i, (campaign_t)-1)) ? TYPE_PREMIER : TYPE_BO,
                   (campaign_t)-1,
                    XI(m_i)
               );
    // éventuels outages sur le premier noeud
    noeuds_eventuels_outage( *m_noeuds.back(), (week_t) 0, XI(m_i), XI(m_i) );

    // Parcours des faisceaux d'arcs et création de nouveaux noeuds
    for ( week_t h=(week_t) 0; h<H(); ++h)
        // pour chaque semaine possédant des noeuds
        if ( m_map_noeuds.find(h)!=m_map_noeuds.end())
            for ( campaign_levels_t::iterator ik=m_map_noeuds[h].begin(); ik!=m_map_noeuds[h].end(); ik++)
                for (type_levels_t::iterator it=ik->second.begin(); it!=ik->second.end(); it++)
                    for( niveaux_t::iterator in=it->second.begin(); in!=it->second.end(); in++)
                        // pour chaque noeud
                        if (in->second!=NULL)
#ifdef VERBOSE_NB_NOEUDS
                        {
                            std::cout << "faisceau (centrale " << m_i << ", sem." << h << ", campagne " << ik->first << ", type " << TYPE(it->first)
                                    << ", niveau " << in->first << ") : " << m_noeuds.size() << " noeuds avant" << std::endl;
                            std::cout.flush();
#endif
                            faisceau(*in->second) ;
#ifdef VERBOSE_NB_NOEUDS
                            std::cout << "faisceau (centrale " << m_i << ", sem." << h << ", campagne " << ik->first << ", type " << TYPE(it->first)
                                    << ", niveau " << in->first << ") : " << m_noeuds.size() << " noeuds après." << std::endl;
                        }
#endif

    //
    // Création du dernier noeud
    //
    creer_noeud( H(), TYPE_DERNIER, (campaign_t)-1, 0 );
}


/**
 * création des noeuds auxquels un noeud de départ peut se brancher
 * 
 * @param depart 
 */
void CGraphe::faisceau( const CNoeud & depart )
{
    // début de la production
    week_t debut_production = depart.m_date_campaign ;
    // niveaux de fuel possibles au départ de la production :
    double _smax_d ;
    double _smin_d ;
    if (depart.m_type == TYPE_OUTAGE) {
        refuel( m_i, depart.m_campaign, depart.m_fuel_avant_refuel, _smin_d, _smax_d) ;
    } else { // TYPE_BO et TYPE_PREMIER
        _smin_d = _smax_d = depart.m_fuel ;
    }

    if (depart.m_type == TYPE_BO) {
        faisceau_suite_BO( depart, debut_production, _smin_d, _smax_d ) ;
    } else if (depart.m_type == TYPE_PREMIER) {
        faisceau_suite_modul( depart, (week_t) 1, _smin_d, _smax_d ) ;
    } else { // depart.m_type == TYPE_OUTAGE
        faisceau_suite_modul( depart, debut_production, _smin_d, _smax_d ) ;
    }
}


void CGraphe::faisceau_suite_BO( const CNoeud & depart, week_t h, double smin, double smax )
{
    // production de la semaine
    consommation_semaine_BO( m_i, depart.m_campaign, h, smin, smax ) ;
    // création éventuelle de noeuds
    ++h ;
    // test sur la date
    if (h>=H()) return ;
    noeuds_eventuels_outage( depart, h, smin, smax ) ;
    // suite
    faisceau_suite_BO( depart, h, smin, smax ) ;
}

#define _EMAX(h0,h1) (m_somme_EMAX[h1]-m_somme_EMAX[h0])
void CGraphe::faisceau_suite_modul( const CNoeud & depart, week_t h, double smin_d, double smax_d )
{
    if (h>H()) return ;
    // stock restant min et max :
    double smin = smin_d - _EMAX(depart.m_date_campaign, h) ;
    double smax = smax_d - std::max( (double) 0,  _EMAX(depart.m_date_campaign, h) - MMAX(m_i, depart.m_campaign) ) ;
#ifdef DEBUG_CREATION_NOEUDS
    std::cout << "h=" << h << "\tsmin=" << smin << "\tsmax=" << smax << "\tEMAX=" << _EMAX(h+(wlen_t)(-1), h) << std::endl ;
#endif
    // création de noeuds de type BO
    if (smin<BO(m_i, depart.m_campaign))
    {
        double smin_hmoins1_vers_BO = std::max(
                                               smin_d - _EMAX(depart.m_date_campaign,h-1),
                                               BO(m_i,depart.m_campaign)
                                      );
        double smax_hmoins1_vers_BO = std::min(
                                               smax_d - std::max( (double) 0,  _EMAX(depart.m_date_campaign,h-1) - MMAX(m_i, depart.m_campaign) ),
                                               BO(m_i,depart.m_campaign) + EMAX(m_i,h+(wlen_t)(-1)));

        noeuds_eventuels_BO( depart, (week_t)(h-1), smin_hmoins1_vers_BO, smax_hmoins1_vers_BO);
    }
    // test sur la date et la quantité de stock restante
    if (h>=H() or smax<BO(m_i,depart.m_campaign)) return ;

    // création de noeuds de type Outage
    noeuds_eventuels_outage( depart, h, std::max(smin, BO(m_i,depart.m_campaign)), smax ) ;

    // suite
    faisceau_suite_modul( depart, ++h, smin_d, smax_d ) ;
}



void CGraphe::noeuds_eventuels_outage( const CNoeud & depart, week_t h, double smin, double smax )
{
    // création de noeuds de type OUTAGE
    for (vector<campaign_t>::iterator kit=m_campaigns[h].begin(); kit!=m_campaigns[h].end(); kit++)
    {
        if ( *kit > depart.m_campaign ) // l'indice de campagne k doit croître strictement d'un outage à l'autre
        {
            creer_noeuds_intervalle( h, TYPE_OUTAGE, *kit,
                                     smin, std::min(smax, AMAX(m_i,*kit)),
                                             smin==0 or smin==BO(m_i,depart.m_campaign),
                                             smax==0 or smax==BO(m_i,depart.m_campaign) or smax>=AMAX(m_i,*kit) );
        }
    }
}

/**
 * Création éventuelle de noeuds une semaine donnée.
 * La création dépend des bornes du fuel de la semaine courante et précédente (pour les noeuds BO)
 * et des dates d'outage autorisées (pour les noeuds OUTAGE)
 * @param depart noeud de départ
 * @param h semaine courante
 * @param smin_p min fuel de la semaine précédente
 * @param smax_p max fuel de la semaine précédente
 * @param smin min fuel de la semaine courante
 * @param smax max fuel de la semaine courante
 */
void CGraphe::noeuds_eventuels_BO( const CNoeud & depart, week_t h, double smin, double smax )
{
    if (smax<smin) return ;
    // création de noeuds de type BO
    creer_noeuds_intervalle( h, TYPE_BO, depart.m_campaign, smin, smax,
                             smin==0 or smin==BO(m_i,depart.m_campaign),
                                     smax==0 or smax==BO(m_i,depart.m_campaign));
}

struct DansIntervalle {
    double m_min, m_max;
    DansIntervalle(double min, double max) {m_min=min; m_max=max;}
    bool operator()(const double & x) {return (x >= m_min) and (x <= m_max);}
};
struct PairDansIntervalle {
    double m_min, m_max;
    PairDansIntervalle(double min, double max) {m_min=min; m_max=max;}
    bool operator()(const pair<double, CNoeud *> & x) {return (x.first >= m_min) and (x.first <= m_max);}
};
struct HorsIntervalle {
    double m_min, m_max;
    HorsIntervalle(double min, double max) {m_min=min; m_max=max;}
    bool operator()(const double & x) {return (x > m_min) or (x > m_max);}
};
template <class T1, class T2>
T1 First( const pair<T1,T2> & p ) { return p.first; }

/**
 * Création de noeuds une semaine donnée à l'intérieur d'un intervalle donné de fuel
 * @param h semaine de la création
 * @param type type des noeuds à créer
 * @param k campagne de production des noeuds à créer
 * @param min fuel min
 * @param max fuel max
 * @param force_min forcer la création d'un noeud au niveau de fuel min
 * @param force_max forcer la création d'un noeud au niveau de fuel max
 */
void CGraphe::creer_noeuds_intervalle( week_t h, noeud_t type, campaign_t k, double min, double max, bool force_min, bool force_max)
{
    if (min>max) return ;

    // création des noeuds imposés
    if (force_min) creer_noeud( h, type, k, min );
    if (force_max) creer_noeud( h, type, k, max );

    // création d'un noeud entre min et max si aucun des deux n'est imposé et qu'il n'en existe aucun
    if (not force_min and not force_max and
        count_if(m_map_noeuds[h][k][type].begin(), m_map_noeuds[h][k][type].end(), PairDansIntervalle(min,max))==0)
        creer_noeud( h, type, k, (min+max)/2 );

    // liste des niveaux existants compris entre min et max
    vector<double> niveaux_existants ;
    for (std::map<double,CNoeud *>::iterator it=m_map_noeuds[h][k][type].begin(); it!=m_map_noeuds[h][k][type].end(); it++) {
        if (DansIntervalle(min,max)(it->first))
            niveaux_existants.push_back(it->first) ;
    }

    // création de noeuds intermédiaires entre min et le premier existant
    double premier = niveaux_existants.front() ;
    for (premier -= pavage(h); premier>=min; premier -= pavage(h)) {
        creer_noeud( h, type, k, premier ) ;
    }
    // création de noeuds intermédiaires entre le dernier existant et max
    double dernier = niveaux_existants.back() ;
    for (dernier += pavage(h); dernier<=max; dernier += pavage(h)) {
        creer_noeud( h, type, k, dernier ) ;
    }
    // création de noeuds intermédiaires entre les noeuds existants
    if (niveaux_existants.size()>=2) {
        premier = niveaux_existants.front() ;
        for ( vector<double>::iterator it = niveaux_existants.begin()+1; it!=niveaux_existants.end(); it++ )
        {
            double second = *it ;
            int nb_noeuds = (int)((second-premier) / pavage(h)) ;
            double pas = (second-premier) / (nb_noeuds+1) ;
            for ( premier+=pas; premier<second ; premier += pas ) {
                creer_noeud( h, type, k, premier ) ;
            }
            premier = *it ;
        }
    }
}

/**
 * pas de pavage : distance maximale entre deux niveaux de stocks
 * @param h semaine
 * @return 
 */
double CGraphe::pavage(week_t h)
{
    return EMAX(m_i,h)==0 ? 1 : EMAX(m_i,h) * FACTEUR_DE_PAVAGE ;
}

/**
 * Création d'un noeud
 * @param h semaine
 * @param type type de noeud
 * @param k campagne
 * @param fuel niveau de fuel
 */
void CGraphe::creer_noeud( week_t h, noeud_t type, campaign_t k, double fuel )
{
    // sortie si le noeud existe déjà
    if (m_map_noeuds.find(h) != m_map_noeuds.end() and
        m_map_noeuds[h].find(k) != m_map_noeuds[h].end() and
        m_map_noeuds[h][k].find(type) != m_map_noeuds[h][k].end() and
        m_map_noeuds[h][k][type].find(fuel) != m_map_noeuds[h][k][type].end())
        return;

    // création d'un nouveau noeud
    CNoeud * noeud_ptr = new CNoeud;
    CNoeud & noeud(*noeud_ptr);
    noeud.m_type        = type ;
    noeud.m_campaign    = k ;
    noeud.m_fuel        = fuel ;
    noeud.m_fuel_avant_refuel = (type==TYPE_OUTAGE ? fuel - perte_refuel(m_i,k,fuel) : fuel) ;
    noeud.m_date        = h ;
    switch (type) {
        case TYPE_BO :
            noeud.m_date_outage     = (week_t) -1 ;
            noeud.m_date_campaign   = h ;
            break;
        case TYPE_OUTAGE :
            noeud.m_date_outage     = h ;
            noeud.m_date_campaign   = h + DA(m_i,k) ;
            break;
        case TYPE_PREMIER :
            if (h!=0) throw string("Mauvaise date pour le premier noeud");
            noeud.m_date_outage     = (week_t) -1 ;
            noeud.m_date_campaign   = h ;
            break ;
        case TYPE_DERNIER :
            if (h!=H()) throw string("Mauvaise date pour le dernier noeud");
            noeud.m_date_outage     = (week_t) -1 ;
            noeud.m_date_campaign   = (week_t) h ;
            break ;
    }
    m_noeuds.push_back( noeud_ptr ) ;
    m_map_noeuds[h][k][type][fuel] = noeud_ptr ;
#ifdef VERBOSE_CREATION_NOEUD
    std::cout  << "\tsem. " << noeud_ptr->m_date << "\tk=" << noeud_ptr->m_campaign << "\t" << TYPE(noeud_ptr->m_type)
            << "\ts="  << noeud_ptr->m_fuel << "\ts'=" << noeud_ptr->m_fuel_avant_refuel << std::endl ;
#endif
}




/**
 *
 * creerArcs : création des arcs
 *
 */

void CGraphe::creerArcs()
{
#ifdef VERBOSE_CREATION_ARCS
    std::cout << "Tri des noeuds. Nb noeuds = " << m_noeuds.size() << std::endl ;
#endif
    // tri des noeuds
    std::sort(m_noeuds.begin(), m_noeuds.end(), OrdreCNoeud());

    // initialisation de la liste des arcs
    m_arcs.resize(0,CArc(NULL,-1,NULL,-1,0,0));

    // parcours des noeuds de départ

#ifdef VERBOSE_CREATION_ARCS
//     unsigned int completion  = -1 ;
    size_t nb_arcs_precedent = -1 ;
#endif
    for (size_t dep=0; dep!=m_noeuds.size()-1; dep++)
    {
#ifdef VERBOSE_CREATION_ARCS
/*        if ( completion != (dep * 100) / m_noeuds.size() ) {
            completion = (dep * 100) / m_noeuds.size() ;*/
            std::cout << "Parcours des " << m_noeuds.size() << " noeuds :\t" << dep << " (" << m_arcs.size() << " arcs : + " << m_arcs.size() - nb_arcs_precedent << ")" << std::endl ;
//         }
            if (m_arcs.size()==nb_arcs_precedent) throw string("Aucun arc créé partant du noeud précédent");
            nb_arcs_precedent = m_arcs.size() ;
#endif
        CNoeud & depart = *m_noeuds[dep];

        // niveaux de fuel possibles au départ de la production
        double _smax_d ;
        double _smin_d ;
        if (depart.m_type == TYPE_OUTAGE) {
            refuel( m_i, depart.m_campaign, depart.m_fuel_avant_refuel, _smin_d, _smax_d) ;
        } else { // TYPE_BO et TYPE_PREMIER
            _smin_d = _smax_d = depart.m_fuel ;
        }

        if (depart.m_type == TYPE_BO)
        {
            double _smax, _smin ;
            _smax = _smin = _smin_d ;
            // semaine du noeud d'arrivée courant
            week_t _semaine_atteinte = depart.m_date_campaign ;
            for (size_t arr=dep+1; arr<m_noeuds.size()-1; arr++)
            {
                CNoeud & arrivee = *m_noeuds[arr];

                // vérification des types et des numéros de campagne
                if ( (arrivee.m_type!=TYPE_OUTAGE) or (depart.m_campaign >= arrivee.m_campaign) )
                        continue ;

                // productions possibles sur la période :
                for (week_t h = _semaine_atteinte ; h < arrivee.m_date; ++ h ) {
                    consommation_semaine_BO( m_i, depart.m_campaign, h, _smin, _smax );
                }
                _semaine_atteinte = arrivee.m_date ;

                // ajout de l'arc, s'il est réalisable
                if (_smax >= arrivee.m_fuel and _smin <= arrivee.m_fuel)
                {
                    m_arcs.push_back( CArc(m_noeuds[dep],dep, m_noeuds[arr], arr, 0, 0) ) ;
                }
            }
            // dernier_noeud
            m_arcs.push_back( CArc(m_noeuds[dep], dep, m_noeuds.back(), m_noeuds.size()-1, 0, 0) ) ;
        }
        else // types TYPE_PREMIER et TYPE_OUTAGE
        {
            // parcours des noeuds d'arrivée
            for (size_t arr=dep+1; arr<m_noeuds.size()-1; arr++)
            {
                CNoeud & arrivee = *m_noeuds[arr];

                // vérification de la date d'arrivée
                if ( (depart.m_type==TYPE_OUTAGE and depart.m_date_campaign>arrivee.m_date) )
                    continue ;

                // vérification des numéros de campagne
                if ( (arrivee.m_type==TYPE_BO and depart.m_campaign != arrivee.m_campaign)
                    or (arrivee.m_type==TYPE_OUTAGE and depart.m_campaign >= arrivee.m_campaign) )
                            continue ;

                // productions possibles sur la période :
                double _smin = _smin_d - _EMAX(depart.m_date_campaign, arrivee.m_date) ;
                double _smax = _smax_d - std::max( (double) 0,  _EMAX(depart.m_date_campaign, arrivee.m_date) - MMAX(m_i, depart.m_campaign) ) ;

                // arrêt du parcours des noeuds d'arrivée si on est passé sous BO
                if (_smax < BO(m_i,depart.m_campaign))
                    break ;

                // ajout de l'arc, s'il est réalisable
                if (_smax >= arrivee.m_fuel and _smin <= arrivee.m_fuel)
                {
                    m_arcs.push_back( CArc(m_noeuds[dep], dep, m_noeuds[arr], arr, 0, 0) ) ;
                }
            }
            // dernier noeud
            double _smax = _smax_d
                        - std::max( (double) 0,
                                     _EMAX(std::min(H(),depart.m_date_campaign), H())
                                             - MMAX(m_i, depart.m_campaign) ) ;

            if (_smax>=BO(m_i,depart.m_campaign))
            {
                m_arcs.push_back( CArc(m_noeuds[dep], dep, m_noeuds.back(), m_noeuds.size()-1, 0, 0) ) ;
            }
        }
    }

    // tri des arcs
    std::sort(m_arcs.begin(), m_arcs.end(), OrdreArc());

    // création des vecteurs de taille nb_arcs
    m_cout          .resize(m_arcs.size(), 0 );
    m_noeud_debut   .resize(m_arcs.size(), -1);
    m_noeud_fin     .resize(m_arcs.size(), -1);
    std::transform(m_arcs.begin(), m_arcs.end(), m_noeud_debut  .begin(), CArc::de);
    std::transform(m_arcs.begin(), m_arcs.end(), m_noeud_fin    .begin(), CArc::vers);

    // création du vecteur m_premier_arc
    m_premier_arc.resize( m_noeuds.size(), -1 );
    for (size_t index_noeud=-1, index_arc=0; index_arc<m_arcs.size(); index_arc++)
    {
        if (m_noeud_debut[index_arc]!=index_noeud) {
            do {
                m_premier_arc[++index_noeud] = index_arc ;
            } while (index_noeud<m_noeud_debut[index_arc]);
        }
    }
    m_premier_arc[m_noeuds.size()-1] = m_arcs.size() ; // arc fictif, adapté au parcours dans Viterbi

#ifdef VERBOSE_NB_ARCS
    printf("nombre d'arcs pour la centrale %d : %u\n", m_i, (unsigned int)m_noeud_debut.size());
#endif
#ifdef __DUMP_GRAPHE__
    dump() ;
#endif
}



void CGraphe::dump()
{
    std::cout << "********* dumping graph " << m_i << " ***"<<std::endl ;

    for (size_t i=0; i<m_noeuds.size(); i++)
    {
        CNoeud & n = *(m_noeuds[i]);
        std::cout << "\tnoeud " << i << " (" << TYPE(n.m_type) << "," << n.m_date << "," << n.m_fuel << "," << n.m_fuel_avant_refuel << ")\n" ;
    }
    
    for (size_t n=0; n<m_arcs.size(); n++)
    {
        CArc & a = m_arcs[n];
        std::cout << "\tarc "<< n << " : de " << a.ide << " ("<<TYPE(a.m_de->m_type)<<","<<a.m_de->m_date<<")"
                << " vers " << a.ivers << " ("<<TYPE(a.m_vers->m_type)<<","<<a.m_vers->m_date<<")" << std::endl;
    }
    std::cout << "********* fin dump *********" << std::endl << std::endl ;
}



void CGraphe::creer_colonnes( CGraphe * graphe )
{
    graphe->creer_colonne() ;
}

void CGraphe::creer_colonne()
{
           /* calculer les coûts de chaque arc */
            calcule_couts_des_arcs();

            /* viterbi pour extraire le meilleur chemin */
            viterbise();

            /* création d'une nouvelle colonne en fonction du meilleur chemin */
            // plan de production
            plan_de_production_apres_viterbi();
}








/**
 * calcule_couts_des_arcs
 */

// création d'un vecteur ordonné des u_i

struct CreerUi {
    CreerUi( const std::vector<double> * _u ) : u(_u) {
        n = 0 ;
    }
    Ui operator()() {
        return Ui(u->at(n), (week_t) n++);
    }
    const vector<double> * u ;
    size_t n ;
};

bool ordreUi(const Ui & a, const Ui & b) { return a.u < b.u ; }

void CGraphe::prepare_viterbi(const vector<double> * u)
{
    // conservation d'une référence à ui,vi
    m_u = u ;

    // création du vecteur ordonné des Ui
    U.resize( H() );
    std::generate(U.begin(), U.end(), CreerUi(u));
    std::sort(U.begin(), U.end(), ordreUi);

    // création de somme_ui_EMAX
    m_somme_ui_EMAX.resize( H() + 1) ;
    m_somme_ui_EMAX[0] = 0 ;
    double somme = 0 ;
    for (week_t h=(week_t)1; h<=H(); ++h)
    {
        somme += u->at(h-1) ;
        m_somme_ui_EMAX[h] = somme ;
    }
}

void CGraphe::calcule_couts_des_arcs( )
{

    for (size_t arc=0; arc<m_noeud_fin.size(); arc++)
    {
        m_cout[arc] = calcule_cout_arc<false>( arc ) ;
    }
}

/**
 * appel à Viterbi
 */

void CGraphe::viterbise()
{
    m_meilleur_chemin = Viterbi().minimise(m_premier_arc, m_noeud_fin, m_cout,
                                &m_cout_minimal);
}


/**
 * plan_de_production_apres_viterbi
 */

void CGraphe::plan_de_production_apres_viterbi()
{
    // création d'un plan de production vide
    if (m_plan) delete m_plan ; // m_plan vaut NULL si le plan a été accepté quelque part
    m_plan = new CPlan(m_i, H(),I());

    // appel à la fonction calcule_cout_arc pour les campagnes de production
    m_cout_minimal = 0 ;
    for (std::list<size_t>::const_iterator it_arc = m_meilleur_chemin.begin();
        it_arc!=m_meilleur_chemin.end(); it_arc++)
    {
        calcule_cout_arc<true>(*it_arc);
    }
}


/**
 *
 * calcul du coût d'un arc
 * effet de bord : remplit le plan de production m_plan du graphe
 *
 */

struct EstDans{
    EstDans(size_t i, size_t j) {m_i=i; m_j=j;}
    bool operator()(const Ui & ui) { return (size_t)ui.i>=m_i and (size_t)ui.i<m_j;}
    size_t m_i; size_t m_j;
} ;


template<bool enregistrer_plan>
double CGraphe::calcule_cout_arc( size_t arc_idx )
{
#ifdef DEBUG_COUT_ARCS
    std::cout << "centrale n°" << m_i << "\tcout arc n° " << arc_idx << std::endl ;
#endif
    // noeud de départ
    const CNoeud & depart = *m_arcs[arc_idx].m_de ;
    const CNoeud & arrivee = *m_arcs[arc_idx].m_vers ;
    const double Cur = depart.m_campaign == -1 ? 0 : C(m_i, depart.m_campaign) ;
    double _cout = 0 ;
    double _r = 0; // recharge
    double _fuel = depart.m_fuel_avant_refuel ;

    /**************************************/
    /*   cas facile du type BO            */
    /**************************************/

    if ( depart.m_type == TYPE_BO )
    {
            // facile : la production est imposée
            for (week_t h=depart.m_date_campaign; h<arrivee.m_date ; ++h)
            {
                if (!enregistrer_plan and _fuel<=0) break ;
                if (enregistrer_plan and _fuel<=0)
                    m_plan->set_prod(h,0);
                else {
                    double prod = PB( m_i, depart.m_campaign, _fuel, h, (wlen_t) 1 );
                    _fuel -= prod ;
                    _cout -= (m_u->at(h)) * prod ;
                    if (enregistrer_plan) {
                        m_plan->set_prod(h,prod) ;
                    }
                }
            }
            _cout += Cur * _r ;
            if (enregistrer_plan) {
                m_plan->push_back(depart, arrivee, _r, _cout, Cur*_r, arrivee.m_fuel);
            }
            return _cout ;
    }

    /*****************************************************/
    /*    cas des campagnes de productions non BO        */
    /*****************************************************/

    // Prédicat pour parcourir les Ui triés appartenant à la période de production courante */
    EstDans estDans(depart.m_date_campaign, arrivee.m_date) ;

    // stock min et max
    double _smin_d, _smax_d ;
    if (depart.m_type == TYPE_OUTAGE) {
        refuel( m_i, depart.m_campaign, depart.m_fuel_avant_refuel, _smin_d, _smax_d) ;
    } else { // TYPE_PREMIER
        _smin_d = _smax_d = depart.m_fuel ;
    }

    if (enregistrer_plan) {
        if (depart.m_type == TYPE_OUTAGE) {
            for (week_t h=depart.m_date; h<depart.m_date_campaign and h<H(); ++h)
                m_plan->set_prod(h,0) ;
        }
    }
    /*****************************************************/
    /*    cas des campagnes de productions non finales   */
    /*****************************************************/
    // calcul des productions et du coût
    if ( arrivee.m_type != TYPE_DERNIER )
    {
        // somme des u_i*EMAX sur la durée
        _cout = m_somme_ui_EMAX[arrivee.m_date] - m_somme_ui_EMAX[depart.m_date_campaign] ;

        // prod min et max sur la durée
        double _emax = _EMAX(depart.m_date_campaign, arrivee.m_date) ;
        //double _emin = std::max( (double) 0,  _emax - MMAX(m_i, depart.m_campaign) ) ;
        // prod min et max pour l'arc
        double _emax_arc = _smax_d - arrivee.m_fuel ; // forcément positif car l'arc a été créé
        double _emin_arc = std::max( (double) 0, _smin_d - arrivee.m_fuel ) ;
        // modulation min et max
        double _mmax = std::min(MMAX(m_i,depart.m_campaign), _emax - _emin_arc) ;
        double _mmin = std::max((double)0, _emax - _emax_arc) ;

        // modulation accomplie
        double _modul = 0 ;
        vector<Ui>::iterator itui = find_if( U.begin(), U.end(), estDans);
        // modulation tant qu'on y est forcé
        double modul = 0;
        while (_modul<_mmin) {
            modul = std::min( EMAX(m_i, itui->i), _mmin-_modul );
            _modul += modul ;
            _cout += (itui->u) * modul ;
            if (enregistrer_plan) {
                m_plan->set_prod(itui->i, EMAX(m_i, itui->i)-modul) ;
            }
            if (_modul<_mmin)
                itui = find_if( ++itui, U.end(), estDans ) ;
        }
        // annulation de la modulation courante si elle est réimposée dans
        // la boucle suivante
        if (Cur-itui->u > 0) {
            _modul -= modul;
            _cout -= (itui->u) * modul ;
        }

        // encore modulation sur la suite si on y a intérêt
        while (itui!=U.end() and _modul<_mmax and Cur-itui->u > 0) {
            modul = std::min( EMAX(m_i, itui->i), _mmax-_modul );
            _modul += modul ;
            _cout += (itui->u) * modul ;
            if (enregistrer_plan) {
                m_plan->set_prod(itui->i, EMAX(m_i, itui->i)-modul) ;
            }
            itui = find_if( ++itui, U.end(), estDans ) ;
        }

        double _prod = _emax - _modul ;
        _r = arrivee.m_fuel + _prod - depart.m_fuel_avant_refuel ;
        _cout += Cur * _r ;
        if (enregistrer_plan) {
            m_plan->push_back(depart, arrivee, _r, _cout, Cur*_r, arrivee.m_fuel);
        }
        return _cout ;
    }

    /****************************************************************/
    /*    cas des campagnes de productions finales non BO           */
    /****************************************************************/

    double _seuil_augmentation_revente; // seuil de modulation au-delà duquel
                                        // on fait monter la revente, plutôt
                                        // que baisser l'achat

    // prod min et max pour l'arc
    double _emax = _EMAX(std::min(depart.m_date_campaign,H()), arrivee.m_date) ;
    double _emax_arc = _smax_d - BO(m_i,depart.m_campaign) ; // forcément positif si l'arc a été créé
    double _emin_arc = 0 ;
    double _e_augmentation_revente = std::max((double)0, _smin_d - BO(m_i, depart.m_campaign)) ;

    // modulation min et max
    double _mmax = std::min(MMAX(m_i,depart.m_campaign), _emax - _emin_arc) ;
    double _mmin = std::max((double)0, _emax - _emax_arc) ;

    if (depart.m_type == TYPE_PREMIER) {
        // pas de recharge, tout ce qu'on consomme est pris sur le prix de revente
        _r = 0 ;
        _seuil_augmentation_revente = _mmin ;
    } else {
        if (CiTp1(m_i)>Cur) {
            // on charge au maximum, tout ce qu'on consomme est pris sur le prix de revente
            _r = std::min( RMAX(m_i, depart.m_campaign), SMAX(m_i, depart.m_campaign)-depart.m_fuel_avant_refuel );
            _seuil_augmentation_revente = _mmin ;
        } else {
            // on charge au minimum, on commence par faire baisser la quantité de fuel achetée,
            // puis on augmente la quantité revendue
            _r = RMIN(m_i, depart.m_campaign);
            _seuil_augmentation_revente = std::max(_mmin, std::min(_mmax, _emax - _e_augmentation_revente));
        }
    }

    // modulation accomplie
    double _modul = 0 ;
    vector<Ui>::iterator itui = find_if( U.begin(), U.end(), estDans);
    // modulation tant qu'on y est forcé
    double modul = 0;
    while (_modul<_mmin) {
        modul = std::min( EMAX(m_i, itui->i), _mmin-_modul );
        _modul += modul ;
        _cout += (itui->u) * modul ;
        if (enregistrer_plan) {
            m_plan->set_prod(itui->i, EMAX(m_i, itui->i)-modul) ;
        }
        if (_modul<_mmin)
            itui = find_if( ++itui, U.end(), estDans ) ;
    }
    // annulation de la modulation courante si elle est acceptée dans la boucle suivante
    if (Cur-itui->u > 0) {
        _modul -= modul;
        _cout -= (itui->u) * modul ;
    }
    // encore modulation sur la suite si on y a intérêt
    while (_modul<_seuil_augmentation_revente and Cur-itui->u > 0) {
        modul = std::min( EMAX(m_i, itui->i), _seuil_augmentation_revente-_modul );
        _modul += modul ;
        _cout += (itui->u) * modul ;
        if (enregistrer_plan) {
            m_plan->set_prod(itui->i, EMAX(m_i, itui->i)-modul) ;
        }
        if (_modul<_seuil_augmentation_revente)
            itui = find_if( ++itui, U.end(), estDans ) ;
    }
    // annulation de la modulation courante si elle est acceptée dans la boucle suivante
    if (CiTp1(m_i)-itui->u > 0) {
        _modul -= modul;
        _cout -= (itui->u) * modul ;
    }
    // encore modulation sur la suite si on y a intérêt
    while (itui!=U.end() and _modul<_mmax and CiTp1(m_i)-itui->u > 0) {
        modul = std::min( EMAX(m_i, itui->i), _mmax-_modul );
        _modul += modul ;
        _cout += (itui->u) * modul ;
        if (enregistrer_plan) {
            m_plan->set_prod(itui->i, EMAX(m_i, itui->i)-modul) ;
        }
        itui = find_if( ++itui, U.end(), estDans ) ;
    }


    // calcul de la recharge
    _r -= std::max( (double)0, _seuil_augmentation_revente - _modul ) ;
    // calcul de la revente
    double _revente = depart.m_fuel_avant_refuel + _r - _emax + _modul ;

    // calcul du coût de l'arc
    _cout += (Cur * _r) // prix de la recharge
             - CiTp1(m_i) * _revente ; // prix de la revente

    if (enregistrer_plan) {
        m_plan->push_back(depart, arrivee, _r, _cout,
                              (Cur * _r),
                               _revente);
    }
    return _cout;
}