#include "DA.h"
// #include <algorithm>
#ifndef EPSILON
#define EPSILON 0.001
#endif
// #include <cmath>

#include <cinstance.h>
#include <cinstancereduite.h>
#include <algorithm>

static double s_infini = INFINI ;
using namespace std;

DualAscent::DualAscent(vector<kWh>  dem, powerplants2_t nb_centrales)
{
    b = dem ;
    b_EMAX = dem;
        // Ajout Affectation sur b et b_EMAX
    b_EMAX.resize( dem.size() + nb_centrales , (kWh)1 );
    b.resize( dem.size() + nb_centrales , (kWh)1 );

    euros coutTherm = (euros) 0;
    euros b_PMAX_coutTherm = (euros) 0;

    for (week_t w=(week_t) 0;w<H();w=(week_t) (w+1) )
    {
            // Remplissage de b_EMAX
        kWh sumEMAXweek = (kWh) 0;
        for (powerplants2_t i = (powerplants2_t) 0 ; i < nb_centrales ; i = (powerplants2_t) (i+1) )
            sumEMAXweek += CInstanceReduite::EMAX(i,w);
            // std::cout << "sumPMAXweek = " << sumPMAXweek << std::endl;
        b_EMAX[w] =sumEMAXweek;

            // remplissage de coutTherm et de b_PMAX_coutTherm
        coutTherm += CInstanceReduite::ThCout(w) * b[w];
            // std::cout<<"w = "<<w << " & b_EMAX.size() = "<<b_EMAX.size() << std::endl;
        b_PMAX_coutTherm += CInstanceReduite::ThCout(w) * b_EMAX[w];
    }

    m_lambda.resize( H()+I() );
    m_u.resize( H()+I() );

    std::vector<CPlan::CArc> vide;
#ifdef PMAX_POUR_DEM
    ajoutcolonne(&b_EMAX, b_PMAX_coutTherm, vide, 0);
#else
    ajoutcolonne(&b, coutTherm, vide, (kWh)0);
#endif
}


void DualAscent::ajoutcolonne( vector<kWh> * colonne, euros cou, std::vector<CPlan::CArc>& arcs, kWh resteCol )
{
    A.push_back(colonne);
    cost.push_back(cou);
    vector<kWh> quantiteRefuelee;
    for ( size_t i=0 ; i < arcs.size() ; i++ )
    {
        if ( arcs[i].m_type == TYPE_OUTAGE )
        {
            quantiteRefuelee.push_back(arcs[i].m_recharge);
        }
    }
    recharge.push_back(quantiteRefuelee);
    reste.push_back(resteCol);
}

void DualAscent::ajoutcolonne(CPlan* cplan)
{
    // ajout du plan à la liste des plans
    m_plans.push_back(cplan);

    cout << "Cout réduit de la nouvelle colonne : "<<cplan->m_cout_reduit<<" pour la centrale "<< cplan->m_i << endl;
    ajoutcolonne(&cplan->m_colonne, cplan->m_cout, cplan->m_arcs, cplan->m_reste);
}

/// retrait d'un plan de production nucléaire
void DualAscent::retraitColonne(const CPlan* plan)
{
    size_t indice = std::find(m_plans.begin(), m_plans.end(), plan) - m_plans.begin() + 1 ;
    m_plans.erase(m_plans.begin()+indice-1);
    A.erase(A.begin()+indice);
    cost.erase(cost.begin()+indice);
    recharge.erase(recharge.begin()+indice);
    reste.erase(reste.begin()+indice);
}

