Skip to content

Commit

Permalink
Added files for do_multimfit_tests.sh; added missing files for multim…
Browse files Browse the repository at this point in the history 
…fit compilation.
  • Loading branch information
@perwin
perwin committed Jun 23, 2024
1 parent 0644cc1 commit 5176524
Show file tree
Hide file tree
Showing 39 changed files with 1,732 additions and 52 deletions.
15 changes: 9 additions & 6 deletions SConstruct
View file
Original file line number Diff line number Diff line change
Expand Up @@ -666,9 +666,8 @@ modelobject_sources = [name + ".cpp" for name in modelobject_objs]
functionobject_obj_string = """function_object func_gaussian func_exp func_gen-exp
func_sersic func_gen-sersic func_core-sersic func_broken-exp
func_broken-exp2d func_moffat func_flatsky func_tilted-sky-plane
func_flatbar func_gaussian-ring
func_gaussian-ring2side func_gaussian-ring-az func_edge-on-disk_n4762
func_edge-on-disk_n4762v2 func_edge-on-ring func_edge-on-ring2side
func_flatbar func_gaussian-ring func_gaussian-ring2side func_gaussian-ring-az
func_edge-on-ring func_edge-on-ring2side
func_king func_king2 func_ferrersbar2d func_peanut_dattathri
helper_funcs helper_funcs_3d psf_interpolators"""
#if useGSL:
Expand Down Expand Up @@ -705,6 +704,8 @@ if useExtraFuncs:
functionobject_obj_string += " func_logspiral3"
functionobject_obj_string += " func_logspiral_gauss"
functionobject_obj_string += " func_logspiral_arc"
func_edge-on-disk_n4762
func_edge-on-disk_n4762v2
functionobject_obj_string += " func_polynomial_d1"
functionobject_obj_string += " func_triaxbar3d"
functionobject_obj_string += " func_triaxbar3d_sq"
Expand Down Expand Up @@ -821,6 +822,7 @@ if scanBuild:
env_debug["ENV"].update(x for x in os.environ.items() if x[0].startswith("CCC_"))



# *** Finally, define the actual targets for building
# specify ".do" as the suffix for "full-debug" object code
imfit_dbg_objlist = [ env_debug.Object(obj + ".do", src) for (obj,src) in zip(imfit_objs, imfit_sources) ]
Expand All @@ -839,8 +841,8 @@ env.Program("imfit-mcmc", mcmc_sources)
multi_objlist = [ env.Object(obj, src) for (obj,src) in zip(makemultimages_objs, makemultimages_sources) ]
env.Program("makemultimages", multi_objlist)

# multimfit_objlist = [ env_multi.Object(obj, src) for (obj,src) in zip(multimfit_objs, multimfit_sources) ]
# env_multi.Program("multimfit", multimfit_objlist)
multimfit_objlist = [ env.Object(obj, src) for (obj,src) in zip(multimfit_objs, multimfit_sources) ]
env.Program("multimfit", multimfit_objlist)


