vmc.cpp

/*

vmc.cpp ... Variational Monte Carlo.
Copyright (C) 2010  KennyTM~ <kennytm@gmail.com>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <cmath>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_rng.h>
#include <ctime>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <cstring>
#include <algorithm>
#include <cassert>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_siman.h>

using namespace std;

int N = 30000;

// These are functions to map a number from [x, y] to (-inf, inf), so the minimizer can respect the bounds.
static double forward_transform(double x, double low, double high);
static double backward_transform(double hx, double low, double high);
static double move_by_bounded(double& x, double delta, double low, double high);

struct Parameters {
	enum { Count = 2 };
	
	double params[Count];
	
	void read(const gsl_vector* v) {
		memcpy(params, v->data, sizeof params);
	}
	void read(const gsl_vector* v, const Parameters& low, const Parameters& high) {
		read(v);
		for (int i = 0; i < Count; ++ i)
			params[i] = backward_transform(params[i], low.params[i], high.params[i]);
	}
	void write(gsl_vector* v) const {
		memcpy(v->data, params, sizeof params);
	}
	void write(gsl_vector* v, const Parameters& low, const Parameters& high) const {
		for (int i = 0; i < Count; ++ i)
			gsl_vector_set(v, i, forward_transform(params[i], low.params[i], high.params[i]));
	}

	void print() const {
		printf("\ta->%g, b->%g", params[0], params[1]);
	}
	
/*
	void random_walk(const gsl_rng* rng, double sigma, const Parameters& low, const Parameters& high) {
		move_by_bounded(a, (gsl_rng_uniform(rng)*2-1)*sigma, low.a, high.a);
		move_by_bounded(n, (gsl_rng_uniform(rng)*2-1)*sigma, low.n, high.n);
	}
	
	double distance(const Parameters& other, const Parameters& low, const Parameters& high) const {
		double da = forward_transform(a, low.a, high.a) - forward_transform(other.a, low.a, high.a);
		double dn = forward_transform(n, low.n, high.n) - forward_transform(other.n, low.n, high.n);
		return hypot(da, dn);
	}
*/
};

struct Point {
	double x, y, z;
	
	enum { Dimensions = 1 };
	
//	Point(double x_ = 0) : x(x_) {}

	void fill(double* t) const {
		t[0] = x;
		/*
		t[1] = y;
		t[2] = z;
		*/
	}
};


double V(const Point& p) {
	if (fabs(p.x) < 1)
		return -1;
	else
		return 0;
}

/*
double laplacian_V(const Point& p) {
	return 0;
}
*/

double psi(const Point& p, const Parameters& a_) {
	double a = a_.params[0], b = a_.params[1];
	
	double xx = fabs(p.x);
	if (xx >= 1)
		return exp(-a * xx);
	else
		return exp(-a) * cos(b * xx) / cos(b);
}


double laplacian_psi(const Point& p, const Parameters& a_) {
	double a = a_.params[0], b = a_.params[1];
	
	double xx = fabs(p.x);
	if (xx >= 1)
		return a * a * exp(-a * xx);
	else
		return -b * b  * cos(b * xx) * exp(-a) / cos(b);
}


/*
double biharmonic_psi(const Point& p, const Parameters& params) {
	if (p.x < 1)
		return -params.n * (params.n-1) * (params.n-2) * (params.n-3) * pow(p.x, params.n-4);
}
*/

//------------------------------------------------------------------------------

static double forward_transform(double x, double low, double high) {
	double nx = (x - low) / (high - low);
	nx = nx * 2 - 1;
	double hx = nx / (1 - nx * nx);

	return hx;
}

static double backward_transform(double hx, double low, double high) {
	double nx;
	hx *= 2;
	if (fabs(hx) > 1) {
		if (isinf(hx))
			nx = signbit(hx);
		else
			nx = (hypot(hx, 1) - 1) / hx;
	} else {
		hx = 1/hx;
		nx = hypot(1, hx) - hx;
	}
	nx = (nx + 1)/2;
	return nx*(high - low) + low;
}

