# Difference between revisions of "Tutorials:Code-06 MonteCarloHOWTO"

Work in progress ... The tutorial is neither finished nor corrected ...

The plan for this page is to explain how to write a Monte Carlo program from scratch using the ALPS libraries.

## Introduction

While several Monte Carlo (MC) programs are already available in the ALPS applications, clearly there is a need for helping developers to program their own MC program in the easiest and quickest possible way. Indeed there are several common tasks performed by all (Markov chain) MC algorithms which do not need to be re-programmed each time. The goal of this article is to convince that the ALPS libraries propose a simple framework to do so, and to demonstrate with examples how to use them. This includes both classical Monte Carlo as well as Quantum Monte Carlo algorithms.

### What are the advantages of the ALPS libraries over one's own tools for developing a new MC application ?

• Automatic calculation of error bars and autocorrelation times
• Automatic parallelization
• Automatic output of results (in XML) with the corresponding extraction tools
• Automatic management of the sign (for Quantum Monte Carlo simulations with a sign problem)
• Easy checkpointing & Easy restarting of simulations
• Easy management of input parameters
• Easy way to add new measurements
• Easy management of random numbers
• (Optional) Easy management of lattices
• (Optional) Easy management of quantum Hamiltonians (for Quantum Monte Carlo)
• ...

All this comes with only a small amount of programming/modification that we describe now.

## Starting your own ALPS MC program

First of all, using the ALPS libraries implies that you use C++ as the programming language. From now on, let us assume that you want to program a new MC application and that you already know what are going to be your internal data structure and the algorithm you will use for a single Monte Carlo step.

To use the ALPS libraries, you will need to create your own C++ class, that will be derived from an internal ALPS class. Let's call this class MyMonteCarlo and create a MyMonteCarlo.hpp header file.

``` #ifndef MYMC_HPP
#define MYMC_HPP
#include <alps/scheduler.h>
#include <alps/alea.h>
#include <alps/scheduler/montecarlo.h>
#include <alps/osiris.h>
#include <alps/osiris/dump.h>
#include <alps/expression.h>
using namespace std;
class MyMonteCarlo : public alps::scheduler::MCRun{
public :
MyMonteCarlo(const alps::ProcessList &,const alps::Parameters &,int);
void save(alps::ODump &) const;
void dostep();
bool is_thermalized() const;
double work_done() const;
private :
// your own internal data here ...
};
#endif
```

We now see what all these lines mean.

• First of all you need to include a few ALPS headers (some more will go here for other ALPS functionalities).
• Then we will define our MyMonteCarlo class out by deriving it from a ALPS class. The above header file implies that you will use the MCRun ALPS class, which is the simplest one. If you want to enjoy the lattice functionalities of ALPS, replace the class definition line by:
``` class MyMonteCarlo : public alps::scheduler::LatticeMCRun<>{
```

``` #include <alps/lattice.h>
```

If you moreover want to use the Model library (useful for lattice quantum models), use instead:

``` class MyMonteCarlo : public alps::scheduler::LatticeModelMCRun<>{
```
• In the public part, MyMonteCarlo(const alps::ProcessList&,const alps::Parameters&,int) is the constructor of your class, and will be needed to initialize all parameters.
• print_copyright(std::ostream&) is a simple function to output any useful information that you want to be output at the beginning of each simulation.
• The save and load functions are very useful functions where you will describe what is needed to save (load) on disk to restart the simulations after each checkpoint
• The dostep function is your main MC function: this is the function that is going to be executed every MC step.
• The is_thermalized function will tell the ALPS libraries when the thermalization part of your simulation is finished (that is, when can the measurement series start)
• The work_done function will tell the ALPS libraries what is the percentage of the simulations that is already done.

The rest of the job is now to correctly interface these functions with your program. This will be done in a file called, say, MyMonteCarlo.cpp

## Building your own ALPS MC program

``` #include "MyMonteCarlo.hpp";
/************************************ ALPS functions **********************************************/
{
out << "My own ALPS Monte Carlo program v. 0.1\n"
<< "  copyright (c) 2006 by Myself,\n"
<< "  available from the author on request\n\n";
}
```

### Management of Monte Carlo steps

Let us now concentrate on the management of MC steps. A typical situation is the following: one would like to perform a fixed number of steps for thermalization, then a fixed number of steps for the measurement part. Very often, one would like also to perform a certain amount of Monte Carlo steps between each measurement. Therefore it could be useful to define in your internal data structure the following variables (in MyMonteCarlo.hpp):

``` private :
// your own internal data here ...
int Nb_Steps;
int Nb_Thermalisation_Steps;
int Each_Measurement;
int Steps_Done_Total;
int Measurements_Done;
void do_update();
void do_measurements();
// the rest of your own internal data here ...
};
```

