#include "i3_model.h"
#include "i3_model_tools.h"
#include "i3_load_data.h"
#include "i3_data_set.h"
#include "i3_great.h"
#include "i3_mcmc.h"
#include "i3_math.h"
#include "i3_minimizer.h"
#include "i3_fisher.h"
#include "i3_psf.h"
#include "i3_image_fits.h"
#include "i3_model_tools.h"
#include "i3_utils.h"
#include "i3_options.h"
#include "i3_analyze.h"
#include <time.h>
#include "extern/levmar/levmar.h"
#define BEERMAT_MIN_R 0.2    // THIS VALUE SHOULD BE TESTED/OPTIMIZED
#include <stdlib.h>

#ifdef I3_USE_DOUBLE
#define levmar_bc_dif dlevmar_bc_dif
#define levmar_bc_der dlevmar_bc_der
#else
#define levmar_bc_dif slevmar_bc_dif
#define levmar_bc_der slevmar_bc_der
#endif

#define NBC_IMAGE_INFO_SZ 9

i3_flt measure_entropy(i3_image * image){

	
	i3_flt entropy = 0;

	// normalise
	i3_flt Z = i3_array_normalise( image->data, image->n ); 

	for(int i=0;i<image->n;i++){
		
		double l = image->data[i] * log(image->data[i]);
		if(isnan(l)) l = 0;
		if(isinf(l)) l = 0;
		entropy -= l;
	} 
		
	// scale back
	i3_image_scale( image, Z );

	return entropy;
	 

}

i3_flt measure_kldivergence(i3_image * image1,i3_image * image2){

	
	double kl = 0;
	
	// normalise
	
	i3_flt Z1 = i3_array_normalise( image1->data, image1->n );
	i3_flt Z2 = i3_array_normalise( image2->data, image2->n ); 

	if(image1->n != image2->n) I3_FATAL("Image sizes don't match",1);

	for(int i=0;i<image1->n;i++){
		
		double l = image1->data[i] * log( (double)image1->data[i] / (double)image2->data[i] );
		if(isnan(l)) l = 0;
		if(isinf(l)) l = 0;
		kl += l;
	} 

	i3_image_scale( image1, Z1 );
	i3_image_scale( image2, Z2 );
		
	return (i3_flt)kl;
	 

}

void minimizer_wrapper_levmar(i3_flt *p, i3_flt *hx, int m, int n, void *adata){
	
	i3_data_set * dataset = (i3_data_set*)adata;
	int n_pix = dataset->stamp_size;
	i3_image * model_image = i3_image_create(n_pix,n_pix);
	
	i3_sersics_parameter_set * params = malloc(sizeof(i3_sersics_parameter_set));
	memcpy( params,dataset->start_params,sizeof(i3_sersics_parameter_set));

	//printf("wrapper org bi di %f %f \n",dataset->start_params->bulge_index,dataset->start_params->disc_index);
	//printf("wrapper cpy bi di %f %f \n",params->bulge_index,params->disc_index);


	params->x0=p[0];
	params->y0=p[1];
	params->e1=p[2];
	params->e2=p[3];
	params->radius=p[4];
	params->bulge_A=p[5];
	params->disc_A=p[6];
	
	i3_sersics_likelihood(model_image,params,dataset);
	
	i3_array_copy_into(hx,model_image->data,n);
	
	dataset->signal_to_noise = i3_signal_to_noise_model_fitted(dataset->image,model_image);
	//dataset->signal_to_noise = i3_signal_to_noise_model_fitted(dataset->image,model_image);
	
	i3_image_destroy(model_image);
	free(params);		
}

// Following finite difference function taken from levmar - misc_core.c
static void levmar_fdif_cent_jac_approx(void (*func)(i3_flt *p, i3_flt *hx, int m, int n, void *adata), i3_flt *p, i3_flt *hxm, i3_flt *hxp, i3_flt delta, i3_flt *jac, int m, int n, void *adata){

register int i, j;
i3_flt tmp;
register i3_flt d;

 i3_data_set * dataset = (i3_data_set*)adata;
 int n_pix = dataset->stamp_size;
 i3_image * galaxy = i3_image_create(n_pix,n_pix);
 i3_image * jacobian = i3_image_create(n_pix,n_pix);


  for(j=0; j<m; ++j){
    // MH: commented this out from the original finite difference code as done in levmar/misc_core.c
    /* determine d=max(1E-04*|p[j]|, delta), see HZ */

    d=1E-04*p[j]; // force evaluation
    d=FABS(d);
    if(d<delta)
      d=delta;

    //d=delta;

    tmp=p[j];
    //p[j]-=d;
    (*func)(p, hxm, m, n, adata);
    
    p[j]=tmp+d;
    (*func)(p, hxp, m, n, adata);
    p[j]=tmp; /* restore */

    //d=0.5/d; /* invert so that divisions can be carried out faster as multiplications */
    d=1.0/d; /* invert so that divisions can be carried out faster as multiplications */
    for(i=0; i<n; ++i){
      jac[i*m+j]=(hxp[i]-hxm[i])*d;    
    }
  }
}

