Monte-carlo exercise

homework/Monte-Carlo-integration/Makefile

+CFLAGS = -O1 -std=gnu11
+CFLAGS += $(shell /usr/bin/gsl-config --cflags)
+LDLIBS += $(shell /usr/bin/gsl-config --libs)
+LDLIBS += -lm
+CC = gcc
+
+.PHONEY: default
+default: main run out.gnuplot_error.png
+
+out.gnuplot_error.png: data.txt Makefile
+	echo '\
+	set terminal pngcairo dashed size 1500,600;\
+	set output "$@";\
+	set key top center;\
+	set tics out;\
+	set xlabel "Nr of Samples in each dimension";\
+	set xlabel font "Helvetica,20";\
+	set xzeroaxis linetype -1 linewidth 2.5;\
+	set yzeroaxis linetype -1 linewidth 2.5;\
+	set ylabel "Err. estimate";\
+	set ylabel font "Helvetica,20";\
+	set title "Comparison of errors";\
+	set title font "Helvetica,20";\
+	set grid ;\
+	set xrange [0:1500];\
+	set yrange [0:1];\
+	plot \
+	 "data.txt" using 1:2 with line linewidth 2 title "Pseudo random"\
+	,"data.txt" using 1:3 with line linewidth 2 title "Quasi random"\
+	' | gnuplot
+
+
+.PHONEY: run
+run:
+	./main > out.txt
+
+.PHONEY: clean
+clean:
+	rm main *.o *.txt *.png
+
+main: main.o functions.o
+
+functions.o: functions.c functions.h
+	$(CC) $(CFLAGS) -c functions.c $(LDLIBS) 
+
+main.o: main.c functions.h
+	$(CC) $(CFLAGS) -c main.c $(LDLIBS) 
+
+

homework/Monte-Carlo-integration/functions.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include "functions.h"
+
+
+double f(double* x)
+{
+    return 1.0/pow(M_PI,3)*1.0/(1-cos(x[0])*cos(x[1])*cos(x[2]));  
+}
+double g(double* x) 
+{
+    return 1.0/pow(x[0],2)+x[1]+1.0/pow(x[2],2);
+}
+
+void rnd_xs(double* a, double* b, double* x, int dim)
+{
+    for(int i = 0; i < dim; i++)
+    {
+        x[i] = a[i]+RND*(b[i]-a[i]);
+    }
+}
+
+double volume(double* a, double* b, int dim)
+{
+    double V = 1.0;
+    for(int i = 0; i < dim; i++)
+    {
+        V *= b[i]-a[i];
+    }
+    return V;
+}
+
+double plain_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err)
+{
+    double V = volume(a,b,dim);
+    double sum1 = 0.0, sum2 = 0.0, x[dim];
+    for(int i = 0; i < N; i++)
+    {
+        rnd_xs(a,b,x,dim);
+        sum1 += f(x);
+        sum2 += f(x)*f(x);
+    }
+    double variance = sqrt(sum2/((double) N)-pow(sum1/((double) N),2));
+    *err += V*variance/sqrt(N);
+    double result = (V/N)*sum1;
+    return result;
+}
+
+double corput_number(int base, int number)
+{
+    double q = 0, bk = (double)1/base;
+    while(number>0)
+    {
+        q += (number % base)*bk;
+        number /= base;
+        bk /= base;
+    }
+    return q;
+}
+
+double quasi_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err)
+{
+    double V = volume(a,b,dim);
+    double sum1 = 0.0, sum2 = 0.0, x[dim], x2[dim];
+    int primes[20] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71}; 
+    assert(dim < 11 && "The number of dimensions exceeds the number of primes given\n");
+    for(int i = 1; i < N+1; i++)
+    {
+        for(int j = 0; j < dim; j++)
+        {
+            x[j] = a[j]+corput_number(primes[j],i)*(b[j]-a[j]);
+        }
+        sum1 += f(x);
+    }
+    double result1 = (V/N)*sum1;
+    for(int i = 1; i < N + 1 ; i++)
+    {
+        for(int j = 0; j < dim; j++)
+        {
+            x2[j] = a[j]+corput_number(primes[j+dim],i)*(b[j]-a[j]);
+        }
+        sum2 += f(x2);
+    }
+    double result2 = (V/N)*sum2;
+    *err += fabs(result1-result2);
+    return result2;
+}
+
+double stratified_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err)
+{
+    double nmin = 5000;
+    if(N<nmin) {return plain_monte_carlo(a,b,f,dim,N,err);}
+    double x[dim], midway[dim], midway2[dim], subdiv_vars[2*dim];
+    for(int d = 0; d < dim; d++) {
+        double subdiv1_sum1=0, subdiv1_sum2=0, subdiv2_sum1=0, subdiv2_sum2=0;
+        for(int i=0;i<dim;i++) {midway[i]=b[i];}
+        midway[d] -= (b[d]-a[d])/2;
+        for(int i = 0; i<nmin/2; i++) {
+            rnd_xs(a,midway,x,dim);
+            subdiv1_sum1 += f(x);
+            subdiv1_sum2 += f(x)*f(x);
+        }
+        for(int i=0;i<dim;i++) {midway[i]=a[i];}
+        midway[d] += (b[d]-a[d])/2;
+        for(int i = 0; i<nmin/2; i++) {
+            rnd_xs(midway,b,x,dim);
+            subdiv2_sum1 += f(x);
+            subdiv2_sum2 += f(x)*f(x);
+        }
+        subdiv_vars[2*d] = sqrt(subdiv1_sum2/((double)(nmin/2))-pow(subdiv1_sum1/((double)(nmin/2)),2));
+        subdiv_vars[2*d+1] = sqrt(subdiv2_sum2/((double)(nmin/2))-pow(subdiv2_sum1/((double)(nmin/2)),2));
+    }
+    int max_index = 0;
+    for(int i=1;i<2*dim;i++) {
+        if(subdiv_vars[i]>subdiv_vars[max_index]) {max_index=i;}
+    }
+    int best_dim = (int)(max_index/2);
+    int N1, N2;
+    if(max_index%2==0) {N2=(int)((subdiv_vars[max_index+1]/(subdiv_vars[max_index]+subdiv_vars[max_index+1]))*N);N1 = N-N2;}
+    if(max_index%2!=0) {N1=(int)((subdiv_vars[max_index-1]/(subdiv_vars[max_index]+subdiv_vars[max_index-1]))*N);N2 = N-N1;}
+    for(int i = 0; i < dim; i++) {
+        midway[i]=b[i];
+        midway2[i]=a[i];
+    }
+    midway[best_dim] -= (b[best_dim]-a[best_dim])/2;
+    midway2[best_dim] += (b[best_dim]-a[best_dim])/2;
+    return stratified_monte_carlo(a,midway,f,dim,N1,err)+stratified_monte_carlo(midway2,b,f,dim,N2,err);
+}

