#include "cgraphe.h"

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

#include <algorithm>
#include <functional>
#include <sstream>
#include <cmath>

using std::pair ;

const euros_par_kWh CGraphe::V(powerplants2_t i)
{ return (*m_u)[H()+i]; }

/**
 * 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 kWh perte_refuel( powerplants2_t i, campaign_t k, kWh x )
{
    return BO(i,(campaign_t)(k-1)) * ((Q(i,k)-(coeff_t)1)/Q(i,k)) - BO(i,k) + x / Q(i,k) ;
}

inline kWh fuel_depart( powerplants2_t i, campaign_t k, kWh fuel_avant_refuel )
{
    return  BO(i,(campaign_t)(k-1)) + (fuel_avant_refuel - BO(i,k)) * (Q(i,k)/(Q(i,k)-(coeff_t)1)) ;
}

/**
 * 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, kWh scourant, kWh & smin, kWh & smax)
{
    smin = scourant + RMIN(i, k) ;
    smax = std::min(scourant + RMAX(i, k), SMAX(i,k));
    if (smin > smax)
        throw ExceptionRefuelImpossible(i,k,scourant,smin,smax,RMIN(i,k),RMAX(i,k),SMAX(i,k)) ;
}

ExceptionRefuelImpossible::ExceptionRefuelImpossible(powerplants2_t i, campaign_t k, kWh s, kWh smin, kWh smax,
        kWh rmin, kWh rmax, kWh ssmax)
{
    std::ostringstream oss;
    oss << "ExceptionRefuelImpossible (i=" << (int)i << ", k=" << (int)k << ", s=" << s << ", rmin=" << rmin << ", rmax=" << rmax
            << ", smin=" << smin << ", smax=" << smax << ", smax(i,k)=" << ssmax << ")" ;
    append(oss.str()) ;
}

/**
 * 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 CIntervallesBO consommation_semaine_BO( powerplants2_t i, campaign_t k, week_t h,
                                     CIntervallesBO intervalles)
{
    CIntervallesBO reponse ;
    for (timestep_t t=week2timestep(h); t<week2timestep(h+(wlen_t)1); ++t)
    {
        reponse = type2(i).PBlimits(k,t,intervalles) ;
        intervalles = reponse ;
    }
    return reponse ;
}


/**
 *                  classe CGraph
 */


/* attributs de classe */

const std::vector<euros_par_kWh> *  CGraphe::m_u ;
std::vector<Ui>                     CGraphe::U ;
bool                                CGraphe::s_contraintes_cliques = false ;


/* 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 ;
    m_plan = NULL ;
//     std::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, (kWh)0);
    kWh somme = (kWh)0 ;
    for (week_t h=(week_t)0; h<H(); ++h) {
        somme+=EMAX(m_i,h);
        m_somme_EMAX[h+1] = somme;
    }
    m_pavage = ((somme==(kWh)0) ? (kWh) INFINI : ((somme / (coeff_t) H()) * FACTEUR_DE_PAVAGE ) ) ;
}

/**
 * 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) and h<H(); ++h)
        {
            m_campaigns[h].push_back(k) ;
        }
    }

    //
    // CREATION du premier noeud
    //

    creer_noeud((week_t)0, TYPE_PREMIER, (campaign_t)-1, XI(m_i) );

    // éventuel noeud BO en semaine 0
    if (XI(m_i)-EMAX(m_i,(week_t)0)<=BO(m_i, (campaign_t)-1))
        creer_noeud( (week_t)0, TYPE_BO, (campaign_t)-1, XI(m_i) );

    // éventuels outages en semaine 0
    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, (kWh)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 :
    kWh _smax_d ;
    kWh _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, CIntervallesBO((kWh)-1,_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 ) ;
    }
}


/**
 * faisceau_suite_BO
 * @param depart
 * @param h
 * @param smin
 * @param smax
 */
void CGraphe::faisceau_suite_BO( const CNoeud & depart, week_t h, CIntervallesBO limits )
{
    // production de la semaine
    // on évite les croisements, même temporaires, pour que les arcs soient gérables dans le filtre
    // aval avec la méthode CType2::back_level
    limits = consommation_semaine_BO( m_i, depart.m_campaign, h, limits );
    // création éventuelle de noeuds
    ++h ;
    // test sur la date
    if (h>=H()) return ;
    noeuds_eventuels_outage( depart, h, (kWh)0, limits.f0 ) ;
    noeuds_eventuels_outage( depart, h, limits.f1, limits.f2 ) ;
    // suite
    faisceau_suite_BO( depart, h, limits ) ;
}

