#ifndef _H_I3_MODEL_TOOLS
#define _H_I3_MODEL_TOOLS

#include "i3_image.h"

void i3_add_real_space_gaussian(i3_image * image, i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0);

void i3_add_real_space_king_profile(i3_image * image, i3_flt rt, i3_flt rc, i3_flt e1, i3_flt e2, i3_flt k, i3_flt x0, i3_flt y0);

void i3_add_real_space_beta_profile(i3_flt rc, i3_flt e1, i3_flt e2, i3_flt A, i3_flt beta, i3_flt x0, i3_flt y0, i3_image * image);


void i3_add_real_space_sersic(i3_image * image, i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0, i3_flt sersic_n);

void i3_add_fourier_space_sersic(i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_fourier * fourier_image);

void i3_add_fourier_space_exponential(i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0, i3_fourier * fourier);


void i3_compute_standard_precision_matrix(i3_flt e1, i3_flt e2, i3_flt ab, i3_flt P[4]);

i3_flt i3_moffat(i3_flt r2_in, i3_flt rd2_in, i3_flt rc2_in, i3_flt beta_in);
i3_flt i3_moffat_rd2(i3_flt fwhm_in, i3_flt beta_in);

/**
 * Add a sersic profile truncated at some radius (the other function truncates at some value instead).  
 * The radius is specified in units of r_0 = sqrt(ab).
*/
void i3_add_real_space_sersic_truncated_radius(i3_image * image, i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_flt truncation_factor);

/**
 * Add a sersic profile truncated at some radius (the other function truncates at some value instead).
 * The radius is specified in units of r_0 = sqrt(ab).
 * n_central_pixels_to_upsample is the number of pixels around central pixel, for which the average value will be computed of n_central_upsampling subixels, evenly spaced. 
*/
void i3_add_real_space_sersic_truncated_radius_upsample_central(i3_image * image, i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_flt truncation_factor, int n_central_pixels_to_upsample, int n_central_upsampling);

/** 
 * This variant does do central pixel upsampling but is not truncated
*/
void i3_add_real_space_sersic_upsample_central(i3_image * image, i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, int n_central_pixels_to_upsample, int n_central_upsampling);


/**
 * Compute Jacobians
 */
void i3_add_real_space_sersic_truncated_radius_upsample_central_jac(i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_flt truncation_factor, int n_central_pixels_to_upsample, int n_central_upsampling,  i3_image ** image);

void i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact(i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_flt truncation_factor, int n_central_pixels_to_upsample, int n_central_upsampling,  i3_image ** image);

void i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact_vec(i3_flt ab,i3_flt e1,i3_flt e2,i3_flt A,i3_flt x0,i3_flt y0,i3_flt sersic_n, i3_flt truncation_factor, int n_central_pixels_to_upsample, int n_central_upsampling,  i3_image ** image);

/**
Calculate the chisq between model image and observed image - both real space
 */ 
i3_flt i3_chi2_white_noise(i3_image * model_image, i3_image * observed_image, i3_flt sigma);

i3_flt i3_chi2_wn_reduced(i3_image * model_image, i3_image * observed_image, i3_flt noise_std);

/**
Calculate the chisq between model image and observed image - both in fourier space
 */
i3_flt i3_chi2_white_noise_fourier(i3_fourier * model_fourier, i3_fourier * observed_fourier, i3_flt sigma);

/**
Calculate the logL of image given model with marginalised amplitude.
 */
i3_flt i3_logL_white_noise_marg_A(i3_image * model_image, i3_image * observed_image, i3_flt noise_std);

/**
NOT WORKING YET Calculate the logL of image given model with marginalised amplitude (analytically) and x-y centre position (numerically) */
i3_flt i3_logL_white_noise_marg_A_xy0(i3_image * model_image, i3_image * observed_image, i3_flt noise_std, int n_additional_pix, int n_sub_pix);

/**
 * Calculate the logL of image given x0 and y0 with marginalised amplitude (analytically). Outupts log-pdf in 2 parameters: x0 and y0. Uses the convolved subpixel grid to calculate it.
 */
i3_flt * i3_logL_in_xy0(i3_image * model_image, i3_image * observed_image, i3_flt noise_std, int n_additional_pix, int n_sub_pix);


/**
 * Calculate the chisq between a set of model images observed image - both in real space. Uses reduced model set, which contain only models with positive e1 and e2. For other values the chis2 is marked as 1e4.
 */
i3_flt * i3_chi2_with_models_white_noise_reduced(i3_image ** all_models, int n_models, i3_image * observed_image, i3_flt noise_std);

/**
 * Calculates chi2 vector ot the observed_image with the models
 */
i3_flt * i3_chi2_with_models_white_noise(i3_image ** all_models, int n_models, i3_image * observed_image, i3_flt noise_std);
i3_flt * i3_logL_with_models_white_noise_marg_A(i3_image ** all_models, int n_models, i3_image * observed_image, i3_flt noise_std);


/**
 * Calculates the likelihood of models and observed_image. Uses the information about true e1 and e2 to avoid searching in regions of parameter space where pdf = 0. It tests only models, which abs(e) - e_true (<) than cut_off 
 * \param all_models Contains the array of models
 * \param n_models Number of models. n_models = n_grid_e*n_grid_e*n_grid_r
 * \param models_e1 Values of e1 for all_models
 * \param models_e2 Values of e1 for all_models
 * \param observed_image Noisy image
 * \param e1 True e1 of the noisy image
 * \param e2 True e2 of the noisy image
 * \param noise_std Std of the noise
 */
