tmontecarlo.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) 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!\n";
           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 }