#include "DA.h"
#include <limits>
static double s_infini = std::numeric_limits<double>::infinity() ;

using namespace std;

// constructeur

//entree::entree()
//{};

vector<double> DualAscent::calculePoids() const
{
    // dimension du probleme
    int n,m;
    n= cost.size();
    m= b.size();

    vector<double> q (n);
    double qj;

    for (int j =0;j<n;j++)
    {
        qj = 0;
        for(int i =0;i<m;i++)
        {
            qj+=a(j,i);
        }
        q[j] = qj;
    }
    return q;
}

//methode de résolution
vector<double> DualAscent::DAcalc(const vector<double> &lambda, vector<double> &L,const vector<double> &q ) const
{
    // dimension du probleme
    int n=(int) cost.size();
    int m=(int) b.size();

    // Calcul des Lj
    for (int j =0;j<n;j++)
    {
        double Lj = 0;
//        cout << "Calcul de L(aj) pour j = " << j << endl;
        for(int i =0;i<m;i++)
        {
            Lj+=a(j,i)*lambda[i];
//            cout << "lambda_i = " << lambda[i] << endl;
//            cout << "a_ij = " << a(j,i) << endl;
//            cout << "Lj = " << Lj << endl;
        }
        L[j] = Lj;
    }

    // Calcul du Dual Ascent
    vector <double> w (m,0);
//    double min_ ;

//    cout << " Lambda(aj) " ;
//    for (int i =0;i< L.size();i++)
//    cout << L[i] << "   ";
//    cout << endl;

//    cout << " qj " ;
//    for (int i =0;i< q.size();i++)
//    cout << q[i] << "   ";
//    cout << endl;

    for (int i =0;i<m;i++)
    {
        double min_ = s_infini;
//        cout << "m = " << m << endl;
        for (int j=0;j<n;j++)
        {
            if (a(j,i)>0)
            {
                min_ = min(min_,(cost[j] - L[j])/q[j]);
//                cout << "  c"<<j<<"-L"<<j<<" = " << cost[j] - L[j] << endl;
//                cout << " q" << j << "=" << q[j] << endl;
//                cout << " wij "<< (cost[j] - L[j])/q[j]<< endl;
//                cout << " min "<< j << " = " << min_<< endl;
            }
        }
//        cout << "min sur i " << min_ << endl;
        w[i] = min_ + lambda[i];
    }

//    cout << "wi calcules : " << endl;
//    for(int i =0;i<(int) w.size();i++)
//        cout << w[i] << endl;

//    cout << endl;

    return w;
}

vector <int> DualAscent::EcritureCplex(const vector <double> &w, const char *fichier) const
{
    // Ecriture modele Cplex
    int n = cost.size();
    int m = b.size();
    // Reperage des contraintes saturees
    vector <int> Q;
    double som_;
    for (int j=0;j<n;j++)
    {
        som_= 0 ;
        for (int i =0;i<m;i++)
        som_+=a(j,i)*w[i];
        cout << "somme" << j << " = " << som_ << " =? " << cost[j] << endl;

        if (epsiqual(som_, cost[j]))//(som_==cost[j])
        Q.push_back(j);
    }

    fstream fileout;
    char path[80];

    strcpy(path,"SortieCplex/");
    strcat(path,fichier);
    fileout.open(path,ios::out);

    fileout << "# Modele lineaire pour calcul de sous-gradient" << endl<< endl;
    fileout << "Minimize" << endl;
    fileout << "   Obj : " ;
    for (int i=0;i<m-1;i++)
    { fileout << "s" << i << " + " ;}
    fileout << "s" << m-1 << endl;
    fileout << "Subject To" << endl;
    for (int k=0;k<m;k++)
    {
        fileout << "   c"<< k <<" : ";
        for (int j=0;j<(int)Q.size();j++)
        {
            fileout << a(Q[j],k) << "*X" << k << " + ";
        }
        fileout << "s" << k << " = " << b[k] << endl;
    }
    fileout << "Bounds" << endl;
    for (int k=0;k<m;k++)
    {
        fileout << "   0 <= X" << k << " <=1" << endl;
    }
    fileout << "    All other variables are >= 0." << endl;
    fileout << "End" << endl;

    fileout.close();

    return Q;
}

double DualAscent::updateLambda(const vector<double> &S, vector<double> & lambda,const vector<double> &w, double pas) const
{
    int m=(int) b.size();
    double denominateur =0;
    double Zd=0;


    for (int i=0;i<m;i++)
    {
        denominateur+=S[i]*S[i];
        Zd+=b[i]*w[i];
    }
//    cout << "Zd " << Zd << endl;
    for (int i=0;i<m;i++)
    {
        lambda[i]+= -pas*0.5*Zd*(S[i])/denominateur;
    }
    return Zd;
}


