Merge pull request #13 from DCM-UPB/mpi_wrapper

Mpi wrapper

config_template.sh

@@ -6,6 +6,9 @@ LIBNAME="mci"
 # C++ compiler
 CC="g++"
 
+# MPI compiler wrapper
+MPICC="mpic++"
+
 # C++ flags (std=c++11 is necessary)
 FLAGS="-std=c++11 -Wall -Werror"
 

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,87 @@
+#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);
+    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);
+    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);
+    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);
+    assert( abs(average[0]-CORRECT_RESULT) < 2.*error[0] );
+
+
+    delete pdf;
+    delete obs;
+    delete mci;
+    delete [] x;
+    delete average;
+    delete error;
+
+    return 0;
+}

doc/user_manual.tex

@@ -568,6 +568,34 @@ \subsection{Multidimensional estimations} % (fold)
 Data (\verb+x+) are supposed to be organised as a $\verb+n+ \times \verb+ndim+$ matrix.
 % subsection multidimensional_estimations (end)
 
+\section{MPI-MCI} % (fold)
+\label{sec:mpimci}
+
+It is possible to use Message Passing Interface (MPI) for evaluating the MC integrals
+in parallel. The library provides a namespace \verb+MPIMCI+ that contains a few
+simple methods to use for creating a parallel MC program. Although the MCI library can be
+compiled without any MPI implementation present, to use the MPIMCI methods in
+executable code you need OpenMPI or an alternative implementation.
+\\\\Then just use MPIMCI in your code as follows and compile\&run the executable with your system's MPI wrappers (e.g mpic++ and mpirun).
+
+\begin{verbatim}
+#include "MPIMCI.hpp"
+// start of main
+const int myrank = MPIMCI::init();
+// ...
+MPIMCI::integrate(mci, Nmc, average, error);
+// ...
+if (myrank==0) { /* printout etc. */ }
+// ...
+MPIMCI::finalize();
+// end of main
+\end{verbatim}
+Remember not to encapsulate any code within \verb+if (myrank==0){}+ clauses,
+that is actually relevant to all threads. Typically this means only file and
+console output belongs there.
+\\\\For further information on options and use, check out \verb+ex2+ and the \verb+MPIMCI.hpp+ header file.
+
+% section mpivmc (end)
 
 
 \printindex

examples/README.md

@@ -5,3 +5,8 @@
 ## Example 1
 
 `ex1/`: integration of a 1-dimensional quadratic function, without and with a sampling function.
+
+
+## Example 2
+
+`ex2/`: like ex2, but using the MPI wrapper for parallel integration. Pass the wanted number of threads to run.sh.

examples/ex2/main.cpp

