#include "ctype2.h"
#include "lecture.h"
#include "indices.h"
#include "cinstance.h"

#include <string>
#include <iterator>
#include <sstream>
#include <algorithm>

using namespace std ;
using instanceNs::lire_valeur;
using instanceNs::lire_liste ;

void CType2::lire( timestep_t T )
{
    if (                    lire_valeur<string>("begin")!="powerplant")
        throw string("section powerplant manquante");

    m_name =                lire_valeur<string>("name");
    if (                    lire_valeur<int>("type") != 2)
        throw string("centrale de type 2 attendue");
    i = (powerplants2_t)    lire_valeur<int>("index") ;
    XI =                    lire_valeur<kWh>("stock") ;
    K  =       (campaign_t) lire_valeur<campaign_t>("campaigns") ;
    // préparation du stockage des contraintes de type 13
    m_TO.resize(K, (week_t) CConstraint13::inconnu); // inconnu : pas de contrainte lue
    m_TA.resize(K, (week_t) CConstraint13::inconnu);

    //
    m_DA.resize(K);         lire_liste("durations",                  m_DA) ;
    m_CMMAX =               lire_valeur<kWh>("current_campaign_max_modulus");
    m_MMAX.resize(K);       lire_liste("max_modulus",                m_MMAX);
    m_RMAX.resize(K);       lire_liste("max_refuel" ,                m_RMAX);
    m_RMIN.resize(K);       lire_liste("min_refuel" ,                m_RMIN);
    m_Q   .resize(K);       lire_liste<coeff_t>("refuel_ratio",      m_Q);
    m_CBO  =                lire_valeur<kWh>("current_campaign_stock_threshold");
    m_BO  .resize(K+1);     lire_liste("stock_threshold",            m_BO);
    m_PMAX.resize(T);       lire_liste("pmax",                       m_PMAX) ;
    m_EMAX.resize(T);
    for (timestep_t t=(timestep_t)0; t<T; ++t) EMAX(t) = PMAX(t)*D(t) ;
    m_AMAX.resize(K);       lire_liste("max_stock_before_refueling", m_AMAX);
    m_SMAX.resize(K);       lire_liste("max_stock_after_refueling",  m_SMAX);
    m_C   .resize(K);       lire_liste("refueling_cost",             m_C);
#warning pseudo-lecture des coûts de déplacement des arrêts
    kWh _emax;
    for (timestep_t t; t<T; ++t) _emax += EMAX(t) ;
    double _emax_par_semaine = _emax.v / H().v;
    for (size_t k=0; k<K; ++k) {
        // déplacer d'une semaine coûte comme 0.2 semaines de production
	  m_C_post.push_back(euros_par_semaine(m_C[k].v*_emax_par_semaine * 0.2 ));
	  m_C_pre .push_back(euros_par_semaine(m_C[k].v*_emax_par_semaine * 0.2 ));
    }
    
    CiTp1  =                lire_valeur<euros_par_kWh>("fuel_price");
    m_CPB.lire(true,m_CBO) ;
    m_PB.resize(K) ;
    for (size_t k=0; (signed)k<K; k++)
    {
        CProfile & profile = m_PB[k] ;
        profile.lire(false,m_BO[k]);
        if ( profile.k != (signed)k )
            throw string("Indice du profile ne correspond pas");
    }

    if (lire_valeur<string>("end")!="powerplant")
        throw string("end powerplant manquant");
}



void CProfile::lire( bool current, kWh bo )
{
    if (current)
    {
        if (lire_valeur<string>("begin")!="current_campaign_profile")
            throw string("section current_campaign_profile manquante");
        k = (campaign_t) -1 ;
    } else {
        if (lire_valeur<string>("begin")!="profile")
            throw string("section profile manquante");
        k = (campaign_t) lire_valeur<int>("campaign_profile");
    }

    Np = (points_t) lire_valeur<int>("profile_points");
    vector<double> points = lire_liste<double>("decrease_profile", 0, 2*Np-1) ;
    CPoint dernier((kWh)-1,(coeff_t)-1);
    for (size_t p=0; (signed)p<Np; p++)
    {
        CPoint pt((kWh)points[2*p], (coeff_t)points[2*p+1]);
        if (pt.f!=dernier.f or pt.c!=dernier.c)
            m_points.push_back( pt );
        dernier = pt;
    }
    // Insertion éventuelle d'un premier point (f0,1)
    if (m_points.front().f<bo)
        m_points.insert(m_points.begin(), CPoint(bo, (coeff_t)1));


    if (current)
    {
        if (lire_valeur<string>("end")!="current_campaign_profile")
            throw string("end current_campaign_profile manquant");
    } else {
        if (lire_valeur<string>("end")!="profile")
            throw string("end profile manquant");
    }
}