static double move_by_bounded(double& x, double delta, double low, double high) {
	x = backward_transform(forward_transform(x, low, high) + delta, low, high);
}


struct Result {
	double energy;
	double d_energy;
};

struct MeanAndParam {
	double mean;
	Parameters params;
};

double denom_function(double* x, size_t dim, void* params_) {
	MeanAndParam& params = *static_cast<MeanAndParam*>(params_);
	Point p;
	memcpy(&p, x, sizeof p);
	double y = psi(p, params.params);
	return y*y;
}

double energy_function(double* x, size_t dim, void* params_) {
	MeanAndParam& params = *static_cast<MeanAndParam*>(params_);
	Point p;
	memcpy(&p, x, sizeof p);
	double y = psi(p, params.params);
	return y * (-0.5 * laplacian_psi(p, params.params) + (V(p) - params.mean)*y);
}

/*
double variance_function(double* x, size_t dim, void* params_) {
#if 0
	MeanAndParam& meanAndParam = *static_cast<MeanAndParam*>(params_);
	Point p(x);
	double y = psi(p, meanAndParam.params);
	double v = V(p) - meanAndParam.mean;
	double h1 = 0.25 * biharmonic_psi(p, meanAndParam.params);
	double h2 = -0.5 * laplacian_V(p) * y;
	double h3 = v * laplacian_psi(p, meanAndParam.params);
	double h4 = v * v * y;
	return y * (h1 + h2 + h3 + h4);
#endif
	double e = energy_function(x, dim, params_);
	return fabs(e);
}
*/

Result expected_energy(const Parameters& params, const Point& lower, const Point& upper, gsl_monte_vegas_state* vegas, gsl_rng* rng) {
	double denom, denom_err, ene, ene_err, var, var_err;
	
	MeanAndParam meanAndParam = {0.0, params};
	
	gsl_monte_function f;
	f.f = denom_function;
	f.dim = Point::Dimensions;
	f.params = &meanAndParam;
	
	double l[Point::Dimensions], u[Point::Dimensions];
	lower.fill(l);
	upper.fill(u);
	
//	gsl_monte_vegas_init(vegas);
	gsl_monte_vegas_integrate(&f, l, u, Point::Dimensions, N, rng, vegas, &denom, &denom_err);
	
	f.f = energy_function;
//	gsl_monte_vegas_init(vegas);
	gsl_monte_vegas_integrate(&f, l, u, Point::Dimensions, N, rng, vegas, &ene, &ene_err);
	
	Result res;
	denom = 1 / denom;
	res.energy = ene * denom;
	res.d_energy = hypot(ene_err * denom, denom_err * denom * denom);

/*
	meanAndParam.mean = res.energy;
	f.f = variance_function;
	gsl_monte_vegas_integrate(&f, l, u, Point::Dimensions, N, rng, vegas, &var, &var_err);	
	
	res.variance = var * denom;
*/
	return res;
}

struct EEParams {
	gsl_monte_vegas_state* vegas;
	gsl_rng* rng;
	const Parameters& params_low;
	const Parameters& params_high;
	const Point& lower;
	const Point& upper;
	double err;
};

double energy_extracter(const gsl_vector* x, void* p) {
	EEParams& eep = *static_cast<EEParams*>(p);
	Parameters params;
	params.read(x, eep.params_low, eep.params_high);
	Result res = expected_energy(params, eep.lower, eep.upper, eep.vegas, eep.rng);
	eep.err = res.d_energy;
	return res.energy;
}

