#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"

i3_data_set * load_data_set(char * filename, i3_options * options){
	int n = options->stamp_size;
	
	i3_data_set * data_set = malloc(sizeof(i3_data_set));
	i3_image * image = i3_read_fits_image(filename);
	data_set->image = i3_image_copy_part(image,0,0,n,n);
	data_set->PSF=i3_image_create(n,n);
	for (int p=0;p<n;p++) data_set->PSF->data[p]=0.0;
	data_set->PSF->row[0][0]=1.0;
//	data_set->PSF_FT=i3_image_fft(data_set->PSF);
	data_set->noise=1.0e-3;
	i3_image_destroy(image);
	return data_set; 
	
}

i3_data_set * stack_file(char * filename, i3_options * options){
	int images_per_side=100;
	int n = options->stamp_size;
	if (n!=40){fprintf(stderr,"stamp size should = 40\n");exit(9);}

	i3_data_set * data_set = malloc(sizeof(i3_data_set));
	data_set->image = i3_image_create(n,n);
	i3_image * full_image = i3_read_fits_image(filename);

	for(int p=0;p<n*n;p++) data_set->image->data[p]=0;
	
	for(int x=0;x<images_per_side*n;x++){
		i3_flt * stack_row=data_set->image->row[x%n];
		i3_flt * row=full_image->row[x];
		for(int y=0;y<images_per_side*n;y++){
			stack_row[y%n]+=row[y];
		}		
	}
	for(int p=0;p<n*n;p++) data_set->image->data[p]/=(images_per_side*images_per_side);
	return data_set;	
}


i3_options * make_options(char * filename){
	i3_options * options = malloc(sizeof(i3_options));
	i3_options_set_defaults(options);
	i3_options_read(options, filename);
	return options;
	
}

i3_mcmc * make_mcmc(i3_model * model, i3_options * options, i3_parameter_set * width){
	i3_mcmc * mcmc = i3_mcmc_create(model, options);
	i3_model_copy_parameters(model,mcmc->width,width);
	return mcmc;
}

i3_minimizer * make_minimizer(i3_model * model, i3_options * options){
	i3_minimizer * m = i3_minimizer_create(model, options);
	free(m->step_size);
	m->step_size = i3_model_option_widths(model,options);
	return m;
}


int main(int argc, char * argv[]){
	i3_init_rng();
	char param_string[i3_max_string_length];
	
	if(argc<2){fprintf(stderr,"Usage: %s parameter_filename\n\n",argv[0]);exit(1);}
	i3_options * options = make_options(argv[1]);
	
	i3_model * model = i3_model_create(options->model_name,options);
		
	i3_parameter_set * initial_parameters = i3_model_option_starts(model,options);
	i3_model_parameter_pretty_string(model, initial_parameters, param_string, i3_max_string_length);
	printf("Start position =\n{%s}\n",param_string);
	
	
	
	i3_data_set * data_set = load_data_set("../great08/Stars/set0001.fit",options);
	i3_parameter_set * widths = i3_model_option_widths(model,options);

	i3_image * model_image = i3_image_create(options->stamp_size,options->stamp_size);

	i3_model_parameter_pretty_string(model, widths, param_string, i3_max_string_length);
	printf("Width =\n{%s}\n",param_string);

	i3_minimizer * minimizer = make_minimizer(model,options);
	i3_minimizer_result status = i3_minimizer_run(minimizer, initial_parameters,data_set);
	printf("Status = %d\n",status);
	i3_parameter_set * best_fit = i3_minimizer_result(minimizer);
	printf("Like = %.2f\n\n",i3_model_posterior(model,model_image,best_fit,data_set) );
	i3_model_parameter_pretty_string(model, best_fit, param_string, i3_max_string_length);
	
	printf("Best fit =\n{%s}\n",param_string);

	i3_parameter_set * fisher = i3_fisher_compute_diagonal(model, data_set, best_fit, widths);

	i3_model_parameter_pretty_string(model, fisher, param_string, i3_max_string_length);

	printf("Fisher diagonal = {\n%s}\n",param_string);
	
	i3_mcmc * mcmc = make_mcmc(model,options,fisher);
	i3_mcmc_run(mcmc,best_fit,data_set,options);
	i3_mcmc_save_results(mcmc,options->output_filename);
	
	int best_index = i3_mcmc_best_index(mcmc);
	i3_parameter_set * best_parameters = i3_mcmc_sample_number(mcmc,best_index);
	i3_model_parameter_pretty_string(model, best_parameters, param_string, i3_max_string_length);
	printf("New best fit = {\n%s}\n",param_string);
	
	printf("Best Like = %.2f\n",mcmc->likelihoods[best_index]);
	printf("Acceptance = %.1f%%\n\n",100*i3_mcmc_acceptance_rate(mcmc));
	
	
	
	return 0;
	
}