#define _EMAX(h0,h1) (m_somme_EMAX[h1]-m_somme_EMAX[h0])
/**
 * faisceau_suite_modul
 * @param depart
 * @param h
 * @param smin_d
 * @param smax_d
 */
void CGraphe::faisceau_suite_modul( const CNoeud & depart, week_t h, kWh smin_d, kWh smax_d )
{
    if (h>H()) return ;
    // contrainte 13bis
    if ((campaign_t) (depart.m_campaign+1)<K() and h>TA(m_i,(campaign_t)(depart.m_campaign+1))) return ;

    // stock restant min et max :
    kWh smin = smin_d - _EMAX(depart.m_date_campaign, h) ;
    kWh smax = smax_d - std::max( (kWh) 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))
    {
        kWh smin_hmoins1_vers_BO = std::max( smin_d - _EMAX(depart.m_date_campaign,h-1),
                                             inf(BO(m_i,depart.m_campaign)) + EMAX(m_i,week2timestep(h)-(tlen_t)1) );

        kWh smax_hmoins1_vers_BO = std::min( smax_d - std::max( (kWh) 0,  _EMAX(depart.m_date_campaign,h-1) - MMAX(m_i, depart.m_campaign) ),
                                             inf(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, sup(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, kWh smin, kWh smax )
{
    kWh _smax_avant_refuel ;
    kWh _smax_depart ;
    kWh _smax1 ;
    kWh _smax2 ;
    kWh _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
#ifdef SAUT_CAMPAGNE_AUTORISE
              >
#else
              == 1 +
#endif
              depart.m_campaign ) // l'indice de campagne k doit croître strictement d'un outage à l'autre
        {
            _smax_avant_refuel = SMAX(m_i,*kit)-RMIN(m_i,*kit) ;
            _smax_depart = fuel_depart(m_i, *kit, _smax_avant_refuel);
            _smax1 = smax > AMAX(m_i,*kit) ? AMAX(m_i,*kit) : smax;
            _smax2 = _smax1 > _smax_depart ? _smax_depart : _smax1 ;
            _smax = _smax2 ;
//             std::cout << _smax_avant_refuel << smax << _smax1 << _smax2 << _smax << std::endl ;
            //(type==TYPE_OUTAGE ? fuel - perte_refuel(m_i,k,fuel) : fuel)

            creer_noeuds_intervalle(    h, TYPE_OUTAGE, *kit,
                                        smin, _smax,
                                        smin==(kWh)0 or smin==sup(BO(m_i,depart.m_campaign)),
                                        _smax==(kWh)0
                                                    or _smax==AMAX(m_i,*kit)
                                                    or _smax==_smax_depart);
        }
    }
}

/**
 * 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, kWh smin, kWh smax )
{
    if (smax<smin) return ;
    // création de noeuds de type BO
    creer_noeuds_intervalle( h, TYPE_BO, depart.m_campaign, smin, smax, false, false );
}

struct DansIntervalle {
    kWh m_min, m_max;
    DansIntervalle(kWh min, kWh max) {m_min=min; m_max=max;}
    bool operator()(const kWh & x) {return (x >= m_min) and (x <= m_max);}
};
struct PairDansIntervalle {
    kWh m_min, m_max;
    PairDansIntervalle(kWh min, kWh max) {m_min=min; m_max=max;}
    bool operator()(const pair<kWh, CNoeud *> & x) {return (x.first >= m_min) and (x.first <= m_max);}
};
struct HorsIntervalle {
    kWh m_min, m_max;
    HorsIntervalle(kWh min, kWh max) {m_min=min; m_max=max;}
    bool operator()(const kWh & 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, kWh min, kWh 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)/(coeff_t)2 );

    // liste des niveaux existants compris entre min et max
    vector<kWh> niveaux_existants ;
    for (std::map<kWh,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
    kWh 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
    kWh 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<kWh>::iterator it = niveaux_existants.begin()+1; it!=niveaux_existants.end(); it++ )
        {
            kWh second = *it ;
            coeff_t nb_noeuds = (coeff_t) floor(((second-premier) / pavage(h)).v) ;
            kWh pas = (second-premier) / (nb_noeuds+(coeff_t)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
 */
kWh CGraphe::pavage(week_t)
{
    return m_pavage ;
//     return EMAX(m_i,h)==(kWh)0 ? (kWh)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, kWh 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 ;
}




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

void CGraphe::creerArcs()
{
#ifdef VERBOSE_CREATION_NOEUD
    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());

#ifdef VERBOSE_CREATION_NOEUD
    for (size_t i=0; i<m_noeuds.size(); i++)
    {
        CNoeud * noeud_ptr = m_noeuds[i];
        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

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

    // parcours des noeuds de départ

#ifdef VERBOSE_CREATION_ARCS
    size_t nb_arcs_precedent = -1 ;
#endif
    for (std::vector<CNoeud *>::const_iterator dep=m_noeuds.begin(); dep!=m_noeuds.end() and (*dep)->m_type!=TYPE_DERNIER; 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 = **dep;

        // niveaux de fuel possibles au départ de la production
        kWh _smax_d ;
        kWh _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)
        {
            CIntervallesBO limits((kWh)-1, _smin_d, _smax_d ) ;
            // semaine du noeud d'arrivée courant
            week_t _semaine_atteinte = depart.m_date_campaign ;
            for (std::vector<CNoeud *>::const_iterator arr=dep+1; arr!=m_noeuds.end() and (*arr)->m_type!=TYPE_DERNIER; arr++)
            {
                CNoeud & arrivee = **arr;

                // vérification des types et des numéros de campagne
                if ( (arrivee.m_type!=TYPE_OUTAGE) or not (arrivee.m_campaign
#ifdef SAUT_CAMPAGNE_AUTORISE
                      >
#else
                      == 1 +
#endif
                      depart.m_campaign) )
                        continue ;

                // on ne branche pas un noeud BO vers un noeud de fuel>=BO
                if (arrivee.m_fuel>=BO(m_i,depart.m_campaign))
                    continue ;

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

                // ajout de l'arc, s'il est réalisable
                if ( limits.contient(arrivee.m_fuel) )
                {
                    // vérification
                    m_arcs.push_back( CArc( *dep, dep-m_noeuds.begin(), *arr, arr-m_noeuds.begin() ) ) ;
                }
            }
            // dernier_noeud
            if (
                // contrainte 13bis
                (depart.m_campaign+1==K() or (depart.m_campaign+1<K() and H()<=TA(m_i,(campaign_t)(depart.m_campaign+1)))))
                m_arcs.push_back( CArc( *dep, dep-m_noeuds.begin(), m_noeuds.back(), m_noeuds.size()-1) ) ;
        }
        else if (depart.m_type == TYPE_OUTAGE)// type TYPE_OUTAGE
        {
            // parcours des noeuds d'arrivée
            for (std::vector<CNoeud *>::const_iterator arr=dep+1; arr!=m_noeuds.end() and (*arr)->m_type!=TYPE_DERNIER; arr++)
            {
                CNoeud & arrivee = **arr;

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

                // vérification du stock d'arrivée
                if ( arrivee.m_fuel<sup(BO(m_i,depart.m_campaign)))
                    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 not (arrivee.m_campaign
#ifdef SAUT_CAMPAGNE_AUTORISE
                      >
#else
                      == 1 +
#endif
                      depart.m_campaign)) )
                            continue ;

                // stocks possibles à l'arrivée :
                kWh _smin = _smin_d - _EMAX(depart.m_date_campaign, arrivee.m_date) ;
                kWh _smax = _smax_d - std::max( (kWh) 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 < sup(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(*dep, dep-m_noeuds.begin(), *arr, arr-m_noeuds.begin()) ) ;
                }
            }
            // dernier noeud
            kWh _smax = _smax_d
                        - std::max( (kWh) 0,
                                     _EMAX(std::min(H(),depart.m_date_campaign), H())
                                             - MMAX(m_i, depart.m_campaign) ) ;

            if (_smax>=sup(BO(m_i,depart.m_campaign)) and
                // contrainte 13bis
                (depart.m_campaign+1==K() or (depart.m_campaign+1<K() and H()<=TA(m_i,(campaign_t)(depart.m_campaign+1)))))
            {
                m_arcs.push_back( CArc(*dep, dep-m_noeuds.begin(), m_noeuds.back(), m_noeuds.size()-1) ) ;
            }
        }
        else // type TYPE_PREMIER
        {
            // parcours des noeuds d'arrivée
            for (std::vector<CNoeud *>::const_iterator arr=dep+1; arr!=m_noeuds.end() and (*arr)->m_type!=TYPE_DERNIER; arr++)
            {
                CNoeud & arrivee = **arr;

                // vérification du stock d'arrivée
                if ( arrivee.m_date>0 and arrivee.m_fuel<sup(BO(m_i,depart.m_campaign)))
                    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 not (arrivee.m_campaign
#ifdef SAUT_CAMPAGNE_AUTORISE
                      >
#else
                      == 1 +
#endif
                      depart.m_campaign)) )
                    continue ;

                // stocks possibles à l'arrivée :
                kWh _smin = _smin_d - _EMAX(depart.m_date_campaign, arrivee.m_date) ;
                kWh _smax = _smax_d - std::max( (kWh) 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 < sup(BO(m_i,depart.m_campaign)) and arrivee.m_date>0)
                    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(*dep, dep-m_noeuds.begin(), *arr, arr-m_noeuds.begin()) ) ;
                }
            }
            // dernier noeud
            kWh _smax = _smax_d
                    - std::max( (kWh) 0,
                                 _EMAX(std::min(H(),depart.m_date_campaign), H())
                                         - MMAX(m_i, depart.m_campaign) ) ;

            if (_smax>=sup(BO(m_i,depart.m_campaign)) and
                // contrainte 13bis
                (depart.m_campaign+1==K() or (depart.m_campaign+1<K() and H()<=TA(m_i,(campaign_t)(depart.m_campaign+1)))))
            {
                m_arcs.push_back( CArc(*dep, dep-m_noeuds.begin(), m_noeuds.back(), m_noeuds.size()-1) ) ;
            }
        }
    }

    // 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());
    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", (int)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()
{
    // création de somme_ui_EMAX
    m_somme_ui_EMAX.resize ( H() + 1 ) ;
    m_somme_ui_EMAX[0] = (euros)0 ;
    euros somme = (euros)0 ;

    for ( week_t h = ( week_t ) 1; h <= H(); ++h )
    {
        somme += u( h - (wlen_t)1 ) * EMAX ( m_i, (week_t) (h - 1) ) ;
        m_somme_ui_EMAX[h] = somme ;
    }

#ifdef _MODUL_REPARTIE_
    m_somme_ui_EMAX_cout_positif = vector<vector<euros> >( K()+2,
            vector<euros>(H()+1, (euros)0 ) );
    m_somme_EMAX_cout_positif = vector<vector<kWh> >( K()+2,
            vector<kWh>(H()+1, (kWh)0 ) );
    for (campaign_t k = (campaign_t)(-1); k<K()+1; ++k)
    {
        euros_par_kWh Cur ;
        if (k<K()) Cur = C(m_i,k);
        else Cur = CiTp1(m_i);

        m_somme_ui_EMAX_cout_positif[k+1].resize( H() + 1 );
        m_somme_EMAX_cout_positif   [k+1].resize( H() + 1 );
        m_somme_ui_EMAX_cout_positif[k+1][0] = (euros)0 ;
        m_somme_EMAX_cout_positif   [k+1][0] = (kWh)0 ;

        euros   somme_ui_EMAX   = (euros) 0 ;
        kWh     somme_EMAX      = (kWh)   0 ;
        for ( week_t h = ( week_t ) 1; h <= H(); ++h )
        {
            if (Cur - (*m_u)[h-1] >= (euros_par_kWh) 0) {
                somme_ui_EMAX   += u( h - (wlen_t) 1 ) * EMAX ( m_i, h - (wlen_t)1 ) ;
                somme_EMAX      +=                       EMAX ( m_i, h - (wlen_t)1 ) ;
            }
            m_somme_ui_EMAX_cout_positif[k+1][h] = somme_ui_EMAX ;
            m_somme_EMAX_cout_positif   [k+1][h] = somme_EMAX ;
        }
    }
#endif

    /* 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<euros_par_kWh> * _u ) : u(_u) {
        n = 0 ;
    }
    Ui operator()() {
        return Ui( (*u)[n], (week_t) n++);
    }
    const vector<euros_par_kWh> * u ;
    size_t n ;
};

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

void CGraphe::prepare_viterbi(const vector<euros_par_kWh> * u, bool contraintes_cliques)
{
    s_contraintes_cliques = contraintes_cliques ;

    // 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);
}

void CGraphe::calcule_couts_des_arcs( )
{
    for (size_t noeud=0; noeud<m_noeuds.size(); noeud++)
    {
        calcule_cout_contraintes( m_noeuds[noeud] );
    }

    for (size_t arc=0; arc<m_noeud_fin.size(); arc++)
    {
#ifdef DEBUG_COUT_ARCS
        std::cout << "centrale n°" << m_i << "\tcout arc n° " << arc ;
#endif
        m_cout[arc] = calcule_cout_arc<false>( arc ) ;
#ifdef DEBUG_COUT_ARCS
    std::cout << " = " << m_cout[arc].contraintes << ", " << m_cout[arc].prix << std::endl ;
#endif
    }
}

void CGraphe::calcule_cout_contraintes( CNoeud * noeud )
{
    if (s_contraintes_cliques and noeud->m_type == TYPE_OUTAGE) {
        noeud->m_cout_contraintes = CConstraint::cout_blocage_contraintes_13_a_21(m_i, noeud);
    } else {
        noeud->m_cout_contraintes = 0;
    }
}

/**
 * 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 = cout_t() ;
    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;
} ;


#ifndef _MODUL_REPARTIE_

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 euros_par_kWh Cur = depart.m_campaign == -1 ? (euros_par_kWh)0 : C(m_i, depart.m_campaign) ;
    euros _cout = 0 ;
    kWh _r = 0; // recharge
    kWh _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 {
                    kWh prod = PB( m_i, depart.m_campaign, _fuel, h, (wlen_t) 1 );
                    _fuel -= prod ;
                    _cout -= (*m_u)[h] * prod ; // cout reduit <- (-somme(u_h*p_h))
                    if (enregistrer_plan) {
                        m_plan->set_prod(h,prod) ;
                    }
                }
            }
            _cout += Cur * _r ; // cout reduit <- (-somme(u_h*p_h)) + Cur * R
            if (enregistrer_plan) {
                m_plan->push_back(depart, arrivee,
                                  _r, arrivee.m_type!=TYPE_DERNIER ? arrivee.m_fuel : _fuel,
                                  _cout, Cur*_r);
            }
            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
    kWh _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 )
    {
                                                                                // cout reduit <- -(somme(u(h)*EMAX(h)))
        _cout = - (m_somme_ui_EMAX[arrivee.m_date] - m_somme_ui_EMAX[depart.m_date_campaign]) ;

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

        // modulation accomplie
        kWh _modul = 0 ;
        vector<Ui>::iterator itui = find_if( U.begin(), U.end(), estDans);
        // modulation tant qu'on y est forcé
        kWh modul = 0;
        while (_modul<_mmin) {
            modul = std::min( EMAX(m_i, itui->i), _mmin-_modul );
            _modul += modul ;
                                                 // cout reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
            _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
            _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
            _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 ) ;
        }

        kWh _prod = _emax - _modul ;
        _r = arrivee.m_fuel + _prod - depart.m_fuel_avant_refuel ;
                                                 // cout reduit <- Cu * R - (somme(u(h)*prod(h)))
        _cout += Cur * _r ;
        if (enregistrer_plan) {
            m_plan->push_back(depart, arrivee,
                              _r, arrivee.m_fuel,
                              _cout, Cur*_r);
        }
        return _cout ;
    }

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

    kWh _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
    kWh _emax = _EMAX(std::min(depart.m_date_campaign,H()), arrivee.m_date) ;
    kWh _emax_arc = _smax_d - BO(m_i,depart.m_campaign) ; // forcément positif si l'arc a été créé
    kWh _emin_arc = 0 ;
    kWh _e_augmentation_revente = std::max((kWh)0, _smin_d - BO(m_i, depart.m_campaign)) ;

    // modulation min et max
    kWh _mmax = std::min(MMAX(m_i,depart.m_campaign), _emax - _emin_arc) ;
    kWh _mmin = std::max((kWh)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));
        }
    }

                                                                                // cout reduit <- -(somme(u(h)*EMAX(h)))
    if (arrivee.m_date > depart.m_date_campaign)
        _cout = - (m_somme_ui_EMAX[arrivee.m_date] - m_somme_ui_EMAX[depart.m_date_campaign]) ;
    // modulation accomplie
    kWh _modul = 0 ;
    vector<Ui>::iterator itui = find_if( U.begin(), U.end(), estDans);
    // modulation tant qu'on y est forcé
    kWh modul = 0;
    while (_modul<_mmin) {
        modul = std::min( EMAX(m_i, itui->i), _mmin-_modul );
        _modul += modul ;
                                                 // cout reduit <-  -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
        _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
        _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
        _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
        _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 reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
        _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( (kWh)0, _seuil_augmentation_revente - _modul ) ;
    // calcul de la revente
    kWh _revente = depart.m_fuel_avant_refuel + _r - _emax + _modul ;

    // calcul du coût de l'arc
                                                 /* cout reduit <- -(somme(u(h)*EMAX(h))) + somme(u(h)*(EMAX(h)-prod(h))
                                                         + Cu * R - C(t+1) * reste
                                                    = Cu * R - C(t+1) * reste - somme(u(h)*prod(h))
                                                 */
    _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, _revente,
                          _cout, (Cur * _r));
    }
    return _cout;
}