// Following finite difference function taken from levmar - misc_core.c
static void levmar_fdif_cent_jac_approx_j(void (*func)(i3_flt *p, i3_flt *hx, int m, int n, void *adata), i3_flt *p, i3_flt *hxm, i3_flt *hxp, i3_flt delta, i3_flt *jac, int m, int n, int j, void *adata){

  register int i;
  i3_flt tmp;
  register i3_flt d;

 i3_data_set * dataset = (i3_data_set*)adata;
 int n_pix = dataset->stamp_size;
 i3_image * galaxy = i3_image_create(n_pix,n_pix);
 i3_image * jacobian = i3_image_create(n_pix,n_pix);

 d=1E-04*p[j]; // force evaluation
 d=FABS(d);
 if(d<delta)
   d=delta;

 //d=delta;
 
 tmp=p[j];
 //p[j]-=d;
 (*func)(p, hxm, m, n, adata);

 p[j]=tmp+d;
 (*func)(p, hxp, m, n, adata);
 p[j]=tmp; /* restore */
      
 //d=0.5/d; /* invert so that divisions can be carried out faster as multiplications */
 d=1.0/d; /* invert so that divisions can be carried out faster as multiplications */
 for(i=0; i<n; ++i){
   jac[i*m+j]=(hxp[i]-hxm[i])*d;
 }
}

static void minimizer_wrapper_levmar_jacobian(i3_flt *p, i3_flt *jac, int m, int n, void *adata){

	i3_data_set * dataset = (i3_data_set*)adata;
	int n_pix = dataset->stamp_size;
	
	i3_sersics_parameter_set * params = malloc(sizeof(i3_sersics_parameter_set));
	memcpy( params,dataset->start_params,sizeof(i3_sersics_parameter_set));

	params->x0=p[0];
	params->y0=p[1];
	params->e1=p[2];
	params->e2=p[3];
	params->radius=p[4];
	params->bulge_A=p[5];
	params->disc_A=p[6];

	// prepare the beermatted param set
	i3_sersics_parameter_set * paramset_beermat = i3_sersics_beermat_params(params);
        
	int sanity_check = 0;
	if(sanity_check){	  
	    // Sanity check: compute Jacobian with finite differences, should give same results like 
	    // calling levmar without Jacobian function
	    i3_flt delta = 1e-6;
	    i3_flt * hxm = malloc(sizeof(i3_flt)*n);
	    i3_flt * hxp = malloc(sizeof(i3_flt)*n);

	    levmar_fdif_cent_jac_approx(minimizer_wrapper_levmar, p, hxm, hxp, delta, jac, m, n, adata); 

	    free(hxm);
	    free(hxp);
	  }
	else{
	  // get the jacobian of the model image                
	  i3_sersics_model_jacobian(paramset_beermat, dataset, jac);

	  /*
	  i3_flt delta = 1e-6;
	  i3_flt * hxm = malloc(sizeof(i3_flt)*n);
	  i3_flt * hxp = malloc(sizeof(i3_flt)*n);
	  
	  levmar_fdif_cent_jac_approx_j(minimizer_wrapper_levmar, p, hxm, hxp, delta, jac, m, n, 2, adata); 
	  levmar_fdif_cent_jac_approx_j(minimizer_wrapper_levmar, p, hxm, hxp, delta, jac, m, n, 3, adata); 
	  
	  free(hxm);
	  free(hxp);
	  */
	}

	//Free memory
	free(params);		

}


