#include "i3_image.h"
#include "i3_model_tools.h"
#include "i3_mcmc.h"
#include "i3_minimizer.h"
#include "i3_psf.h"

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;
	
}


void log_callback(i3_mcmc * mcmc ,i3_data_set * data_set,i3_options *options){
	i3_flt L = i3_mcmc_last_like(mcmc);
	fprintf(stdout,"%.12e\t",L);
	i3_model_fprintf_parameters(mcmc->model,stdout,i3_mcmc_last_sample(mcmc));
	fprintf(stdout,"\n");
}


void adapt_log(i3_mcmc * mcmc, i3_data_set * data_set, i3_options * options)
{	
	i3_mcmc_adaptive_callback(mcmc,data_set,options);
	log_callback(mcmc,data_set,options);
}


i3_mcmc * make_mcmc(i3_model * model, i3_options * options, i3_parameter_set * width){
	i3_mcmc * mcmc = i3_mcmc_create(model, options);
	if (width!=NULL){
		i3_model_copy_parameters(model,mcmc->width,width);
	}
	else{
		mcmc->width = i3_model_option_widths(model,options);
	}
	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;
}

void report_and_save_model(i3_model * model, int npix, i3_parameter_set * p, i3_data_set * data_set, char * name){
	i3_image * model_image = i3_image_create(npix,npix);
	i3_flt like = i3_model_posterior(model,model_image,p,data_set);
	fprintf(stderr,"parameters %s: {\n",name);
	i3_model_pretty_fprintf_parameters(model,stderr,p);
	fprintf(stderr,"}\n");
	fprintf(stderr,"Like = %f\n\n",like);
	char filename[i3_max_string_length];
	snprintf(filename,i3_max_string_length,"%s_image.txt",name);
	i3_image_save_text(model_image,filename);
	i3_image_destroy(model_image);
	
	if (like==BAD_LIKELIHOOD||like!=like){
		fprintf(stderr,"BAD LIKE at %s :\n",name);
//		i3_model_fprintf_prior_violations(model, p, stderr);
		exit(1);
	}
}



void milestone_process_image(i3_image * postage_stamp, i3_fourier * kernel, i3_flt noise, i3_model * model,i3_options * options)
{


	/* Create a data set from the image.  Use a fake noise level. */
	i3_data_set * data_set = malloc(sizeof(i3_data_set));
//	int n = options->stamp_size;
	data_set->image = postage_stamp;
	data_set->psf_downsampler_kernel = kernel;
	data_set->upsampling=options->upsampling;

	data_set->noise=noise;
//	i3_parameter_set * width = i3_model_option_widths(model, options);

	i3_parameter_set * estimated_starting_parameters;
	estimated_starting_parameters = i3_model_option_starts(model,options);
	i3_model_starting_position(model, estimated_starting_parameters, data_set);

//	report_and_save_model(model,options->stamp_size,estimated_starting_parameters,data_set,"estimated_start");

	i3_minimizer * minimizer = make_minimizer(model,options);
	i3_parameter_set * best_fit_minimizer = i3_minimizer_run_praxis(minimizer,estimated_starting_parameters,data_set);
//	report_and_save_model(model,options->stamp_size,best_fit_minimizer,data_set,"minimzed");

/*	i3_mcmc * mcmc = make_mcmc(model,options,width);
	i3_mcmc_run_callback(mcmc,best_fit_minimizer,data_set,options,adapt_log);
	i3_parameter_set * best_mcmc = i3_mcmc_best_sample(mcmc);
	if (best_mcmc!=NULL){
		report_and_save_model(model,options->stamp_size,best_mcmc,data_set,"mcmc");
		i3_model_fprintf_parameters(model,stderr,best_mcmc);
		fprintf(stderr,"\n");
	}
	else{
		fprintf(stderr,"-------------\n");
		
	}
*/	
	i3_milestone_parameter_set * q = (i3_milestone_parameter_set*)best_fit_minimizer;
	printf("%lf\t%lf\n",q->e1,q->e2);
	
	free(best_fit_minimizer);
//	i3_mcmc_destroy(mcmc);
	i3_minimizer_destroy(minimizer);
	free(data_set);
}


i3_image * generate_milestone_image(int N, i3_flt e1, i3_flt e2, int upsampling, i3_fourier * kernel){
	i3_flt x0 = N/2.;
	i3_flt y0 = N/2.;
	// i3_flt bulge_amplitude = 1.0;
	i3_flt disc_amplitude = 2.0;
	i3_flt bulge_radius = 2.0;
//	i3_flt disc_radius = MILESTONE_BULGE_DISC_SCALE*bulge_radius;
	
//	i3_flt bulge_sersic_index = 4.0;
//	i3_flt disc_sersic_index = 1.0;
	
	i3_milestone_parameter_set p;
	i3_data_set data_set;
	
//	p.B=disc_amplitude    ;
	p.x0=x0   ;
	p.y0=y0   ;
	p.e1=e1   ;
	p.e2=e2  ;
	p.radius=bulge_radius	  ;
	
	i3_image * image = i3_image_create(N,N);
	i3_image * empty_image = i3_image_create(N,N);
	data_set.image=empty_image;
	data_set.noise=1.0;
	data_set.upsampling=upsampling;
	data_set.psf_downsampler_kernel = kernel;
	
	i3_milestone_likelihood(image, &p, &data_set);

	i3_image_destroy(empty_image);

	i3_image_moments mom;
	i3_image_compute_moments(image,&mom);
//	fprintf(stderr,"e1,e2 = %f, %f\tx0,y0 = %f, %f\t sum = %f,  qxx,qxy,qyy = %f,%f,%f\n",mom.e1,mom.e2,mom.x0,mom.y0,mom.sum,mom.qxx,mom.qxy,mom.qyy);
	return image;
	
}


i3_fourier * generate_downsampler(int N, int upsampling){
	int M = N+1;
//	int H = M*upsampling;
	i3_fourier * downsampling_kernel = i3_fourier_downsampling_psf(M, M, upsampling, M);
	return downsampling_kernel;
}



int main(int argc, char * argv[]){
	i3_fftw_load_wisdom();
	atexit(i3_fftw_save_wisdom);
	i3_init_rng();
	
	if (argc==1){
		fprintf(stderr,"Syntax: %s [parameter_filename]\n\n",argv[0]);
		exit(1);
		
	}
	i3_options * options = make_options(argv[1]);


	int upsampling = options->upsampling;
	int N=options->stamp_size;
	// i3_image * high_res_image = generate_test_image((N+1)*upsampling,e1,e2);
	// printf("Have image\n");
	

	i3_fourier * kernel = generate_downsampler(N,upsampling);


	
		
		
	i3_model * model = i3_model_create(options->model_name, options);
	if (!model) I3_FATAL("Model not found",1);
	
	
	for (int i=0;i<100;i++){
		i3_flt e1 = 1.1, e2=1.1;
		while(e1*e1+e2*e2>1){
			e1 = i3_gsl_random_normal_float() * 0.2;
			e2 = i3_gsl_random_normal_float() * 0.2;
		}
		i3_image * image  = generate_milestone_image(N, e1, e2, upsampling, kernel);
//		i3_image_save_text(image,"target.txt");
		i3_flt noise = 1e-3;
		printf("%f\t%f\t",e1,e2);
		milestone_process_image(image, kernel, noise, model, options);
		i3_image_destroy(image);
	}
}