#include "i3_mcmc.h"
#include "i3_load_data.h"
#include "i3_math.h"
#include "i3_model.h"
#include "i3_minimizer.h"
#include "i3_fisher.h"
#include "i3_mcmc.h"
#include "i3_image_fits.h"
#ifdef MPI
#include "mpi.h"
#endif


#include "i3_analyze.h"
#include "i3_model_tools.h"
#include "i3_data_set.h"
#include "i3_great.h"
#include "i3_psf.h"
#include "i3_image_fits.h"
#include "i3_utils.h"
#include "i3_options.h"
#include <time.h>



static void print_result(i3_model * model, i3_parameter_set * best, i3_flt like, i3_data_set * data_set)
{
	//Print the likelihood, followed by the model parameters, all on one line.
	if (like==BAD_LIKELIHOOD) like = 0.0/0.0;
	printf("% 1.4e\t",like);
	i3_model_printf_parameters(model, best);
	//printf("\n"); i3_model_pretty_printf_parameters(model,best);	printf("\n");
	printf ("\t%d",data_set->n_logL_evals);
	printf ("\t% 1.4e",data_set->signal_to_noise);
	printf ("\t% 1.4e",data_set->amplitude1);
	printf ("\t% 1.4e",data_set->amplitude2);
	printf ("\t% 1.4e",data_set->amplitude1/(data_set->amplitude2 + data_set->amplitude1));
	printf("\n");
	fflush(stdout);
}



int main(int argc, char * argv[]){
	
	//Initialize FFTW
	i3_init_rng_environment();
//		i3_init_rng();
	
	i3_fftw_load_wisdom();
	atexit(i3_fftw_save_wisdom);

	int N=48;
	
	i3_flt e1 = 0.60;
	i3_flt e2 = 0.05;
	i3_flt x0 = N/2.;
	i3_flt y0 = N/2.;
	i3_flt SN_scale = 5.e3;
	i3_flt bulge_amplitude = 1.0;
	i3_flt disc_amplitude = 0.5;
	i3_flt bulge_radius = 3.0;
	i3_flt disc_radius = 3.75;  // Matches fixed assumption

	i3_flt noise_sigma = 1.;	
	i3_image * image = i3_image_create(N,N);
	i3_image * model_image = i3_image_create(N,N);
	i3_image * weight = i3_image_create(N,N); 
	i3_image_fill(weight, 1. / noise_sigma / noise_sigma);


	// Build/read options
	// Use defaults and change if command line options specified
	i3_options * options = i3_options_default();

	options->stamp_size = N;
	options->minimizer_loops = 1;
	options->milestone_radius_ratio_max = 10.;

	// Read in the options if given on command line in file/multiple option form!
	if (argc>1){
		i3_options_read(argv[1], options);
	}

	if (argc>2){
		i3_options_read_command_line(options, argc, 2, argv);
	}


	//Create basic information and space for the chosen model (e.g. milestone)
	i3_model * model = i3_model_create("milestone", options);


	//Set up the basic PSF
	i3_moffat_psf psf_form;
	psf_form.beta = 3.;
	psf_form.e1 = 0.;
	psf_form.e2 = 0.;
	psf_form.fwhm = 4.;
        
	i3_data_set * dataset = i3_build_dataset(options, 0, image, weight, &psf_form, GREAT10_PSF_TYPE_MOFFAT);


	// Build the model parameters
	i3_milestone_parameter_set * params_gal = i3_model_option_starts(model, options);
	params_gal->radius = disc_radius;
	params_gal->radius_ratio = bulge_radius / disc_radius;
	params_gal->flux_ratio = bulge_amplitude / (bulge_amplitude + disc_amplitude);
	params_gal->bulge_index = 4.0;
	params_gal->disc_index = 1.0;
	params_gal->x0 = x0;
	params_gal->y0 = y0;
	params_gal->e1 = e1;
	params_gal->e2 = e2;

	
	i3_image_zero(image);
	// Build the image
	i3_milestone_model_image(params_gal, dataset, image);
	i3_image_scale(image, SN_scale);
	// Add noise
	i3_image_add_white_noise(image, noise_sigma);
	i3_image_save_fits(image, "optest_tmp.fits");

	// Analyse the test data
	// i3_analyze_dataset_method(dataset, model, options, stdout, 1, i3_minimizer_method_praxis);
	i3_flt like_test;
	i3_parameter_set * params_fitted;

	printf("input |e1| = %f\n", params_gal->e1);
	like_test = BAD_LIKELIHOOD;
	params_fitted = i3_analyze_dataset_method(dataset, model, options, &like_test, model_image, i3_minimizer_method_levmar);
	print_result(model, params_fitted, like_test, dataset);

	i3_image_zero(image);
	// Build the image
	i3_milestone_model_image(params_fitted, dataset, image);
        i3_image_save_fits(image, "optest_model.fits");
	exit(0);

	// Build the image again with a larger |e| each time
	i3_image_zero(image);
        dataset->n_logL_evals = 0;
	params_gal->e1 += 0.2;
	i3_milestone_model_image(params_gal, dataset, image);
	i3_image_scale(image, SN_scale / i3_image_sum(image));
	// Add noise
	i3_image_add_white_noise(image, noise_sigma);
	i3_image_save_fits(image, "optest_tmp.fits");

	printf("input |e1| = %f\n", params_gal->e1);	
	like_test = BAD_LIKELIHOOD;
	params_fitted = i3_analyze_dataset_method(dataset, model, options, &like_test, model_image, i3_minimizer_method_levmar);
	print_result(model, params_fitted, like_test, dataset);

	// Build the image again with a larger e
	i3_image_zero(image);
        dataset->n_logL_evals = 0;
	params_gal->e1 += 0.2;
	i3_milestone_model_image(params_gal, dataset, image);
	i3_image_scale(image, SN_scale);
	// Add noise
	i3_image_add_white_noise(image, noise_sigma);
	i3_image_save_fits(image, "optest_tmp.fits");

	printf("input |e1| = %f\n", params_gal->e1);	
	like_test = BAD_LIKELIHOOD;
	params_fitted = i3_analyze_dataset_method(dataset, model, options, &like_test, model_image, i3_minimizer_method_levmar);
	print_result(model, params_fitted, like_test, dataset);

	// Build the image again with a larger e
	i3_image_zero(image);
        dataset->n_logL_evals = 0;
	params_gal->e1 += 0.2;
	i3_milestone_model_image(params_gal, dataset, image);
	i3_image_scale(image, SN_scale);
	// Add noise
	i3_image_add_white_noise(image, noise_sigma);
	i3_image_save_fits(image, "optest_tmp.fits");




	printf("input |e1| = %f\n", params_gal->e1);	
	like_test = BAD_LIKELIHOOD;
	params_fitted = i3_analyze_dataset_method(dataset, model, options, &like_test, model_image, i3_minimizer_method_levmar);
	print_result(model, params_fitted, like_test, dataset);

}