//
//  mainMPI.cpp
//  Tests Pohlig-Hellman pseuodrandom number generator randomph()
//
//  Created by Paul Beale on 12/15/14.
//  Copyright (c) 2014 Paul Beale. All rights reserved.
//
// Performs multiple processor statistical tests of randomph() using MPI 
//

#include "randomph.h"
#include <iostream>
#include <iomanip>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <mpi.h>
#define MPI

using namespace std;


// Return the one-sided p-value from a chi square test for a n-bin histogram
// from a uniform distribution with a large number of degreesoffreedom=n-1. 
// In the large n limit, the distribution approaches a gaussian with 
// average = degreesoffreedom and variance 2*degreesoffreedom.
// Uses complementary error function erfc() to determine p-value.
double ChiSquareLargeNPValue(int* histogram,int n)
{
    double degreesoffreedom=n-1;
    double average=0.0,chisqr=0.0,pvalue,sigma;
    for (int i=0;i<n;i++) average += histogram[i];
    average /= n;
    for (int i=0;i<n;i++) chisqr += pow(histogram[i]-average,2)/average;
    sigma = (chisqr-degreesoffreedom)/sqrt(2.*degreesoffreedom);
    pvalue=0.5*erfc(fabs(sigma)/sqrt(2.0));
    return pvalue;
}



int main(int argc, char ** argv)
{
    // number of processes (numprocesses) and process number (processid) for each process
    // 0<= processid < numprocesses
    int numprocesses=1;
    int processid=0;
    
#ifdef MPI
    MPI_Status status;
    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &processid);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocesses);
#endif


    long ntests = 20000000; //number of pseudorandom numbers used in
    //   3D, 4D, and 5D chi square tests, twice as many are used in 2D test, 
    //   and five times as are used many used in 1D test
    long nrounds=10; //number of rounds of tests

    double time=0.0,timelast=0.0;
    time =((double) clock())/((double) CLOCKS_PER_SEC); timelast=time;

    if(processid==0)
    {
        cout << "Pohlig-Hellman pseudorandom number generator based on "<<endl;
        cout << "'A new class of scalable parallel pseudorandom number generators"<<endl;
        cout<< "based on Pohlig-Hellman exponentiation ciphers,'"<<endl;
        cout << "Paul D. Beale, (2014)."<<endl;
        cout << "http://arxiv.org/abs/1411.2484"<<endl;
        cout << endl;


        cout << "Intrastream tests of sequences of "<< 5*ntests*nrounds <<" pseudorandom numbers per process"<<endl;
        cout << "using 1D, 2D, 3D, 4D, and 5D histogram tests."<<endl;
        cout << "Number of processes = "<<numprocesses<<endl;
        cout << "With these parameters, the run time for each process should be roughly " << 5*ntests*nrounds/8e6/60. << " minutes."<<endl<<endl;

    }
    
    //histograms for 1D, 2D, 3D, 4D, and 5D correlations tests
    int * hist1d = new int[1048576];
    int * hist2d = new int[1048576];
    int * hist3d = new int[1000000];
    int * hist4d = new int[1048576];
    int * hist5d = new int[1048576];
    
    //countpvalues[0]= will be total number of number of pvalue tests in one process with
    //countpvalues[1]=number with pvalue < 1e-1
    //countpvalues[2]=number with pvalue < 1e-2
    //countpvalues[3]=number with pvalue < 1e-3
    //countpvalues[4]=number with pvalue < 1e-4
    //countpvalues[5]=number with pvalue < 1e-5
    //countpvalues[6]=number with pvalue < 1e-6
    //countpvalues[7]=number with pvalue < 1e-7
    //countpvalues[8]=number with pvalue < 1e-8
    long countpvalues[9]={0};
    
    // array to hold state of the generator
    unsigned long integerparameters[7];
    unsigned long m0, s0, c0, n0, e0, p0, a0; // initial values of state parameters
    
    // p-value and base 10 log from chi square tests
    double pvalue,log10pvalue;
    
    // first call to pseudorandom number generator initializes the state
    // determine and store the initial state in each process
    randomph();
    getparameters_randomph(integerparameters);
    cout << "initial state of the generator in process " <<processid <<": m, s, c, n, e, p, a" << endl<<endl;
    for (int j=0;j<7;j++) cout << integerparameters[j] << " "; cout << endl;
    m0=integerparameters[0];
    s0=integerparameters[1];
    c0=integerparameters[2];
    n0=integerparameters[3];
    e0=integerparameters[4];
    p0=integerparameters[5];
    a0=integerparameters[6];

    for (int irounds=0;irounds<nrounds;irounds++)
    {
    
    //initialize the histograms
    for (int i=0;i<1048576;i++) hist1d[i]=hist2d[i]=hist4d[i]=hist5d[i]=0;
    for (int i=0;i<1000000;i++) hist3d[i]=0;
    
    //five pseudorandom numbers per pass 
    double r1,r2,r3,r4,r5;
    
    for (long itest=0;itest<ntests;itest++)
    {
        r1=randomph();
        r2=randomph();
        r3=randomph();
        r4=randomph();
        r5=randomph();
        
        // increment 1D, 2D, 3D, 4D, and 5D histograms
        hist1d[(int)(1048576.0*r1)]++;
        hist1d[(int)(1048576.0*r2)]++;
        hist1d[(int)(1048576.0*r3)]++;
        hist1d[(int)(1048576.0*r4)]++;
        hist1d[(int)(1048576.0*r5)]++;
        hist2d[1024*(int) (1024.0*r2) +(int)(1024*r1)]++;
        hist2d[1024*(int) (1024.0*r4) +(int)(1024*r3)]++;
        hist3d[10000*(int) (100.0*r3)+100*(int) (100.0*r2) +(int)(100.0*r1)]++;
        hist4d[32768*(int) (32.0*r4)+1024*(int) (32.0*r3)+32*(int) (32.0*r2) +(int)(32.0*r1)]++;
        hist5d[65536*(int)(16.0*r5)+4096*(int) (16.0*r4)+256*(int) (16.0*r3)+16*(int) (16.0*r2) +(int)(16.0*r1)]++;
    }
    
    //analyze the chi square values from the histograms and increment p-value counts
    pvalue=ChiSquareLargeNPValue(hist1d,1048576);
    log10pvalue = log10(pvalue);
    countpvalues[0]++;
    for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;
    
    pvalue=ChiSquareLargeNPValue(hist2d,1048576);
    log10pvalue = log10(pvalue);
    countpvalues[0]++;
    for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;
 
    pvalue = ChiSquareLargeNPValue(hist3d,1000000);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;

    pvalue = ChiSquareLargeNPValue(hist4d,1048576);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;

    pvalue = ChiSquareLargeNPValue(hist5d,1048576);
        log10pvalue = log10(pvalue);
        countpvalues[0]++;
        for (int ip=1;ip<=8;ip++) if (-log10pvalue > ip) countpvalues[ip]++;


    }
    
    // get final state of the generator in each process and 
    // compare to the predicted state given in eqn. (7) in the reference above
    getparameters_randomph(integerparameters);
    unsigned long mf, sf, cf, q, r;
    q=(5*ntests*nrounds)/(p0-1);
    r=(5*ntests*nrounds)%(p0-1);
    sf=s0;
    mf=(p0*(p0-1)/2)%n0;
    mf= (q*mf)%n0;
    mf=(mf+m0)%n0;  // m_{q (p-1)} = (m_0 + q p(p-1)/2 ) mod n
    // step forward r steps
    for (int j=1;j<=r;j++)
    {
        sf=(a0*sf)%p0;
        mf=(mf+sf)%n0;
    }
    cf=powermod(mf,e0,n0);
    cout << "final state of the generator in process " <<processid << ": m, s, c, n, e, p, a" << endl;
    for (int j=0;j<7;j++) cout << integerparameters[j] << " "; cout << endl;
    
    cout << "expected final state of the generator in process " <<processid << ": m, s, c, n, e, p, a" <<endl; 
    cout << mf << " " << sf << " " << cf << " " << n0 << " " << e0 << " " << p0 << " " << a0 << endl<<endl;

    