// Function to check Jacobian with finite difference approximation
void check_jacobian(
    void (*func)(i3_flt *p, i3_flt *hx, int m, int n, void *adata),
    void (*jacf)(i3_flt *p, i3_flt *j,  int m, int n, void *adata),
    i3_flt *p, int m, int n, void *adata, i3_flt delta, i3_flt *err){

  i3_data_set * dataset = (i3_data_set*)adata;
  int n_pix = dataset->stamp_size;

  // Allocate memory
  i3_flt * jac_approx = malloc(sizeof(i3_flt)*m*n);
  i3_flt * jac_exact  = malloc(sizeof(i3_flt)*m*n);

  i3_flt * hxm = malloc(sizeof(i3_flt)*n);
  i3_flt * hxp = malloc(sizeof(i3_flt)*n);

  // set values for e1 and e2 manually to determine breakdown points

  //p[2] = 0.20000009999999;
  //p[3] = +0.1;
  //p[2] = -0.2000001;
  //p[3] = -0.1;
  // e1=0.265616 and e2=-0.176652
  // e1=0.323314 and e2=0.381255
  // e1=-0.068910 and e2=-0.076020
  // e1=-0.002733 and e2=0.011096
  //p[2] = -0.002733;//0.1*(i3_random_uniform()-0.5)*2;
  //p[3] = 0.011096;//0.1*(i3_random_uniform()-0.5)*2;
  p[2] = (i3_random_uniform()-0.5)*2;
  p[3] = (i3_random_uniform()-0.5)*2;
  printf("\n e1=%f and e2=%f\n",p[2],p[3]);

  // Compute finite difference approximation	
  levmar_fdif_cent_jac_approx(minimizer_wrapper_levmar, p, hxm, hxp, delta, jac_approx, m, n, adata); 

  // Compute analytic jacobain
  minimizer_wrapper_levmar_jacobian(p, jac_exact, m, n, adata);
  
  // =====================================================================
  // Set variable save to true if you want to save Jacobian images to disk
  int save = true;
  if (save == 1){
    
      i3_image * approx = i3_image_create( n_pix, n_pix ); i3_image_zero( approx );
      i3_image * exact  = i3_image_create( n_pix, n_pix ); i3_image_zero( exact  );

      char fname_approx[20];
      char fname_exact[20];

      register int i;
      register int j;
      register int k;
      register int prm; // choose which component to plot
      for (prm=0; prm<m; ++prm){
	k = 0;
	for (i=0; i<n_pix; ++i){
	  for (j=0; j<n_pix; ++j){
	    approx->row[i][j] = jac_approx[k*m+prm];
	    exact->row[i][j]  = jac_exact[k*m+prm];
	    k += 1;
	  }
	}
	sprintf(fname_approx, "jac_approx_%d.fits",prm);
	sprintf(fname_exact, "jac_exact_%d.fits",prm );
	i3_image_save_fits( approx, fname_approx);
	i3_image_save_fits( exact, fname_exact);
      }

      i3_image_destroy(approx);
      i3_image_destroy(exact);
  }
  // =====================================================================
  // Compute difference between numerical and analytical Jacobians
  register int i;
  for(i=0; i<n*m; ++i){
    err[i]=jac_approx[i] - jac_exact[i];
  }	
  
  // Free memory 
  free(jac_approx);
  free(jac_exact);
  free(hxm);
  free(hxp);
}

i3_sersics_parameter_set * minimizer_run_levmar(i3_data_set * dataset, i3_options * options){
	
        i3_flt levmar_options[LM_OPTS_SZ];
        levmar_options[0]= dataset->options->levmar_LM_INIT_MU*LM_INIT_MU;
	levmar_options[1]= dataset->options->levmar_eps1;
        levmar_options[2]= dataset->options->levmar_eps2;
	levmar_options[3]= dataset->options->levmar_eps3;
	levmar_options[4]= dataset->options->levmar_eps4;
        levmar_options[5]= dataset->options->levmar_tau;
	levmar_options[6]= dataset->options->minimizer_verbosity;
	
	i3_flt levmar_info[LM_INFO_SZ];
		
	int n_params = 7;
	int n_pix = dataset->image->n;
	
	i3_flt * p = malloc(sizeof(i3_flt)*n_params);
	i3_sersics_parameter_set * params_start = (i3_sersics_parameter_set*) dataset->start_params;
	p[0] = params_start->x0;      // 10*i3_random_uniform();//params_start->x0;
	p[1] = params_start->y0;      // 10*i3_random_uniform();//params_start->y0;
	p[2] = params_start->e1;      // 10*i3_random_uniform();//params_start->e1;
	p[3] = params_start->e2;      // 10*i3_random_uniform();//params_start->e2;
	p[4] = params_start->radius;  // 10*i3_random_uniform();//params_start->radius;
	p[5] = params_start->bulge_A; // 10*i3_random_uniform();//params_start->bulge_A;
	p[6] = params_start->disc_A;  // 10*i3_random_uniform();//params_start->disc_A;
	
        i3_flt * covmat = malloc(n_params*n_params*sizeof(i3_flt));

	//i3_image * model_image = i3_image_like( dataset->image );
	
	int check_jac = true;
	if(check_jac)
	{
	  // MH: Check computation of Jacobian with levmar builtin function
	  printf("========================================\n");
      	  printf("Checking Jacobian\n");
     	  i3_flt * err = malloc(sizeof(i3_flt)*n_params*n_pix);
     	  i3_flt delta = 1e-6;
     	  check_jacobian(minimizer_wrapper_levmar, minimizer_wrapper_levmar_jacobian, p, n_params, n_pix, dataset, delta, err);	
     
     	  i3_flt err_tot=0.0;
     	  i3_flt errs[n_params];
     	  register int i,j;
     	   
     	  for(j=0; j<n_params; ++j){
     	    errs[j] = 0.0;
     	    for(i=0; i<n_pix; ++i){
	      // Compute max. error
	      if(i3_fabs(err[i*n_params+j] - errs[j]) > 0){
		errs[j] = i3_fabs(err[i*n_params+j]);
	      }
	      // Compute total error
	      err_tot += i3_fabs(err[i*n_params+j]);
     	    }
     	  }
     	   
     	  printf("Max error in x0  	 %.10f\n", errs[0]); 
	  printf("Max error in y0      	 %.10f\n", errs[1]); 
	  printf("Max error in e1      	 %.10f\n", errs[2]); 
	  printf("Max error in e2      	 %.10f\n", errs[3]); 
	  printf("Max error in radius 	 %.10f\n", errs[4]); 
	  printf("Max error in bulge_A 	 %.10f\n", errs[5]); 
	  printf("Max error in disc_A  	 %.10f\n", errs[6]); 
            
     	  printf("----------------------------------------\n");
     	  printf("Total error %f\n", err_tot);
     	  printf("========================================\n");
     	   
     	  exit(0);
     	  return 0;
	}

	if(dataset->options->levmar_use_analytic_gradients){
	  levmar_bc_der(minimizer_wrapper_levmar, minimizer_wrapper_levmar_jacobian, p, dataset->image->data, n_params, n_pix, NULL, NULL, options->minimizer_max_iterations, levmar_options, levmar_info, NULL, covmat, dataset);
	}
	else{
	  levmar_bc_dif(minimizer_wrapper_levmar, p, dataset->image->data, n_params, n_pix, NULL, NULL, options->minimizer_max_iterations, levmar_options, levmar_info, NULL, covmat, dataset);
	}


	
	i3_sersics_parameter_set * params_ml = malloc(sizeof(i3_sersics_parameter_set));
	memcpy( params_ml,dataset->start_params,sizeof(i3_sersics_parameter_set));

	params_ml->x0 = p[0];
	params_ml->y0 = p[1];
	params_ml->e1 = p[2];
	params_ml->e2 = p[3];
	params_ml->radius = p[4];
	params_ml->bulge_A = p[5];
	params_ml->disc_A = p[6];
	
	i3_array_copy_into(dataset->levmar_info,levmar_info,LM_INFO_SZ);
	
	
	
	return params_ml;
	
	
}