vector<euros_par_kWh> DualAscent::calculeDuals() const
{
//	lambda = vector<double>(lambda.size(),0);
	int nbCol = cost.size();
	int nbLig = b.size();
	vector<kWh> q = calculePoids();
	vector<euros> L(nbCol, (euros) 0);
	coeff_t pas = (coeff_t) 1;
	euros meilleurZd = (euros) - s_infini;
    vector<euros_par_kWh> meilleursU(m_u);
    vector<euros_par_kWh> meilleursLambda(m_lambda);

#ifdef VERBOSE_COUT_COLONNES
    cout << "Cout des colonnes : ";
	for (int c = 0; c < cost.size(); c++)
        cout << cost[c] << "\t" ;
    cout << endl;
#endif

	// Recherche d'un borne sup
#ifndef BORNE_SUP_VERSION_2
	if (UB < 1)
	{
		UB = (double) testUB();
	}

	if (UB>0)
	{
		if (UB > 1)
		{
			UB = min(UB,recherche_UB());
		}
		else
		{
			UB = recherche_UB();
		}
	}
#else
	if (m_UB > (euros) 0)
        m_UB = min(m_UB,recherche_UB2(q));
	else
        m_UB = recherche_UB2(q);
#endif
	int nbItMax = 3000;
	int nbItMaxSansAmelioration = 100;
	int derniereItAmeliorant = 0;
	int nbIt_ = 0;
	for ( int nbIt=0 ; (nbIt < nbItMax) && (derniereItAmeliorant > nbIt-nbItMaxSansAmelioration) ; nbIt++ )
	{
		nbIt_=nbIt;
        m_u = DAcalc(m_lambda, L, q);
		vector<coeff_t> x(nbCol, (coeff_t)0);
        vector<kWh> s = sousgrad(x, m_u, m_lambda, L, q);

		//Calcul de Zd
		euros zd = (euros)0;
        kWh_carre normeS = (kWh_carre)0;
		for ( int i=0 ; i < nbLig ; i++ )
		{
#ifdef _DEBUG_LAMBDA
			if (nbIt ==0)
			{
				cout <<"lambda["<<i<<"] = " <<lambda[i]<<" u[i] = "<<u[i] << " b[i] = " << b[i] <<endl;
			}
#endif
            zd += ( m_u[i] + m_lambda[i] ) * b[i];
			normeS += s[i]*s[i];
		}

		//Eventuelle mémorisation si meilleur score
		if ( zd > meilleurZd )
		{
			meilleurZd = zd;
            meilleursU = m_u;
            meilleursLambda = m_lambda;
			derniereItAmeliorant = nbIt;
		}

		//Maj lambda
		if ( (kWh_carre) EPSILON > normeS )
		{
			normeS = (kWh_carre) 1;
		}
        if (m_UB>(euros)0)
		{
			for ( int i=0 ; i < nbLig ; i++ )
			{
                m_lambda[i] -= ((m_UB -meilleurZd) / (normeS / s[i])) * pas;
			}
		}
		else
		{
			for ( int i=0 ; i < nbLig ; i++ )
			{
                m_lambda[i] -= ((cost[0] - meilleurZd) / (normeS/s[i])) * pas ;
//					lambda[i] -= pas * s[i] * 0.1 * zd / normeS;
			}
		}
#ifdef _DEBUG_ZD
		cout<<"zd = "<<zd<<" à l'itération "<<nbIt<<endl;
#endif
		pas *= (coeff_t) 0.99;
	}

#ifdef _DEBUG_BORNESUP
    if (m_UB >(euros)0)
        cout << "Borne Sup utilisée = "<< m_UB <<endl;
	else
		cout << "Borne Sup utilisée = cost[0] = "<<cost[0] <<endl;
	cout << "Meilleur Zd obtenu pour la passe : "<<meilleurZd<<endl;
	cout << "Nb itérations "<<nbIt_<<endl;
#endif
    m_u = meilleursU;
    m_lambda = meilleursLambda;
    return m_u;
}

vector<kWh> DualAscent::calculePoids() const
{
	vector<kWh> q(cost.size(), (kWh)0);
	for ( size_t j=0 ; j < cost.size() ; j++ )
	{
		for ( size_t i=0 ; i < b.size() ; i++ )
		{
			q[j] += a(j, i);
		}
	}
	return q;
}