homework/Monte-Carlo-integration/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include "functions.h"
+#define RND (double)rand()/RAND_MAX
+
+void rnd_x(double* a, double* b, double* x, int dim);
+double volume(double* a, double* b, int dim);
+double f(double* x);
+double g(double* x);
+double corput_number(int base, int number);
+double plain_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err);
+double quasi_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err);
+double stratified_monte_carlo(double* a, double* b, double f(double* x), int dim, int N, double* err);
+
+
+
+
+#endif
\ No newline at end of file

homework/Monte-Carlo-integration/main.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include "functions.h"
+
+
+
+int main()
+{
+FILE* data = fopen("data.txt","w");
+double a[3] = {0.0,0.0,0.0}, b[3] = {M_PI,M_PI,M_PI};
+double err=0, err1=0;
+double result = plain_monte_carlo(a,b,f,3,100000,&err);
+printf(" 1/pi^3*[1-cos(x)cos(y)cos(z)]^(-1) dxdydz all from 0-PI with plain mc  = %g, with err = %g\n",result,err);
+err=0; err1=0;
+result = quasi_monte_carlo(a,b,f,3,100000,&err);
+printf(" 1/pi^3*[1-cos(x)cos(y)cos(z)]^(-1) dxdydz all from 0-PI with quasi rnd  = %g, with err = %g\n",result,err);
+err=0; err1=0;
+
+for(int i = 0; i < 1500; i++)
+{
+    result = plain_monte_carlo(a,b,f,3,i,&err);
+    double result1 = quasi_monte_carlo(a,b,f,3,i,&err1);
+    fprintf(data,"%10d %10g %10g\n",i,err,err1);
+    err=0; err1=0;
+}
+fclose(data);
+
+double a_test[3] = {1.0,1.0,1.0}, b_test[3] = {2.0,2.0,2.0};
+result = stratified_monte_carlo(a_test,b_test,g,3,50000,&err);
+printf("Less singular function chosen for stratisfied\n");
+printf(" 1/x^2+y+1(z^2) dxdydz all from 1-2 with stratisfied = %g, with err = %g\n",result,err);
+printf("Answer should be 2.5\n");
+err=0; err1=0;
+}

homework/Monte-Carlo-integration/out.txt

+ 1/pi^3*[1-cos(x)cos(y)cos(z)]^(-1) dxdydz all from 0-PI with plain mc  = 1.41279, with err = 0.0333786
+ 1/pi^3*[1-cos(x)cos(y)cos(z)]^(-1) dxdydz all from 0-PI with quasi rnd  = 1.39537, with err = 0.0161598
+Less singular function chosen for stratisfied
+ 1/x^2+y+1(z^2) dxdydz all from 1-2 with stratisfied = 2.50207, with err = 0.00632601
+Answer should be 2.5