Adapt MPIMCI interface to changes and add unittest 3

Author: Ithanil

debug/unittests/README.md

@@ -9,3 +9,7 @@
 ## Unit Test 2
 
 `ut2/`: Like ut1, but with one one-dimensional and one three-dimensional observable.
+
+## Unit Test 3
+
+`ut3/`: Like ut1, but checking with fixed number of findMRT2step/decorrelation steps and fixed blocking technique

debug/unittests/ut3/main.cpp

@@ -0,0 +1,91 @@
+#include "MCIntegrator.hpp"
+#include "MCISamplingFunctionInterface.hpp"
+
+#include <iostream>
+#include <math.h>
+#include <assert.h>
+
+using namespace std;
+
+
+
+class ThreeDimGaussianPDF: public MCISamplingFunctionInterface{
+public:
+    ThreeDimGaussianPDF(): MCISamplingFunctionInterface(3, 1){}
+    ~ThreeDimGaussianPDF(){}
+
+    void samplingFunction(const double *in, double * protovalues){
+        protovalues[0] = (in[0]*in[0]) + (in[1]*in[1]) + (in[2]*in[2]);
+    }
+
+
+    double getAcceptance(const double * protoold, const double * protonew){
+        return exp(-protonew[0]+protoold[0]);
+    }
+
+};
+
+
+class XSquared: public MCIObservableFunctionInterface{
+public:
+    XSquared(): MCIObservableFunctionInterface(3, 1){}
+    ~XSquared(){}
+
+    void observableFunction(const double * in, double * out){
+        out[0] = in[0] * in[0];
+    }
+};
+
+
+
+int main(){
+    const long NMC = 10000;
+    const double CORRECT_RESULT = 0.5;
+
+    ThreeDimGaussianPDF * pdf = new ThreeDimGaussianPDF();
+    XSquared * obs = new XSquared();
+
+    MCI * mci = new MCI(3);
+    mci->setSeed(5649871);
+    mci->addSamplingFunction(pdf);
+    mci->addObservable(obs);
+    // the integral should provide 0.5 as answer!
+
+    double * x = new double[3];
+    x[0] = 5.; x[1] = -5.; x[2] = 10.;
+
+    double * average = new double;
+    double * error = new double;
+
+    // this integral will give a wrong answer! This is because the starting point is very bad and initialDecorrelation is skipped (as well as the MRT2step automatic setting)
+    mci->setX(x);
+    mci->integrate(NMC, average, error, 0, 0);
+    std::cout << average[0] << " " << error[0] << std::endl;
+    assert( abs(average[0]-CORRECT_RESULT) > 2.*error[0] );
+
+    // this integral, instead, will provide the right answer
+    mci->setX(x);
+    mci->integrate(NMC, average, error, 10, 1000);
+    std::cout << average[0] << " " << error[0] << std::endl;
+    assert( abs(average[0]-CORRECT_RESULT) < 2.*error[0] );
+
+    // now, doing an integral without finding again the MRT2step and doing the initialDecorrelation will also result in a correct result
+    mci->integrate(NMC, average, error, 0, 0);
+    std::cout << average[0] << " " << error[0] << std::endl;
+    assert( abs(average[0]-CORRECT_RESULT) < 2.*error[0] );
+
+    // and using fixed blocking also gives the same result
+    mci->integrate(NMC, average, error, 0, 0, 15);
+    std::cout << average[0] << " " << error[0] << std::endl;
+    assert( abs(average[0]-CORRECT_RESULT) < 2.*error[0] );
+
+
+    delete pdf;
+    delete obs;
+    delete mci;
+    delete [] x;
+    delete average;
+    delete error;
+
+    return 0;
+}

src/MPIMCI.hpp

@@ -65,7 +65,13 @@ namespace MPIMCI
 
 
     // integrate in parallel and accumulate results
-    void integrate(MCI * const mci, const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation=true, bool use_mpi=true) // by setting use_mpi to false you can use this without requiring MPI
+
+    void integrate(MCI * const mci, const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation=true, size_t nblocks=0, bool use_mpi=true) // auto-mode-wrapper
+    {
+        integrate(mci, Nmc, average, error, findMRT2step ? -1 : 0, initialdecorrelation ? -1 : 0, nblocks, use_mpi);
+    }
+
+    void integrate(MCI * const mci, const long &Nmc, double * average, double * error, int NfindMRT2stepIterations, int NdecorrelationSteps, size_t nblocks=0, bool use_mpi=true) // by setting use_mpi to false you can use this without requiring MPI
     {
         if (use_mpi) {
             // make sure the user has MPI in the correct state
@@ -81,7 +87,7 @@ namespace MPIMCI
             double myAverage[mci->getNObsDim()];
             double myError[mci->getNObsDim()];
 
-            mci->integrate(Nmc, myAverage, myError, findMRT2step, initialdecorrelation);
+            mci->integrate(Nmc, myAverage, myError, NfindMRT2stepIterations, NdecorrelationSteps, nblocks);
 
             for (int i=0; i<mci->getNObsDim(); ++i) {
                 MPI_Allreduce(&myAverage[i], &average[i], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
@@ -94,7 +100,7 @@ namespace MPIMCI
             }
         }
         else {
-            mci->integrate(Nmc, average, error, findMRT2step, initialdecorrelation); // regular single thread call
+            mci->integrate(Nmc, average, error, NfindMRT2stepIterations, NdecorrelationSteps, nblocks); // regular single thread call
         }
     }