#include "DA.h"
#include <algorithm>
#ifndef EPSILON
#define EPSILON 0.001
#endif
#include <cmath>
static double s_infini = INFINI ;
using namespace std;

void DualAscent::ajoutcolonne( vector<double> * colonne, double cou, std::vector<CPlan::CArc>& arcs, double resteCol )
{
    A.push_back(colonne);
    cost.push_back(cou);
    vector<double> 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)
{
    Clique::ajoutcolonne(cplan, A.size());
    ajoutcolonne(&cplan->m_colonne, cplan->m_cout, cplan->m_arcs, cplan->m_reste);
    m_plans.push_back(cplan);
}

DualAscent::~DualAscent()
{
    std::for_each(m_plans.begin()   , m_plans.end() , Delete());
    // les colonnes font partie des objets CPlan. Pas besoin de détruire les vecteurs
    // de A.
}

void DualAscent::boucle(vector<double>& u, vector <double> & lambda, double& UB, const char*x)  const
{
//	lambda = vector<double>(lambda.size(),0);
	int nbCol = cost.size();
	int nbLig = b.size();
	vector<double> q = calculePoids();
	vector<double> L(nbCol, 0);
	double pas = 1;
	double meilleurZd = - s_infini;
	vector<double> meilleursU(u);
	vector<double> meilleursLambda(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 (UB > 0)
		UB = min(UB,recherche_UB2(q));	
	else
		UB = recherche_UB2(q);
#endif
	int nbItMax = 3000;
	int nbItMaxSansAmelioration = 100;
	int derniereItAmeliorant = 0;
	for ( int nbIt=0 ; (nbIt < nbItMax) && (derniereItAmeliorant > nbIt-nbItMaxSansAmelioration) ; nbIt++ )
	{
		u = DAcalc(lambda, L, q);
		vector<double> x(nbCol, 0);
		vector<double> s = sousgrad(x, u, lambda, L, q);
		
		//Calcul de Zd
		double zd = 0;
		double normeS = 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 += ( u[i] + lambda[i] ) * b[i];
			normeS += s[i]*s[i];
		}
		
		//Eventuelle mémorisation si meilleur score
		if ( zd > meilleurZd )
		{
			meilleurZd = zd;
			meilleursU = u;
			meilleursLambda = lambda;
			derniereItAmeliorant = nbIt;
		}
		
		//Maj lambda
		if ( normeS > EPSILON )
		{
			if (UB>0)
			{
				for ( int i=0 ; i < nbLig ; i++ )
				{
					lambda[i] -= pas * s[i] * (UB -meilleurZd) / normeS;
				}
			}
			else
			{
				for ( int i=0 ; i < nbLig ; i++ )
				{
 					lambda[i] -= pas * s[i] * (cost[0] -meilleurZd) / normeS;
//					lambda[i] -= pas * s[i] * 0.1 * zd / normeS;
				}
			}
		}
#ifdef _DEBUG_ZD
		cout<<"zd = "<<zd<<" à l'itération "<<nbIt<<endl;
#endif
		pas *= 0.99;
	}
	
#ifdef _DEBUG_BORNESUP
	if (UB >0)
		cout << "Borne Sup utilisée = "<< UB <<endl;
	else
		cout << "Borne Sup utilisée = cost[0] = "<<cost[0] <<endl;
	cout << "Meilleur Zd obtenu pour la passe : "<<meilleurZd<<endl;	
#endif
	u = meilleursU;
	lambda = meilleursLambda;
}

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

vector<double> DualAscent::DAcalc(const vector<double> &lambda, vector<double> &L,const vector<double> &q ) const
{
	vector<double> u(b.size(), 0);
	
	//Calcul des L[j]
	for ( int j=0 ; j < L.size() ; j++ )
	{
		L[j] = 0;
		for ( int i=0 ; i < b.size() ; i++ )
		{
			L[j] += a(j, i) * lambda[i];
		}
	}
	
	//Calcul des Ui
	for ( int i=0 ; i < b.size() ; i++ )
	{
		double wmin = s_infini;
		for ( int j=0 ; j < L.size() ; j++ )
		{
			if ( a(j, i) > 0 )
			{
				double w = (cost[j]-L[j])/q[j];
				wmin = (wmin < w) ? wmin : w;
			}
		}
		u[i] = wmin + lambda[i];
	}
	return u;
}