i3_flt * i3_logL_with_models_white_noise_fast(i3_image ** all_models, int n_models, i3_flt * models_e1, i3_flt * models_e2, i3_image * observed_image, i3_flt e1, i3_flt e2, i3_flt cut_off , i3_flt noise_std, int metric);

/**
 * Calculates the likelihood of models and observed_image. Uses the information about true e1 and e2 to avoid searching in regions of parameter space where pdf = 0. It tests only models, which abs(e) - e_true (<) than cut_off 
 * \param all_models Contains the array of models
 * \param n_models Number of models. n_models = n_grid_e*n_grid_e*n_grid_r
 * \param models_e1 Values of e1 for all_models
 * \param models_e2 Values of e1 for all_models
 * \param observed_image Noisy image
 * \param e1 True e1 of the noisy image
 * \param e2 True e2 of the noisy image
 * \param noise_std Std of the noise
 */
i3_flt *  i3_marginalise_r(i3_flt * pdf_full, int n_pdf_full, int n_pdf_marg);


void i3_covariance_matrix_proposal(i3_mcmc_ptr mcmc, i3_parameter_set * current, i3_parameter_set * proposed);
i3_flt i3_fwhm_to_ab(i3_flt fwhm, i3_flt sersic_n);
i3_flt i3_r_to_fwhm(i3_flt r, i3_flt sersic_n);


/**
 * Calculate the standard deviation of the noise from signal to noise ratio for a given image
 */
i3_flt i3_snr_to_noise_std(i3_flt snr, i3_image * image);

/**
 * Returns an output array which contains exponent of input array. To avoid precision issues, the input array is scaled, so that the output array value of maximum value in the input array is 1. The input array is not overwritten.
 * \param array_log Input array to exponentiate
 * \param n Size of the input array
*/

i3_flt * i3_pdf_from_logL(i3_flt * logL, int n);

/** 
 * Calculate the chi2 between model and data images given a weight map that describes the varying (though still uncorrelated) noise.
 * The weight map should be the inverse of the noise variance on the pixel
*/
i3_flt i3_chi2_weight_map(i3_image * model_image, i3_image * observed_image, i3_image * weight_image);


/** 
 * Calculate the likelihood of model and data images given a weight map that describes the varying (though still uncorrelated) noise.
 * The weight map should be the inverse of the noise variance on the pixel.  This variant marginalizes over the amplitude of the image.
*/
i3_flt i3_logL_weight_map_marginalize_amplitude(i3_image * model_image, i3_image * observed_image, i3_image * weight_image);

/** 
 * Calculate the likelihood of model and data images given a weight map that describes the varying (though still uncorrelated) noise.
 * The weight map should be the inverse of the noise variance on the pixel.  This variant maximises over the amplitude of two components of an image.
*/
i3_flt i3_logL_weight_map_maximise_amplitude_2components(i3_image * model_image, i3_image * model_image2, i3_image * observed_image, i3_image * weight_image, i3_flt * amplitude1, i3_flt * amplitude2);

/** 
 * Calculate the likelihood of model and data images given a weight map that describes the varying (though still uncorrelated) noise.
 * The weight map should be the inverse of the noise variance on the pixel.  This variant marginalizes over the amplitude of two components of an image with the constraint that they have to be positive
*/
i3_flt i3_logL_weight_map_maximise_amplitude_2components_positive(i3_image * model_image1, i3_image * model_image2, i3_image * observed_image, i3_image * weight_image, i3_flt * amplitude1, i3_flt * amplitude2);

/**
 * Use the fitting formula From Voigt & Bridle 2010 (via Capaccioli 1989) 
 * to get the kappa value for a given sersic index n
 * which means you can use the sensible sersic radius r_0.
 * 
 * This function is inline, which can be faster for smaller functions
 * (basically means the code is copied+pasted into other functions by)
 * the compiler instead of called as a separate function with those overheads
*/
static
inline i3_flt i3_sersic_kappa(i3_flt n) {return 1.9992*n - 0.3271;};
i3_flt i3_sersic_total_flux(i3_flt sersic_index, i3_flt amplitude, i3_flt ab);
i3_flt i3_sersic_flux_to_radius(i3_flt n, i3_flt A, i3_flt ab, i3_flt r);

void i3_precision_matrix_jacobian(i3_flt e1, i3_flt e2, i3_flt ab, i3_flt * J);

i3_flt i3_estimate_noise_rms_observed_model(i3_image * observed, i3_image * model, i3_image * mask);

i3_flt i3_signal_to_noise_fit(i3_image * model, i3_flt noise, i3_image * mask);
i3_flt i3_signal_to_noise_model_fitted(i3_image * observed, i3_image * model, i3_image * mask);
i3_flt i3_signal_to_noise_with_weight(i3_image * model, i3_image * weight);
i3_flt i3_signal_to_noise_with_weight_and_data(i3_image * model, i3_image * data, i3_image * weight);

void i3_ellipticity_eta_to_e(i3_flt eta1, i3_flt eta2, i3_flt * e1,   i3_flt * e2);
void i3_ellipticity_e_to_eta(i3_flt e1,   i3_flt e2,   i3_flt * eta1, i3_flt * eta2);

i3_flt i3_sersic_log_flux_to_amplitude(i3_flt log_flux, i3_flt sersic_index, i3_flt ab);

i3_flt i3_poisson_like(i3_image * model_image, i3_image * observed_image);

i3_flt i3_poisson_like_exposure_background(i3_image * model_image, i3_image * observed_image, 
	i3_image * exposure_map, i3_flt background);

i3_flt i3_pixel_shift_image(i3_image * image, i3_image * shifted, int dx, int dy);

#endif