#else

template<bool enregistrer_plan>
        cout_t CGraphe::calcule_cout_arc( size_t arc_idx )
{
    // noeud de départ
    const CNoeud & depart = *m_arcs[arc_idx].m_de ;
    const CNoeud & arrivee = *m_arcs[arc_idx].m_vers ;
    const euros_par_kWh Cur = depart.m_campaign == -1 ? (euros_par_kWh) 0 : C(m_i, depart.m_campaign) ;
    euros _ui_prod = (euros) 0 ;
    kWh _r = (kWh)0; // recharge
    kWh _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<=(kWh) 0) break ;
            if (enregistrer_plan and _fuel<=(kWh) 0)
                m_plan->set_prod(h,(kWh)0);
            else {
                kWh prod = PB( m_i, depart.m_campaign, _fuel, h, (wlen_t) 1 );
                _fuel -= prod ;
                _ui_prod += u(h) * prod ; // cout reduit <- (-somme(u_h*p_h))
                if (enregistrer_plan) {
                    m_plan->set_prod(h,prod) ;
                }
            }
        }
        if (enregistrer_plan) {
            m_plan->push_back(depart, arrivee,
                              _r, arrivee.m_type!=TYPE_DERNIER ? arrivee.m_fuel : _fuel,
                              _ui_prod, Cur*_r);
        }
        return cout_t(depart.m_cout_contraintes, Cur*_r - _ui_prod) ;
    }

    /*****************************************************/
    /*    cas des campagnes de productions non BO        */
    /*****************************************************/
    if (arrivee.m_type!=TYPE_DERNIER)
        if ( arrivee.m_fuel < BO(m_i, depart.m_campaign)
             and not (depart.m_type==TYPE_PREMIER and arrivee.m_date==0) )
            throw string("Arrivée de l'arc avec un fuel inférieur à 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
    kWh _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,(kWh)0) ;
        }
    }
    /*****************************************************/
    /*    cas des campagnes de productions non finales   */
    /*****************************************************/
    // calcul des productions et du coût
    if ( arrivee.m_type != TYPE_DERNIER )
    {
        // prod min et max sur la durée
        kWh _emax = _EMAX(depart.m_date_campaign, arrivee.m_date) ;
        //kWh _emin = std::max( (kWh) 0,  _emax - MMAX(m_i, depart.m_campaign) ) ;
        // prod min et max pour l'arc
        kWh _emax_arc = _smax_d - arrivee.m_fuel ; // forcément positif car l'arc a été créé
        kWh _emin_arc = std::max( (kWh) 0, _smin_d - arrivee.m_fuel ) ;
        // modulation min et max
        kWh _mmax = std::min(MMAX(m_i,depart.m_campaign), _emax - _emin_arc) ;
        kWh _mmin = std::max((kWh)0, _emax - _emax_arc) ;

        // somme(u(h)*EMAX(h) :
        euros _somme_ui_EMAX = (m_somme_ui_EMAX[arrivee.m_date] - m_somme_ui_EMAX[depart.m_date_campaign]) ;

        // energie disponible sur les périodes de cout réduit positif
        kWh _e_cout_reduit_positif = somme_EMAX_cout_positif(depart.m_campaign, depart.m_date_campaign, arrivee.m_date) ;
        // coefficient à appliquer pour rammener cette énergie à _mmax
        coeff_t _coeff_cout_reduit_positif = (_e_cout_reduit_positif==(kWh) 0 ? (coeff_t)1 : std::min((coeff_t) 1, _mmax / _e_cout_reduit_positif)) ;
        // coût réduit de la modulation sur les périodes de coût réduit positif
        euros _cout_reduit_modul_positif = somme_ui_EMAX_cout_positif(depart.m_campaign, depart.m_date_campaign, arrivee.m_date)
                                    * _coeff_cout_reduit_positif;

        // modulation restante à répartir sur les périodes de coût réduit négatif
        kWh _modul_restante = std::max((kWh)0, _mmin - _e_cout_reduit_positif);
        // energie disponible sur les périodes de cout réduit négatif
        kWh _e_cout_reduit_negatif = _emax - _e_cout_reduit_positif ;
        // coefficient à appliquer sur ces périodes pour atteindre la modulation restante
        coeff_t _coeff_cout_reduit_negatif = (_e_cout_reduit_negatif==(kWh) 0 ? (coeff_t)1 : _modul_restante / _e_cout_reduit_negatif) ;
        // coût réduit disponible sur les périodes de coût réduit négatif
        euros _cout_reduit_dispo_negatif = _somme_ui_EMAX
                    - somme_ui_EMAX_cout_positif(depart.m_campaign, depart.m_date_campaign, arrivee.m_date) ;
        // coût réduit de la modulation sur les périodes de coût réduit négatif
        euros _cout_reduit_modul_negatif = _cout_reduit_dispo_negatif * _coeff_cout_reduit_negatif ;

        // production sur l'arc
        kWh _prod = _e_cout_reduit_positif * ((coeff_t)1-_coeff_cout_reduit_positif)
                + _e_cout_reduit_negatif * ((coeff_t)1-_coeff_cout_reduit_negatif) ;
        // recharge
        _r = arrivee.m_fuel + _prod - depart.m_fuel_avant_refuel ;

        // coût réduit de l'arc
        _ui_prod = (_somme_ui_EMAX - (_cout_reduit_modul_positif + _cout_reduit_modul_negatif) );

        if (enregistrer_plan) {
            for (week_t h=depart.m_date_campaign; h<arrivee.m_date; ++h)
            {
                m_plan->set_prod(h, EMAX(m_i, h) * ((coeff_t)1 - (Cur-u(h) < (euros_par_kWh) 0 ? _coeff_cout_reduit_negatif : _coeff_cout_reduit_positif ) ) ) ;
            }
            m_plan->push_back(depart, arrivee,
                              _r, arrivee.m_fuel,
                              _ui_prod, Cur*_r);
        }
        return cout_t(depart.m_cout_contraintes, Cur*_r - _ui_prod) ;
    }

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

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

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

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

    // prix unitaire de la revente
    euros_par_kWh C_tp1 = CiTp1(m_i) ;

    if (depart.m_type == TYPE_PREMIER) {
        // pas de recharge, tout ce qu'on consomme est pris sur le prix de revente
        _r = (kWh)0 ;
        _mseuil = _mmin ;
    } else {
        if (Cur<C_tp1) {
            // 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 );
            _mseuil = _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);
            _mseuil = std::max(_mmin, std::min(_mmax, _emax - _e_augmentation_revente));
        }
    }

                                                                                // cout reduit <- -(somme(u(h)*EMAX(h)))
    week_t date_depart = std::min(arrivee.m_date, depart.m_date_campaign);

    // somme(u(h)*EMAX(h) :
    euros _somme_ui_EMAX = (m_somme_ui_EMAX[arrivee.m_date] - m_somme_ui_EMAX[date_depart]) ;

    // energie disponible sur les périodes de cout réduit sur Ctp1 positif
    kWh _e_cout_reduit_positif_tp1 = somme_EMAX_cout_positif(K_Tp1, date_depart, arrivee.m_date) ;
    // coefficient à appliquer pour rammener cette énergie à _mmax
    coeff_t _coeff_cout_reduit_positif_tp1 = (_e_cout_reduit_positif_tp1==(kWh)0 ? (coeff_t) 1 : std::min((coeff_t) 1, _mmax / _e_cout_reduit_positif_tp1)) ;
    // coût réduit disponible sur les périodes de coût réduit sur Ctp1 positif
    euros _cout_reduit_dispo_positif_tp1 = somme_ui_EMAX_cout_positif(K_Tp1, date_depart, arrivee.m_date);
    // coût réduit de la modulation sur les périodes de coût réduit sur CTp1 positif
    euros _cout_reduit_modul_positif_tp1 = _cout_reduit_dispo_positif_tp1 * _coeff_cout_reduit_positif_tp1;

    // modulation restante à répartir sur les périodes de coût réduit sur cur positif et sur ctp1 négatif
    kWh _modul_restante_cur = std::max((kWh)0, _mseuil - _e_cout_reduit_positif_tp1);
    // energie disponible sur les périodes de cout réduit sur cur positif et sur ctp1 négatif
    kWh _e_cout_reduit_positif_cur = std::max((kWh)0, somme_EMAX_cout_positif(depart.m_campaign, date_depart, arrivee.m_date) - _e_cout_reduit_positif_tp1) ;
    // coefficient à appliquer sur ces périodes pour atteindre la modulation restante
    coeff_t _coeff_cout_reduit_positif_cur = (_e_cout_reduit_positif_cur==(kWh)0 ? (coeff_t)1 : _modul_restante_cur / _e_cout_reduit_positif_cur) ;
    // coût réduit disponible sur les périodes de coût réduit sur cur positif et sur ctp1 négatif
    euros _cout_reduit_dispo_positif_cur = std::max((euros)0, somme_ui_EMAX_cout_positif(depart.m_campaign, date_depart, arrivee.m_date) - _cout_reduit_dispo_positif_tp1 ) ;
    // coût réduit de la modulation sur les périodes de coût réduit sur cur positif et sur ctp1 négatif
    euros _cout_reduit_modul_positif_cur = _cout_reduit_dispo_positif_cur * _coeff_cout_reduit_positif_cur ;

    // modulation restante à répartir sur les périodes de coût réduit négatif
    kWh _modul_restante = std::max((kWh)0, _mmin - _e_cout_reduit_positif_cur - _e_cout_reduit_positif_tp1);
    // energie disponible sur les périodes de cout réduit négatif
    kWh _e_cout_reduit_negatif = _emax - _e_cout_reduit_positif_cur - _e_cout_reduit_positif_tp1 ;
    // coefficient à appliquer sur ces périodes pour atteindre la modulation restante
    coeff_t _coeff_cout_reduit_negatif = (_e_cout_reduit_negatif==(kWh)0 ? (coeff_t)1 : _modul_restante / _e_cout_reduit_negatif) ;
    // coût réduit disponible sur les périodes de coût réduit négatif
    euros _cout_reduit_dispo_negatif = _somme_ui_EMAX
                - _cout_reduit_dispo_positif_cur
                - _cout_reduit_dispo_positif_tp1 ;
    // coût réduit de la modulation sur les périodes de coût réduit négatif
    euros _cout_reduit_modul_negatif = _cout_reduit_dispo_negatif * _coeff_cout_reduit_negatif ;

    // production sur l'arc
    kWh _prod    = _e_cout_reduit_positif_cur * ((coeff_t)1-_coeff_cout_reduit_positif_cur)
                 + _e_cout_reduit_positif_tp1 * ((coeff_t)1-_coeff_cout_reduit_positif_tp1)
                 + _e_cout_reduit_negatif     * ((coeff_t)1-_coeff_cout_reduit_negatif    ) ;
    // calcul de la recharge
    if (Cur<C_tp1) {
            // 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 );
    } else {
            // on charge au minimum, on commence par faire baisser la quantité de fuel achetée,
            // puis on augmente la quantité revendue
        _r = std::max( RMIN(m_i, depart.m_campaign),
                       sup(BO(m_i, depart.m_campaign)) + _prod - depart.m_fuel_avant_refuel );
    }
    // calcul de la revente
    kWh _revente = depart.m_fuel_avant_refuel + _r - _prod ;

    // coût réduit de l'arc
    _ui_prod = (_somme_ui_EMAX - (_cout_reduit_modul_positif_cur + _cout_reduit_modul_positif_tp1 + _cout_reduit_modul_negatif) );

    // enregistrement du plan
    if ( enregistrer_plan )
    {
        for ( week_t h = date_depart; h < arrivee.m_date; ++h )
        {
            m_plan->set_prod ( h, EMAX ( m_i, h ) * ( (coeff_t)1 - (
                                   C_tp1 - (*m_u)[h] >= (euros_par_kWh) 0 ?
                                   _coeff_cout_reduit_positif_tp1
                                   : Cur - (*m_u)[h] >= (euros_par_kWh) 0 ?
                                   _coeff_cout_reduit_positif_cur
                                   : _coeff_cout_reduit_negatif ) ) ) ;
        }

        m_plan->push_back ( depart, arrivee,
                            _r, _revente,
                            _ui_prod, Cur*_r );
    }
    return cout_t(depart.m_cout_contraintes, Cur*_r - C_tp1 * _revente - _ui_prod);
}

#endif