#include
#include
#include
#include
using namespace std;
unsigned int MonteCarlo(double (*)(double), double, double, double*, double*, double*); //Monte Carlo method
unsigned int MonteCarloCauchy(double (*)(double), double, double, double*, double*, double*); //Monte Carlo method using Cauchy distribution
void Log(double*, unsigned int); //Applies the common logrithim to the contents of an array except for the values of 0.
inline double Const(double); //Returns 1
inline double Linear(double); //Returns x
inline double Quatric(double); //Returns x^3
inline double Cauchy(double); //Returns a Cauchy distribution of x_0=.5, gamma=.01
inline double CauchyPDF(double);//Returns the value of the Cauchy
inline double RandC(); //Returns a Cauchy(.5,.1) random number between 0 and 1
inline double Rand(); //Returns a uniformally distributed random number between 0 and 1
void Plot(double*, double*, double*, unsigned int, char*); //Plot the values generated by the MonteCarlo Method using ROOT
const unsigned int MAX_Array = 38500;
const unsigned int Sample_Rate = 40000;
const double Gamma = .01; //Spread parameter of Cauchy Random number generator
const double Mode = .5; //Center of Cauchy Random number generator
int main()
{
srand(time(NULL)); //Seed the random number generator
double* Value;
double* RMS;
double* Error;
unsigned int Number;
#pragma omp parallel num_threads(4) private(Value, RMS, Error, Number) //Multi-threading will cause race conditions for output
switch(omp_get_thread_num())
{
case 0:
Value = new double[MAX_Array]; //Assigns a new block of array different from the other Value pointers in the other threads
RMS = new double [MAX_Array];
Error = new double [MAX_Array];
Number = MonteCarlo(Const, 0., 1., Value, RMS, Error);
cout << "Const Value = " << Value[Number] << endl; //Output information for use by the user.
cout << "Const RMS = " << RMS[Number] << endl;
cout << "Const Error = " << Error[Number] << endl;
cout << "Number = " << Number << endl;
Plot(Value, RMS, Error, Number, "Const.pdf");
Number = MonteCarlo(Cauchy, 0., 1., Value, RMS, Error);
cout << "Cauchy Value = " << Value[Number] << endl; //Output information for use by the user.
cout << "Cauchy RMS = " << RMS[Number] << endl;
cout << "Cauchy Error = " << Error[Number] << endl;
cout << "Number = " << Number << endl;
Log(Error, Number);
Log(RMS, Number); //It makes it more readible
Plot(Value, RMS, Error, Number, "Cauchy.pdf");
break;
case 1:
Value = new double[MAX_Array]; //Assigns a new block of array different from the other Value pointers in the other threads
RMS = new double [MAX_Array];
Error = new double [MAX_Array];
Number = MonteCarlo(Linear, 0., 1., Value, RMS, Error);
cout << "Linear Value = " << Value[Number] << endl; //Output information for use by the user.
cout << "Linear RMS = " << RMS[Number] << endl;
cout << "Linear Error = " << Error[Number] << endl;
cout << "Number = " << Number << endl;
Log(Error, Number); //It makes it more readible
Log(RMS, Number); //It makes it more readible
Plot(Value, RMS, Error, Number, "Linear.pdf");
break;
case 2:
Value = new double[MAX_Array]; //Assigns a new block of array different from the other Value pointers in the other threads
RMS = new double [MAX_Array];
Error = new double [MAX_Array];
Number = MonteCarlo(Quatric, 0., 1., Value, RMS, Error);
cout << "Quatric Value = " << Value[Number] << endl; //Output information for use by the user.
cout << "Quatric RMS = " << RMS[Number] << endl;
cout << "Quatric Error = " << Error[Number] << endl;
cout << "Number = " << Number << endl;
Log(Error, Number);
Log(RMS, Number); //It makes it more readible
Plot(Value, RMS, Error, Number, "Quatric.pdf");
break;
case 3:
Value = new double[MAX_Array]; //Assigns a new block of array different from the other Value pointers in the other threads
RMS = new double [MAX_Array];
Error = new double [MAX_Array];
Number = MonteCarloCauchy(Cauchy, 0., 1., Value, RMS, Error);
cout << "Cauchy Value = " << Value[Number] << endl; //Output information for use by the user.
cout << "Cauchy RMS = " << RMS[Number] << endl;
cout << "Cauchy Error = " << Error[Number] << endl;
cout << "Number = " << Number << endl;
Log(Error, Number);
Log(RMS, Number); //It makes it more readible
Plot(Value, RMS, Error, Number, "CauchyCauchy.pdf");
break;
default:
cout << "Something has gone wrong" << endl;
}
return(0);
}
unsigned int MonteCarlo(double (*Func)(double), double a, double b, double* Value, double* RMS, double* Error) //Monte Carlo method
{
unsigned int i;
double Entry; //Measurment
double First = 0; //Expectation of the first moment
double Second = 0; //Expectation of the second moment
double Random; //Random number holder
for(i = 0; i 2*Sample_Rate+1 && Error[i/Sample_Rate] < .0001) //Early break conditions, but must sample at least 300 elements before breaking for the graph
break;
}
}
Value[i/Sample_Rate] = First; //To make sure that the last value has been placed
RMS[i/Sample_Rate] = sqrt((Second-pow(First,2)));
Error[i/Sample_Rate] = sqrt((Second-pow(First,2))/double(i))/abs(First);
Error[0] = Error[1]; //Resolve a divide by zero proplem at i==0
return(i/Sample_Rate);
}
unsigned int MonteCarloCauchy(double (*Func)(double), double a, double b, double* Value, double* RMS, double* Error) //Monte Carlo method
{
unsigned int i;
double Entry; //Measurment
double First = 0; //Expectation of the first moment
double Second = 0; //Expectation of the second moment
double Random; //Random number holder
for(i = 0; i= a && Random <= b) //Normal order
{
Entry = Func(Random)/CauchyPDF(Random); //Find the value of F(x) and normalize for the width of the integral and the distribution of random numbers
First = (double(i)*First+Entry)/(i+1.); //Calculate the running average of the first moment
Second = (double(i)*Second+pow(Entry,2))/(i+1.);//Calculate the running average of the second moment
}
else if(Random >= b && Random <= a) //Reverse order
{
Entry = -Func(Random)/CauchyPDF(Random); //Find the value of F(x) and normalize for the width of the integral and the distribution of random numbers
First = (double(i)*First+Entry)/(i+1.); //Calculate the running average of the first moment
Second = (double(i)*Second+pow(Entry,2))/(i+1.);//Calculate the running average of the second moment
}
else //Average in a zero for out of range
{
First *= double(i)/(i+1.); //Calculate the running average of the first moment
Second *= double(i)/(i+1.); //Calculate the running average of the second moment
}
if(i%Sample_Rate == 0) //Every Sample_Rate samples record and check the statistics of the run to date
{
Value[i/Sample_Rate] = First; //Running average
RMS[i/Sample_Rate] = sqrt((Second-pow(First,2))); //Running RMS of the data
Error[i/Sample_Rate] = sqrt((Second-pow(First,2))/double(i))/abs(First); //Running uncertainity of the data
if(i > 2*Sample_Rate+1 && Error[i/Sample_Rate] < .0001) //Early break conditions, but must sample at least 300 elements before breaking for the graph
break;
}
}
Value[i/Sample_Rate] = First; //To make sure that the last value has been placed
RMS[i/Sample_Rate] = sqrt((Second-pow(First,2)));
Error[i/Sample_Rate] = sqrt((Second-pow(First,2))/double(i))/abs(First);
Error[0] = Error[1]; //Resolve a divide by zero proplem at i==0
return(i/Sample_Rate);
}
void Log(double* Array, unsigned int N)
{
for(unsigned int i = 0; i < N; i++)
if(Array[i] > 0) //If the content positive, take the common log of it. if not positive leave it as doing anything else is not well defined and will only cause problms later
Array[i] = log10(Array[i]);
return;
}
inline double RandC()
{
return(Mode+Gamma*tan(M_PI*(Rand()-.5)));
}
inline double Const(double Arg)
{
return(1.);
}
inline double Linear(double Arg)
{
return(Arg);
}
inline double Quatric(double Arg)
{
return(pow(Arg,3));
}
inline double Cauchy(double Arg)
{
return(.01/(M_PI*(pow(Arg-.5,2)+.0001)));
}
inline double CauchyPDF(double Arg)
{
return(Gamma/(M_PI*(pow(Arg-Mode,2)+pow(Gamma,2))));
}
inline double Rand()
{
return(rand()/double(RAND_MAX));
}