#ifndef __SDP_H_
#define __SDP_H_

#include <vector>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define pi 3.141592653589793238462643383279502884197169399375105

/*------Constants for rnd_uni()--------------------------------------------*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-07
#define RNMX (1.0-EPS)

using namespace std;

int randVal, // save random integer
    randVal2;  // save random interger (used in SDP15)

int nobj;
int
nobj_l=3,
nobj_u=3;

int nvar;
int
nvar_l=10,
nvar_u=20;

int seed = 177;
long rnd_uni_init;

double   rnd_rt, // random number between 0 and 1;
Tt,		// time instant;
sdp_y[256]; // preserve PS of SDP1 every time step;

double rnd_uni(long *idum);
double minPeaksOfSDP7(double x, double pt);
double SafeCot(double x);
double SafeAcos(double x);
int changeEnvironments(string strTestInstance, int gen, int T0, int nt, int taut);

int changeEnvironments(string strTestInstance, int gen, int T0, int nt, int taut)
{
    seed = (seed + 23) % 1377;
    rnd_uni_init = -(long)seed;
    
//    if (strTestInstance=="SDP13"){ nobj_l=2; nobj_u=5;} else nobj_l=nobj_u=3;
    
    nobj = nobj_l;
    nvar = nvar_l;
    
    randVal2=0;
    
    int g = (gen - T0)>0 ? (gen - T0) : 0;
    if (g > 0)
        Tt = 1.0 / nt*(floor(1.0*(g - 1) / taut) + 1);
    else
    {
        Tt = 0.0;
        rnd_rt = 0.0;
        
        // set number of objectives and variables, change here if needed
    }
    
    if (g%taut == 1)
    {
        rnd_rt = rnd_uni(&rnd_uni_init);
        
        if ((strTestInstance== "SDP1")){
            for (int i = nobj - 1; i < nvar; i++)
            {
                rnd_rt = rnd_uni(&rnd_uni_init);
                sdp_y[i] += 5*(rnd_rt - 0.5)*sin(0.5*pi*Tt);
                if ((sdp_y[i]>1.0) || (sdp_y[i] < 0.0))
                    sdp_y[i] = rnd_rt;
            }
        }
        
        else if (strTestInstance== "SDP12"){
            nvar = nvar_l + floor(0.5+rnd_rt*(nvar_u - nvar_l));
            randVal=nvar;
            return nvar;
        }
        else if (strTestInstance=="SDP13"){
            nobj = nobj_l + floor(0.5+rnd_rt*(nobj_u - nobj_l));
            randVal=nobj;
            return nobj;
        }
        else if (strTestInstance=="SDP15"){
            randVal=1 + floor(0.5+rnd_rt*(nobj - 2));
            randVal2=1 + floor(0.5+rnd_uni(&rnd_uni_init)*(nobj - 2));
            return nobj;
        }
    }
    else
    {
        if ((strTestInstance== "SDP1"))
        {
            for (int i = nobj; i < nvar; i++)
            {
                sdp_y[i] = 1.0 * (i + 1) / nvar;
            }
        }
        
        if ((strTestInstance== "SDP13"))
        {
            nobj = 2;
        }
        
        if ((strTestInstance== "SDP15"))
        {
            randVal=nobj-1; randVal2=0;
        }

    }
    return 1;
}


double SafeAcos(double x)
{
    if (x < -1.0) x = -1.0;
    else if (x > 1.0) x = 1.0;
    return acos(x);
}

double SafeCot(double x)
{
    if (sin(x) == 0)
        return INFINITY;
    else
        return cos(x) / sin(x);
}

double minPeaksOfSDP7(double x, double pt)
{
    if (pt<1) pt=1;
    if (pt>5) pt=5;
    double min = INFINITY;
    double tmp;
    for (int i = 1; i < 6; i++)
    {
        if (i == pt)
            tmp = 0.5 + 10 * pow(10 * x - 2 * (i - 1), 2.0);
        else
            tmp = i + 10 * pow(10 * x - 2 * (i - 1), 2.0);
        if (tmp <= min)
            min = tmp;
    }
    return min;
}


// the objective vector of SDP.
void objective(const char* strTestInstance,vector<double> &x_var, vector<double> &y_obj)
{
    double G = sin(0.5*pi*Tt);
    // SDP PROBLEMS
    
    if (strcmp(strTestInstance,"SDP1")==0) //SDP1
    {
        double g = 0;
        
        // ps functions
        for (int i = nobj; i < nvar; i++)
        {
            g += pow(x_var[i] - sdp_y[i],2.0);
        }
        
        // pf functions
        double prd = 1;
        for (int j = 0; j < nobj; j++){
            double xx = 3 * x_var[j] + 1;
            prd *= xx;
        }
        for (int j = 0; j < nobj; j++)
        {
            double xx = 3 * x_var[j] + 1;
            y_obj[j] = (1 + g)*xx / pow(prd / xx, 1.0 / (nobj-1));
        }
        return;
    }
    
    if (strcmp(strTestInstance,"SDP2")==0) //SDP2
    {
        double g = 0;
        
        double x1 = 3 * x_var[0] + 1;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double xx = 2 * (x_var[i] - 0.5) - cos(Tt+2*x1); 
            g += xx*xx;
        }
        g*=sin(pi*x1/8);
        // pf functions
        double sm = 0;
        for (int j = 0; j < nobj - 1; j++)
        {
            double xx = 3 * x_var[j] + 1;
            sm += xx;
        }
        sm += 1;
        
        for (int j = 0; j < nobj - 1; j++)
        {
            double xx = 3 * x_var[j] + 1;
            y_obj[j] = (1 + g)*(sm - xx+Tt) / xx;
        }
        y_obj[nobj - 1] = (1 + g)*(sm-1)/(1+Tt);
        return;
    }
    
    if (strcmp(strTestInstance,"SDP3")==0) //SDP3
    {
        double g = 0;
        int pt = floor(5 * fabs(sin(pi*Tt)));
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = 2 * (x_var[i] - 0.5) - cos(Tt); 
            g += 4*pow(y, 2.0) -cos(2 * pt*pi*y) + 1;
        }
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj-1; j++)
        {
            y_obj[j] = (1 + g)*(1 - x_var[j] + 0.05*sin(6 * pi*x_var[j]))*pd;
            pd *= x_var[j] + 0.05*sin(6 * pi*x_var[j]);
        }
        y_obj[nobj - 1] = (1 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP4")==0) //SDP4
    {
        double w = (rnd_rt - 0.5); //----------------------
        if (w > 0)
            w = 1.0;
        else if (w < 0)
            w = -1.0;
        w = w*floor(6 * fabs(G));
        
        randVal = w; // save this random value
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = 2 * (x_var[i] - 0.5) - cos(Tt+x_var[i-1]+x_var[0]); 
            g += y*y;
        }
        
        // pf functions
        double sm = 0;
        for (int j = 0; j < nobj - 2; j++)
        {
            y_obj[j] = x_var[j];
            sm += x_var[j];
        }
        sm += x_var[nobj - 2];
        sm = sm / (nobj - 1);
        
        y_obj[nobj - 2] = (1 + g)*(sm + 0.05*sin(w*pi*sm));
        y_obj[nobj - 1] = (1 + g)*(1.0 - sm + 0.05*sin(w*pi*sm));
        return;
    }
    
    if (strcmp(strTestInstance,"SDP5")==0) //SDP5
    {
        double Gt = fabs(G);
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = x_var[i] - 0.5*Gt*x_var[0]; 
            g += y*y;
        }
        g += Gt;
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            double y = pi*Gt / 6 + (pi / 2 - pi*Gt / 3)*x_var[j];
            y_obj[j] = (1 + g)*sin(y)*pd;
            pd = cos(y)*pd;
        }
        y_obj[nobj - 1] = (1 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP6")==0) //SDP6
    {
        double at = 0.5*fabs(sin(pi*Tt)); 
        double pt = 10 * cos(2.5*pi*Tt);
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = x_var[i] - 0.5;
            g += y*y *(1 + fabs(cos(8 * pi*x_var[i])));
        }
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[nobj - 1 - j] = (1 + g)*sin(0.5*pi*x_var[j])*pd;
            pd *= cos(0.5*pi*x_var[j]);
        }
        y_obj[0] = (1 + g)*pd;
        
        
        if (x_var[0]<at)
            y_obj[nobj - 1] = (1 + g)*fabs(pt*(cos(0.5*pi*x_var[0]) - cos(0.5*pi*at)) + sin(0.5*pi*at)); 
        //		else
        //            y_obj[nobj - 1] = (1 + g)*sin(0.5*pi*x_var[j])
        return;
    }
    
    if (strcmp(strTestInstance,"SDP7")==0) //SDP7
    {
        double at = 0.5*sin(pi*Tt);
        int pt = 1+floor(5 * rnd_rt);
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            g += minPeaksOfSDP7(x_var[i], pt);
        }
        g /= (nvar - nobj + 1);
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[j] = (0.5 + g)*(1 - x_var[j])*pd;
            pd *= x_var[j];
        }
        y_obj[nobj - 1] = (0.5 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP8")==0) //SDP8
    {
        double Gt = fabs(G);
        int pt = floor(6 * Gt);
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = x_var[i] - fabs(atan(SafeCot(3 * pi*Tt*Tt))) / pi;
            g += y*y;
        }
        g += Gt;
        
        // pf functions
        double sm = 0;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[j] = (1 + g)*pow(cos(0.5*pi*x_var[j]), 2.0) + Gt;
            sm += pow(sin(0.5*pi*x_var[j]), 2.0) + sin(0.5*pi*x_var[j])*pow(cos(pt*pi*x_var[j]), 2.0);
        }
        y_obj[nobj - 1] = sm + Gt;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP9")==0) //SDP9
    {
        double kt = 10 * sin(pi*Tt);
        
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
            pd *= sin(floor(kt*(2 * x_var[j] - 1.0))*pi / 2);
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = 2 * (x_var[i] - 0.5) - sin(Tt*x_var[0]);
            g += y*y;
        }
        g += fabs(pd);
        
        // pf functions
        pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[nobj - 1 - j] = (1 + g)*sin(0.5*pi*x_var[j])*pd;
            pd *= cos(0.5*pi*x_var[j]);
        }
        y_obj[0] = (1 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP10")==0) //SDP10
    {
        int r = floor(10 * fabs(G));
        
        double sm = 0;
        for (int j = 0; j < nobj - 1; j++)
            sm += pow(x_var[j], 2.0);
        sm /= (nobj - 1);
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = 2 * (x_var[i] - 0.5) - sin(0.5*pi*Tt+x_var[0]); 
            g += y*y;
        }
        
        // pf functions
        for (int j = 0; j < nobj - 1; j++)
            y_obj[j] = x_var[j];
        
        y_obj[nobj - 1] = (1 + g)*(2 - sm - pow(sm, 0.5)*pow(-sin(2.5*pi*sm), r));
        return;
    }
    
    if (strcmp(strTestInstance,"SDP11")==0) //SDP11
    {
        double at = 3 * Tt - floor(3 * Tt);
        double bt = 3 * Tt + 0.2 - floor(3 * Tt + 0.2);
        
        //double sm = 0;
        /*for (int j = 0; j < nobj - 1; j++)
         sm += x_var[j];
         sm /= (nobj - 1);*/
        
        
        double g = 0;
        double ps=0;
        for (int i=0; i<nobj-1; i++)
            ps+=x_var[i];
        
        if (ps >= at && ps < bt){
            for (int i = nobj - 1; i < nvar; i++) {
                double p = x_var[i] - fabs(G);
                g += -0.9 * p * p + pow(fabs(p), 0.6);
            }
        }
        else{
            for (int i = nobj - 1; i < nvar; i++) {
                double p = x_var[i] - fabs(G); 
                g += p*p;
            }
        }
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            double yj = 0.5*pi*x_var[j];
            y_obj[j] = (1 + g)*sin(yj)*pd;
            pd *= cos(yj);
        }
        y_obj[nobj - 1] = (1 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP12")==0) //SDP12 -change variables
    {
        //int nvar = nvar_l + floor(rnd_rt*(nvar_u - nvar_l));
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = 2 * (x_var[i] - 0.5) - sin(Tt)*sin(2 * pi*x_var[1]);
            g += y*y;
        }
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[j] = (1.0 + g)*(1 - x_var[j])*pd;
            pd *= x_var[j];
        }
        y_obj[nobj - 1] = (1.0 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP13")==0) //SDP13-change objectives
    {
        // *** this is for many objectives only.
        //nobj = nobj_l + floor(rnd_rt*(nobj_u - nobj_l));
        
        randVal= nobj;  // save this random value.
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = x_var[i] - (i*Tt)/(nobj+i*Tt);
            g += pow(y, 2.0);
        }
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj; j++)
        {
            double yj = pi*(x_var[j] + 1) / 6;
            y_obj[j] = (1 + g)*sin(yj)*pd;
            pd *= cos(yj);
        }
        
        for (int j=nobj; j<y_obj.size(); j++)
            y_obj[j]=0;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP14")==0) //SDP14 -degenerate
    {
        double pt = 1 + floor(fabs((nobj - 2)*cos(0.5*pi*Tt)));
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            g += pow(x_var[i] - 0.5, 2.0);
        }
        
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            double yj = x_var[j];
            if (j >= pt) yj = 0.5 + x_var[j] * g*fabs(G);
            y_obj[j] = (1 + g)*(1+g - yj)*pd;
            pd *= yj;
        }
        y_obj[nobj - 1] = (1.0 + g)*pd;
        return;
    }
    
    if (strcmp(strTestInstance,"SDP15")==0) //SDP15-degenerate
    {
        int dt = randVal;
        int pt = randVal2;
        
        double g = 0;
        for (int i = nobj - 1; i < nvar; i++)
        {
            double y = x_var[i] - dt/ (nobj - 1);
            g += pow(y, 2.0);
        }
        
        vector<double> yk(nobj-1,0);
        for (int i=1; i<nobj; i++)
        {
            int k=(pt+i-1)%(nobj-1);
            if (k<=dt) yk[k]=0.5*pi*x_var[k];
            else yk[k]=SafeAcos(1.0 / (pow(2.0, 0.5)*(1.0 + x_var[k] * g*fabs(G))));
        }
        
        // pf functions
        double pd = 1;
        for (int j = 0; j < nobj - 1; j++)
        {
            y_obj[j] = pow(1 + g, j + 1)*sin(yk[j])*pd;
            pd *= cos(yk[j]);
        }
        y_obj[nobj - 1] = pow(1 + g, nobj)*pd;
        return;
    }
    
}