void set_start_params(i3_model * model, i3_data_set * dataset, i3_options * options, i3_sersics_parameter_set * start, i3_sersics_parameter_set * truth){

	i3_model_copy_parameters(model, start, truth);
	//i3_sersics_start(start,dataset,options);

	//start->e1 = 0.0;
        //start->e2 = 0.0;

	start->e1 = truth->e1;
        start->e2 = truth->e2;;
		
        //if(truth->e1 > 0.3) start->e1 = 0.3;
	//if(truth->e2 > 0.3) start->e2 = 0.3;
	
	//if(truth->e1 > 0.6) start->e1 = 0.6;
	//if(truth->e2 > 0.6) start->e2 = 0.6;

	//if(truth->e1 < -0.3) start->e1 = -0.3;
	//if(truth->e2 < -0.3) start->e2 = -0.3;
	
	//if(truth->e1 < -0.6) start->e1 = -0.6;
	//if(truth->e2 < -0.6) start->e2 = -0.6;
	
		
	start->bulge_A = truth->bulge_A;
        start->disc_A  = truth->disc_A;
	start->x0 = truth->x0;
	start->y0 = truth->y0;
	start->radius = truth->radius;
	
	i3_flt change_r  = 0.8;
	i3_flt change_xy0 = 0.1;
	i3_flt change_A = 0.8;
		
        start->radius *=	(1.+	change_r	*(i3_random_uniform()-0.));
        start->x0 *= 		(1.+	change_xy0	*(i3_random_uniform()-0.5));
        start->y0 *= 		(1.+	change_xy0	*(i3_random_uniform()-0.5));
        start->bulge_A *=	(1.+	change_A	*(i3_random_uniform()-0.5));
        start->disc_A  *= 	(1.+	change_A	*(i3_random_uniform()-0.5));
	
   
        
        }
        
static void save_result(i3_model * model, i3_sersics_parameter_set * hp, FILE * output_file, int identifier, i3_data_set * dataset, i3_flt * image_info){

	
	i3_flt r1 = hp->radius;
	i3_flt r2 = hp->radius/hp->radius_ratio;

	//printf("bi di %f %f \n",hp->bulge_index,hp->disc_index);

	i3_flt f1 = i3_sersic_total_flux(hp->bulge_index,hp->bulge_A,r1*r1);
	i3_flt f2 = i3_sersic_total_flux(hp->disc_index ,hp->disc_A ,r2*r2);
        i3_flt fr = f1/(f1+f2); 

        fprintf(output_file, "%d",identifier);  
        fprintf(output_file, "\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e",hp->x0,hp->y0,hp->e1,hp->e2,hp->radius,hp->bulge_A,hp->disc_A); ; 
	fprintf(output_file, "\t% 2.3e\t% 2.3e",f1,f2); 
	for(int i=0;i<NBC_IMAGE_INFO_SZ;i++) fprintf(output_file, "\t% 2.4e",image_info[i]); 	
	fprintf(output_file, "\t% 2.2e",dataset->signal_to_noise);
        fprintf(output_file, "\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e",dataset->levmar_info[0],dataset->levmar_info[1],dataset->levmar_info[2],dataset->levmar_info[3],dataset->levmar_info[4],dataset->levmar_info[5]);
	fprintf(output_file, "\t% d\t% d\t% d\t% d\t% d",(int)dataset->levmar_info[6],(int)dataset->levmar_info[7],(int)dataset->levmar_info[8],(int)dataset->levmar_info[9],(int)dataset->levmar_info[10]);
        fprintf(output_file, "\n");
        fflush(output_file);
	
}