vector<euros_par_kWh> DualAscent::DAcalc(const vector<euros_par_kWh> &lambda, vector<euros> &L,const vector<kWh> &q ) const
{
    vector<euros_par_kWh> u(b.size(), (euros_par_kWh)0);

	//Calcul des L[j]
	for ( size_t j=0 ; j < L.size() ; j++ )
	{
		L[j] = (euros) 0;
		for ( size_t i=0 ; i < b.size() ; i++ )
		{
            L[j] += lambda[i] * a(j, i) ;
		}
	}

	//Calcul des Ui
	for ( size_t i=0 ; i < b.size() ; i++ )
	{
        euros_par_kWh wmin = (euros_par_kWh) s_infini;
		for ( size_t j=0 ; j < L.size() ; j++ )
		{
			if ( a(j, i) > (kWh) 0 )
			{
                euros_par_kWh w = (cost[j]-L[j])/q[j];
				wmin = (wmin < w) ? wmin : w;
			}
		}
		u[i] = wmin + lambda[i];
	}
	return u;
}

vector<kWh> DualAscent::sousgrad(vector<coeff_t> &x,
                                    const vector<euros_par_kWh> &u,
                                    const vector<euros_par_kWh> &lambda,
                                    const vector<euros> &L,
                                    const vector<kWh> &q) const
{
	//Calcul de l'indexMin
	vector<int> indexMin(b.size(), -1);
	for ( size_t i=0 ; i < b.size() ; i++ )
	{
        euros_par_kWh mini = (euros_par_kWh) s_infini;
		for ( size_t j=0 ; j < cost.size() ; j++ )
		{
			if ( a(j, i) > (kWh)0 )
			{
                euros_par_kWh w = (cost[j]-L[j])/q[j];
				if ( w < mini )
				{
					mini = w;
					indexMin[i] = j;
				}
			}
		}
	}

	//Calcul de x
		//initialisation
    vector<kWh> somme_b(x.size(), (kWh) 0);
	/* Créé NULL juste avant sousgrad
	for ( size_t j=0 ; j < cost.size() ; j++ )
	{
		x[j] = (coeff_t) 0;
	}
	//*/
		//Remplissage
    for ( size_t i=0 ; i < indexMin.size() ; i++)
	{
        somme_b[indexMin[i]] += b[i];
	}
		//Division
    for ( size_t j=0 ; j < cost.size() ; j++ )
	{
        x[j] = somme_b[j] / q[j];
	}

	//Calcul de s
	vector<kWh> s(b.size(), (kWh) 0);
    for ( size_t i=0 ; i < s.size() ; i++ )
	{
        for ( size_t j=0 ; j < cost.size() ; j++ )
		{
            s[i] += a(j, i) * x[j];
		}
		s[i] -= b[i];
	}

	return s;
}

void DualAscent::vireColonneThermique()
{
    therm_A         = A.front() ;
    therm_recharge  = recharge.front();
    therm_cost      = cost.front();
    therm_reste     = reste.front();

    A.erase(A.begin());
    recharge.erase(recharge.begin());
    cost.erase(cost.begin());
    reste.erase(reste.begin());
}

void DualAscent::remetColonneThermique()
{
    A.insert        ( A.begin(),        therm_A );
    recharge.insert ( recharge.begin(), therm_recharge );
    cost.insert     ( cost.begin(),     therm_cost );
    reste.insert    ( reste.begin(),    therm_reste);
}


// Test Borne Sup

bool DualAscent::testUB () const
{
    bool out = false;
    bool STOP = false;
    size_t compteur = 0;
    kWh SumProd;
    while (!STOP)
    {
        SumProd = (kWh) 0;
        for (size_t l =1 ; l < cost.size() ; l++)
        {
         	SumProd+=a(l,compteur);
        }
#ifdef PMAX_POUR_DEM
        out = (SumProd>=b_EMAX[compteur]);
#else
        out = (SumProd>=b[compteur]);
#endif
        STOP = (!out);

        compteur++;

        if (compteur==b.size())
        {
            STOP = true;
        }
    }
#ifdef VERBOSE_SUMPROD
    if (!(out))
    {
#ifdef PMAX_POUR_DEM
    	cout << "Sumprod : "<< SumProd<< " et b_EMAX["<<(compteur-1)<<"] = "<<b_EMAX[compteur]<<endl;
#else
    	cout << "Sumprod : "<< SumProd<< " et b["<<(compteur-1)<<"] = "<<b[compteur]<<endl;
#endif
    }
#endif
    return out;
}