Here Nb_Thermalisation_Steps will be the number of thermalization steps that you ask for and Nb_Steps the number of steps after thermalization that you ask for. Each_Measurement will be the number of steps between each measurement. These three numbers will be initialized later in the constructor. Steps_Done_Total will store the current number of MC finished steps (including those of thermalization). Finally Measurements_Done will store the intermediate number of steps done between each measurement. The do_update() and do_measurements() function (that you will have to define yourself) will perform respectively one single MC step and one series of measurement.

All these definitions lead to the following simple definitions of your class member functions in MyMonteCarlo.cpp

``` bool MyMonteCarlo::is_thermalized() const
{  return (Steps_Done_Total >= Nb_Thermalisation_Steps); }
```

This function will indeed return 1 if the current number of MC steps is superior to the total number of thermalization steps asked for , 0 otherwise. The second function

``` double MyMonteCarlo::work_done() const
{ return (is_thermalized() ? (Steps_Done_Total-Nb_Thermalisation_Steps)/double(Nb_Steps) :0.); }
```

will return 0 if the simulation is not thermalized, and a number between 0 and 1 corresponding to the percentage of asked measurement steps already performed. Finally the dostep() function will look like this

``` void MyMonteCarlo::dostep()
{ do_update(); // you'll have to define what this function does later
++Steps_Done_Total; // increment the number of steps done
if (is_thermalized()) // do measurements only if simulation thermalized
{ if (++Measurements_Done==Each_Measurement) // do a measurement every Each_Measurement
{ Measurements_Done=0;
do_measurements(); // you'll have to define what this function does later
}
}
}
```

At some point you have of course to define what your program really does during one MC step, and this will be done in the do_update() function. I can't help you with that one! The measurements have been grouped in our example in the do_measurements() function that will be described below.

### The save and load functions

ALPS allows an easy treatment of the internal data to be checkpointed to disk. Let us imagine that in your internal data structure you need to checkpoint a few integers and an array of double to be able to restart your simulations

``` private :
// the rest of your own internal data here ...
int Number_of_Spins;
int MyOwnVariable;
std::vector<double> SpinArray;
...
};
```

This can be accomplished very easily by writing

``` void MyMonteCarlo::save(alps::ODump& dump) const
{ dump <<  Number_of_Spins << MyOwnVariable << SpinArray;}
```

and the load function is easily guessed to be

``` void MyMonteCarlo::load(alps::IDump& dump) const
{ dump >>  Number_of_Spins >> MyOwnVariable >> SpinArray;}
```

Note that you can dump this way most of the usual types (int,double,bool etc) and most standard containers are also available (such as vector or set).

``` struct Vertex
{ int vertex_type;
std::vector<double> coordinates;
int SomeOtherVariable; }
```

and you want to save an instantiation thereof in your Monte Carlo class.

``` private :
// the rest of your own internal data here ...
Vertex MyVertex;
...
};
```

Well, you just have to teach ALPS what to save and load in your structure

``` alps::ODump& operator<<(alps::ODump& dump, const Vertex& v)
{ return dump << v.vertex_type << v.coordinates; }
```
``` alps::IDump& operator>>(alps::IDump& dump, Vertex& v)
{ return dump >> v.vertex_type >> v.coordinates;}
```

and then you'll easily be able to add MyVertex in the main ALPS save and load functions:

``` void MyMonteCarlo::save(alps::ODump& dump) const
{ dump <<  Number_of_Spins << MyOwnVariable << SpinArray << MyVertex;}
```

If you have used the management of Monte Carlo steps described above, you also want to add this in your save/load functions:

``` void MyMonteCarlo::save(alps::ODump& dump) const
{ dump <<  Number_of_Spins << MyOwnVariable << SpinArray << MyVertex;
dump <<  Nb_Steps << Measurements_Done; }
```

### The do_measurements() function

In this function you will perform your measurements, and give them to ALPS to be treated. This is pretty simple.

``` void MyMonteCarlo::do_measurements()
{ double Energy; std::valarray<double> Correl(L);
// do the measurements in your code (update the Energy and Correl variable) ...
...
```
```     // give them to ALPS
measurements["Energy"] << Energy;
measurements["Spin Correlations"] << Correl;
}
```

and that's it? Now you might wonder: how does ALPS know what "Energy" or "Spin Correlations" is ? How does it distinguish between scalar measurements (such as Energy here) and vector ones (such as Correl) ? These two questions will be addressed below, in the Constructor of your class.

Another question might be: why did you use a std::valarray and not a std::vector for the vector observable ? This is for internal ALPS reasons, but just remember that for vector observables, ALPS need the vectors to *always* be of the same size for each different measurements (in this example, measurements["Spin Correlations"] will fail if the Correl object has not always the same size -L here-).