static void save_starts(i3_model * model, i3_sersics_parameter_set * hp, FILE * output_file, int identifier){

	
	i3_flt r1 = hp->radius;
	i3_flt r2 = hp->radius/hp->radius_ratio;

	//printf("bi di %f %f \n",hp->bulge_index,hp->disc_index);

	i3_flt f1 = i3_sersic_total_flux(hp->bulge_index,hp->bulge_A,r1*r1);
	i3_flt f2 = i3_sersic_total_flux(hp->disc_index ,hp->disc_A ,r2*r2);
        i3_flt fr = f1/(f1+f2); 

        fprintf(output_file, "%d",identifier);  
        fprintf(output_file, "\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e",hp->x0,hp->y0,hp->e1,hp->e2,hp->radius,hp->bulge_A,hp->disc_A); ; 
	fprintf(output_file, "\t% 2.3e\t% 2.3e",f1,f2); 
	fprintf(output_file, "\n");
        fflush(output_file);
	
}

static void save_truth(i3_model * model, i3_sersics_parameter_set * hp, FILE * output_file, int identifier){

	
	i3_flt r1 = hp->radius;
	i3_flt r2 = hp->radius/hp->radius_ratio;

	//printf("bi di %f %f \n",hp->bulge_index,hp->disc_index);

	i3_flt f1 = i3_sersic_total_flux(hp->bulge_index,hp->bulge_A,r1*r1);
	i3_flt f2 = i3_sersic_total_flux(hp->disc_index ,hp->disc_A ,r2*r2);
        i3_flt fr = f1/(f1+f2); 

        fprintf(output_file, "%d",identifier);  
        fprintf(output_file, "\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e",hp->x0,hp->y0,hp->e1,hp->e2,hp->radius,hp->bulge_A,hp->disc_A); ; 
	fprintf(output_file, "\t% 2.3e\t% 2.3e",f1,f2); 
	fprintf(output_file, "\n");
        fflush(output_file);
	
}

i3_flt measure_fwxm(i3_image * image, i3_flt fwxm){

	// this will be a function soon
	int nx = image->nx;
	i3_flt xy0 = nx/2;
	i3_flt f1 = 0;
	int x1 = 0;
	i3_flt f2 = 0;
	int x2 = 0;

	int s = nx/2+1;

	i3_flt * profile = malloc(sizeof(i3_flt)*s);
	i3_flt * diff = malloc(sizeof(i3_flt)*s);
	for(int i=0;i<s;i++) profile[i] = image->row[nx/2][i];	

	int max_ind = nx/2;
	i3_flt max_val = profile[max_ind];
	i3_flt f3 = max_val*fwxm;

	for(int i=0;i<nx/2;i++)	diff[i] = i3_fabs(profile[i] - f3); 

	//printf("max_val fwxm f3 %2.4f %2.4f %2.4f \n",max_val,fwxm,f3);

	x1=i3_array_min_index( diff, nx/2 );
	f1=profile[x1];
	if( f1 < f3 ) 	x2 = x1+1;
	else 		x2 = x1-1;
 	f2 = profile[x2];



	i3_flt a = (f1-f2)/((i3_flt)x1 - (i3_flt)x2);
	i3_flt b = f1 - a*(i3_flt)x1; 
	i3_flt x3 = (f3 - b)/a;

	i3_flt fwm = 2.*i3_fabs(xy0-x3); 

	//printf("max_val f1 f2 f3 %2.6f %2.6f %2.6f %2.6f \n",max_val,f1,f2,f3);
	//printf("max_ind x1 x2 x3 %d %d %d %2.6f\n",max_ind,x1,x2,x3);
	//printf("xy0 fwm %2.4f %2.4f \n",xy0,fwm);


	return fwm;

}