euros DualAscent::recherche_UB () const
{
    euros UB = (euros) 0;

    size_t n = cost.size();
    size_t m = b.size();

    vector<bool> choix (n,false);
    choix[0] = true ; // la première colonne est la colonne thermique

    vector<kWh> SumProd (m, (kWh) 0);

    for (int k =  b.size()-1;k>=0;k--)
    {
#ifdef PMAX_POUR_DEM
        if (SumProd[k]>=b_EMAX[k]) {;}
#else
        if (SumProd[k]>=b[k]) {;}
#endif
        else
        {
            bool STOP = false;
            size_t  compteur = 0 ;
            int     index_min=-1;

            euros cost_min = (euros) s_infini;

            while (!STOP)
            {
                if (!choix[compteur])
                {
                // Colonne non déjà sélectionnée
                    if (a(compteur,k) > (kWh) 0)
                    {
                    	// Colonne à production positive
                        if (cost[compteur] < cost_min)
                        {
                            index_min=compteur;
                            cost_min=cost[compteur];
                        }
                    }
                }

                compteur++;

                if (compteur==n)
                {
                    if (index_min>=0)
                    {
                        choix[index_min]=true;
                        UB += cost_min;
                        for (size_t aa=0;aa<m;aa++)
                        {
                            SumProd[aa]+=a(index_min,aa);
                        }
                        compteur = 0 , index_min=-1;
                        cost_min = (euros) s_infini;
                    }
                    else
                    {
                        break;
                    }
#ifdef PMAX_POUR_DEM
                    STOP = (SumProd[k]>=b_EMAX[k]);
#else
                    STOP = (SumProd[k]>=b[k]);
#endif
                }
            }
        }
    }
#ifdef VERBOSE_BORNE_SUP
    cout << "Calcul de borne Sup : "<<UB<<endl;
/*    for (int t=0;t<n;t++)
        cout << choix[t] << " ";cout << endl;*/
#endif
return UB;
}


euros DualAscent::recherche_UB2 (const vector<kWh> &q) const
{
    euros UB = (euros) 0;
    size_t n = cost.size();
    size_t m = b.size();
    if ( n > (unsigned) (3*I()) )
    {
	//	vector<bool> choix (n,false);
		// choix[0] = true ; // la première colonne est la colonne thermique

		vector<kWh> SumProd (m, (kWh) 0);
		vector<int> index_min(I(), -1); // Ne pas sélectionner les premières colonnes (sans production)
		vector<euros_par_kWh> ratio_min (I(), (euros_par_kWh) s_infini) ;

        for (size_t k = H(); k<m; k++)
		{
		    /// Choix des colonnes
// #warning : le +1 est-il pour la colonne thermique ? A modifier si on en prend plusieurs
		    for ( size_t i = I()+1 ; i < n ; i++)
		    {
//		    	cout << "i = " << i << endl;
	    	    if (a(i,k) > (kWh) 0) // ==1 car partie affectation
	            {
	            	int id_cent = (i-1) % I();
		            // Colonne à production positive
		            euros_par_kWh ratio = cost[i]/q[i];
	                if (ratio < ratio_min[id_cent])
	                {
	                    index_min[id_cent]=i;
	                    ratio_min[id_cent]=ratio ; //cost[i];
	                }
	            }
	        }
        }

       	/// Somme des productions nucléaires
		for (int i=0;i<I();i++)
		{
			UB += cost[index_min[i]];
			for( int w = 0;w<H();w++)
			{
		    	SumProd[w]+=a(index_min[i],w);
		    }
		}

		/// Calcul de la différence et production thermique
		for (week_t w = (week_t) 0 ; w < H() ; w = (week_t) (w+1) )
		{
			#ifdef PMAX_POUR_DEM
		    kWh diff = max((kWh) 0, b_EMAX[w]-SumProd[w]);
			#else
            kWh diff = max((kWh) 0, b[w]-SumProd[w]);
			#endif
			UB += CInstanceReduite::ThCout( w )*diff;
		}
	#ifdef VERBOSE_BORNE_SUP
		cout << "Calcul de borne Sup : "<<UB<<endl;
/*		for (int t=0;t<I();t++)
		    cout << index_min[t] << " ";*/
// 		cout << endl;
	#endif
	}
	else
	{
		UB = cost[0];
	}
	return UB;
}