void CType2::corrigeTO( vector<CConstraint13> & liste )
{
    
    week_t w = (week_t) 0;
    for ( campaign_t k = (campaign_t) 0; k<K; k = (campaign_t)( k+1 ) )
    {
        // si TO est indéfini ou inconnu, c'est qu'il n'a pas de limite inférieure
        // autre que le précédent TO + DA
        if (TO(k)==CConstraint13::inconnu or TO(k)==CConstraint13::indefini)
            TO(k) = w;
        w = (week_t) (TO(k) + DA(k)) ;
    }

    w = H() ;
    for ( campaign_t k = (campaign_t)(K-1) ; k>=(campaign_t)0; k = (campaign_t)( k-1 ) )
    {
        // si TA est indéfini ou inconnu, on le met après la dernière semaine
        if (TA(k)==CConstraint13::inconnu or TA(k)==CConstraint13::indefini
           )
            TA(k) = w;
        // si TA est différent de H, c'est que l'outage k doit être programmé.
        // les outages précédents doivent donc aussi être programmés
        if (TA(k)!=H() and k>(campaign_t)0)
                w = TA(k) - DA((campaign_t)(k-1)) ;
    }
    // ajout des contraintes implicites à la liste des contraintes
    for (campaign_t k = (campaign_t) 0; k<K; ++k )
    {
        if (std::find_if(m_constraints.begin(), m_constraints.end(), [k](CConstraint * c){return c->getType()==C13 and ((CConstraint13 *)c)->k==k;})==m_constraints.end())
        {
            CConstraint13 nouvelle_contrainte;
            nouvelle_contrainte.interpole(i, k, TO(k), TA(k));
            liste.push_back(nouvelle_contrainte);
            m_constraints.push_back(&(liste.back()));
        }
    }
}


/// production sur BO au temps t
kWh CType2::PBmin        ( campaign_t k, timestep_t t, kWh fuel)
{
    if ( fuel >= BO(k) )
        return EMAX(t);
    else
        return std::min(fuel, PB(k,t,fuel)*((coeff_t)1-epsilon()));
}

kWh CType2::PBmax        ( campaign_t k, timestep_t t, kWh fuel)
{
    if ( fuel >= BO(k) )
        return EMAX(t);
    else
        return std::min(fuel, PB(k,t,fuel)*((coeff_t)1+epsilon()));
}

CIntervallesBO CType2::PBlimits        ( campaign_t k, timestep_t t, kWh cut, kWh in_f1, kWh in_f2 )
{
    kWh f0,f1,f2 ;
    // faut-il couper la sortie en deux ?
    if ( cut > in_f1 and cut < in_f2) {
        f1 = in_f1 ; // [f1,cut] se projette en [f1,cut]
        f2 = cut ;
        f0 = in_f2 - PBmin( k,t,in_f2) ; // [cut,f2] se projette en [0, f2-PBmin(f2)]
            // si les deux images sont sécantes, il faut les rassembler
        if (f1<=f0) {
            f2 = std::max(f0,f2) ;
            f1 = (kWh)0 ;
            f0 = (kWh)-1 ;
        }
    } else {
        f0 = (kWh)-1 ;
        if (cut<=in_f1) {
            f1 = in_f1 - PBmax( k, t, in_f1 );
            f2 = in_f2 - PBmin( k, t, in_f2 );
        } else {
            f1 = in_f1 ;
            f2 = in_f2 ;
        }
    }
    return CIntervallesBO(f0,f1,f2);
}

kWh CType2::PBinverse0( campaign_t k, timestep_t t )
{
    kWh cut ;
    size_t ptindex2 = 1 ;
    vector<CProfile::CPoint> & prof = PB(k).m_points ;
    if (EMAX(t)==(kWh)0) return (kWh)INFINI;

    while (ptindex2<prof.size()) {
        CProfile::CPoint & p1 = prof[ptindex2-1];
        CProfile::CPoint & p2 = prof[ptindex2];
        cut = EMAX(t)* (( (p2.f-p1.f)*p1.c - p1.f*(p2.c-p1.c) ) / ((p2.f-p1.f)-EMAX(t)*(p2.c-p1.c))) ;
        if (p2.f<cut) return cut ;
        ptindex2++;
    }
    throw string("Inverse de 0 non trouvé dans CType2::PBinverse0");
}

