#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <vector>
#include <fstream>
#include <algorithm>
#include <sstream>
#include <string>
using namespace std;
using std::string;

#define ITMAX 10000
#define EPS 1.0e-10
#define EPSn 1.0e-8
/*************************nn is sample size; 
*will be MODIFIED ***/
const int nn = 100;
/***** pp is number of biomarkers; 
*will be MODIFIED ***/
#define pp 10 

#define GAMM 10.0    //gamma in sigmoid 

#define R 0.61803399  //golden section
#define C (1.0-R)     //the other part
#define SHFT2(a,b,c) (a)=(b);(b)=(c);
#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);


ofstream testout("y_and_xdata_golden.txt");
ofstream testout2("estimated_alpha.txt");
ofstream testout4("result.txt");

struct readdata {
	double **data;
};

struct paradata {
    double lambda;
    double gamma;
    double **myxydata;
    int mn;
    int p;
};


//Sigmoid function
double sigmoid(double az, double gamma){    
    return 1.0/(1.0+exp(-gamma*az));
}

//Derivative of Sigmoid function
double derisigm(double az, double gamma){
    double x;
    x=exp(-gamma*az);
    return gamma*x/((1+x)*(1+x));
}

//Regular AUC with regularization 
double reguauc(double *alpha, struct paradata mydata) {
    double xx, az, pnty, lambda0, gamma0, **zz;
    int i,j,p0, mn0;
    p0 = mydata.p;
    mn0 = mydata.mn;
    lambda0 = mydata.lambda;
    gamma0 = mydata.gamma;
    zz= mydata.myxydata;
    
    xx=pnty=0.0;
    for(i=0;i<mn0;i++) {
        az=0.0;
        for(j=0;j<p0;j++) 
				az += (zz[i][j]) * (alpha[j]);
				xx += sigmoid(az, gamma0);
    }
    for(j=0;j<p0;j++) pnty += (alpha[j])*(alpha[j]);  

    xx = 1.0 - xx/((double) mn0) + pnty * lambda0;

    return xx;
}


//Regular AUC without regularization
double reguaucwol(double *alpha, int p0, int mn, int gamma0,double **zz) {
    double xx, az;
    int i,j;
	int mn0=mn;
    
    xx=0.0;
    for(i=0;i<mn0;i++) {
        az=0.0;
        for(j=0;j<p0;j++) 
				az += (zz[i][j]) * (alpha[j]);
				xx += sigmoid(az, gamma0);
    }  

    xx = 1.0 - xx/((double) mn0);

    return xx;
}


//Derivative function 
void drauc(double *alpha, double *gdv, struct paradata mydata){
    double xx, lambda0, gamma0, *az, **zz;
    int i,j, pp0, mn0;
    
    pp0 = mydata.p;
    mn0 = mydata.mn;
    lambda0 = mydata.lambda;
    gamma0 = mydata.gamma;
    zz= mydata.myxydata;
    
    az=(double *)malloc(mn0*sizeof(double));

    for(i=0;i<mn0;i++) {
        xx=0.0;
        for(j=0;j<pp0;j++) xx += (zz[i][j]) * (alpha[j]);
		az[i]=xx;
    }
    
    for(j=0;j<pp0;j++) {
        xx=0.0;
        for(i=0;i<mn0;i++) xx += ( derisigm(az[i], gamma0) * (zz[i][j]));
        gdv[j] = - xx/((double) mn0) + 2*lambda0*alpha[j];
    }
    
    free(az);
}


// Empirical AUC function
double empauc(double *alpha, int p0, int mn, double **zz) {
    double xx, az;
    int i,j;
    xx=0.0;
    for(i=0;i<mn;i++) {
        az=0.0;
        for(j=0;j<p0;j++) az += (zz[i][j]) * (alpha[j]);
        if(az>0.0) xx +=1.0;
        else if (az==0.0) xx +=0.5;
    }    
    xx = xx/((double) mn) ;
    
    return xx;
}

//TGDR-AUC algorithm
//Given a starting point alpha[1..p], Threshold Generalized Gradient Descent minimization is 
// performed on a function func, using its gradient as calculated by a function dfunc. 
// The threshold is set as tau.
// The infinitesimal steps is input as dnu.  Returned quantities are the aug min alpha, 
// and gret, the minimum value of the function. 
void tggdmn(double *alpha, double *gret, double tau, double dnu, 
            struct paradata mydata, double(*func)(double *, struct paradata), 
            void (*dfunc)(double *, double *, struct paradata))
{
    double gg,fp, xmax;
    double *gdn, *devr;
    
    int its, j, pp0;
    
    pp0 = mydata.p;
    gdn = (double *)malloc(pp0*sizeof(double)); 
    devr = (double *)malloc(pp0*sizeof(double));

    fp=(*func)(alpha, mydata);
    for (its=0;its<ITMAX;its++) {
        (*dfunc)(alpha, devr, mydata);

        xmax=fabs(devr[0]);
        gg=0.0;
        for (j=0;j<pp0;j++) {
			gg += (devr[j])*(devr[j]);
            if (fabs(devr[j])>xmax) xmax=fabs(devr[j]);
        }
        if (gg < EPSn) break;

        for (j=0;j<pp0;j++) {
            if(fabs(devr[j]) < tau*xmax) gdn[j]=0.0;
            else gdn[j]=devr[j];
            alpha[j] -= dnu * (gdn[j]);
            if(alpha[j] < 0.0) alpha[j]=0.0;   // positive constraints.
        }
        
        (*gret)=(*func)(alpha, mydata);
        if (fabs( (*gret)-fp ) < EPS) break;
        fp=(*gret);
        
        printf("fget=%f fp=%f gg=%f \n", *gret, fp, gg);
    }
    if (its>=ITMAX) printf("\nUsed all ietration.\n");

    free(gdn);
    free(devr);
}