# Run tests
Expand All @@ -853,7 +855,6 @@ env.Command("alltests", None,




# Build Imfit library
base_for_lib_objstring = """mp_enorm statistics mersenne_twister utilities
config_file_parser add_functions bootstrap_errors count_cpu_cores"""
Expand Down Expand Up @@ -926,6 +927,8 @@ psfconvolve1d_sources = [name + ".cpp" for name in psfconvolve1d_objs]





# source+obj lists for older or less-used programs:
# readimage: put all the object and source-code lists together
readimage_sources = ["readimage_main.cpp", "image_io.cpp"]
Expand Down
170 changes: 128 additions & 42 deletions core/bootstrap_errors.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
* nonlinfit (imfit's conceptual predecessor), so yay for reuse!
*/

// Copyright 2013-2024 by Peter Erwin.
// Copyright 2013-2019 by Peter Erwin.
//
// This file is part of Imfit.
//
Expand Down Expand Up @@ -38,6 +38,7 @@
#include <tuple>

#include "definitions.h"
#include "definitions_multimage.h"
#include "model_object.h"
#include "levmar_fit.h"
#ifndef NO_NLOPT
Expand All @@ -56,33 +57,40 @@ using namespace std;
const int MIN_ITERATIONS_FOR_STATISTICS = 3;
const int PROGRESS_BAR_WIDTH = 80;

const vector<string> imageParamLabels = {"PIXEL_SCALE", "ROTATION", "FLUX_SCALE",
"X0", "Y0"};


/* ------------------- Function Prototypes ----------------------------- */

int BootstrapErrorsBase( const double *bestfitParams, vector<mp_par> parameterLimits,
const bool paramLimitsExist, ModelObject *theModel, const double ftol,
const int nIterations, const int nFreeParams, const int whichStatistic,
double **outputParamArray, FILE *outputFile_ptr, unsigned long rngSeed=0 );
void PrintBootstrapSummary( const double *bestfitParams, double **outputParamArray, int nParams,
int nIterations, vector<mp_par> parameterLimits, ModelObject *theModel );
void PrintBootstrapSummaryMultimfit( const double *bestfitParams, double **outputParamArray,
int nParams, int nIterations, vector<mp_par> parameterLimits,
ModelObject *theModel );




/* ---------------- FUNCTION: BootstrapErrors -------------------------- */
/// Primary wrapper function (meant to be called from main() in imfit_main.cpp, etc.).
/// If saving of all best-fit parameters to file is requested, then outputFile_ptr
/// should be non-nullptr (i.e., should point to a file object opened for writing, possibly
/// should be non-NULL (i.e., should point to a file object opened for writing, possibly
/// with header information already written).
/// Returns the number of (successful) bootstrap iterations, or returns -1 if error
/// encountered.
int BootstrapErrors( const double *bestfitParams, vector<mp_par> parameterLimits,
const bool paramLimitsExist, ModelObject *theModel, const double ftol,
const int nIterations, const int nFreeParams, const int whichStatistic,
FILE *outputFile_ptr, unsigned long rngSeed )
FILE *outputFile_ptr, unsigned long rngSeed, bool multimfitMode )
{
double *paramSigmas;
double *bestfitParams_offsetCorrected, *paramOffsets;
double **outputParamArray;
double lower, upper, plus, minus, halfwidth;
// double lower, upper, plus, minus, halfwidth;
int i, nSuccessfulIterations;
int nParams = theModel->GetNParams();

Expand All @@ -96,11 +104,11 @@ int BootstrapErrors( const double *bestfitParams, vector<mp_par> parameterLimits
for (i = 0; i < nParams; i++)
outputParamArray[i] = (double *)calloc( (size_t)nIterations, sizeof(double) );

paramOffsets = (double *) calloc(nParams, sizeof(double));
bestfitParams_offsetCorrected = (double *) calloc(nParams, sizeof(double));
// paramOffsets = (double *) calloc(nParams, sizeof(double));
// bestfitParams_offsetCorrected = (double *) calloc(nParams, sizeof(double));

// write column header info to file, if user requested saving to file
if (outputFile_ptr != nullptr) {
if (outputFile_ptr != NULL) {
string headerLine = theModel->GetParamHeader();
fprintf(outputFile_ptr, "#\n# Bootstrap resampling output (%d iterations requested):\n%s\n",
nIterations, headerLine.c_str());
Expand All @@ -117,44 +125,25 @@ int BootstrapErrors( const double *bestfitParams, vector<mp_par> parameterLimits
else {
// Apply image-offset corrections to best-fit parameters, so their values
// get printed corrected
theModel->GetImageOffsets(paramOffsets);
for (i = 0; i < nParams; i++)
bestfitParams_offsetCorrected[i] = bestfitParams[i] + paramOffsets[i];
// theModel->GetImageOffsets(paramOffsets);
// for (i = 0; i < nParams; i++)
// bestfitParams_offsetCorrected[i] = bestfitParams[i] + paramOffsets[i];

// Calculate sigmas and 68% confidence intervals for the parameters
// vector to hold estimated sigmas for each parameter
paramSigmas = (double *)calloc( (size_t)nParams, sizeof(double) );
for (i = 0; i < nParams; i++)
paramSigmas[i] = StandardDeviation(outputParamArray[i], nSuccessfulIterations);
// Print parameter values + standard deviations, for non-fixed parameters
// (note that calling ConfidenceInterval() sorts the vectors in place!)
printf("\nStatistics for parameter values from bootstrap resampling");
printf(" (%d successful iterations):\n", nSuccessfulIterations);
printf("Best-fit\t\t Bootstrap [68%% conf.int., half-width]; (mean +/- standard deviation)\n");
for (i = 0; i < nParams; i++) {
if (parameterLimits[i].fixed == 0) {
std::tie(lower, upper) = ConfidenceInterval(outputParamArray[i], nSuccessfulIterations);
plus = upper - bestfitParams_offsetCorrected[i];
minus = bestfitParams_offsetCorrected[i] - lower;
halfwidth = (upper - lower)/2.0;
printf("%s = %g +%g, -%g [%g -- %g, %g]; (%g +/- %g)\n",
theModel->GetParameterName(i).c_str(),
bestfitParams_offsetCorrected[i], plus, minus, lower, upper, halfwidth,
Mean(outputParamArray[i], nSuccessfulIterations), paramSigmas[i]);
}
else {
printf("%s = %g [fixed parameter]\n", theModel->GetParameterName(i).c_str(),
bestfitParams_offsetCorrected[i]);
}
}
free(paramSigmas);
if (multimfitMode)
PrintBootstrapSummaryMultimfit(bestfitParams, outputParamArray, nParams, nIterations,
parameterLimits, theModel);
else
PrintBootstrapSummary(bestfitParams, outputParamArray, nParams, nIterations,
parameterLimits, theModel);
}

for (i = 0; i < nParams; i++)
free(outputParamArray[i]);
free(outputParamArray);
free(paramOffsets);
free(bestfitParams_offsetCorrected);
// free(paramOffsets);
// free(bestfitParams_offsetCorrected);

return nSuccessfulIterations;
}
Expand All @@ -164,7 +153,7 @@ int BootstrapErrors( const double *bestfitParams, vector<mp_par> parameterLimits
/* ---------------- FUNCTION: BootstrapErrorsBase ---------------------- */
/// Base function called by the wrapper functions (above), which does the main work
/// of overseeing the bootstrap resampling.
/// Saving individual best-fit vales to file is done *if* outputFile_ptr != nullptr.
/// Saving individual best-fit vales to file is done *if* outputFile_ptr != NULL.
/// Returns the number of successful iterations performed (-1 if an error was
/// encountered)
int BootstrapErrorsBase( const double *bestfitParams, vector<mp_par> parameterLimits,
Expand All @@ -180,13 +169,13 @@ int BootstrapErrorsBase( const double *bestfitParams, vector<mp_par> parameterLi
bool saveToFile = false;
string outputLine, iterTemplate;

if (outputFile_ptr != nullptr)
if (outputFile_ptr != NULL)
saveToFile = true;

if (rngSeed > 0)
init_genrand(rngSeed);
else
init_genrand((unsigned long)time((time_t *)nullptr));
init_genrand((unsigned long)time((time_t *)NULL));

paramsVect = (double *) calloc(nParams, sizeof(double));
paramOffsets = (double *) calloc(nParams, sizeof(double));
Expand Down Expand Up @@ -282,7 +271,7 @@ int BootstrapErrorsArrayOnly( const double *bestfitParams, vector<mp_par> parame
if (rngSeed > 0)
init_genrand(rngSeed);
else
init_genrand((unsigned long)time((time_t *)nullptr));
init_genrand((unsigned long)time((time_t *)NULL));

paramsVect = (double *) calloc(nParams, sizeof(double));
paramOffsets = (double *) calloc(nParams, sizeof(double));
Expand Down Expand Up @@ -347,5 +336,102 @@ int BootstrapErrorsArrayOnly( const double *bestfitParams, vector<mp_par> parame
}


void PrintBootstrapSummary( const double *bestfitParams, double **outputParamArray, int nParams,
int nIterations, vector<mp_par> parameterLimits, ModelObject *theModel )
{
double lower, upper, plus, minus, halfwidth;
double *paramSigmas;
int i;

paramSigmas = (double *)calloc( (size_t)nParams, sizeof(double) );
for (i = 0; i < nParams; i++)
paramSigmas[i] = StandardDeviation(outputParamArray[i], nIterations);
// Print parameter values + standard deviations, for non-fixed parameters
// (note that calling ConfidenceInterval() sorts the vectors in place!)
printf("\nStatistics for parameter values from bootstrap resampling");
printf(" (%d successful iterations):\n", nIterations);
printf("Best-fit\t\t Bootstrap [68%% conf.int., half-width]; (mean +/- standard deviation)\n");
for (i = 0; i < nParams; i++) {
if (parameterLimits[i].fixed == 0) {
std::tie(lower, upper) = ConfidenceInterval(outputParamArray[i], nIterations);
plus = upper - bestfitParams[i];
minus = bestfitParams[i] - lower;
halfwidth = (upper - lower)/2.0;
printf("%s = %g +%g, -%g [%g -- %g, %g]; (%g +/- %g)\n",
theModel->GetParameterName(i).c_str(),
bestfitParams[i], plus, minus, lower, upper, halfwidth,
Mean(outputParamArray[i], nIterations), paramSigmas[i]);
}
else {
printf("%s = %g [fixed parameter]\n", theModel->GetParameterName(i).c_str(),
bestfitParams[i]);
}
}
free(paramSigmas);
}



void PrintBootstrapSummaryMultimfit( const double *bestfitParams, double **outputParamArray,
int nParams, int nIterations, vector<mp_par> parameterLimits,
ModelObject *theModel )
{
double lower, upper, plus, minus, halfwidth;
double *paramSigmas;
int i, i_last, nImages;

nImages = theModel->GetNImages();

paramSigmas = (double *)calloc( (size_t)nParams, sizeof(double) );
for (i = 0; i < nParams; i++)
paramSigmas[i] = StandardDeviation(outputParamArray[i], nIterations);
// Print parameter values + standard deviations, for non-fixed parameters
// (note that calling ConfidenceInterval() sorts the vectors in place!)
printf("\nStatistics for parameter values from bootstrap resampling");
printf(" (%d successful iterations):\n", nIterations);
printf("Best-fit\t\t Bootstrap [68%% conf.int., half-width]; (mean +/- standard deviation)\n");

i = -1;
// print results for image-description parameters (for 2nd and subsequent images)
for (int n = 1; n < nImages; n++) {
printf("# Image-description parameters for image %d:\n", n + 1);
for (int j = 0; j < N_IMAGE_PARAMS; j++) {
i += 1;
if (parameterLimits[i].fixed == 0) {
std::tie(lower, upper) = ConfidenceInterval(outputParamArray[i], nIterations);
plus = upper - bestfitParams[i];
minus = bestfitParams[i] - lower;
halfwidth = (upper - lower)/2.0;
printf("%s = %g +%g, -%g [%g -- %g, %g]; (%g +/- %g)\n", imageParamLabels[i].c_str(),
bestfitParams[i], plus, minus, lower, upper, halfwidth,
Mean(outputParamArray[i], nIterations), paramSigmas[i]);
}
else
printf("%s = %g [fixed parameter]\n", imageParamLabels[i].c_str(),
bestfitParams[i]);
}
}

// now print results for model parameters
i_last = i;
printf("# Model parameters:\n");
for (i = i_last; i < nParams; i++) {
if (parameterLimits[i].fixed == 0) {
std::tie(lower, upper) = ConfidenceInterval(outputParamArray[i], nIterations);
plus = upper - bestfitParams[i];
minus = bestfitParams[i] - lower;
halfwidth = (upper - lower)/2.0;
printf("%s = %g +%g, -%g [%g -- %g, %g]; (%g +/- %g)\n",
theModel->GetParameterName(i).c_str(), bestfitParams[i], plus, minus,
lower, upper, halfwidth, Mean(outputParamArray[i], nIterations), paramSigmas[i]);
}
else
printf("%s = %g [fixed parameter]\n", theModel->GetParameterName(i).c_str(),
bestfitParams[i]);
}

free(paramSigmas);
}


/* END OF FILE: bootstrap_errors.cpp ----------------------------------- */
2 changes: 1 addition & 1 deletion core/bootstrap_errors.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
int BootstrapErrors( const double *bestfitParams, vector<mp_par> parameterLimits,
const bool paramLimitsExist, ModelObject *theModel, const double ftol,
const int nIterations, const int nFreeParams, const int whichStatistic,
FILE *outputFile_ptr, unsigned long rngSeed=0 );
FILE *outputFile_ptr, unsigned long rngSeed=0, bool multimfitMode=false );


// NOTE: The following function is used in PyImfit
Expand Down
Loading

0 comments on commit 5176524

Please sign in to comment.