Work in progress with Metropolis-Hastings sampling function

$CC -o zio metropolis_gibbs.c  ${MC_INSTDIR}/../../gr_code/clock_random.c -I ${HEADERS_GSL}   -I ${MC_INSTDIR}/../../include/ -I ${HEADERS_MPI} -I ${HEADERS_CFITSIO}  -L ${LIB_GSL}  -lm -lgsl -lgslcblas

$CC -o zio metropolis_gibbs.c  ${MC_INSTDIR}/../../gr_code/clock_random.c ${MC_INSTDIR}/../../gr_code/synchrotron/sync_distribution.c -I ${HEADERS_GSL}   -I ${MC_INSTDIR}/../../include/ -I ${HEADERS_MPI} -I ${HEADERS_CFITSIO}  -L ${LIB_GSL}  -lm -lgsl -lgslcblas

#include <math.h>

#include <time.h>

#include <local_header.h>

#include <synchrotron.h>

#define NUM_SAMPLES 400000

#define BURN_IN 1000  // Numero di iterazioni di "bruciatura"

#define INITIAL_X 0.1  // Valore iniziale di x (deve essere positivo)

#define INITIAL_Y 0.1  // Valore iniziale di y (deve essere positivo)

#define INITIAL_X 0.29  // Valore iniziale di x (deve essere positivo)

#define INITIAL_Y 0  // Valore iniziale di y (deve essere positivo)

#define MIN_X 0.01

#define MIN_X 0.1

#define MAX_X 10.0  // Massimo valore di x

#define MIN_Y -5 // Massimo valore di y

#define MAX_Y 5  // Massimo valore di y

double z_function(double x, double y);

double acceptance_ratio(double x_new, double x_old, double y_new, double y_old);

void metropolis_gibbs_sampling(double samples[NUM_SAMPLES][2]);

double z_function(double x, double y) 

{

void metropolis_gibbs_sampling(double* x, double* y, int Nburns);

    double f1;

    double f2;

    double bessel_arg;

    double value;

    bessel_arg=x / 2.0 * pow(1 + y * y, 1.5);

    f1=x *  pow(1 + y * y, 2)     * pow(gsl_sf_bessel_Knu(2.0/3.0, bessel_arg),2);

    f2=x *  pow(y,2) *(1 +y*y)    * pow(gsl_sf_bessel_Knu(1.0/3.0, bessel_arg),2);

    return f1+f2;  // Z(x, y) = x * BesselK^(2/3)(...)

}

/*====================================================================================*/

double acceptance_ratio(double x_new, double x_old, double y_new, double y_old) 

{

    double z_new;

    double z_old;

    double FlagValue;

    z_new = Psync_total(x_new, y_new);

    z_old = Psync_total(x_old, y_old);

    if (x_new <= MIN_X ||  x_new > MAX_X || x_old <= MIN_X ||  x_old > MAX_X 

    		|| y_new > MAX_Y || y_old > MAX_Y)

    if (z_new > z_old || clock_random() < z_new / z_old)

    {

    	return 0;

    	FlagValue=1;

    }

    else

    {

    	FlagValue=0;

    }

      z_new = z_function(x_new, y_new);

      z_old = z_function(x_old, y_old);

        if (z_new > z_old || clock_random() < z_new / z_old)

    return FlagValue;

            return 1.0;

        else

            return 0;	

}

}

/*====================================================================================*/

void metropolis_gibbs_sampling(double* x_sample, double* y_sample, int Nburns) {

void metropolis_gibbs_sampling(double samples[NUM_SAMPLES][2]) {

    // Inizializza il generatore di numeri casuali

    srand(time(NULL));

	int ii;

    double x = INITIAL_X;

    double y = INITIAL_Y;

    // Valori iniziali

    double x = INITIAL_X;  // Inizializza con un valore positivo

    double y = INITIAL_Y;  // Inizializza con un valore positivo

    double x_new;

    double y_new;

    double acceptance_x;

    double acceptance_y;

    for (int i = 0; i < NUM_SAMPLES; i++) 

    for (ii = 0; ii <= Nburns; ii++)

    {

    	/*===================================================*/

        // Campiona x utilizzando il Metropolis step

        double x_new;

    	/*===================================================*/

        do 

        {

            x_new = x + (clock_random() - 0.5);  // Propone un nuovo valore di x

        } while (x_new < MIN_X || x_new > MAX_X);  // Assicura che x sia positivo e rientri nel dominio

        }

        while (x_new < MIN_X || x_new > MAX_X);

        double acceptance_x = acceptance_ratio(x_new, x, y, y);

        acceptance_x = acceptance_ratio(x_new, x, y, y);

        if (clock_random() < acceptance_x)

        if (acceptance_x==1)

            x = x_new;

        /*===================================================*/

        // Campiona y utilizzando il Metropolis step

        double y_new;

        /*===================================================*/

        do 

        {

            y_new = y + (clock_random() - 0.5);  // Propose un nuovo valore di y

        } while (y_new < MIN_Y || y_new > MAX_Y);  // Assicura che y sia positivo e rientri nel dominio

        }

        while (y_new < MIN_Y || y_new > MAX_Y);

        double acceptance_y = acceptance_ratio(x, x_new, y_new, y);

        acceptance_y = acceptance_ratio(x, x_new, y_new, y);

        if (clock_random() < acceptance_y)

        if (acceptance_y==1)

            y = y_new;

        /*===================================================*/

        // Salva il campione dopo il periodo di "bruciatura"

        if (i >= BURN_IN) {

            samples[i - BURN_IN][0] = x;

            samples[i - BURN_IN][1] = y;

            printf("%lf %lf  %5.4e\n", x, y, z_function(x, y));

        /*===================================================*/

        if (ii == Nburns)

        {

        	*x_sample = x;

        	*y_sample = y;

        }

    }

}

int main() {

    double samples[NUM_SAMPLES][2];

    // Esegui il campionamento utilizzando Metropolis-Gibbs

    metropolis_gibbs_sampling(samples);

int main()

{

    // Stampa i campioni

    //printf("Campioni dalla distribuzione Z(x, y) = x * BesselK^(2/3)(...):\n");

	int ii;

	double energy=0;

	double csi=0;

  /*  for (int i = 0; i < NUM_SAMPLES; i++) 

	for (ii=0; ii< 150000; ii++)

	{

        printf("(%5.4e, %5.4e)\n", samples[i][0], samples[i][1]);

		metropolis_gibbs_sampling(&energy, &csi, 50);

		printf("%lf %lf  %5.4e\n", energy, csi, Psync_total(energy, csi));

	}

*/

    return 0;

}