void get_image_info(i3_sersics_parameter_set * params, i3_data_set * dataset, i3_flt * image_args, i3_flt * image_info){

// get the image parameters
 
	int original_dataset_upsampling=dataset->upsampling ;
	dataset->upsampling = 9;
        int n_sub = dataset->upsampling;
	int n_pix = dataset->image->nx;
        int n_pad = dataset->padding;
        int n_all = n_pix+n_pad;

// prepare the circular psf
	
        i3_moffat_psf psf_form_circular;
        psf_form_circular.beta = dataset->psf_form->beta;
        psf_form_circular.fwhm = dataset->psf_form->fwhm;
	psf_form_circular.e1 = 0.;
       	psf_form_circular.e2 = 0.;
        
// prepare the parameter set

	i3_sersics_parameter_set * params_circular = malloc(sizeof(i3_sersics_parameter_set));
	memcpy(params_circular,params,sizeof(i3_sersics_parameter_set));
	params_circular->x0 = n_pix/2;
	params_circular->y0 = n_pix/2;
	params_circular->e1 = 0;
	params_circular->e2 = 0;
	
// create image space

	//i3_image * image_npsf = i3_image_create( n_all*n_sub, n_all*n_sub ); i3_image_zero( image_npsf );
	i3_image * image_wpsf = i3_image_create( n_all*n_sub, n_all*n_sub ); i3_image_zero( image_wpsf );
	i3_fourier * fourier_ker = dataset->psf_downsampler_kernel;
	i3_fourier * kernel_circular = i3_build_great10_kernel(&psf_form_circular, n_pix, n_sub, n_pad, GREAT10_PSF_TYPE_MOFFAT);

	//i3_fourier * kernel_elliptic = i3_build_great10_kernel( dataset->psf_form, n_pix, n_sub, n_pad, GREAT10_PSF_TYPE_MOFFAT);
	

// get the unconvolved image
	
	//dataset->psf_downsampler_kernel = NULL;
	//i3_sersics_model_image(params_circular,dataset,image_npsf);
	//if(dataset->options->save_images) i3_image_save_fits( image_npsf, "image_test_npsf.fits" );

// get the convolved image

	dataset->psf_downsampler_kernel = kernel_circular;
        i3_sersics_model_image(params_circular,dataset,image_wpsf);
	if(dataset->options->save_images) i3_image_save_fits( image_wpsf, "image_test_wpsf.fits" );

// measure what needs to be measured
	
	image_info[0] = measure_fwxm(image_wpsf, image_args[0])/(i3_flt)n_sub ;
	image_info[1] = measure_fwxm(image_wpsf, image_args[1])/(i3_flt)n_sub ;
	image_info[2] = measure_fwxm(image_wpsf, image_args[2])/(i3_flt)n_sub ;
	//image_info[3] = measure_entropy(image_npsf);
	//image_info[4] = measure_entropy(image_wpsf);
	//image_info[5] = measure_kldivergence(image_npsf,image_wpsf);
	//image_info[6] = measure_kldivergence(image_wpsf,image_npsf);
	image_info[3] = 66.;
	image_info[4] = 66.;
	image_info[5] = 66.;
	image_info[6] = 66.;

// do the same with elliptical profiles

	//dataset->psf_downsampler_kernel = NULL;
	//i3_sersics_model_image(params,dataset,image_npsf);
	//dataset->psf_downsampler_kernel = kernel_elliptic;
	//i3_sersics_model_image(params,dataset,image_wpsf);

	//image_info[7] = measure_kldivergence(image_npsf,image_wpsf);
	//image_info[8] = measure_kldivergence(image_wpsf,image_npsf);

	image_info[7] = 66.;
	image_info[8] = 66.;

// clean up

	dataset->psf_downsampler_kernel = fourier_ker;
	dataset->upsampling = original_dataset_upsampling;

	i3_image_destroy( image_wpsf );
	//i3_image_destroy( image_npsf );
	i3_fourier_destroy( kernel_circular );
	//i3_fourier_destroy( kernel_elliptic );
	free(params_circular);
	
}


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

// init 

        //Initialize FFTW
        i3_fftw_load_wisdom();
        atexit(i3_fftw_save_wisdom);

// check the arguments
        
        int n_args = 16+1;
        if(argc !=n_args) {
        
                printf("You supplied %d of %d args.\n",argc,n_args);
                printf("1  filename_ini - name of the ini file\n");
                printf("2  dir_out - output directory \n");
                printf("3  proc_id - task id [int]\n");
                printf("4  params_set_id - id of the parameter set id [int]\n");
                printf("5  n_reals - number of noise realisation per core [int] (if 0, then just save the model images and quit)\n");
                printf("6  psf_type - 0 for Moffat, 1 for Airy\n");
                printf("7  psf_beta\n");
                printf("8  psf_fwhm\n");
                printf("9  psf_e1\n");
                printf("10 psf_e2\n");
                printf("11 gal_r - r_bulge, r_disc will be calculated\n");
                printf("12 gal_r_ratio - r_bulge/r_disc (overwrites the vals in .ini file)\n");
                printf("13 gal_f_ratio - flux ratio defined as flux_bulge/(flux_bulge+flux_disc)\n");
                printf("14 gal_e1\n");
                printf("15 gal_e2\n");
                printf("16 SNR - signal to noise ratio\n");

                return 0;
        }
        
                
// read the arguments

        char * filename_ini = argv[1];
        char * dir_out = argv[2];
        int proc_id = atoi(argv[3]);
        int params_set_id = atoi(argv[4]);
        int n_reals = atoi(argv[5]);
        int psf_type = atoi(argv[6]);
        i3_flt psf_beta = atof(argv[7]);
        i3_flt psf_fwhm = atof(argv[8]);
        i3_flt psf_e1 = atof(argv[9]);
        i3_flt psf_e2 = atof(argv[10]);
        i3_flt gal_r = atof(argv[11]) - BEERMAT_MIN_R;
        i3_flt gal_r_ratio = atof(argv[12]);
        i3_flt gal_f_ratio = atof(argv[13]);            
        i3_flt gal_e1 = atof(argv[14]);
        i3_flt gal_e2 = atof(argv[15]);
        i3_flt snr = 0.;//atof(argv[16]);
        

