#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" #include "i3_utils.h" #include "i3_great.h" #include "i3_image_fits.h" #include "i3_analyze.h" #include "i3_version.h" #include "extern/levmar/lm.h" #define MAX_MINIMIZER_LOOPS_NAN 1 i3_data_set * i3_build_dataset_truncated_psf(i3_options * options, gal_id identifier, i3_image * image, i3_image * weight, i3_moffat_psf * psf_form, int psf_type, i3_flt truncation){ i3_data_set * data_set = (i3_data_set*)malloc(sizeof(i3_data_set)); data_set->identifier=identifier; data_set->image = image; data_set->psf_downsampler_kernel = i3_fourier_conv_kernel_moffat(image->nx,options->upsampling, options->padding, psf_form->beta, psf_form->fwhm, psf_form->e1, psf_form->e2, truncation, psf_type); data_set->noise=options->noise_sigma; data_set->psf_form = psf_form; data_set->psf_fwhm = psf_form->fwhm; data_set->central_pixel_upsampling=options->central_pixel_upsampling; data_set->n_central_pixel_upsampling=options->n_central_pixel_upsampling; data_set->n_central_pixels_to_upsample=options->n_central_pixels_to_upsample; data_set->weight=weight; data_set->amplitude1 = 0.0; data_set->amplitude2 = 0.0; data_set->stamp_size=options->stamp_size; data_set->upsampling=options->upsampling; data_set->padding=options->padding; data_set->start_params = NULL; data_set->options=options; i3_dataset_reset_counters(data_set); data_set->levmar_info = malloc(sizeof(i3_flt)*LM_INFO_SZ); data_set->covariance_estimate = NULL; return data_set; } i3_data_set * i3_build_dataset(i3_options * options, gal_id identifier, i3_image * image, i3_image * weight, i3_moffat_psf * psf_form, int psf_type){ i3_flt truncation = FLT_MAX; return i3_build_dataset_truncated_psf(options, identifier, image, weight, psf_form, psf_type, truncation); } i3_data_set * i3_build_dataset_with_psf_image(i3_options * options, gal_id identifier, i3_image * image, i3_image * weight, i3_image * psf_image){ i3_data_set * data_set = (i3_data_set*)malloc(sizeof(i3_data_set)); data_set->identifier=identifier; data_set->image = image; data_set->psf_form = NULL; data_set->psf_fwhm = i3_image_standard_fwhm(image); if (options->use_hankel_transform){ // fft high res psf image, then downsample (to match model fourier image) to create low res fourier psf image } else data_set->psf_downsampler_kernel = i3_fourier_conv_kernel(image->nx,options->upsampling, options->padding, psf_image, options->perform_pixel_integration); data_set->noise=options->noise_sigma; data_set->central_pixel_upsampling=options->central_pixel_upsampling; data_set->n_central_pixel_upsampling=options->n_central_pixel_upsampling; data_set->n_central_pixels_to_upsample=options->n_central_pixels_to_upsample; data_set->weight=weight; data_set->stamp_size=options->stamp_size; data_set->upsampling=options->upsampling; data_set->padding=options->padding; data_set->start_params = NULL; data_set->options=options; i3_dataset_reset_counters(data_set); data_set->levmar_info = malloc(sizeof(i3_flt)*LM_INFO_SZ); data_set->covariance_estimate = NULL; return data_set; } // static void save_result(i3_model * model, i3_parameter_set * best, i3_image * model_image, i3_image * data, i3_flt like, char * output_dir, FILE * output_file, int identifier, i3_data_set * data_set){ // // if (like==BAD_LIKELIHOOD) like = 0.0/0.0; // //Print the likelihood, followed by the model parameters, all on one line. // fprintf(output_file,"%d\t%15.10e\t",identifier,like); // i3_model_fprintf_parameters(model,output_file,best); // fprintf (output_file, "\t%d",data_set->n_logL_evals); // fprintf (output_file, "\t%f",data_set->signal_to_noise); // fprintf (output_file, "\t%f",data_set->amplitude1); // fprintf (output_file, "\t%f",data_set->amplitude2); // fprintf(output_file,"\n"); // fflush(output_file); // // if (output_dir){ // //Now save the image and the model to file for later perusal. // char filename[i3_max_string_length]; // snprintf(filename,i3_max_string_length,"%s/model.%d.txt",output_dir,identifier); // i3_image_save_text(model_image,filename); // // snprintf(filename,i3_max_string_length,"%s/image.%d.txt",output_dir,identifier); // i3_image_save_text(data,filename); // } // // } const char * i3_analyze_sersics_header_line() { return "# ID\tcatalog_x\tcatalog_y\tLikelihood\tx0\ty0\te1\te2\tradius\tbulge_A\tdisc_A\tflux_ratio\tfwhm\tsnr\tmin_residuals\tmax_residuals\tmodel_min\tmodel_max\tpsf_e1\tpsf_e2\tpsf_fwhm\tpsf_beta\tRgpp_Rp\tlevmar_e0\tlevmar_e\tlevmar_J\tlevmar_Dp\tlevmar_mu\tlevmar_De\tlevmar_nit\tlevmar_reason\tlevmar_neval\tlevmar_njac\tlevmar_nlin\ttime_taken\n#imshape-version: " HG_VERSION_STRING "\n"; } // currently this is just a copy of save_result_nbc // This should migrate into a standard im3shape output format eventually // Need to rewrite this output. // remove levmar // think about rgp/rp void copy_result_im3shape_to_struct(i3_options * options, i3_model * model, i3_parameter_set * params, i3_flt likelihood, i3_result * results, gal_id identifier, i3_data_set * dataset, i3_flt min_residuals, i3_flt max_residuals, i3_flt model_min, i3_flt model_max, double time_taken, i3_flt * image_info){ // get the fluxeas as if they were not truncated untruncated if (00==strcmp(model->name,"unknown_model")){ results->fluxes[0] = -1.0; results->fluxes[1] = -1.0; results->flux_ratio = -1.0; } #ifdef HAVE_SERSICS_MODEL else if (0==strcmp(model->name,"sersics")){ i3_sersics_parameter_set * sersic_params = (i3_sersics_parameter_set * ) params; i3_flt bulge_radius = sersic_params->radius; i3_flt disc_radius = sersic_params->radius/sersic_params->radius_ratio; i3_flt bulge_flux = i3_sersic_total_flux( sersic_params->bulge_index, sersic_params->bulge_A, bulge_radius); i3_flt disc_flux = i3_sersic_total_flux( sersic_params->disc_index, sersic_params->disc_A, disc_radius); results->fluxes[0] = bulge_flux; results->fluxes[1] = disc_flux; results->flux_ratio = bulge_flux/(bulge_flux + disc_flux); } #endif #ifdef HAVE_MULTISERSICS_MODEL else if (0==strcmp(model->name,"multisersics")){ i3_multisersics_parameter_set * sersic_params = (i3_multisersics_parameter_set * ) params; i3_flt bulge_radius = sersic_params->radius; i3_flt disc_radius = sersic_params->radius/sersic_params->radius_ratio; i3_flt bulge_flux = i3_sersic_total_flux( sersic_params->bulge_index, sersic_params->bulge_A, bulge_radius); i3_flt disc_flux = i3_sersic_total_flux( sersic_params->disc_index, sersic_params->disc_A, disc_radius); results->fluxes[0] = bulge_flux; results->fluxes[1] = disc_flux; results->flux_ratio = bulge_flux/(bulge_flux + disc_flux); } #endif #ifdef HAVE_LOGSERSICS_MODEL else if (0==strcmp(model->name,"logsersics")){ i3_logsersics_parameter_set * sersic_params = (i3_logsersics_parameter_set * ) params; i3_flt bulge_flux = exp(sersic_params->bulge_A); i3_flt disc_flux = exp(sersic_params->disc_A); results->fluxes[0] = bulge_flux; results->fluxes[1] = disc_flux; results->flux_ratio = bulge_flux/(bulge_flux + disc_flux); } #endif else{ results->fluxes[0] = -1.0; results->fluxes[1] = -1.0; results->flux_ratio = -1.0; } results->identifier = identifier; results->time_taken = time_taken; i3_model_copy_parameters(model, results->params, params); // Image statistics results->stats_signal_to_noise = dataset->signal_to_noise; results->stats_old_signal_to_noise = dataset->old_signal_to_noise; results->stats_min_residuals = min_residuals; results->stats_max_residuals = max_residuals; results->stats_model_min = model_min; results->stats_model_max = model_max; /* * O: information regarding the minimization. Set to NULL if don't care * info[0]= ||e||_2 at initial p. * info[1-5]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, mu/max[J^T J]_ii ], D||e||_2, all computed at estimated p. * info[6]= # iterations, * info[7]=reason for terminating: 1 - stopped by small gradient J^T e * 2 - stopped by small Dp * 3 - stopped by itmax * 4 - singular matrix. Restart from current p with increased mu * 5 - no further error reduction is possible. Restart with increased mu * 6 - stopped by small ||e||_2 * 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error * 8 - stopped by small D||e||_2 * info[8]= # function evaluations * info[9]= # Jacobian evaluations * info[10]= # linear systems solved, i.e. # attempts for reducing error * printf("Minimizer stopped because: %s\n",message); * i3_flt norm_start = i3_array_norm(start,minimizer->nparam); * printf("||J^T e||_inf = % 2.4e (% 2.4e)\n",info[2],levmar_options[1]); * printf("||Dp||_2 = % 2.4e (% 2.4e)\n",info[3]*norm_start,levmar_options[2]); * printf("||e||_2 = % 2.4e (% 2.4e)\n",info[1],levmar_options[3]); * printf("D||e||_2 = % 2.4e (% 2.4e)\n",info[4],levmar_options[4]); * printf("mu/max[J^T J]_ii = % 2.4e\n",info[5]); * printf("niter = % 1.0f\n",info[6]); */ // Info on minimization results->likelihood = likelihood; results->levmar_residual_error_at_start = dataset->levmar_info[0]; results->levmar_residual_error_at_end = dataset->levmar_info[1]; results->levmar_maximal_gradient_projected_residual = dataset->levmar_info[2]; results->levmar_difference_in_parameter_vector = dataset->levmar_info[3]; results->levmar_difference_in_residual_error = dataset->levmar_info[4]; results->levmar_maximal_gradient_component = dataset->levmar_info[5]; results->levmar_number_of_iterations = (int)dataset->levmar_info[6]; // int results->levmar_reason_of_termination = (int)dataset->levmar_info[7]; // int results->levmar_number_of_likelihood_evaluations = (int)dataset->levmar_info[8]; // int results->levmar_number_of_gradient_evaluations = (int)dataset->levmar_info[9]; // int results->levmar_number_of_linear_systems_solved = (int)dataset->levmar_info[10]; // int results->number_varied_params = i3_model_number_varied_params(model); if (dataset->covariance_estimate){ int np = results->number_varied_params; if (np*np>MAX_NP2){ I3_FATAL("Too many varied parameters in model - more than limit of 32", np*np); } memcpy(results->covariance_estimate, dataset->covariance_estimate, np*np*sizeof(i3_flt)); } } void save_result_sersics(i3_options * options, i3_model * model, i3_parameter_set * params_input, i3_flt likelihood, FILE * output_file, i3_catalog_row * row, i3_data_set * dataset, i3_flt min_residuals, i3_flt max_residuals, i3_flt model_min, i3_flt model_max, double time_taken, i3_flt * image_info) { #if NO_SERSICS I3_FATAL("You should never see this error!",666); #else i3_sersics_parameter_set * params = (i3_sersics_parameter_set * ) params_input; i3_flt bulge_radius = params->radius; i3_flt disc_radius = params->radius/params->radius_ratio; // get the fluxeas as if they were not truncated untruncated i3_flt bulge_flux = i3_sersic_total_flux( params->bulge_index, params->bulge_A, pow(bulge_radius,2)); i3_flt disc_flux = i3_sersic_total_flux( params->disc_index, params->disc_A, pow(disc_radius, 2)); //i3_flt flux_ratio = params->bulge_A/(params->bulge_A + params->disc_A); i3_flt flux_ratio = bulge_flux/(bulge_flux + disc_flux); i3_flt fwhm = image_info[0]; fprintf(output_file, "%ld\t% 2.8e\t% 2.8e \t%e",row->id,row->x,row->y,likelihood); //i3_model_fprintf_parameters( model, output_file, params); fprintf(output_file, "\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e",params->x0,params->y0,params->e1,params->e2,params->radius,params->bulge_A,params->disc_A); fprintf(output_file, "\t% 2.8f",flux_ratio); fprintf(output_file, "\t% 2.2e",fwhm); fprintf(output_file, "\t% 2.4f",dataset->signal_to_noise); fprintf(output_file, "\t% 2.4e",min_residuals); fprintf(output_file, "\t% 2.4e",max_residuals); fprintf(output_file, "\t% 2.4e",model_min); fprintf(output_file, "\t% 2.4e",model_max); if (dataset->psf_form){ fprintf(output_file, "\t% 2.6e",dataset->psf_form->e1); fprintf(output_file, "\t% 2.6e",dataset->psf_form->e2); fprintf(output_file, "\t% 2.6e",dataset->psf_form->fwhm); fprintf(output_file, "\t% 2.6e",dataset->psf_form->beta); fprintf(output_file, "\t% 2.6e",fwhm / dataset->psf_form->fwhm); } else{ fprintf(output_file, "\t% 2.6e",0.0); fprintf(output_file, "\t% 2.6e",0.0); fprintf(output_file, "\t% 2.6e",dataset->psf_fwhm); fprintf(output_file, "\t% 2.6e",0.0); fprintf(output_file, "\t% 2.6e",fwhm / dataset->psf_fwhm); } fprintf(output_file, "\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% d\t% d\t% d\t% d\t% d", dataset->levmar_info[0], dataset->levmar_info[1], dataset->levmar_info[2], dataset->levmar_info[3], dataset->levmar_info[4], dataset->levmar_info[5], (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, "\t% 2.4f", time_taken); fprintf(output_file, "\n"); fflush(output_file); #endif } void save_result_logsersics(i3_options * options, i3_model * model, i3_parameter_set * params_input, i3_flt likelihood, FILE * output_file, i3_catalog_row * row, i3_data_set * dataset, i3_flt min_residuals, i3_flt max_residuals, i3_flt model_min, i3_flt model_max, double time_taken, i3_flt * image_info) { #if NO_LOGSERSICS I3_FATAL("You should never see this error!",666); #else i3_logsersics_parameter_set * params = (i3_logsersics_parameter_set * ) params_input; i3_flt bulge_radius = params->radius; i3_flt disc_radius = params->radius/params->radius_ratio; // get the fluxeas as if they were not truncated untruncated i3_flt bulge_flux = i3_exp(params->bulge_A); i3_flt disc_flux = i3_exp(params->disc_A); // i3_sersic_log_flux_to_amplitude(i3_flt log_flux, i3_flt sersic_index, i3_flt ab){ i3_flt bulge_amp = i3_sersic_log_flux_to_amplitude(params->bulge_A, params->bulge_index, bulge_radius*bulge_radius); i3_flt disc_amp = i3_sersic_log_flux_to_amplitude(params->disc_A, params->disc_index, disc_radius*disc_radius); //i3_flt flux_ratio = params->bulge_A/(params->bulge_A + params->disc_A); i3_flt flux_ratio = bulge_flux/(bulge_flux + disc_flux); i3_flt fwhm = image_info[0]; fprintf(output_file, "%ld\t% 2.8e\t% 2.8e \t%e",row->id,row->x,row->y,likelihood); //i3_model_fprintf_parameters( model, output_file, params); fprintf(output_file, "\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e\t% 2.8e",params->x0,params->y0,params->e1,params->e2,params->radius,bulge_amp,disc_amp); fprintf(output_file, "\t% 2.8f",flux_ratio); fprintf(output_file, "\t% 2.2e",fwhm); fprintf(output_file, "\t% 2.4f",dataset->signal_to_noise); fprintf(output_file, "\t% 2.4e",min_residuals); fprintf(output_file, "\t% 2.4e",max_residuals); fprintf(output_file, "\t% 2.4e",model_min); fprintf(output_file, "\t% 2.4e",model_max); fprintf(output_file, "\t% 2.6e",dataset->psf_form->e1); fprintf(output_file, "\t% 2.6e",dataset->psf_form->e2); fprintf(output_file, "\t% 2.6e",dataset->psf_form->fwhm); fprintf(output_file, "\t% 2.6e",dataset->psf_form->beta); fprintf(output_file, "\t% 2.6e",fwhm / dataset->psf_form->fwhm); fprintf(output_file, "\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% d\t% d\t% d\t% d\t% d", dataset->levmar_info[0], dataset->levmar_info[1], dataset->levmar_info[2], dataset->levmar_info[3], dataset->levmar_info[4], dataset->levmar_info[5], (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, "\t% 2.4f", time_taken); fprintf(output_file, "\n"); fflush(output_file); #endif } void save_result_im3shape(i3_options * options, i3_model * model, i3_parameter_set * params_input, i3_flt likelihood, FILE * output_file, i3_catalog_row * row, i3_data_set * dataset, i3_flt min_residuals, i3_flt max_residuals, i3_flt model_min, i3_flt model_max, double time_taken, i3_flt * image_info){ //This function is only valid for sersic models. if (!strcmp(model->name,"sersics")) { save_result_sersics(options, model, params_input, likelihood, output_file, row, dataset, min_residuals, max_residuals, model_min, model_max, time_taken, image_info); return; } if (!strcmp(model->name,"logsersics")){ save_result_logsersics(options, model, params_input, likelihood, output_file, row, dataset, min_residuals, max_residuals, model_min, model_max, time_taken, image_info); return; } // For general models fprintf(output_file, "%ld %le ", row->id, likelihood); i3_model_fprintf_parameters(model, output_file, params_input); fprintf(output_file, "\n"); fflush(output_file); } void save_model_data_images(i3_image * model_image, i3_image * data, char * output_dir, char * filetail){ //Save the image and the model to file for later perusal. char filename[i3_max_string_length]; snprintf(filename,i3_max_string_length,"%s/model%s.fits",output_dir,filetail); i3_image_save_fits(model_image,filename); snprintf(filename,i3_max_string_length,"%s/image%s.fits",output_dir,filetail); i3_image_save_fits(data,filename); } static void print_result(i3_model * model, i3_parameter_set * best, i3_flt like, i3_data_set * data_set){ //Print the likelihood, followed by the model parameters, all on one line. if (like==BAD_LIKELIHOOD) like = 0.0/0.0; printf("% 1.4e\t",like); i3_model_printf_parameters(model,best); //printf("\n"); i3_model_pretty_printf_parameters(model,best); printf("\n"); printf ("\t%d",data_set->n_logL_evals); printf ("\t% 1.4e",data_set->signal_to_noise); printf ("\t% 1.4e",data_set->amplitude1); printf ("\t% 1.4e",data_set->amplitude2); printf ("\t% 1.4e",data_set->amplitude1/(data_set->amplitude2 + data_set->amplitude1)); printf("\n"); fflush(stdout); } i3_parameter_set * i3_analyze_dataset(i3_data_set * data_set, i3_model * model, i3_options * options, i3_flt * best_like, i3_image * best_image){ return i3_analyze_dataset_method(data_set, model, options, best_like, best_image, i3_minimizer_method_levmar); } // This monster function is the core im3shape routine // It launches the minimization and possibly MCMC of the method. i3_parameter_set * i3_analyze_dataset_method(i3_data_set * data_set, i3_model * model, i3_options * options, i3_flt * best_like, i3_image * best_image, i3_minimizer_method method){ i3_dataset_reset_counters(data_set); i3_image * image = data_set->image; //Choose the point from which to start the minimization. //We do this first so we can check for any prior violation bugs before doing any heavy lifting i3_parameter_set * current_best = NULL; //If we have some pre-set starting params, start from those if(data_set->start_params!=NULL){ current_best = malloc(model->nbytes); i3_model_copy_parameters(model,current_best,data_set->start_params); i3_print(i3_verb_noisy, "Starting using PRESET parameters."); } else { // Otherwise we can start with whatever the inifile/options say current_best = i3_model_option_starts(model, options); //And then optionally fill in some of the parameters based on the data if (options->use_computed_start) { i3_model_starting_position(model, current_best, data_set, options); i3_print(i3_verb_noisy, "Starting using COMPUTED parameters."); } // Or just use them as-is else{ i3_print(i3_verb_noisy, "Starting using FIXED INI parameters."); } } // If we are at high verbosity we can print out these starting parameters if (options->verbosity>2) { printf("\n"); i3_model_printf_parameters(model, current_best); printf("\n"); } //If the starting parameters are outside our prior than that is an immediate problem! // Give up if so. if ((options->minimizer_loops>0) && i3_model_prior_violations(model, current_best, stderr)) { free(current_best); return NULL; } //Set up a minimizer object to perform the minimization i3_minimizer * minimizer = i3_minimizer_create(model,options); //Some starting likelihood values and similar. i3_flt like = BAD_LIKELIHOOD; data_set->signal_to_noise=0.0/0.0; int accepted_loops=0; // We may do a number of minimizer loops to keep trying to improve for(int i=0 ; iminimizer_loops || (options->minimizer_loops>0 && ibest_like is actually a chi^2 bool any_nan = i3_model_any_nan(model, new_best); i3_flt best_like = -0.5*minimizer->best_like; // If nothing is NaN and this is an improvement then all is well if (best_like>like && (!any_nan) ) { //Report success i3_print(i3_verb_debug, "Next minimizer iteration was an improvement: like %le -> %le (delta=%le)", like, best_like, best_like-like); //Update current "best" values like = best_like; i3_model_copy_parameters(model,current_best,new_best); // Record S/N (this is just a convenient place) data_set->signal_to_noise = i3_signal_to_noise_with_weight_and_data(minimizer->model_image, data_set->image, data_set->weight); data_set->old_signal_to_noise = i3_signal_to_noise_with_weight(minimizer->model_image, data_set->weight); //And keep track of the number of successful improvement loops accepted_loops++; // We might wish if high verbosity to keep track of the separate minimizer run best-fits if ((options->verbosity>=i3_verb_noisy && iminimizer_loops-1) || options->verbosity>=i3_verb_debug){ printf("loop %02d ML parameters: \n",i); print_result(model,current_best,like,data_set); } } // It may be (usually in second and subsequent minimizer runs, or in the event of failure) // that the new result is no better. else{ i3_print(i3_verb_noisy, "Next minimizer iteration no better: like=%le", like); //Ignore this result - it is not better } free(new_best); } if (options->minimizer_loops<=0){ //No minimization at all - we just use the start value (for testing/debugging) i3_minimizer_setup_model_image(minimizer,data_set->image); like = i3_model_posterior(model,minimizer->model_image,current_best,data_set); data_set->signal_to_noise = i3_signal_to_noise_with_weight_and_data(minimizer->model_image, data_set->image, data_set->weight); data_set->old_signal_to_noise = i3_signal_to_noise_with_weight(minimizer->model_image, data_set->weight); i3_print(i3_verb_quiet, "Warning: no minimization"); } // We now have a best fit from the minimizer. // Get its posteriors probability. // There is some complexity in where we use likelihoods and where posteriors - // basically we have to minimize on the likelihood because then we use the full image // But we also want the prior value. // If the prior is non-flat then this is inconsistent. // BUT the main thing this does is fill in model_image with the best-fit // image rather than just the last one that the minimizer happened to try! i3_model_posterior(model,minimizer->model_image,current_best,data_set); //Some more nice output for debugging, mainly if (options->verbosity>=i3_verb_noisy) { printf("i3_analyze ML parameters: \n"); print_result(model,current_best,like,data_set); printf("\n"); } // We run an optional MCMC at the end of the optimizer. // Right now we just take the best value in its chain, // but depending on the outcome of Tomek's work we // may want to change this behaviour. bool violations = false; i3_parameter_set * mcmc_start = NULL; if (options->do_mcmc){ //If we are doing MCMC we need to check for prior violation //issues mcmc_start = malloc(model->nbytes); //Because of beermatting the ML can generate a starting position //outside the prior. Ideally we would fix this up properly so the beermatting //is done consistently in the right place but for now we just //remap the start position. The beermat function is *nearly* idempotent. i3_model_map_physical(model, current_best, mcmc_start, options); violations = i3_model_prior_violations( model, mcmc_start, stdout); if (violations) i3_print(i3_verb_quiet, "Prior violation in max-like; no MCMC"); } if (options->do_mcmc && (!violations)){ i3_mcmc * mcmc = i3_mcmc_create(model,options); // Some minimizers provide estimates of the covariance that we // can use to construct a good proposal if (minimizer->has_covariance_estimate){ // But - we need to check that the number of parameters is the same in the two systems // They may not be if some "width" parameters in the options are zero if (minimizer->nparam==mcmc->nParam){ // Just by allocating the covariance matrix in the mcmc // we force it to be used when we run below. int nb = minimizer->nparam*minimizer->nparam*sizeof(i3_flt); mcmc->cov = malloc(nb); // copy in the minimizer covmat into the mcmc memcpy(mcmc->cov, minimizer->covariance_estimate, nb); } else{ // If there is a difference in nparam, we give the user a single warning // it would be annoying to get this every time. static int have_done_warning=0; if (have_done_warning==0){ i3_print(i3_verb_quiet, "WARNING: Failed to use minimizer covariance for mcmc - widths wrong (%d,%d)", minimizer->nparam, mcmc->nParam); have_done_warning=1; } } } // This is the main MCMC run. // The details are in the options int mcmc_status = i3_mcmc_run(mcmc, mcmc_start, data_set, options); free(mcmc_start); if (mcmc_status) i3_print(i3_verb_standard,"MCMC failure! Status %d", mcmc_status); // It has proven convenient to put the acceptance on stderr for fiddling and debugging. if (options->verbosity>2)fprintf(stderr,"Acceptance = %f\n",i3_mcmc_acceptance_rate(mcmc)); // There are three ways we can extract a single sample from the MCMC to put in the catalog. if (options->use_mcmc_mean){ // WAY 1 - use the main of the MCMC samples. // Compute the mean i3_mcmc_compute_statistics(mcmc,current_best,NULL); // And then fill in the model image with the image from the mean i3_model_posterior(model,mcmc->model_image, current_best, data_set); data_set->signal_to_noise = i3_signal_to_noise_with_weight_and_data(mcmc->model_image, data_set->image, data_set->weight); data_set->old_signal_to_noise = i3_signal_to_noise_with_weight(mcmc->model_image, data_set->weight); } else if (options->use_mcmc_sample) { // WAY 2 - Just take the last MCMC element as a sample i3_model_copy_parameters(model, current_best, i3_mcmc_last_sample(mcmc)); like = i3_mcmc_last_like(mcmc); // Again, fill in the MCMC i3_model_posterior(model,mcmc->model_image, current_best, data_set); data_set->signal_to_noise = i3_signal_to_noise_with_weight_and_data(mcmc->model_image, data_set->image, data_set->weight); data_set->old_signal_to_noise = i3_signal_to_noise_with_weight(mcmc->model_image, data_set->weight); } else{ // WAY 3 - Use the MCMC like an exploratory minimizer. // This is not usually very good as that is not what MCMCs are for, especially in higher dimensions // First find the best sample int best_mcmc_sample_index = i3_mcmc_best_index(mcmc); // We may or may not find any good samples, if we are unlucky. if (best_mcmc_sample_index!=-1){ //This indicates that there is a reasonable likelihood in the //chain somewhere (hopefully this is usually true!) i3_flt best_mcmc_like = 0.5*i3_mcmc_like_number(mcmc, best_mcmc_sample_index); // We need to check if it is actually better than the minimizer sample if (best_mcmc_like > like){ // If so, celebrate and record new best parameters i3_print(i3_verb_noisy, "MCMC better (peak=%lf, mcmc=%lf, diff=%lf)\n\n\n",like, best_mcmc_like, best_mcmc_like-like); like = best_mcmc_like; i3_model_copy_parameters(model, current_best, i3_mcmc_sample_number(mcmc,best_mcmc_sample_index)); i3_model_posterior(model,mcmc->model_image, current_best, data_set); data_set->signal_to_noise = i3_signal_to_noise_with_weight_and_data(mcmc->model_image, data_set->image, data_set->weight); data_set->old_signal_to_noise = i3_signal_to_noise_with_weight(mcmc->model_image, data_set->weight); } else{ //Otherwise, just report the sad news. i3_print(i3_verb_noisy, "MCMC no better (peak=%lf, mcmc=%lf, diff=%lf)\n\n\n",like, best_mcmc_like, best_mcmc_like-like); } } } // End of the different ways of getting an MCMC sample // At last, clean up the MCMC i3_mcmc_destroy(mcmc); } // Now either from the minimizer or possibly the MCMC, we have a sample value. // We need to decide what to do with it. // These tests are really a placeholder for a proper model-independent // mapping from numerical parameters to physical ones. #ifdef HAVE_SERSICS_MODEL if (0==strcmp(model->name,"sersics")){ i3_sersics_parameter_set * beermat = i3_sersics_beermat_params(current_best,options); free(current_best); current_best = beermat; } #endif #ifdef HAVE_MULTISERSICS_MODEL if (0==strcmp(model->name,"multisersics")){ i3_multisersics_parameter_set * beermat = i3_multisersics_beermat_params(current_best, options); free(current_best); current_best=beermat; } #endif #ifdef HAVE_LOGSERSICS_MODEL if (0==strcmp(model->name,"logsersics")){ i3_logsersics_parameter_set * beermat = i3_logsersics_beermat_params(current_best,options); free(current_best); current_best = beermat; } #endif // Copy covmat if (minimizer->has_covariance_estimate){ int nb = sizeof(i3_flt) * minimizer->nparam * minimizer->nparam; data_set->covariance_estimate = malloc(nb); memcpy(data_set->covariance_estimate, minimizer->covariance_estimate, nb); } //Free up memory that we allocated. i3_image_copy_into(minimizer->model_image, best_image); i3_minimizer_destroy(minimizer); // And return two things, the parameters of the best fit, // and their likelihood *best_like = like; return current_best; } void i3_mcmc_callback_save_samples(i3_mcmc * mcmc ,i3_data_set * data_set,i3_options *options){ i3_flt L = i3_mcmc_last_like(mcmc); printf("%.12e\t",L); i3_model_printf_parameters(mcmc->model,i3_mcmc_last_sample(mcmc)); printf("\n"); }