//main function and input data files
int main(int argc, char* argv[])
{

    int i = nn;
	string line;

	/***********Input function name is simup10n50.txt; need to be MODIFIED*/
	ifstream in("simup10n50.txt");

		double **allsp;
		allsp = (double **) malloc(nn*sizeof(double *));
		for(i=0;i<nn;i++)
			allsp[i]=new double[(pp+1)];

		int num_casesp=0;
		int num_controlsp=0;
		if (in.is_open()) {
		
		int j;
		for (i=0; i<nn; i++) {
			for (j=0; j<(pp+1); j++) {
				if (!in.eof()) {
					in>>allsp[i][j];
				}
			}
			if (allsp[i][0]==1.0)
				num_casesp++;
			else
				num_controlsp++;
		
		}
	}
	
		/***********Optimal Tau; will be MODIFIED*/
		double taus = 1;
		/***********Optimal Lambda; will be MODIFIED*/
		double lamda=0.000385795;
	
		double **ydata, **xdata, **zdata;
		double *alpha, *dgval, gval, gval2,xval, dnv;

		int  j,k, m,n,mn; 
		

	 struct paradata alldata;



			ydata = (double **) malloc(num_casesp*sizeof(double *));
			for (i=0; i<num_casesp; i++) 
				ydata[i] = new double[pp];
			xdata = (double **) malloc(num_controlsp*sizeof(double *));
			for (i=0; i<num_controlsp; i++) 
				xdata[i] = new double[pp];

			vector<int> yIndex, xIndex;
			for(i = 0; i < num_casesp; i++){
				yIndex.push_back(i);
			}
			for(i = 0; i < num_controlsp; i++){
				xIndex.push_back(i);
			}


			m=n=0;
			for(i=0;i<nn;i++) {
					if(allsp[i][0]==-1)						
						n++;
					
					else
						m++;
					
			}

						for(i=0;i<num_controlsp;i++){
				for(j=0;j<pp;j++){
					xdata[i][j]=allsp[xIndex[i]][j+1];
				}
			}
			
			for(i=0;i<num_casesp;i++){
				for(j=0;j<pp;j++){
					ydata[i][j]=allsp[(num_controlsp+yIndex[i])][j+1];
				}
			}

			
			mn=m*n;
			
			
	
				for (i=0; i<m; i++) {
					for (int j=0; j<pp; j++)
						testout<<ydata[i][j]<<" ";
					testout<<endl;
				}
				testout<<endl<<endl<<endl;
				for (i=0; i<n; i++) {
					for (int j=0; j<pp; j++)
						testout<<xdata[i][j]<<" ";
						testout<<endl;
				}
			testout<<endl<<endl<<endl;
			


			zdata = (double **) malloc(m*n*sizeof(double *));
			for (i=0; i<m*n; i++) 
				{
				zdata[i]= new double[pp];
			}	

			mn=0;
			for(i=0;i<m;i++){
				 for(j=0;j<n;j++) {
						xval=0.0;
						for(k=0;k<pp;k++) {
							zdata[mn][k]=ydata[i][k]-xdata[j][k];
							xval+=(zdata[mn][k])*(zdata[mn][k]);
						}
					if(xval!=0) for(k=0;k<pp;k++) zdata[mn][k] /= sqrt(xval);
					mn++;
				}
			}

			alpha = (double *) malloc(pp*sizeof(double));

			dgval = (double *) malloc(pp*sizeof(double));
		
			for(k=0;k<pp;k++) {
				alpha[k]=0.01*k/pp;
				dgval[k]=0.0;
			}

	
			dnv=0.01;
			gval=0.0;
			gval2=0.0;

	
			alldata.gamma = GAMM;    // value of gamma  
	
			alldata.lambda = lamda;   // value of lambda
	
			alldata.myxydata = zdata;   // value of z=y-x 
		
			alldata.mn = mn;   // value of mn  
	
			alldata.p = pp;   // value of p    

			tggdmn(alpha, &gval, taus, dnv, alldata, reguauc, drauc);
    
			for(i=0;i<pp;i++) testout2 <<alpha[i]<<" ";
			testout2<<endl;

			testout4 <<"tau= "<<taus<<endl;
			testout4<<"\n\n penalized AUC ="<< gval;
			testout4<<endl;
    
			gval=empauc(alpha, pp, mn, zdata);
			testout4<< gval<<endl;
					
			gval2=reguaucwol(alpha, pp, mn,GAMM,zdata);
			testout4<< gval2<<endl;

			testout4<<endl;

		for(i = 0; i < num_casesp; i++)
			delete ydata[i];
		free(ydata);
		for(i = 0; i < num_controlsp; i++)
			delete xdata[i];
		free(xdata);
		for(i = 0;i < m*n;i++)
			delete zdata[i];
		free(zdata);

		free(alpha);
		free(dgval);
		for(i = 0; i < nn; i++)
			delete allsp[i];		
		free(allsp);


		testout.close();
		testout2.close();
		testout4.close();
		
	 
    return 0;
}