// load the options
 
        i3_options * options = i3_options_default();
        i3_options_read(options, filename_ini);

        int n_pix = options->stamp_size;
        int verb  = options->verbosity;

        i3_model * model = i3_model_create(options->model_name,options);
                        
	i3_init_rng_multiprocess(proc_id);
	//i3_gsl_init_rng_environment();
	// MH: Initialise random number generator with fixed seed
	//i3_gsl_init_rng_fixedseed(722012);

	printf("=====================================\n");
	printf("Check seed of random number generator\n");
	printf("First random number %f\n", i3_random_uniform());
	printf("=====================================\n");

// ==========================================================================
// whether to use gradients or not, global option in ini file is overwritten!
// ==========================================================================

	options->levmar_use_analytic_gradients=true;

	char * expID;
	if(options->levmar_use_analytic_gradients){
	  printf("Using gradients\n");	  
	  expID="exact";
	}
	else{
	  printf("Not using gradients\n");
	  expID="approx_cent";
	}
	printf("=====================================\n");
	
	char filename_result[256]; snprintf(filename_result,256,"%s/i3_jacobian_test_result_%06d_%06d_%s.dat",dir_out,params_set_id,proc_id,expID);
	FILE * file_result = fopen(filename_result,"w");
	if(!file_result) I3_FATAL("Cannot open nbc_result file for writing",0);

	char filename_starts[256]; snprintf(filename_starts,256,"%s/i3_jacobian_test_starts_%06d_%06d_%s.dat",dir_out,params_set_id,proc_id,expID);
	FILE * file_starts = fopen(filename_starts,"w");
	if(!file_starts) I3_FATAL("Cannot open nbc_starts file for writing",0);

	char filename_truth[256]; snprintf(filename_truth,256,"%s/i3_jacobian_test_truth_%06d_%06d_%s.dat",dir_out,params_set_id,proc_id,expID);
	FILE * file_truth = fopen(filename_truth,"w");
	if(!file_truth) I3_FATAL("Cannot open nbc_starts file for writing",0);
	
	char filename_nbc_params[256]; snprintf(filename_nbc_params,256,"%s/i3_jacobian_test_nbc_params_%06d_%s.dat",dir_out,params_set_id,expID);
	FILE * file_params = fopen(filename_nbc_params,"w");
	if(!file_params) I3_FATAL("Cannot open nbc_params file for writing ",0);

// loop timer
     
        time_t time_start, time_end;
	double cpu_time_used;
        time_start = time(NULL);

// make psf
        
        i3_moffat_psf psf_form;
        psf_form.beta = psf_beta;
        psf_form.e1 = psf_e1;
        psf_form.e2 = psf_e2;
        psf_form.fwhm = psf_fwhm;

// prep image used as noisy

        i3_image * image_noisy = i3_image_create( n_pix, n_pix ); i3_image_zero( image_noisy );
	i3_image * image_weight      = i3_image_create( n_pix, n_pix ); i3_image_zero( image_weight );
        i3_data_set * dataset = i3_build_dataset(options, 0, image_noisy, image_weight, &psf_form, psf_type);
	dataset->start_params = malloc(sizeof(i3_sersics_parameter_set));
        
// prepare the paramters

        i3_sersics_parameter_set * params_gal = i3_model_option_starts(model, options);
        
        params_gal->e1 = gal_e1;
        params_gal->e2 = gal_e2;
        params_gal->radius= gal_r;
        params_gal->radius_ratio = gal_r_ratio;
	
	i3_flt tmp_f1 = gal_f_ratio/(1. - gal_f_ratio);
	i3_flt tmp_f2 = 1.;
	i3_flt f1=tmp_f1/(tmp_f1+tmp_f2);
	i3_flt f2=tmp_f2/(tmp_f1+tmp_f2);
	i3_flt r1= params_gal->radius;
	i3_flt r2= params_gal->radius/params_gal->radius_ratio;

	params_gal->bulge_A = f1/i3_sersic_total_flux(params_gal->bulge_index,1.,r1*r1);
	params_gal->disc_A  = f2/i3_sersic_total_flux(params_gal->disc_index ,1.,r2*r2);

// set up space fot the best found

        i3_image * bestI = i3_image_copy( image_noisy );

// get the image info

	i3_flt image_args[NBC_IMAGE_INFO_SZ] = {0.5, 0.75, 0.9,0,0,0,0,0,0};
	i3_flt image_info[NBC_IMAGE_INFO_SZ] = {0,1,2,3,4,5,6,7,8};

        char * str_fmt; 
	str_fmt = "%d\t%d\t%f\t%f\t%2.3f\t%2.3f\t%2.6f\t%2.6f\t%2.6e\t%2.6e\t%2.6e\t%2.2f";
        fprintf(file_params,str_fmt,params_set_id,psf_type,psf_beta,psf_fwhm,psf_e1,psf_e2,gal_e1,gal_e2,gal_r,gal_r_ratio,gal_f_ratio,snr);

	get_image_info(params_gal, dataset, image_args, image_info );
	for(int i=0;i<NBC_IMAGE_INFO_SZ;i++) fprintf(file_params, "\t% 2.4e",image_info[i]); printf("\n");
        fclose(file_params);