vector<double> DualAscent::sousgrad(vector<double> &Xtest,
                                    const vector<double> &w,
                                    const vector<double> &lambda,
                                    const vector<double> &L,
                                    const vector<double> &q) const
{
    int m = b.size();
    int n = cost.size();
    double min_, minmem;

    vector<double> s (m);
    //vector<vector <int> > index;
    vector <int> index (m);
    vector<double> X (n);
//    vector<int> "jmin";
     int jmin;

    for (int i =0;i<m;i++)
    {
        minmem = s_infini;

        for (int j=0;j<n;j++)
        {
            if (a(j,i)>0)
            {
                min_ = min(minmem,(cost[j] - L[j])/q[j]);
                if (minmem>min_)
                {
                    jmin=j;
//                    jmin.clear();
//                    jmin.push_back(min_);
                    minmem=min_;
                }
//                else if (minmem==min_)
//                {
//                    cout << minmem << endl;
//                    jmin.push_back(min_);
//                }
            }
        }
        index[i]=jmin;// probleme si plusieurs colonnes a même cout
//        index.push_back(jmin);
    }

//    for (int c = 0;c < index.size();c++)
//    { cout << index[c].size() << " ";}cout << endl;

    vector<double> B (n,0), Btest (n,0);
    for (int i = 0;i<m;i++)
    {
//        for (int k=0;k<index[i].size();k++)
//        {
//            B[index[i][k]]+=b[i];
            B[index[i]]+=b[i];
            //Hyp : la demande est toujours superieure a 1
            if (b[i]==1)
            {
//            cout << "index[i] : " << index[i] << endl;
//                Btest[index[i][k]]+=b[i];
                Btest[index[i]]+=b[i];
            }
//        }
    }
//    for (int i =0;i<index.size();i++)
//    cout << index[i] << " " ;
//    cout << endl;
//    vector<double> X ;

    for (int j = 0;j<n;j++)
    {
//        double Qt =0 ;
//        for (int k = 0;k<index[j].size();k++)
//        { Qt += q[index[j][k]];}
        X[j]=(B[j]/q[j]);//Qt);
        Xtest[j] = Btest[j];
    }
    double s_;
    for(int i =0;i<m;i++)
    {
        s_=0;
        for (int j =0;j<n;j++)
        {
            s_+=X[j]*a(j,i);
        }
        s[i]= (s_ - b[i]);
    }

    Xtest = X;
    return s;
}

void DualAscent::boucle(vector<double> &lambda, bool BorneSup, const char* titre) const
{
    int iter = 0,lastchange;
    vector<double> u;

    bool impression = false;
    int Niter = 50000;
    char fichier[20], test[20];
    double pas = 1;
    double Zd, Zdmem;
    vector<double> q, s, Xmem, lambdamem, umem;
    vector<double>Zdmax (Niter);

    vector<double> L (cost.size());
    vector<double> X (cost.size());
    q = calculePoids();

//------------Recherche d'une upper bound ----------

    double UB = 0;
    if (!BorneSup)
    { BorneSup = testUB();}

    if (BorneSup)
    {
        UB = recherche_UB();
        cout << "upper bound = " << UB << endl;
    }

    while (iter < Niter)
    {
        if (impression) {
        cout << "Iteration   " << iter << endl;}

        u = DAcalc(lambda,L,q);
        Zd=0;

        //////////////////// Solution 2 ////

        s = sousgrad(X,u,lambda, L, q);

        double NormS = 0, denom;
        for (size_t i =0;i<b.size();i++)
        {
            Zd += u[i]*b[i];
            NormS += s[i]*s[i];

        if (impression) {
            cout << "S" << i << " = " << s[i] << endl;
            }
        }

        if (NormS>0)
        denom=NormS;
        else denom=1;

        for (size_t i =0;i<b.size();i++)
        {
            if (UB==0) { lambda[i] += -pas*0.1*Zd*s[i]/denom;}
            else { lambda[i] += -pas*(UB-Zd)*s[i]/denom;}//  UB : Upper bound
        }

        Zdmem=max(Zdmem,Zd);
        if (Zd==Zdmem)
        {
//            cout << "changement à iter " << iter << " ; Zd = " <<Zd<<endl;
            Xmem.clear();
            lastchange=iter;
            Xmem = X;
            lambdamem = lambda;
            umem = u;
//            cout << " Xmem : ";
//            for(int i =0;i<Xmem.size();i++)
//                cout << Xmem[i] << " " ;
//            cout << endl;
//            cout << " lamba : ";
//            for(int i =0;i<lambdamem.size();i++)
//                cout << lambdamem[i] << " " ;
//            cout << endl;

        }

        if (impression) {
        cout << "Zd "<< Zd << endl;
        }

        Zdmax[iter]=(Zd);

        if (impression) {

         cout << "lambda" << iter+1 << ": ";
            for(int i =0;i<lambda.size();i++)
            cout << lambda[i] << " ";
            cout << endl << "------------------" << endl ;
        }
            if (iter - lastchange > 20)
                {pas=pas*0.99;
                lastchange = iter;}
        iter+=1;


        /////////////////// solution 1 ////////
/*

        strcpy(fichier,titre);
        int taille = sprintf(test,"%d",iter);
        strcat(fichier,test);

        Q = EcritureCplex(E,w,fichier);

        if (Q.size()==1)
        {
            cout << "colonne dans Q : " << Q[0] << endl;
            ligne = S[Q[0]];

            for (int i = 0;i < ligne.size();i++)
            {cout << ligne[i] << " " ;}
            cout << endl;

            Zd = updateLambda(E,ligne,lambda,w,pas);
            Zdmem=max(Zdmem,Zd);
            Zdmax.push_back(Zd);

            cout << "lambda" << iter+1 << ": ";
            for(int i =0;i<lambda.size();i++)
            cout << lambda[i] << " ";
            cout << endl << "------------------" << endl ;

            iter+=1;
            pas=pas*0.9;
        }
        else
        {
            cout << "#####################################################" << endl;
            cout << " Taille de Q : " << Q.size() << endl;
            cout << " iteration "  << iter << endl;
            cout << " Zd mem " << Zdmem << endl;
            cout << " Quelle colonne choisir ? " << endl;
            int a;
            cin >> a ;
            ligne = S[a];

//            exit(0);

            Zd = updateLambda(E,ligne,lambda,w,1);
            Zdmem = max(Zdmem, Zd);

            Zdmax.push_back(Zd);

            cout << "lambda" << iter+1 << ": ";
            for(int i =0;i<lambda.size();i++)
            cout << lambda[i] << " ";
            cout << endl << "------------------" << endl ;

            iter+=1;
            pas=pas*0.95;

        }
*/
    }
    cout << Zdmem << endl;
    cout << endl <<"---------------------------------------------" << endl;

   // for (int i = 0;i < lambdamem.size();i++)
   // cout << lambdamem[i] << " ";
    cout << endl <<"---------------------------------------------" << endl;
    cout << "les UI " << endl;
    for (int i = 0;i < umem.size();i++)
    cout <<  umem[i] << " ";
    cout << endl <<"---------------------------------------------" << endl;

  //  for (int i = 0;i < Xmem.size();i++)
  //  cout << Xmem[i] << " ";
}