//the random generator in [0,1)
double rnd_uni(long *idum)
{
    long j;
    long k;
    static long idum2 = 123456789;
    static long iy = 0;
    static long iv[NTAB];
    double temp;
    
    if (*idum <= 0)
    {
        if (-(*idum) < 1) *idum = 1;
        else *idum = -(*idum);
        idum2 = (*idum);
        for (j = NTAB + 7; j >= 0; j--)
        {
            k = (*idum) / IQ1;
            *idum = IA1*(*idum - k*IQ1) - k*IR1;
            if (*idum < 0) *idum += IM1;
            if (j < NTAB) iv[j] = *idum;
        }
        iy = iv[0];
    }
    k = (*idum) / IQ1;
    *idum = IA1*(*idum - k*IQ1) - k*IR1;
    if (*idum < 0) *idum += IM1;
    k = idum2 / IQ2;
    idum2 = IA2*(idum2 - k*IQ2) - k*IR2;
    if (idum2 < 0) idum2 += IM2;
    j = iy / NDIV;
    iy = iv[j] - idum2;
    iv[j] = *idum;
    if (iy < 1) iy += IMM1;
    if ((temp = AM*iy) > RNMX) return RNMX;
    else return temp;
    
}/*------End of rnd_uni()--------------------------*/
#endif