// Set the method as levmar

	i3_minimizer_method method;
	method = i3_minimizer_method_levmar;

// if only saving images
	
	if(n_reals<=0){ // save images and exit

        	if(verb > 0) printf("nbc is saving images only \n");
		
		get_image_info(params_gal, dataset, image_args, image_info );	

		for(int i=0;i<NBC_IMAGE_INFO_SZ;i++) fprintf(stdout, "\t% 2.4e",image_info[i]); 
		printf("\n");	
		for(int i=0;i<NBC_IMAGE_INFO_SZ;i++) fprintf(stdout, "\t% 2.4e",image_args[i]);
		printf("\n");

		dataset->options->save_images=true;
	
		params_gal->x0= n_pix/2 ; 
                params_gal->y0= n_pix/2 ; 
                        
                i3_image_zero(image_noisy);
                i3_sersics_model_image(params_gal,dataset,image_noisy);

		i3_image_save_fits( image_noisy, "model_image.fits" );

		
		return 0;
		
	}

// verb 

        if(verb > 0) printf("[%d] running noise bias calibration with n_reals=%d minimizer_loops=%d verb=%d\n",proc_id,n_reals,options->minimizer_loops,verb);

// run the noise realisations loop

        for(int nr=0;nr<n_reals;nr++){

        // create the true image
                
                params_gal->x0= n_pix/2 + (i3_random_uniform() - 0.5); 
                params_gal->y0= n_pix/2 + (i3_random_uniform() - 0.5);
                        
                i3_image_zero(image_noisy);
                i3_sersics_model_image(params_gal,dataset,image_noisy);


	// calculate the noise and create the weight map
                
                i3_flt std_noise = 1e-8;
                if(snr>0) std_noise = i3_snr_to_noise_std(snr, image_noisy);
                options->noise_sigma = std_noise;  
				               
                i3_image_fill(image_weight,pow(std_noise,-2));
                //i3_image_fill(weight,1.);             
                
                dataset->noise = std_noise;
		               
	// add noise to image           
                
                if(snr>0) i3_image_add_white_noise( image_noisy,std_noise );

	               
		if(verb>0){
		  char filename[128];
		  snprintf( filename, 128, "image_noisy_%04d.fits", proc_id );
		  i3_image_save_fits( image_noisy, filename );
		}   

	// set up the start positions

		i3_sersics_parameter_set * params_start = malloc(sizeof(i3_sersics_parameter_set));
		set_start_params(model, dataset, options, params_start, params_gal);
		i3_model_copy_parameters(model,dataset->start_params,params_start); 
 
		//i3_sersics_parameter_set * params_ml = i3_analyze_dataset_method(dataset, model, options, &bestL, bestI, method);
		i3_sersics_parameter_set * params_ml_beermat = minimizer_run_levmar(dataset, options);
	        i3_sersics_parameter_set * params_ml = i3_sersics_beermat_params(params_ml_beermat);


		get_image_info(params_ml, dataset, image_args, image_info );	
					
		if(verb>0){
				
		  	save_starts(model, params_gal,		stdout, params_set_id );
		  	save_starts(model, params_start,	stdout, params_set_id );
		        save_result(model, params_ml,		stdout, params_set_id, dataset,image_info);
		  
		}
		
		save_starts(model,params_start,file_starts,params_set_id);
		save_truth(model,params_gal,file_truth,params_set_id);
		save_result(model,params_ml,file_result,params_set_id,dataset,image_info);
		
		if(options->save_images){
			i3_image * bestfit_image = i3_image_create(n_pix,n_pix);
			i3_sersics_likelihood(bestfit_image,params_ml,dataset);
			i3_image_save_fits(bestfit_image,"image_bestfit.fits");
			i3_image_destroy(bestfit_image);
		}

		
		free(params_ml_beermat);
		free(params_ml);
		free(params_start);
	
		time_end =  time(NULL);
        	cpu_time_used = (double)(time_end - time_start);
		if(verb>-1  ) printf("[% d:% d] % 3d / % 3d\ttime total % 4.2fm\ttime per galaxy % 2.2fs \n",params_set_id,proc_id,nr,n_reals,cpu_time_used/60.,cpu_time_used/(i3_flt)(nr+1)); 
		//if(verb>-1 && nr%10==0 ) printf("[%d:%d] %d / %d\n",params_set_id,proc_id,nr,n_reals);
        }

        i3_image_destroy( image_noisy );
	i3_image_destroy( image_weight );

        // show the time
        fclose(file_result);
 	time_end = time(NULL);
        cpu_time_used = (double)(time_end - time_start);
        if(verb>-1) printf("time total %4.2fm\ttime per galaxy %2.2fs\n",cpu_time_used/60.,cpu_time_used/(i3_flt)n_reals);
        
}