vector<double> DualAscent::sousgrad(vector<double> &x,
                                    const vector<double> &u,
                                    const vector<double> &lambda,
                                    const vector<double> &L,
                                    const vector<double> &q) const
{
	//Calcul de l'indexMin
	vector<int> indexMin(b.size(), -1);
	for ( int i=0 ; i < b.size() ; i++ )
	{
		double mini = s_infini;
		for ( int j=0 ; j < cost.size() ; j++ )
		{
			if ( a(j, i) > 0 )
			{
				double w = (cost[j]-L[j])/q[j];
				if ( mini > w )
				{
					mini = w;
					indexMin[i] = j;
				}
			}
		}
	}
	
	//Calcul de x
		//initialisation
	for ( int j=0 ; j < cost.size() ; j++ )
	{
		x[j] = 0;
	}
		//Remplissage
	for ( int i=0 ; i < indexMin.size() ; i++)
	{
		x[indexMin[i]] += b[i];
	}
		//Division
	for ( int j=0 ; j < cost.size() ; j++ )
	{
		x[j] /= q[j];
	}
	
	//Calcul de s
	vector<double> s(b.size(), 0);
	for ( int i=0 ; i < s.size() ; i++ )
	{
		for ( int j=0 ; j < cost.size() ; j++ )
		{
			s[i] += x[j]*a(j, i);
		}
		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;
    int compteur = 0;
    double SumProd;
    while (!STOP)
    {
        SumProd = 0;
        for (int l =1 ; l < cost.size() ; l++)
        {
         	SumProd+=a(l,compteur);
        }
#ifdef PMAX_POUR_DEM
        out = (SumProd>=b_PMAX[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_PMAX["<<(compteur-1)<<"] = "<<b_PMAX[compteur]<<endl;
#else
    	cout << "Sumprod : "<< SumProd<< " et b["<<(compteur-1)<<"] = "<<b[compteur]<<endl;
#endif
    }
#endif
    return out;
}

double DualAscent::recherche_UB () const
{
    double UB = 0;
    
    int n = cost.size();
    int m = b.size();
    
    vector<bool> choix (n,false);
    choix[0] = true ; // la première colonne est la colonne thermique

    vector<double> SumProd (m,0);

    for (int k =  b.size()-1;k>=0;k--)
    {
#ifdef PMAX_POUR_DEM
        if (SumProd[k]>=b_PMAX[k]) {;}
#else
        if (SumProd[k]>=b[k]) {;}
#endif
        else
        {
            bool STOP = false;
            int compteur = 0 , index_min=-1;
 
            double cost_min = s_infini;
 
            while (!STOP)
            {
                if (!choix[compteur])
                {   
                // Colonne non déjà sélectionnée
                    if (a(compteur,k) > 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 (int aa=0;aa<m;aa++)
                        {
                            SumProd[aa]+=a(index_min,aa);
                        }
                        compteur = 0 , index_min=-1;
                        cost_min = s_infini;
                    }
                    else
                    {
                        break;
                    }
#ifdef PMAX_POUR_DEM
                    STOP = (SumProd[k]>=b_PMAX[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;
}


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

		vector<double> SumProd (m,0);
		vector<int> index_min= vector<int> (I(), -1); // Ne pas sélectionner les premières colonnes (sans production)
		vector<double> ratio_min = vector<double> (I(),s_infini) ;
		    
		for (int k = m-1;k>=H();k--)
		{
		    /// Choix des colonnes
		    for ( int i = I()+1 ; i < n ; i++)
		    {
//		    	cout << "i = " << i << endl;
	    	    if (a(i,k) > 0) // ==1 car partie affectation
	            {
	            	int id_cent = (i-1) % I();
		            // Colonne à production positive
		            double 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
		    double diff = max(0.,b_PMAX[w]-SumProd[w]);
			#else
	        double diff = max(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;
}