GopherProxy

	tmontecarlo.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tmontecarlo.cpp (1949B)
	---
	1 #include <iostream>
	2 #include <vector>
	3 #include <cstdlib> // for random functionality
	4 #include <cmath> // NaN
	5 #include "header.h"
	6 #include "montecarlo.h"
	7
	8 // Constructor
	9 MC::MC(Floattype (*function)(const std::vector<Floattype>),
	10 const std::vector<Floattype> a_in,
	11 const std::vector<Floattype> b_in,
	12 const Lengthtype N_in)
	13 : d(a_in.size()), f(function), a(a_in), b(b_in), x(d), N(N_in)
	14 {
	15 // Check that problem is at least 1D
	16 if (d < 1)
	17 std::cerr << "Error! Problem must be at least 1D\n";
	18
	19 // Check that the user has specified at least 1 random sample
	20 if (N < 1)
	21 std::cerr << "Error! The algorithm should evaluate at least 1 point!…
	22
	23 // Check input data
	24 for (Lengthtype i=0; i<d; ++i)
	25 if (a[i] >= b[i])
	26 std::cerr << "Error! a >= b in dimension " << i << '\n';
	27
	28 // Initialize result and error to NaN
	29 Q = NAN;
	30 err = NAN;
	31 }
	32
	33 // Plain Monte Carlo integrator
	34 void MC::plain()
	35 {
	36
	37 // Set volume of sample space
	38 set_volume();
	39
	40 Floattype sum = 0.0f, sumsum = 0.0f, fx;
	41
	42 for (Lengthtype i=0; i<N; ++i) {
	43 x = random_x(); // Random sample point inside intervals
	44 fx = f(x);
	45 sum += fx;
	46 sumsum += fx*fx;
	47 }
	48
	49 Floattype average = sum/N;
	50 Floattype variance = sumsum/N - average*average;
	51
	52 // Write results
	53 Q = average * V;
	54 err = sqrt(variance/N)*V;
	55 }
	56
	57
	58 // Calculate volume by multiplying interval in each dimension
	59 void MC::set_volume()
	60 {
	61 V = 1.0f;
	62 for (Lengthtype i=0; i<d; ++i)
	63 V *= b[i] - a[i];
	64 }
	65
	66 // Draw pseudo-random position in sample space
	67 std::vector<Floattype> MC::random_x()
	68 {
	69 std::vector<Floattype> pos(d);
	70 for (Lengthtype i=0; i<d; ++i) {
	71 // Random number in [a;b] interval in dimension
	72 pos[i] = (b[i] - a[i]) * ((Floattype)rand()/RAND_MAX) + a[i];
	73 }
	74 return pos;
	75 }
	76
	77
	78 // Print results
	79 void MC::show()
	80 {
	81 std::cout << "Integral Q = " << Q << ", error = " << err << '\n';
	82 }
	83
	84 // Return the error
	85 Floattype MC::error()
	86 {
	87 return err;
	88 }