bool DualAscent::epsiqual(const double &a, const double &b) const
{
    double epsilon = 0.0001 ;
    return (max(a,b)-min(a,b) <= epsilon);
}

/// Recherche d'une borne superieure
double DualAscent::recherche_UB () const
{
    double UB = 0;

    int n = cost.size();
    int m = b.size();
    vector<bool> choix (n,false);

    vector<double> SumProd (m,0);

    for (int k =  b.size()-1;k>=0;k--)
    {
        if (SumProd[k]>b[k]) {;}//cout << "kok : " << k << endl;}
        else
        {
            bool STOP = false;
            int compteur = 0 , index_min=-1;;
            double cost_min = 10000000;
            while (!STOP)
            {
                if (!choix[compteur])
                {   // cout << compteur<<endl;
                    if (a(compteur,k) > 0)
                    {
                        if (cost[compteur] < cost_min)
                        {
                            index_min=compteur;
                            cost_min=cost[compteur];
//                            cout <<  "index_min_temp " << compteur << " de cout " << cost_min<< endl;
                        }
                    }
                }
                compteur++;
               // cout << "---------------"<<compteur<<"-------------------------"<<endl;

                if (compteur==n)
                {
                    if (index_min>=0)
                    {
//                        cout << "Pour ligne " << k << " ajout colonne " << index_min << endl;
                        choix[index_min]=true;
                        UB += cost_min; //cout << "UB : " << UB << endl;
                        for (int aa=0;aa<m;aa++)
                        {
                            SumProd[aa]+=a(index_min,aa);
                        }
                        // cout << "SumProdok" << endl;
                        compteur = 0 , index_min=-1;
                        cost_min = 10000000;
                    }
                    else
                    {
                        break;
                    }
                    STOP = (SumProd[k]>=b[k]);
//                    cout << "stop " <<  STOP << endl;
//                    cout << "compteur " << compteur << endl;
                }
            }
        }
    }
    for (int t=0;t<n;t++)
    cout << choix[t] << " ";cout << endl;
return UB;
}

bool DualAscent::testUB () const
{
    bool out = false;
    bool STOP = false;
    int compteur = 0;
    while (!STOP)
    {
//        cout << "compteur " << compteur << endl;
        double SumProd = 0;
        for (int l =0;l<A.size();l++)
        SumProd+=a(l,compteur);
//        cout << SumProd << endl;
        out = (SumProd>b[compteur]);
        STOP = (!out);
        compteur++;
        if (compteur==b.size())
        {
            STOP = true;
        }
    }
//    cout << out << endl;
    return out;
}