void optimize(const Parameters& init_params, const Parameters& lower_params, const Parameters& upper_params,
	const Point& lower, const Point& upper, gsl_monte_vegas_state* vegas, gsl_rng* rng) {
		gsl_multimin_fminimizer* minimizer = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2, Parameters::Count);
		
		EEParams eep = {vegas, rng, lower_params, upper_params, lower, upper, 0};
		gsl_multimin_function f = {energy_extracter, Parameters::Count, &eep};
		
		gsl_vector* param = gsl_vector_alloc(Parameters::Count);
		init_params.write(param, lower_params, upper_params);
		
		/*
		gsl_vector* ss1 = gsl_vector_alloc(Parameters::Count);
		upper_params.write(ss2, lower_params, upper_params);
		lower_params.write(ss1, lower_params, upper_params);
		gsl_vector_sub(ss2, ss1);
		gsl_vector_free(ss1);
		gsl_vector_scale(ss2, 0.01);
		*/
		gsl_vector* ss2 = gsl_vector_alloc(Parameters::Count);
		gsl_vector_set_all(ss2, 1);
		
		gsl_multimin_fminimizer_set(minimizer, &f, param, ss2);
		gsl_vector_free(param);
		gsl_vector_free(ss2);
		
		double prev = 0;
		for (int i = 0; i < 25; ++ i) {
			gsl_multimin_fminimizer_iterate(minimizer);
			double x = gsl_multimin_fminimizer_minimum(minimizer);
			if (x != prev) {
				printf("%3d: %+.15f @ ", i, x);
				Parameters p;
				p.read(gsl_multimin_fminimizer_x(minimizer), lower_params, upper_params);
				p.print();
				printf("\n");
				prev = x;
			}
			double size = gsl_multimin_fminimizer_size(minimizer);
			if (gsl_multimin_test_size(size, min(eep.err, fabs(x))) == GSL_SUCCESS)
				break;
		}

		printf("end: %+.15f\n", gsl_multimin_fminimizer_minimum(minimizer));
		printf("   +- %.15f\n", eep.err);
		
		Parameters vfin;
		vfin.read(gsl_multimin_fminimizer_x(minimizer), lower_params, upper_params);
		printf(" @ ");
		vfin.print();
		printf("\n");
		
		gsl_multimin_fminimizer_free(minimizer);
}

/*
struct SimParams {
	gsl_monte_vegas_state* vegas;
	gsl_rng* rng;
	const Parameters& params_low;
	const Parameters& params_high;
	const Point& lower;
	const Point& upper;
	Parameters params;
	double err;
};

extern "C" {
	double sim_energy (void* xp) {
		SimParams& x = *static_cast<SimParams*>(xp);
		Result res = expected_energy(x.params, x.lower, x.upper, x.vegas, x.rng);
		x.err = res.d_energy;
		return res.energy;
	}
	void random_move (const gsl_rng* rng, void* xp, double step_size) {
		SimParams& x = *static_cast<SimParams*>(xp);
		x.params.random_walk(rng, step_size, x.params_low, x.params_high);
	}
	double config_dist (void *xp, void *yp) {
		SimParams& x = *static_cast<SimParams*>(xp), &y = *static_cast<SimParams*>(yp);
		return x.params.distance(y.params, x.params_low, x.params_high);
	}
	void params_print(void* xp) {
		SimParams& x = *static_cast<SimParams*>(xp);
		x.params.print();
	}
}

void optimize(const Parameters& init_params, const Parameters& lower_params, const Parameters& upper_params,
	const Point& lower, const Point& upper, gsl_monte_vegas_state* vegas, gsl_rng* rng) {
		SimParams init = {vegas, rng, lower_params, upper_params, lower, upper, init_params, 0};
		
		gsl_siman_params_t p = {200, 100, 10, 1, 0.05, 1.003, 1e-6};
		
		gsl_siman_solve(rng, &init, sim_energy, random_move, config_dist,
						params_print, NULL [copy-func], 
						NULL [copy-cons], NULL [destr],
						sizeof init, p);
		
		
}
*/

int main () {
	Parameters params_init = {{1, 1}};
	Parameters params_low = {{0, 0}};
	Parameters params_high = {{10, 10}};
	Point lower = {-100};
	Point upper = {100};
	
	gsl_rng* rng = gsl_rng_alloc(gsl_rng_taus2);
	gsl_rng_set(rng, time(NULL));
	gsl_monte_vegas_state* vegas = gsl_monte_vegas_alloc(Point::Dimensions);
		
	optimize(params_init, params_low, params_high, lower, upper, vegas, rng);
	
	
	gsl_monte_vegas_free(vegas);
	gsl_rng_free(rng);
	
	return 0;
}