#include "i3_quadratic_test.h"
#include "i3_mcmc.h"
#include "math.h"
#include "i3_math.h"
#include "i3_model.h"

int i3_quadratic_test_number_calls=0;


i3_flt i3_quadratic_test_likelihood(i3_image * model_image, i3_quadratic_test_parameter_set * parameters, i3_data_set * dataset){
	/* Compute the likelihood probability of the supplied parameter set */
	/* This is a placeholder to test the MCMC */
	i3_quadratic_test_number_calls++;
	i3_flt chi2 = 0.0;

	chi2 += i3_pow(parameters->x1 - QUADRATIC_TEST_TRUE_MEAN1, 2)/QUADRATIC_TEST_TRUE_VAR1;
	chi2 += i3_pow(parameters->x2 - QUADRATIC_TEST_TRUE_MEAN2, 2)/QUADRATIC_TEST_TRUE_VAR2;
	chi2 += i3_pow(parameters->x3 - QUADRATIC_TEST_TRUE_MEAN3, 2)/QUADRATIC_TEST_TRUE_VAR3;
	chi2 += i3_pow(parameters->x4 - QUADRATIC_TEST_TRUE_MEAN4, 2)/QUADRATIC_TEST_TRUE_VAR4;
	return -chi2/2;
}

void i3_quadratic_test_propose(i3_mcmc * mcmc, i3_quadratic_test_parameter_set * current, i3_quadratic_test_parameter_set * proposed){

	/* Turn the parameter set into a vector*/
	i3_flt * c = (i3_flt*) current;
	i3_flt * p = (i3_flt*) proposed;
	i3_flt * w = (i3_flt*) mcmc->width;
	int np = mcmc->nParam;
	int i = mcmc->nSamples%np;
	/* If necessary, create a new random rotation */
	if(i==0){
		i3_random_rotation(mcmc->random_rotation,mcmc->nParam);
	}
	/* Jump in the random direction by a random normal distance*/
	i3_flt lambda = mcmc->scaling_parameter * i3_random_normal() / i3_sqrt(mcmc->nParam);
	for(int b=0;b<np;b++) p[b] = c[b] + lambda*w[b]*mcmc->random_rotation[i+np*b];
}

void i3_quadratic_test_tune_proposal(i3_mcmc * mcmc){
	
}


int i3_quadratic_test_parameter_string(i3_quadratic_test_parameter_set * parameter,char * parameter_string,int string_length){
	return snprintf(parameter_string,string_length,"%e %e %e %e",parameter->x1,parameter->x2,parameter->x3,parameter->x4);
}

int i3_quadratic_test_parameter_pretty_string(i3_quadratic_test_parameter_set * parameter,char * parameter_string,int string_length){
	return snprintf(parameter_string,string_length,"x1 = %e\nx2 = %e\nx3 = %e\nx4 = %e\n",parameter->x1,parameter->x2,parameter->x3,parameter->x4);
}