/*
   Monte Carlo simulation
   https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm

   NB requires C++ V11 standard.

   http://www.cplusplus.com/reference/random/uniform_real_distribution/
   http://www.cplusplus.com/reference/random/normal_distribution/
*/

#include <iostream>
#include <iomanip>
#include <random>
#include <map>

using namespace std;

// kboltzmann*T is our energy unit, we set it to one
const double kT=1.0;

// Model here is that of a harmonic spring of equilibrium distance 2.0
const double kspring=5.0;
const double xeq = 2.0;


double Energy(double x)
// Energy of spring with a resting position.
{
    return 0.5*kspring*(x-xeq)*(x-xeq);
}

double ProposeShift(double xnow, double &xnew, double rand)
//
{
    double dx=(rand-0.5)*0.25;    // Transforms [0,1] into [-0.25,0.25]
    xnew = xnow+dx;
}

double Target(double x)
// Unnormalized target distribution
{
    return exp(-Energy(x)/kT);
}

int main()
{
    mt19937 re;
    uniform_real_distribution<double> uniform(0,1);

    const int N = 100000;

    // starting position
    double x=0;
    int accepted=0;

    double M1=0;
    double M2=0;
    double M3=0;
    double M4=0;

    map<int,int> hist;

    for (int i=0 ; i<N ; i++)
    {
        // Propose new state
        double rand1 = uniform(re);
        double xpropose;
        ProposeShift(x, xpropose, rand1);

        // Metropolis-Hastings criterion
        double choice=uniform(re);
        double alpha=Target(xpropose)/Target(x);   // Thus exp(-(E(xproposed)-E(x))/kT)
        if (choice<=alpha)
        {
            x=xpropose;
            accepted++;
        }

        // Sample
        M1+=x;
        M2+=x*x;
        M3+=pow(x,3);
        M4+=pow(x,4);

        hist[int(x*4)]++;
    }

    cout << "Fraction excepted: " << 1.0*accepted/N << endl;

    M1/=N;  M2/=N;   M3/=N;   M4/=N;

    double avg=M1;                 // Average
    double var=M2-M1*M1;           // Variance
    double std=sqrt(var);          // standard deviation

    cout << "Avg=<X>="           << avg << endl;
    cout << "Var=<X^2>-<X>^2="   << var << endl;
    cout << "Std.dev=sqrt(Var)=" << std << endl;
    cout << "<X^2>="             << M2 << endl;
    cout << "<X^3>="             << M3 << endl;
    cout << "<X^4>="             << M4 << endl;

    cout << "A discrete_distribution:" << std::endl;

    for (auto& key_val : hist )
    {
        cout <<  setw(10) << 1.0*key_val.first/4 << ": ";
        for (int j=0 ; j<key_val.second*40/N; j++) cout << "*";
        cout << endl;
    }

    return 0;
}