@@ -0,0 +1,192 @@
+#include "mpi.h"
+#include <iostream>
+#include <cmath>
+#include <math.h>
+#include <fstream>
+
+#include "MCIntegrator.hpp"
+#include "MPIMCI.hpp"
+
+
+// Observable functions
+class Parabola: public MCIObservableFunctionInterface{
+public:
+    Parabola(const int &ndim): MCIObservableFunctionInterface(ndim, 1) {}
+
+    void observableFunction(const double * in, double * out){
+        out[0] = 4.*in[0] - in[0]*in[0];
+    }
+};
+
+class NormalizedParabola: public MCIObservableFunctionInterface{
+public:
+    NormalizedParabola(const int &ndim): MCIObservableFunctionInterface(ndim, 1) {}
+
+    void observableFunction(const double * in, double * out){
+        out[0] = (4. - in[0]) * 5.;
+        if (std::signbit(in[0])) out[0] = -out[0];
+    }
+};
+
+
+
+// Sampling function
+// the 48 is for normalization (even if not strictly necessary)
+class NormalizedLine: public MCISamplingFunctionInterface{
+public:
+    NormalizedLine(const int &ndim): MCISamplingFunctionInterface(ndim, 1) {}
+
+    void samplingFunction(const double * in, double * protovalue){
+        protovalue[0] = 0.2 * abs(in[0]);
+    }
+
+    double getAcceptance(const double * protoold, const double * protonew){
+        return protonew[0] / protoold[0];
+    }
+};
+
+
+
+int main() {
+    using namespace std;
+
+    int myrank = MPIMCI::init(); // first run custom MPI init
+
+    if (myrank == 0) {
+        // intro
+        cout << "We want to compute the integral" << endl;
+        cout << "    Integral[-1:3] dx (4x - x^2)" << endl << endl;
+        cout << "Notice that we can re-write the integral as" << endl;
+        cout << "    Integral[-1:3] dx (5*sign(x)*(4-x) * |x|/5)" << endl;
+        cout << "where g(x)=|x|/5 is a candidate sampling function since it's positive on the domain [-1:3] and its integral on the domain is equal to 1." << endl << endl;
+
+        cout << "We start by initializing MPI and setting the MCI:" << endl;
+    }
+
+    const int ndim = 1;
+    MCI * mci = new MCI(ndim);
+
+    // seed the different MPI threads from a file
+    MPIMCI::setSeed(mci, "rseed.txt", 0); // offset x -> start from the x-th seed in file
+
+    if (myrank == 0) cout << "ndim = " << mci->getNDim() << endl;
+
+    // set the integration range to [-1:3]
+    double ** irange = new double*[ndim];
+    irange[0] = new double[2];
+    irange[0][0] = -1.;
+    irange[0][1] = 3.;
+    mci->setIRange(irange);
+
+    if (myrank == 0) cout << "irange = [ " << mci->getIRange(0, 0) << " ; " << mci->getIRange(0, 1) << " ]" << endl;
+
+    // initial walker position
+    double * initpos = new double[ndim];
+    initpos[0] = -0.5;
+    mci->setX(initpos);
+
+    if (myrank == 0) cout << "initial walker position = " << mci->getX(0) << endl;
+
+    // initial MRT2 step
+    double * step = new double[ndim];
+    step[0] = 0.5;
+    mci->setMRT2Step(step);
+
+    if (myrank == 0) cout << "MRT2 step = " << mci->getMRT2Step(0) << endl;
+
+    // target acceptance rate
+    double * targetacceptrate = new double[1];
+    targetacceptrate[0] = 0.7;
+    mci->setTargetAcceptanceRate(targetacceptrate);
+
+    if (myrank == 0) {
+        cout << "Acceptance rate = " << mci->getTargetAcceptanceRate() << endl;
+        cout << endl << endl;
+
+        // first way of integrating
+        cout << "We first compute the integral without setting a sampling function." << endl;
+        cout << "    f(x) = (4x - x^2) " << endl;
+        cout << "    g(x) = - " << endl << endl;
+    }
+
+    // observable
+    MCIObservableFunctionInterface * obs = new Parabola(ndim);
+    mci->addObservable(obs);
+
+    if (myrank == 0) {
+        cout << "Number of observables set = " << mci->getNObs() << endl;
+        // sampling function
+        cout << "Number of sampling function set = " << mci->getNSampF() << endl;
+    }
+
+    // integrate
+    const long Nmc = 1000000;
+    double * average = new double[mci->getNObsDim()];
+    double * error = new double[mci->getNObsDim()];
+
+    MPIMCI::integrate(mci, Nmc, average, error);
+
+    if (myrank == 0) {
+        cout << "The integral gives as result = " << average[0] << "   +-   " << error[0] << endl;
+        cout << "--------------------------------------------------------" << endl << endl;
+
+        // second way of integrating
+        cout << "Now we compute the integral using a sampling function." << endl;
+        cout << "    f(x) = 5*sign(x)*(4-x) " << endl;
+        cout << "    g(x) = |x|/5 " << endl << endl;
+    }
+
+    // observable
+    delete obs;
+    obs = new NormalizedParabola(ndim);
+    mci->clearObservables();  // we first remove the old observable
+    mci->addObservable(obs);
+
+    if (myrank == 0) cout << "Number of observables set = " << mci->getNObs() << endl;
+
+
+    // sampling function
+    MCISamplingFunctionInterface * sf = new NormalizedLine(ndim);
+    mci->addSamplingFunction(sf);
+
+    if (myrank == 0) cout << "Number of sampling function set = " << mci->getNSampF() << endl;
+
+
+    // integrate
+    MPIMCI::integrate(mci, Nmc, average, error);
+
+    if (myrank == 0) {
+        cout << "The integral gives as result = " << average[0] << "   +-   " << error[0] << endl;
+        cout << "--------------------------------------------------------" << endl << endl;
+
+
+        // final comments
+        cout << "Using a sampling function in this case gives worse performance. In fact, the error bar is larger." << endl;
+        cout << "This implies that the variance of the re-factored f(x) written for introducing a sampling function, is larger than the original f(x)." << endl;
+    }
+
+    // deallocate per-thread allocations
+    delete sf;
+
+    delete obs;
+
+    delete[] targetacceptrate;
+
+    delete[] step;
+
+    delete[] initpos;
+
+    delete[] irange[0];
+    delete[] irange;
+
+    delete[] average;
+    delete[] error;
+
+    delete mci;
+
+    // finalize MPI
+    MPIMCI::finalize();
+
+    // end
+    return 0;
+}

examples/ex2/rseed.txt

@@ -0,0 +1 @@
+../../res/rseed.txt

examples/ex2/run.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+source ../../config.sh
+OS_NAME=$(uname)
+
+# Delete old compiled files
+\rm -f exe
+\rm -f *.o
+
+#runtime dynamic library path
+RPATH="$(pwd)/../.."
+
+FLAGS_TO_USE=$OPTFLAGS
+
+# Build the main executable
+echo "$MPICC $FLAGS $FLAGS_TO_USE -I$(pwd)/../../src/ -c *.cpp"
+$MPICC $FLAGS $FLAGS_TO_USE -Wall -I$(pwd)/../../src/ -c *.cpp
+
+case ${OS_NAME} in
+    "Darwin")
+        echo "$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -o exe *.o -l${LIBNAME}"
+        $MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -o exe *.o -l${LIBNAME}
+        ;;
+    "Linux")
+        echo "$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -Wl,-rpath=${RPATH} -o exe *.o -l${LIBNAME}"
+        $MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -Wl,-rpath=${RPATH} -o exe *.o -l${LIBNAME}
+        ;;
+esac
+
+echo "Rebuilt the executable file"
+echo ""
+echo ""
+
+# Run the debugging executable
+echo "Ready to run!"
+echo ""
+echo "--------------------------------------------------------------------------"
+echo ""
+echo ""
+echo ""
+mpirun -np $1 exe