// Send p-value counts from each processid>0 back to processid=0
#ifdef MPI 
    if (processid != 0) MPI_Send(countpvalues, 9, MPI_LONG, 0, 0, MPI_COMM_WORLD);
#endif


    if  (processid==0)
    {


        long mastercountpvalues[9];
        for (int i=0;i<9;i++) mastercountpvalues[i]=countpvalues[i];

#ifdef MPI
// Collect p-value counts from all processid>0
        for (int event=0;event<numprocesses-1;event++)
        {
            //int MPI_Probe(int source, int tag, MPI_Comm comm, MPI_Status *status)
            MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
            int tag = status.MPI_TAG;
            int process=status.MPI_SOURCE;
            if (tag==0)
            {
                MPI_Recv(countpvalues, 9, MPI_LONG, process, 0, MPI_COMM_WORLD, &status);
                for (int i=0;i<9;i++) mastercountpvalues[i]+=countpvalues[i];
            }
        }
#endif
 
        cout << "total number of chi square tests = "<<mastercountpvalues[0]<<endl;
        for (int ip=1;ip<=8;ip+=1) cout << "number of p<1e-"<<ip<<" events = "<<mastercountpvalues[ip]<<",   expected number = " << 2.0*mastercountpvalues[0]*pow(0.1,ip) << " +/- " << sqrt(2.0*mastercountpvalues[0]*pow(0.1,ip))<<endl;


        cout <<endl;
        time =((double) clock())/((double) CLOCKS_PER_SEC);
        cout << "time for process 0 = "<<time-timelast << " seconds"<<endl;
        cout << "total number of pseudorandom numbers tested = "<< 5L*ntests*nrounds*numprocesses <<endl;
        //cout << "number generated and analyzed per second in process 0 = " << (double) ntests * (double) nrounds * 5.0 / time << endl;
        timelast=time;


    }
    
// Finalize MPI session
#ifdef MPI
	MPI_Finalize();
#endif  
    
    return 0;

    }