CIntervallesBO CType2::PBlimits        ( campaign_t k, timestep_t t, CIntervallesBO in)
{
    kWh f0 ;
    // recherche de l'antécédent de 0
    kWh cut = PBinverse0(k,t);

    // cas 1 : un seul intervalle en entrée
    if (in.f0==(kWh)-1) {
        return PBlimits(k,t,cut,in.f1,in.f2);
    }
    // cas 2 : deux intervalles en entrée
    else {
        CIntervallesBO A = PBlimits(k,t,cut,in.f1 ,in.f2) ;
        CIntervallesBO B = PBlimits(k,t,cut,(kWh)0,in.f0) ;
        // assemblage des résultats
        // union de [0,A.f0] et [0,B.f0] => [0,f0]
        f0 = std::max(A.f0, B.f0==(kWh)-1 ? B.f2 : B.f0);
        // union éventuelle avec [B.f1,B.f2]
        if (B.f1 <= f0) {
            f0 = std::max(f0,B.f2);
        }
        // union éventuelle avec [A.f1,A.f2]
        // l'intervalle haut vient forcément de A
        if (A.f1 <= f0) {
            f0 = std::max(f0,A.f2);
            return CIntervallesBO((kWh)-1,(kWh)0,f0);
        } else {
            return CIntervallesBO(f0,A.f1,A.f2);
        }
    }
}



kWh CType2::PB           ( campaign_t k, timestep_t t, kWh fuel)
{
    vector<CProfile::CPoint> & prof = PB(k).m_points ;
    // trouver sur quelle partie de BO on est :
    if ( fuel > prof.front().f )
        return EMAX(t) ;
    else {
        size_t ptindex2 = 1 ;
        while (prof[ptindex2].f>=fuel and ptindex2+1<prof.size()) ++ptindex2 ;
        // test : est-on passé au-dessous du fuel actuel ?
        if (prof[ptindex2].f>fuel) {
            std::ostringstream oss;
            oss << fuel ;
            throw string("Pas de point dans le profil pour le niveau de fuel actuel : " + oss.str());
        }
        // on a trouvé le second point. Calcul du premier :
        CProfile::CPoint & pt2 = prof[ptindex2];
        CProfile::CPoint & pt1 = prof[ptindex2-1];
        // calcul de la production = fonction affine A VERIFIER !
        kWh prod = EMAX(t) * (pt1.c + ( ((fuel - pt1.f) * (pt2.c - pt1.c))/ (pt2.f - pt1.f))) ;
        return prod > fuel ? (kWh)0 : prod ;
    }
}


CType2::PBentree::PBentree( campaign_t k, kWh fuel, week_t h1, wlen_t duree ) {
    hash = duree.v + H().v * ( h1.v + H().v * ( k.v + ::K().v * fuel.v ) );
}


/// production sur BO sur une durée fixée, partant d'une semaine fixée
kWh CType2::PB           ( campaign_t k, kWh fuel, week_t h1, wlen_t duree )
{
    PBentree entree(k,fuel,h1,duree);
    map< PBentree, kWh >::iterator pte = m_pbentree_map.find(entree);
    if (pte!=m_pbentree_map.end()) {
	//cout << "Trouvé" << endl;
	return pte->second;
    }
    kWh p = (kWh)0 ;
    for (timestep_t t=week2timestep(h1); p<fuel and t<week2timestep(h1+duree); ++t)
    {
        p += PB(k,t,fuel-p);
    }
    m_pbentree_map[entree] = p;
    return p ;
}

/// durée minimale en semaines pour vider la centrale sur BO à partir d'une semaine fixée
// wlen_t CType2::longueur_BO  ( campaign_t k, double fuel, week_t h1 )
// {
//     week_t h ;
//     for (h=h1; fuel>0 and h<H(); ++h)
//     {
//         fuel -= PB( k, fuel, h, (wlen_t) 1 );
//     }
//     return h-h1 ;
// }


/// niveau de fuel permettant d'atteindre le niveau y après le temps t, sur BO
// double CType2::back_level( campaign_t k, timestep_t t, double y)
// {
//     vector<CProfile::CPoint> & prof = PB(k).m_points ;
//     // trouver sur quelle partie de BO on était avant d'atteindre x :
//     size_t ptindex2 = 1 ;
//     while (ptindex2<prof.size()) {
//         CProfile::CPoint & p1 = prof[ptindex2-1];
//         CProfile::CPoint & p2 = prof[ptindex2];
//         double a = EMAX(t)*(p2.c-p1.c)/(p2.f-p1.f);
//         double b = EMAX(t)*p1.c - p1.f*a ;
//         double x = (y+b)/(1-a);
//         if (a==1 or p2.f==p1.f)
//             throw string("Division par zéro dans back_level");
//         if (p2.f<x) return (x<BO(k) ? x : y+EMAX(t));
//         ptindex2++;
//     }
//     throw string("Echec du calcul du back_level");
// }

