#include "i3_model.h"
#include "i3_model_tools.h"
#include <time.h>
#ifdef MPI
#include "mpi.h"
#endif

int main(int argc, char *argv[]){

// set up 

	i3_init_rng();
	FILE * file;
	i3_flt snr; 
	if(argc==1) snr = 20;
	else snr = atof(argv[1]);

	int n_pix = 29;	
//	i3_flt cut_off = 0.4;

	int n_grid_e = 71;
	int n_grid_r = 11;
	int n_grid_ee = n_grid_e*n_grid_e;
	int n_grid_eer = n_grid_ee*n_grid_r;

	char * models_id = "test4";
	
	// for a model with this id the chi2 will be calculated 
	int mid = 39441-1;

	int metric = 2;

// read models into memory

	printf("reading %d models \n",n_grid_eer);
	i3_image ** models_images = malloc(sizeof(i3_image*)*n_grid_eer);
	i3_flt * grid_eer_e1 = malloc(sizeof(i3_flt)*n_grid_eer);
	i3_flt * grid_eer_e2 = malloc(sizeof(i3_flt)*n_grid_eer);
	i3_flt * grid_eer_r = malloc(sizeof(i3_flt)*n_grid_eer);

	char filename_models_params[128]; snprintf(filename_models_params,128,"models/models_%s_params.dat",models_id);
	char filename_models_images[128]; snprintf(filename_models_images,128,"models/models_%s_images.dat",models_id);

	FILE * file_models_params = fopen(filename_models_params,"r"); if(!file_models_params) I3_FATAL("Models params file read error",1);
	FILE * file_models_images = fopen(filename_models_images,"r"); if(!file_models_images) I3_FATAL("Models images file read error",1);
	if(!file_models_images || !file_models_params ){  I3_FATAL("models files not found",1);};
	
	for(int i=0;i<n_grid_eer;i++){
		models_images[i] = i3_image_create(n_pix,n_pix);
		fread(models_images[i]->data,sizeof(i3_flt),n_pix*n_pix,file_models_images);
		fread(&grid_eer_e1[i],sizeof(i3_flt),1,file_models_params);
		fread(&grid_eer_e2[i],sizeof(i3_flt),1,file_models_params);
		fread(&grid_eer_r[i],sizeof(i3_flt),1,file_models_params);
	}
	fclose(file_models_params);
	fclose(file_models_images);

// create the grids
	
	i3_flt * grid_ee_e1 = malloc(sizeof(i3_flt)*n_grid_ee);
	i3_flt * grid_ee_e2 = malloc(sizeof(i3_flt)*n_grid_ee);
	
	for(int ge=0;ge<n_grid_ee;ge++){
		grid_ee_e1[ge] = grid_eer_e1[ge*n_grid_r];
		grid_ee_e2[ge] = grid_eer_e2[ge*n_grid_r];
		printf("%2.3f \t %2.2f \n",grid_ee_e1[ge],grid_ee_e2[ge] );
	}

	i3_flt * grid_e = malloc(sizeof(i3_flt)*n_grid_e);
	for(int ge=0;ge<n_grid_e;ge++){
		grid_e[ge] = grid_ee_e1[ge*n_grid_e];
		printf("%2.3f \t ", grid_e[ge] );	
	}

// choose one of them and calculate chi2
	
	printf("e1_true = %2.2f \t e2_true = %2.2f \t r_true = %2.4f \n", grid_eer_e1[mid],grid_eer_e2[mid]);

	i3_image * observed_image = i3_image_copy( models_images[mid] );

	i3_flt noise_std = i3_snr_to_noise_std( snr, observed_image );

	printf("noise_std = %2.10f\n", noise_std);
	
	i3_image_add_white_noise( observed_image, noise_std );

	printf("saving image\n"); 
	i3_image_save_text( observed_image, "i3_likelihood_test_obs_img.txt" ); 

	printf("calculating logL\n");
	i3_flt * logL;
	if(metric == 1) logL = i3_chi2_with_models_white_noise(models_images, n_grid_eer, observed_image, noise_std);
	if(metric == 2) logL = i3_logL_with_models_white_noise_marg_A(models_images, n_grid_eer, observed_image, noise_std);

	//logL = i3_logL_with_models_white_noise_fast(models_images, n_grid_eer, grid_eer_e1, grid_ee_e2, observed_image, grid_eer_e1[mid], grid_ee_e2[mid], cut_off ,noise_std,metric);
	
	printf("getting pdf from logL\n");
	i3_flt * pdf_full = i3_pdf_from_logL( logL, n_grid_eer );
	
	printf("marginalising r\n");
	i3_flt * pdf_marg = i3_marginalise_r( pdf_full, n_grid_eer , n_grid_ee );
	
	printf("saving\n");
	file = fopen("i3_likelihood_test_logL.dat","w");
	fwrite(logL,sizeof(i3_flt),n_grid_eer,file);
 	fclose(file);

	file = fopen("i3_likelihood_test_pdf_full.dat","w");
	fwrite(pdf_full,sizeof(i3_flt),n_grid_eer,file);
 	fclose(file);

	file = fopen("i3_likelihood_test_pdf_marg.dat","w");
	fwrite(pdf_marg,sizeof(i3_flt),n_grid_ee,file);
 	fclose(file);

	printf("model_id was %d\n",mid+1);
	printf("number of models was %d in %s\n",n_grid_eer,models_id);


//	printf("e_true = \t %f \t %f \t te_be =\t %f \t %f \t bias = \t %f \t %f \n",e1,e2,be_e1,be_e2,e1-be_e1,e2-be_e2);



	
	

	return 0;
}