#include "R.h"
#include <iostream>
#include <fstream>
#include <cstdlib> 
#include <ctime> 
#include <cmath>
#include <sstream>
#include <string>
using namespace std;

extern "C"{
       
double ilogit(double x)
{
       return(exp(x)/(1+exp(x)));
}
  
double runif(int size)
{
       if(size==1)      
       {
              int random_integer; 
              int lowest=1, highest=1000; 
              int range=(highest-lowest)+1; 
              random_integer = lowest+int(range*rand()/(RAND_MAX + 1.0)); 
              double random_number = double(random_integer)/double(range+1); 
              return(random_number);
       }
       
       if(size>1)
       cout << "size is greater than 1!" << endl;
}

/*
int random_int(int lower_bound,int upper_bound)
{
              srand((unsigned)time(0)+rand());
              int random_integer; 
              int lowest=lower_bound, highest=upper_bound; 
              int range=(highest-lowest)+1; 
              random_integer = lowest+int(range*rand()/(RAND_MAX + 1.0)); 
              return(random_integer); 
}
*/

int rbinom(int size,double p)
{
       int random_integer; 
       int lowest=1, highest=1000; 
       int range=(highest-lowest)+1; 
       random_integer = lowest+int(range*rand()/(RAND_MAX + 1.0)); 
       double random_number = double(random_integer)/double(range+1); 
    
       if(size==1)
       {
                  double cut_1=(1-p);
                  double cut_2=p;
                  
                  if(random_number <= cut_1)
                  return(0);
                  else return (1);
       }
       
       if(size==2)
       {
                  double cut_1=(1-p)*(1-p);  
                  double cut_2=2*p*(1-p);
                  double cut_3=p*p;     
                  
                  if(random_number <= cut_1)
                  return (0);
                  else if(random_number <=cut_1+cut_2) 
                  return (1);
                  else return (2);
       }
       
       if(size>2)
       cout << "size is greater than 2!!" << endl;
}

double f(int i,int x,double hat_p1,double hat_p2)
{
	if(i==1)
	return(pow(2.0,double(x==1))*pow(hat_p1,double((x==1)+2*(x==2)))*pow((1-hat_p1),double(2*(x==0)+(x==1))));
	if(i==2)
	return(pow(2.0,double(x==1))*pow(hat_p2,double((x==1)+2*(x==2)))*pow((1-hat_p2),double(2*(x==0)+(x==1))));
}

double mean_2(int array[],int N)
{
       double sum = 0 ;
       for (int i = 1; i <= N; i++)
       sum = sum + double(array [i-1]);
       return sum/double(N);
} 

double var(int array[], int N)
{
       double sum = 0;
       double STD_DEV = 0; 

       for (int i = 1; i <= N; i++)
       {
           sum = sum + double(array [i-1]);
       }
       double average = sum/double(N);
       for(int i = 1; i <= N; i++)
       {
               STD_DEV+=(double(array[i-1])-average)*(double(array[i-1])-average);
       }
       return(STD_DEV/double(N-1));

} 
  
double maximum(double stat[],int nSNP,int SNP0)
{
       double max=-100;
       for(int i=0;i<nSNP;i++)
       {
               if(i!=(SNP0-1))
               {
                          if(stat[i]>max)
                          max = stat[i];    
               }
       }
       return max;
}

double mean(double stat[],int nSNP,int SNP0)
{
       double m=0;
       for(int i=0;i<nSNP;i++)
       {
               if(i!=(SNP0-1))
               {
                          m = m + stat[i];    
               }
       }
       m = m/double(nSNP-1);
       return m;      
}  
  
void MLRT(int *x,int *n_case,int *n_control, int *SNP, int *disease_causing, int *Boot,
                 double *p_value_Max_MLRT, double *p_value_Mean_MLRT, double *MLRT_MAX,
                 double *MLRT_MEAN,double *MLRT_STAT,double *theta,double *theta_b,double *mixing_proportion)
{
     int n = n_case[0];
     int m = n_control[0];     
     int nSNP = SNP[0];
     int B = Boot[0];
     int SNP0 = disease_causing[0];
     
     int g = 2;
     int *data = new int[(n+m)*nSNP];     
     int *true_data = new int[(n+m)*nSNP];                                                    
     
     double * replicate_test_statistics = new double[(B+1)*nSNP];                                            
     double * replicate_test_statistics_z = new double[(B+1)*nSNP];                                          
     double * replicate_test_statistics_CA = new double[(B+1)*nSNP];                                       
   
     double hat_p1_backup;
     double hat_p2_backup;
     double mixing_proportion_backup;
   
     for(int rep=1;rep<=(B+1);rep++)
     {
                 if(rep % (B+1)==1)                                                                                      
                 {
                        for(int i=1;i<=nSNP*(n+m);i++){
                        data[i-1]=x[i-1];   // x should be the transpose of what I usually have
                        true_data[i-1]=x[i-1]; }
                 }
                 
                 if(rep % (B+1)!=1)                                                                                     
                 {
                        int N;
                        srand((unsigned)time(0)+rand());
                        int random_integer; 
                        int lowest=50, highest=n+m; 
                        int range=abs((highest-lowest))+1;
                        double temp= rand()/(RAND_MAX + 1.0);
                        random_integer = lowest+int(range*temp);
                        N = random_integer;  // total number of switches, at least lowest that much
                        for(int i=1;i<=N;i++)
                        {
                                srand((unsigned)time(0)+rand());
                                lowest=1;
                                range=(highest-lowest)+1;
                                temp= rand()/(RAND_MAX + 1.0);
                                random_integer = lowest+int(range*temp);
                                int pair_1=random_integer%(n+m);
                                temp= rand()/(RAND_MAX + 1.0);
                                random_integer = lowest+int(range*temp);
                                int pair_2=random_integer%(n+m);
                                for(int j=0;j<nSNP;j++)
                                {
                                        int temp = data[pair_1*nSNP+j];
                                        data[pair_1*nSNP+j]=data[pair_2*nSNP+j];
                                        data[pair_2*nSNP+j]=temp;        
                                }
                        }
                 }
                 
                 
                 for(int allele=0;allele<nSNP;allele++)
                 {
                          int * data_case = new int[n];                                                 
                          int * data_control = new int[m];                                                
                          int s0=0;
                          int s1=0;
                          int s2=0;
                          int r0=0;
                          int r1=0;
                          int r2=0;
     
                          for(int i=1;i<=n;i++)
                          {                                  
                            data_case[i-1]=data[(i-1)*nSNP+allele];   
                            if(data_case[i-1]==0) r0++;
                            if(data_case[i-1]==1) r1++;
                            if(data_case[i-1]==2) r2++; 
                          }

                          for(int i=n+1;i<=n+m;i++)
                          {
                            data_control[i-n-1]=data[(i-1)*nSNP+allele];
                            if(data_control[i-n-1]==0) s0++;
                            if(data_control[i-n-1]==1) s1++;
                            if(data_control[i-n-1]==2) s2++;
                          }
   
                          double hat_p2 = double(s1+2*s2)/double(2*s0+2*s1+2*s2);
                          // cout << "hat_p2 " << hat_p2 << endl;
                          
                          
                          int starting_values;                		                               
                          starting_values=20;
                          
                          double likelihood_record[starting_values];
                          double hat_mixture_record[starting_values][2];
                          double hat_p1_record[starting_values];
                          double C = 1;
                          
                          for(int rate=1;rate<=starting_values;rate++)
                          {
                                  double hat_p1=0.3;
                                  double hat_mixture [] = {double(rate)/double(3*starting_values),
                                  1-double(rate)/double(3*starting_values)};
                                  double * tao = new double[g*n]; 
                                  
                                  int iteration;                                                        
                                  iteration=60;
                                  						                     
                                  for(int z=1;z<=iteration;z++)
                                  {
                                        for(int j=1;j<=n;j++)
                                        {
	                                         for(int i=1;i<=g;i++)
                                             {
		                                         double denom = 0;
	                      	                     for(int ii=1;ii<=g;ii++)
	                                             {
                                                 denom = denom+hat_mixture[ii-1]*f(ii,data_case[j-1],hat_p1,hat_p2);
                                                 }
                             
                                                 tao[(i-1)*n+(j-1)] = hat_mixture[i-1]*f(i,data_case[j-1],hat_p1,hat_p2)/denom;
                                             }
                                        }
                                   
                                        for(int i=1;i<=g;i++)
                                        {
                                                double temp=0;
                                                for(int j=1;j<=n;j++)
                                                {
                                                        temp+=tao[(i-1)*n+(j-1)];
                                                }
                       
		                                        hat_mixture[i-1] = (temp+C)/(n+C);
                                        }
                                   
                                        double temp_d=0;
                                        double temp=0;
                                        for(int j=1;j<=n;j++)
                                        {
                                             temp+=tao[(j-1)]*((data_case[j-1]==1)+2*(data_case[j-1]==2));
                                             temp_d+=tao[(j-1)];
                                        }
                  
                                        hat_p1 = temp/(2*temp_d);
                                    } // end z loop
                                    
                                    double likelihood = 0;
                                    for(int j=1;j<=n;j++)
                                    {
	                                     double sum_temp = 0;
	                                     for(int i=1;i<=g;i++)
	                                     {
                                             sum_temp = sum_temp+hat_mixture[i-1]*f(i,data_case[j-1],hat_p1,hat_p2);
                                         }
	                                     likelihood = likelihood+log(sum_temp);
                                    }
                                    likelihood = likelihood+C*log(hat_mixture[0]);

                                    likelihood_record[rate-1] = likelihood;
                                    hat_p1_record[rate-1]=hat_p1;
                                    hat_mixture_record[rate-1][0] = hat_mixture[0];
                                    hat_mixture_record[rate-1][1] = hat_mixture[1];
                                    
                                    // if(rep % (B+1)==1){
                                    // cout << "hat_p1 is  " << hat_p1 << endl;
                                    // cout << "likelihood " << likelihood << endl;
                                    // }
                           } // the loop for rate 

                            double max=likelihood_record[0];
                            int which=0;
                            for(int i=0;i<starting_values;i++)
                            {
                                    if(likelihood_record[i]>max)
                                    {
                                    max=likelihood_record[i];
                                    which=i;
                                    }
                            }
                            
                            double likelihood_alternative = max-C*log(hat_mixture_record[which][0]);
                            double likelihood_null=0;
                            for(int j=1;j<=n;j++)
                            {
    	                          likelihood_null= likelihood_null + log(f(2,data_case[j-1],hat_p2,hat_p2));
                            }
                            
                            double test_statistics = 2*(likelihood_alternative-likelihood_null);
                            
                                                       
                            if((rep % (B+1)==1)&&(allele==SNP0-1)) 
                            {
                                 hat_p1_backup=hat_p1_record[which];     
                                 hat_p2_backup=hat_p2;  
                                 mixing_proportion_backup=hat_mixture_record[which][0];                  
                            }
                           
                            
                            double test_statistics_z = (mean_2(data_case,n)-mean_2(data_control,m))/
                                                       sqrt(var(data_case,n)/n+var(data_control,m)/m);
                            
                            double p1_CA = double(r1)/double(n);
                            double p2_CA = double(r2)/double(n);
                            double q1_CA = double(s1)/double(m);
                            double q2_CA = double(s2)/double(m);

                            double optimal_x = 0.5; // # since it is additive model
                            double CA_numerator =   pow((optimal_x*(m*r1-n*s1)+(m*r2-n*s2)),2.0)/pow(double(n+m),2.0);
                            
                            double CA_denominator = n*m/(pow(double(n+m),3.0))*
                                   ((n+m)*(pow(optimal_x,2.0)*(r1+s1)+(r2+s2))-pow((optimal_x*(r1+s1)+(r2+s2)),2.0));

                            double test_statistics_CA = CA_numerator/CA_denominator; 
                            
                            replicate_test_statistics[(rep-1)*nSNP+allele] = test_statistics;
                            replicate_test_statistics_z[(rep-1)*nSNP+allele] = test_statistics_z;
                            replicate_test_statistics_CA[(rep-1)*nSNP+allele] = test_statistics_CA;
                   } // the loop for allele 
  
                   if(rep % (B+1)==0) 
                   {
                          double * replicate_test_statistics_max=new double[B+1];            
                          double * replicate_test_statistics_z_max=new double[B+1]; 
                          double * replicate_test_statistics_CA_max=new double[B+1];
           
                          double * replicate_test_statistics_mean=new double[B+1];                      
                          double * replicate_test_statistics_z_mean=new double[B+1];   
                          double * replicate_test_statistics_CA_mean=new double[B+1];  
           
                          int * replicate_test_statistics_count = new int[nSNP]; 
                          for(int i=1;i<=nSNP;i++)
                          replicate_test_statistics_count[i-1]=0;
                          int replicate_test_statistics_max_count=0;
                          int replicate_test_statistics_mean_count=0;
                          
                          for(int print=1;print<=(1+B);print++)
                          {
                            if(print>1)
                            {
                                       for(int i=0;i<nSNP;i++)                                               
                                       {                                 
                                       if(replicate_test_statistics[(print-1)*nSNP+i]<replicate_test_statistics[i])
                                       replicate_test_statistics_count[i]++;
                                       }                                            
                            }
                                  
                            replicate_test_statistics_max[print-1]=maximum(& replicate_test_statistics[(print-1)*nSNP],nSNP,SNP0); 
                            replicate_test_statistics_mean[print-1]=mean(& replicate_test_statistics[(print-1)*nSNP],nSNP,SNP0);  
                          }
                          
                          for(int final=1;final<=B+1;final++)
                          {
                          if(replicate_test_statistics_max[final-1]<replicate_test_statistics_max[0])
                          replicate_test_statistics_max_count++;
                          }

                          for(int final=1;final<=B+1;final++)
                          {
                          if(replicate_test_statistics_mean[final-1]<replicate_test_statistics_mean[0])
                          replicate_test_statistics_mean_count++;
                          }  

                          //for(int i=1;i<=nSNP;i++)
                          //p_value_MLRT[i-1]=double(replicate_test_statistics_count)/double(B);
                          
                          double Max_MLRT = double(replicate_test_statistics_max_count)/double(B);
                          double Mean_MLRT = double(replicate_test_statistics_mean_count)/double(B);
                          * p_value_Max_MLRT=Max_MLRT;  
                          * p_value_Mean_MLRT=Mean_MLRT; 
                          * MLRT_MAX=replicate_test_statistics_max[0];
                          * MLRT_MEAN=replicate_test_statistics_mean[0];
                          
                          for(int i=0;i<nSNP;i++)
                          MLRT_STAT[i]=replicate_test_statistics[i];
                          
                          * theta=hat_p1_backup;
                          * theta_b=hat_p2_backup;
                          * mixing_proportion=mixing_proportion_backup;
                          
                   } // end if
          } // end rep
} // end the whole function
  
//void foo(int *nin, double *x)   <- this is the sample function
//{
//	int n = nin[0];
//
//	int i;
//
//	for (i=0; i<n; i++)
//		x[i] = x[i] * x[i];
//		
//		x[0] = ilogit(x[0]);
//}

// end of extern "C"
}