kWh CType2::back_level( campaign_t k, timestep_t t, kWh y, CIntervallesBO pre )
{
    vector<CProfile::CPoint> & prof = PB(k).m_points ;

    // était-on au-dessus de BO ?
    if (pre.f1==pre.f2) return pre.f1;

    // normalement, on était donc sur le profil imposé

    // peut-on être arrivé ici sans produire ?
    // i.e. est-on au-dessous du seuil de production précédent (y<antécédent(0))
    // et dans l'un des intervalles de pre ?
    kWh cut = PBinverse0(k,t);
    if (y<cut and pre.contient(y)) return y ;

    // donc on a produit pour arriver jusqu'à y
    // quel est l'intervalle d'où on aurait pu arriver ?
    kWh yxmax,yxmin;

    // trouver sur quelle partie de BO on était avant d'atteindre x :
    size_t ptindex2 = 1 ;
    CProfile::CPoint p1((kWh)-1,(coeff_t)-1), p2((kWh)-1,(coeff_t)-1);
    while (ptindex2<prof.size()) {
        p1 = prof[ptindex2-1];
        p2 = prof[ptindex2];

        coeff_t a = EMAX(t)*(p2.c-p1.c)/(p2.f-p1.f);
        kWh     b = EMAX(t)*p1.c - p1.f*a ;
        if (p2.f==p1.f)
            throw string("Division par zéro dans back_level");
        coeff_t denominateur  = (coeff_t)((coeff_t)1-a*((coeff_t)1+epsilon()));
        kWh     numerateur    = (y + (b * ((coeff_t)1+epsilon())));
        yxmax = numerateur / denominateur ;
        if (yxmax>p2.f)
            break ;
        ptindex2++ ;
    }
    while (ptindex2<prof.size()) {
        p1 = prof[ptindex2-1];
        p2 = prof[ptindex2];

        coeff_t a = EMAX(t)*(p2.c-p1.c)/(p2.f-p1.f);
        kWh     b = EMAX(t)*p1.c - p1.f*a ;
        if (p2.f==p1.f)
            throw string("Division par zéro dans back_level");
        coeff_t denominateur  = ((coeff_t)1-a*((coeff_t)1-epsilon()));
        kWh     numerateur    = (y + (b * ((coeff_t)1-epsilon())));
        yxmin = numerateur / denominateur ;
        if (yxmin<=p1.f)
            break ;
        ptindex2++ ;
    }
    yxmin = std::max(yxmin, cut) ;
    yxmax = std::min(yxmax, BO(k));

    // intersection avec [0,pre.f0]
    kWh yf0xmin = yxmin ;
    kWh yf0xmax = std::min( yxmax, pre.f0 ) ;
    kWh ciblef0 = pre.f0 * (coeff_t)0.5 ;
    kWh solf0 = ciblef0 < yf0xmin ? yf0xmin :
                ciblef0 > yf0xmax ? yf0xmax :
                ciblef0 ;
    // intersection avec [pre.f1,pre.f2]
    kWh yf1xmin = std::max( yxmin, pre.f1 ) ;
    kWh yf2xmax = std::min( yxmax, pre.f2 ) ;
    kWh ciblef1f2 = ( pre.f1+pre.f2 ) / (coeff_t)2 ;
    kWh solf1f2 = ciblef1f2 < yf1xmin ? yf1xmin :
                ciblef1f2 > yf2xmax ? yf2xmax :
                ciblef1f2;
    // choix de la solution en fonction de sa distance aux bornes de pre
    // (on cherche la solution la plus au milieu d'un intervalle)
    if (std::min(solf0, pre.f0-solf0) < std::min(solf1f2-pre.f1, pre.f2-solf1f2))
    {
#ifdef _DEBUG_BO_
        if (not pre.contient(solf1f2))
            throw string("Echec du calcul du back_level");
#endif
        return solf1f2 ;
    } else {
#ifdef _DEBUG_BO_
        if (not pre.contient(solf0))
            throw string("Echec du calcul du back_level");
#endif
        return solf0 ;
    }
}