Homec4science
Diffusion Lenstool-HPC (preparemerge) Browse

R1448/src/wrapper_bayesis.cpp
3ca39a32c40bpreparemerge

Last Change
Raw File

wrapper_bayesis.cpp
View Options

/**
* @file BayeSys3.cu
* @Author Markus Rexroth, EPFL (me@example.com)
* @date July 2015
* @version 0,1
* @brief Run Bayesian parameter inference BaySys3
*
*/
/// Include
///==========================================================================================================
#include<stdio.h>
#include<string>
#include<strings.h>
#include<stdlib.h>
#include<iostream>
#include<fstream>
#include<math.h>
#include<time.h>
#include"wrapper_bayesis.hpp"
#include<typeinfo>
#include<module_readParameters.hpp>
/// Functions
///==========================================================================================================
/** @brief Master function to call BayeSys3. In the first argument, we pass the chi2 function as a pointer for further use!
*
* Master function to call BayeSys3. In the first argument, we pass the chi2 function as a pointer for further use!
* @param *chi2_function: Pointer to chi function used to calculate chi value for actual potential, see chiClassic for more information
* @param FOLDER : Path to where results will be stored
* @param runmode : runmode parameter
* @param potential : Initial Gravitational potentials
* @param frame : Frame and Grid information
* @param potentialsToOptimize: Domain where Bayesis3 will sample to optimise potential
* @param images : Input images
* @param nImagesSet : Number of images per set
* @param BayeSys3Rate :
* @param useGPU : 1 Yes 0 NO
*/
int
runBayeSys3 (void
(*chi2_function) (double*, int*, struct runmode_param *, const struct Potential_SOA*, const struct grid_param *, const int*, galaxy*),
char* FOLDER, runmode_param *runmode, grid_param *frame, Potential_SOA* potentials, potentialoptimization *potentialsToOptimize, galaxy* images,
int * nImagesSet, double BayeSys3Rate, int useGPU)
{
//--------------- Create output files, overwrite them if they already exist
std::string fullPathBurninFile = FOLDER;
int DEBUG = 1;
fullPathBurninFile.append ("/Burnin.dat");
std::string fullPathBayesFile = FOLDER;
fullPathBayesFile.append ("/Bayes.dat");
std::ofstream burninFile;
std::ofstream bayesFile;
burninFile.open (fullPathBurninFile.c_str ()); // Open and close to overwrite them if they already exist
bayesFile.open (fullPathBayesFile.c_str ());
burninFile.close ();
bayesFile.close ();
//--------------- Declare structs required by BayeSys3
CommonStr Common[1]; // Common struct, common to all BayeSys3 objects
ObjectStr *Objects; // Pointer to object struct, we have one struct for each individual object in BayeSys3
UserCommonStr UserCommon[1]; // UserCommon struct, for additional, user-defined data
// common to all objects
BayeSys3_gpuStruct gpuStruct[1]; // Struct with variables to be put on the gpu
int numberDim = 0; // number of dimensions of the object, initialize to 0
//--------------- Parameter Space Dimensionality
BayeSys3_getNumberDim (numberDim, potentialsToOptimize, runmode->nhalos); // Get number of dimension of the object
//--------------- Initialize Common and User Common
Common->MinAtoms = 1; // We fix the number of objects ("atoms") in BayeSys3 at
Common->MaxAtoms = 1; // exactly 1
Common->Alpha = -1.0; // alpha has no effect as max = min number of atoms, but should be set
Common->ENSEMBLE = 10; // Number of trial objects ("ensemble") that BayeSys3 runs in parallel
Common->Method = -1; // Algorithm method, see BayeSys3 manual
Common->Rate = BayeSys3Rate; // Speed of calculation (dimensionless)
if (DEBUG != 1)
{
Common->Iseed = -4321; // A negative value gives a random seed as required
}
if (DEBUG == 1)
{
Common->Iseed = 4321; // A positive value gives the same seed every time,
} // use only for debugging
Common->Ndim = numberDim; // dimension = # coordinates
Common->Valency = 0; // We do not use the Massive Interference extension and disable it (=0)
Common->UserCommon = UserCommon; // Add our UserCommon struct to the Common struct, User common is used to transfer variables to the bayesis3 blackbox
UserCommon->Nsample = 0; // Initialize the number of output ensembles to 0
UserCommon->atoms = 0.0; // Initialize the average number of atoms to 0
UserCommon->Nchi2 = 0; // Initialize number of chi2 counter to 0
UserCommon->lowest_chi2 = 100000.0; // Initialize the lowest chi2 to 100,000.0
UserCommon->sum = 0.; // Initialize the accumulated sum (var/mean/mean) for all parameters to 0
UserCommon->nsets = runmode->nsets; // Initialize the number of sets of multiple images
UserCommon->nhalos = runmode->nhalos; // Initialize number of major halos
UserCommon->n_tot_halos = runmode->n_tot_halos; // Initialize number of total halos
UserCommon->referenceTime = clock (); // Initialize the timer
UserCommon->use_gpu = useGPU; // Initialize the gpu use
UserCommon->outputPath = FOLDER; // Initialize path for output files
UserCommon->runmode = runmode;
UserCommon->images = images;
UserCommon->frame = frame;
UserCommon->nImagesSet = nImagesSet;
// Initialize accumulated sum (x^2) for each parameter
UserCommon->err = (double *) malloc ((int) numberDim * sizeof(double));
// Initialize accumulated sum (x) for each parameter
UserCommon->avg = (double *) malloc ((int) numberDim * sizeof(double));
// Initialize all elements to 0
for (int ii = 0; ii < numberDim; ii++)
{
UserCommon->err[ii] = 0.;
UserCommon->avg[ii] = 0.;
}
// Allocate and copy potentials and optimisation information
UserCommon->potentialsToOptimize = (potentialoptimization *) malloc ((int) runmode->nhalos * sizeof(potentialoptimization));
#ifdef __WITH_GPU
PotentialSOAAllocation_GPU(&UserCommon->potentials, runmode.n_tot_halos);
PotentialSOAAllocation_GPU(&UserCommon->reference_potential, runmode.n_tot_halos);
#else
PotentialSOAAllocation (&UserCommon->potentials, runmode->n_tot_halos);
PotentialSOAAllocation (&UserCommon->reference_potential, runmode->n_tot_halos);
#endif
for (int ii = 0; ii < runmode->nhalos; ii++)
{
UserCommon->potentialsToOptimize[ii] = potentialsToOptimize[ii];
}
UserCommon->potentials = potentials;
UserCommon->reference_potential = potentials; // Initialize the reference potentials as the input potentials
// Pass chi2 function as pointer to UserCommon such that is will be accessible
// in later function calls
UserCommon->chi2_function = *chi2_function;
// Add our gpuStruct to UserCommon
//UserCommon->gpuStruct = gpuStruct; CS: @Gilles Not necessary for the moment, might be useful for memory prelallocation
// Objects points to the object structs of the ensemble
Objects = (ObjectStr *) malloc ((int) Common->ENSEMBLE * sizeof(ObjectStr));
//--------------- Run BayeSys3
BayeSys3 (Common, Objects);
//--------------- Return potential parameters for lowest chi2
BayeSys3_returnBestPotential (potentials, Common);
//--------------- Free the memory
std::cerr << "Freeing Memory" << std::endl;
free (UserCommon->err);
free (UserCommon->avg);
//free (UserCommon->potentials);
free (UserCommon->potentialsToOptimize);
free (Objects);
return 0;
}
/** @brief Testfunction for the chimodule . Only for debug use
* @param *chi2_function: Pointer to chi function used to calculate chi value for actual potential, see chiClassic for more information
* @param FOLDER : Path to where results will be stored
* @param runmode : runmode parameter
* @param potential : Initial Gravitational potentials
* @param frame : Frame and Grid information
* @param potentialsToOptimize: Domain where Bayesis3 will sample to optimise potential
* @param images : Input images
* @param nImagesSet : Number of images per set
* @param BayeSys3Rate :
* @param useGPU : 1 Yes 0 NO
*/
#if 0
int
runBayeSys3_testchi (void
(*chi2_function) (double*, int*, struct runmode_param *, const struct Potential*, const struct grid_param *, const int*, galaxy*),
char* FOLDER, runmode_param *runmode, grid_param *frame, Potential* potentials, potentialoptimization *potentialsToOptimize,
galaxy* images, int * nImagesSet, double BayeSys3Rate, int useGPU)
{
// Create output files, overwrite them if they already exist
std::string fullPathBurninFile = FOLDER;
int DEBUG = 1;
fullPathBurninFile.append ("/Burnin.dat");
std::string fullPathBayesFile = FOLDER;
fullPathBayesFile.append ("/Bayes.dat");
std::ofstream burninFile;
std::ofstream bayesFile;
burninFile.open (fullPathBurninFile.c_str ()); // Open and close to overwrite them if they already exist
bayesFile.open (fullPathBayesFile.c_str ());
burninFile.close ();
bayesFile.close ();
// Declare structs required by BayeSys3
CommonStr Common[1]; // Common struct, common to all BayeSys3 objects
ObjectStr *Objects; // Pointer to object struct, we have one struct for each individual object in BayeSys3
UserCommonStr UserCommon[1]; // UserCommon struct, for additional, user-defined data
// common to all objects
BayeSys3_gpuStruct gpuStruct[1]; // Struct with variables to be put on the gpu
int numberDim = 0; // number of dimensions of the object, initialize to 0
// Initialize data structs
BayeSys3_getNumberDim (numberDim, potentials, potentialsToOptimize, runmode->nhalos);
// Get number of dimension of the object
Common->MinAtoms = 1; // We fix the number of objects ("atoms") in BayeSys3 at
Common->MaxAtoms = 1; // exactly 1
Common->Alpha = -1.0; // alpha has no effect as max = min number of atoms, but should be set
Common->ENSEMBLE = 10; // Number of trial objects ("ensemble") that BayeSys3 runs in parallel
Common->Method = -1; // Algorithm method, see BayeSys3 manual
Common->Rate = BayeSys3Rate; // Speed of calculation (dimensionless)
if (DEBUG != 1)
{
Common->Iseed = -4321; // A negative value gives a random seed as required
}
if (DEBUG == 1)
{
Common->Iseed = 4321; // A positive value gives the same seed every time,
} // use only for debugging
Common->Ndim = numberDim; // dimension = # coordinates
Common->Valency = 0; // We do not use the Massive Interference extension and disable it (=0)
Common->UserCommon = UserCommon; // Add our UserCommon struct to the Common struct
UserCommon->Nsample = 0; // Initialize the number of output ensembles to 0
UserCommon->atoms = 0.0; // Initialize the average number of atoms to 0
UserCommon->Nchi2 = 0; // Initialize number of chi2 counter to 0
UserCommon->lowest_chi2 = 100000.0; // Initialize the lowest chi2 to 100,000.0
UserCommon->sum = 0.; // Initialize the accumulated sum (var/mean/mean) for all parameters to 0
UserCommon->nsets = runmode->nsets; // Initialize the number of sets of multiple images
UserCommon->nhalos = runmode->nhalos; // Initialize number of halos
UserCommon->referenceTime = clock (); // Initialize the timer
UserCommon->use_gpu = useGPU; // Initalize the gpu use
UserCommon->outputPath = FOLDER; // Initialize path for output files
// Initialize accumulated sum (x^2) for each parameter
UserCommon->err = (double *) malloc ((int) numberDim * sizeof(double));
// Initialize accumulated sum (x) for each parameter
UserCommon->avg = (double *) malloc ((int) numberDim * sizeof(double));
// Initialize all elements to 0
for (int i = 0; i < numberDim; i++)
{
UserCommon->err[i] = 0.;
UserCommon->avg[i] = 0.;
}
// Allocate memory for the potential, reference_potential and potentialsToOptimize arrays
UserCommon->potentialsToOptimize = (potentialoptimization *) malloc ((int) runmode->nhalos * sizeof(potentialoptimization));
UserCommon->potentials = (Potential *) malloc ((int) runmode->nhalos * sizeof(Potential));
UserCommon->reference_potential = (Potential *) malloc ((int) runmode->nhalos * sizeof(Potential));
// Add structs potentialsToOptimize[i], reference_potential[i] and potentials[i] to UserCommon such that they are accessible
// in the later function calls
for (int i = 0; i < runmode->nhalos; i++)
{
UserCommon->potentialsToOptimize[i] = potentialsToOptimize[i];
UserCommon->potentials[i] = potentials[i];
UserCommon->reference_potential[i] = potentials[i]; // Initialize the reference potentials as the input potentials
}
// Pass chi2 function as pointer to UserCommon such that is will be accessible
// in later function calls
UserCommon->chi2_function = *chi2_function;
// Add our gpuStruct to UserCommon
UserCommon->gpuStruct = gpuStruct;
// Objects points to the object structs of the ensemble
Objects = (ObjectStr *) malloc ((int) Common->ENSEMBLE * sizeof(ObjectStr));
printf (" %f \n", UserCommon->potentials[0].ellipticity_potential);
// Allocate GPU memory if we use the GPU
// Allocate every variable we might need, e.g.
// allocate all variables for both chi^2 in the
// source and image plane. This only takes a little more
// time but allows us to call the BayeSys with any GPU
// chi^2 function
std::cout << "Allocating memory on the GPU and CPU for Optimisation...\n" << std::endl;
int *nImagesSet_gpu;
struct galaxy *images_gpu;
double *chi2_gpu;
struct runmode_param *runmode_gpu;
struct grid_param *frame_gpu;
struct Potential *lens_gpu;
int *error_gpu;
// Allocate variables on the GPU
//cudasafe(cudaMalloc( (void**)&(gpuStruct), sizeof(BayeSys3_gpuStruct) ),"GpuStruct: ");
cudasafe (cudaMalloc ((void**) &(nImagesSet_gpu), runmode->nsets * sizeof(int)), "Alloc nImagesSet_gpu: ");
cudasafe (cudaMalloc ((void**) &(images_gpu), runmode->nimagestot * sizeof(galaxy)), "Alloc images_gpu: ");
cudasafe (cudaMalloc ((void**) &(chi2_gpu), sizeof(double)), "Alloc chi2_gpu: ");
cudasafe (cudaMalloc ((void**) &(runmode_gpu), sizeof(runmode_param)), "Alloc runmode_gpu: ");
cudasafe (cudaMalloc ((void**) &(frame_gpu), sizeof(grid_param)), "Alloc frame_gpu: ");
cudasafe (cudaMalloc ((void**) &(lens_gpu), runmode->nhalos * sizeof(Potential)), "Alloc lens_gpu: ");
cudasafe (cudaMalloc ((void**) &(error_gpu), sizeof(int)), "Alloc error_gpu: ");
// Copy fixed values to the GPU (only once copied to save computation time in the loop)
cudaMemcpy (nImagesSet_gpu, nImagesSet, runmode->nsets * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy (images_gpu, images, runmode->nimagestot * sizeof(galaxy), cudaMemcpyHostToDevice);
cudaMemcpy (runmode_gpu, runmode, sizeof(runmode_param), cudaMemcpyHostToDevice);
cudaMemcpy (frame_gpu, frame, sizeof(grid_param), cudaMemcpyHostToDevice);
std::cout << "Memory allocated\n" << std::endl;
// Define the number of threads per block the GPU will use
cudaDeviceProp properties_gpu;
cudaGetDeviceProperties (&properties_gpu, 0); // Get properties of 0th GPU in use
if (properties_gpu.maxThreadsDim[0] < threadsPerBlock)
{
fprintf (stderr, "ERROR: The GPU has to support at least %u threads per block.\n", threadsPerBlock);
exit (-1);
}
else
{
UserCommon->nBlocks_gpu = properties_gpu.maxGridSize[0] / threadsPerBlock; // Get the maximum number of blocks with the chosen number of threads
// per Block that the GPU supports
}
/*}
else
{
std::cout << "Running BayeSys3 on the CPU" << std::endl;
}*/
double chi2 = 0; // Initialize chi2 to 0 (otherwise the chi2 function won't work properly)
int error = 0; // To be done: get error value from chi2 function
//UserCommonStr* UserCommon = (UserCommonStr*)Common->UserCommon;
//BayeSys3_gpuStruct* gpuStruct = (BayeSys3_gpuStruct*)UserCommon->gpuStruct;
// Copy the new data from BayeSys3 into the GPU memory
cudaMemcpy (lens_gpu, UserCommon->potentials, (UserCommon->nhalos) * sizeof(Potential), cudaMemcpyHostToDevice);
cudaMemcpy (chi2_gpu, &chi2, sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy (error_gpu, &error, sizeof(int), cudaMemcpyHostToDevice);
// Run the chi2 function on the GPU
// Call the chi2_function specified in the BayeSys3_runBayeSys3(...) function call
// The chi2 function has to return an error value and modify the chi2 passed to it
// negative error values give an error code
UserCommon->chi2_function<<<UserCommon->nBlocks_gpu,threadsPerBlock>>>(chi2_gpu,error_gpu,runmode_gpu,lens_gpu,frame_gpu,nImagesSet_gpu,images_gpu);
// Copy the chi2 and error values form the GPU into the CPU memory
// No need to make the CPU wait for the GPU calculation to finish here,
// because this is automatically done by the cudaMemcpy
cudasafe (cudaMemcpy (&chi2, chi2_gpu, sizeof(double), cudaMemcpyDeviceToHost), "test memcpy chi2");
cudaMemcpy (&error, error_gpu, sizeof(int), cudaMemcpyDeviceToHost);
printf ("error %d, chi2 %f \n", error, chi2);
printf ("error %d, chi2 %f \n", error, chi2);
// Free the memory
free (UserCommon->err);
free (UserCommon->avg);
free (UserCommon->potentials);
free (UserCommon->potentialsToOptimize);
free (Objects);
// Free GPU memory
cudaFree (nImagesSet_gpu);
cudaFree (chi2_gpu);
cudaFree (runmode_gpu);
cudaFree (frame_gpu);
cudaFree (lens_gpu);
cudaFree (error_gpu);
cudaFree (images_gpu);
//cudaFree(gpuStruct);
return 0;
}
#endif
/** @brief Function to get the number of dimensions of the BayeSys3 object. If any additionnal parameter are added in later, this has to be modified
*
* @param numberdim : Empty variable where the result is stored
* @param potentials : Initial Gravitational potentials
* @param potentialsToOptimize: Domain where Bayesis3 will sample to optimise potential
* @param images : Number of halos (runmode_param)
*/
int
BayeSys3_getNumberDim (int &numberDim, potentialoptimization potentialsToOptimize[], int nhalos)
{
for (int i = 0; i < nhalos; i++)
{
if (potentialsToOptimize[i].position.x.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "x coordinate of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].position.y.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "y coordinate of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].weight.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Weight of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].weight.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Weight of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].ellipticity_angle.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Ellipticity angle of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].ellipticity.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Ellipticity of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].ellipticity_potential.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Ellipticity of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].rcore.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Core radius of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].rcut.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Cutoff radius of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].rscale.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Scale radius of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].exponent.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Exponent of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].alpha.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Alpha of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].einasto_kappacritic.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Einasto critical kappa of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
if (potentialsToOptimize[i].z.block == 1) // 1=true, the user wants to optimize this parameter
{
numberDim++;
fprintf (stdout, "Redshift of halo %u/%u will be optimized with BayeSys3\n", i + 1, nhalos);
}
}
fprintf (stdout, "Total number of dimensions for BayeSys3 = %u\n", numberDim);
fflush (stdout);
return 0;
}
/** @brief Function to calculate the chi2
*
* @param chi2 : Chi^2 value used to estimate the fit of the actual potential
* @param Common : Baysesis3 specific variable. Contains UserCommon and gpustruct with the information neccesarry to start CUDA chi² computation
*/
int
BayeSys3_calculateChi2 (double &chi2, CommonStr* Common)
{
chi2 = 0; // Initialize chi2 to 0 (otherwise the chi2 function won't work properly)
int error = 0; // To be done: get error value from chi2 function
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
BayeSys3_gpuStruct* gpuStruct = (BayeSys3_gpuStruct*) UserCommon->gpuStruct;
//Bayesis updates parameters in the rescale function, but does not take into account depended parameter (impact parameter for example)
//calculate potential parameters update these parameters.
printf (" calculate chi2 : ");
Potential_SOA *lenses = UserCommon->potentials;
// Call the chi2_function specified in the BayeSys3_runBayeSys3(...) function call
// The chi2 function has to return an error value and modify the chi2 passed to it
// negative error values give an error code
//std::cerr << "Potential given to chi function: "<< UserCommon->potentials->position_x[0] << UserCommon->potentials->position_y[0] << UserCommon->potentials->ellipticity[0] << std::endl;
//module_readParameters_debug_potential_SOA(1, UserCommon->potentials, UserCommon->runmode->nhalos);
//module_readParameters_debug_image(1, UserCommon->images, UserCommon->nImagesSet, UserCommon->runmode->nsets);
//std::cerr << "Frame xmin : " << UserCommon->frame->xmin << " xmax : " << UserCommon->frame->xmax << std::endl;
//mychi_bruteforce_SOA_CPU_grid_gradient_barycentersource(&chi2, &error, UserCommon->runmode, UserCommon->potentials, UserCommon->frame, UserCommon->nImagesSet, UserCommon->images);
UserCommon->chi2_function(&chi2, &error, UserCommon->runmode, UserCommon->potentials, UserCommon->frame, UserCommon->nImagesSet, UserCommon->images);
printf ("finished running\n");
printf ("Chi2 %f\n",chi2);
//printf ("error %d, chi2 %f UserCommon->lowest_chi2 %f ellipticity %f X %f Y %f rcut %f rcore %f vdisp %f \n", error, chi2, UserCommon->lowest_chi2,
// UserCommon->reference_potential[0].ellipticity, UserCommon->reference_potential[0].position.x, UserCommon->reference_potential[0].position.y,
// UserCommon->reference_potential[0].rcut, UserCommon->reference_potential[0].rcore, UserCommon->reference_potential[0].vdisp);
// if the found chi2 is a new minimum, we set the parameters
// as a future reference
if (chi2 < UserCommon->lowest_chi2)
{
printf ("New Reference error %d, chi2 %f \n", error, chi2);
UserCommon->lowest_chi2 = chi2; // Set new reference
BayeSys3_setReferencePotential (Common);
}
return error;
}
// User defined functions required by BayeSys3
//===========================================================================================================
// Declare UserBuild function
/** @brief User defined functions required by BayeSys3
* */
static int
UserBuild (double*, CommonStr*, ObjectStr*, int, int);
/** @brief User defined functions required by BayeSys3
*
* Purpose: Set Lhood = logLikelihood(coords) for empty sample
// with 0 atoms in Object, and initialise its other information.
//
// Lhood := log L(empty)
//
// I have already put the number 0 in Object->Natoms,
// so you can use a simple "UserBuild" instruction.
* */
int
UserEmpty ( // O >=0, or negative values for error code
double* Lhood, // O loglikelihood
CommonStr* Common, // I O general information
ObjectStr* Object) // I O sample object, output new Lhood
{
UserBuild (Lhood, Common, Object, 0, 1);
return 1;
}
/** @brief User defined functions required by BayeSys3
*
* Purpose: d(logLikelihood(coords)) after supposedly adding one new atom
// to Object, without any side effects on it.
//
// If Valency = 0,
// dLtry := d logL(...,x)
// else mimic the above with cool
// dLtry := (1/cool) log INTEGRAL dPrior(z) dlogL(...,x,z)
//
// I have already put the new atom x after the last location
// of the Atoms list, at Object->Cubes[ Object->Natoms ], so
// if Valency=0 you can use simple "UserBuild" instructions.
* */
int
UserTry1 ( // O +values = OK, 0 = DO NOT USE, -values = error
double* dLtry, // O trial d(logLikelihood) value
CommonStr* Common, // I general information
ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
{
int Natoms = Object->Natoms;
double Lhood = Object->Lhood; // existing Lhood
double Ltry;
int OK;
*dLtry = 0.0;
OK = UserBuild (&Ltry, Common, Object, Natoms + 1, 0); // trial
if (OK > 0)
*dLtry = Ltry - Lhood; // increment
return OK; // OK?
}
/** @brief User defined functions required by BayeSys3
* Purpose: d(logLikelihood(coords)) after supposedly adding one,
// then two new atoms to Object, without any side effects on it.
//
// If Valency = 0,
// dLtry1 := d logL(...,x1)
//
// dLtry2 := d logL(...,x1,x2)
//
// else mimic the above with cool
// dLtry1 := (1/cool) log INTEGRAL dPrior(z1) dlogL(...,x1,z1)
// cool
// dLtry2 := (1/cool) log INTEGRAL dPrior(z1,z2) dlogL(...,x1,z1,x2,z2)
//
// I have already put the new atoms x1,x2 after the last location
// of the Atoms list, at Object->Cubes[ Object->Natoms ] and
// at Object->Cubes[ Object->Natoms + 1 ]
// so if Valency=0 you can use simple "UserBuild" instructions.
//-----------------------------------------------------------------------------
*/
int
UserTry2 ( // O positive values = OK, 0 = DO NOT USE, negative values = error
double* dLtry1, // O trial d(logLikelihood) value for 1st atom
double* dLtry2, // O trial d(logLikelihood) value for both atoms
CommonStr* Common, // I general information
ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood)
{
int Natoms = Object->Natoms;
double Lhood = Object->Lhood; // existing Lhood
double Ltry1, Ltry2;
int OK;
*dLtry1 = *dLtry2 = 0.0;
OK = UserBuild (&Ltry1, Common, Object, Natoms + 1, 0); //trial for 1 more
if (OK > 0)
{
*dLtry1 = Ltry1 - Lhood; // increment for 1
OK = UserBuild (&Ltry2, Common, Object, Natoms + 2, 0); //trial for 2 more
if (OK > 0)
*dLtry2 = Ltry2 - Lhood; // increment for 2
}
return OK; // return OK only if both trials OK
}
/** @brief User defined functions required by BayeSys3
// Purpose: Insert 1 new atom into Object, keeping it up-to-date, and
// set d(loglikelihood(coords)).
//
// If Valency = 0,
// dL := d logL(...,x)
//
// else
// cool
// sample z from Prior(z) L(...,x,z)
//
// and set dL := d logL(...,x,z) at the sampled z.
//
// I have already put the new atom x at the last location of
// the updated Atoms list, at Object->Cubes[ Object->Natoms - 1 ],
//-----------------------------------------------------------------------------
*/
int
UserInsert1 ( // O >=0, or negative values for error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild (&Lnew, Common, Object, Natoms, 1); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
/** @brief User defined functions required by BayeSys3
// Purpose: Insert 2 new atoms into Object, keeping it up-to-date, and
// set d(loglikelihood(fluxes)).
//
// If Valency = 0,
// dL := d logL(...,x1,x2)
//
// else
// cool
// sample z1,z2 from Prior(z1,z2) L(...,x1,z1,x2,z2)
//
// and set dL := d logL(...,x1,z1,x2,z2) at the sampled z1,z2.
//
// I have already put the new atoms x1,x2 at the last location
// of the Atoms list, at Object->Cubes[ Object->Natoms - 2 ] and
// at Object->Cubes[ Object->Natoms - 1 ]
// so if Valency=0 you can use a simple "UserBuild" instruction.
//-----------------------------------------------------------------------------
*/
int
UserInsert2 ( // O >=0, or negative value for error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild (&Lnew, Common, Object, Natoms, 1); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
/** @brief User defined functions required by BayeSys3
// Purpose: Delete 1 old atom from Object, keeping it up-to-date, and
// set d(loglikelihood(fluxes)).
//
// dL := d logL(...)
//
// I have already put the old atom after the last location of
// the updated Atoms list, at Object->Cubes[ Object->Natoms ],
// so if Valency=0 you can use a simple "UserBuild" instruction.
//-----------------------------------------------------------------------------
*/
int
UserDelete1 ( // O >=0, or negative value for error code
double* dL, // O d(loglikelihood)
CommonStr* Common, // I general information
ObjectStr* Object) // I O sample object
{
int Natoms = Object->Natoms; // new number
double Lold = Object->Lhood; // not yet updated
double Lnew;
UserBuild (&Lnew, Common, Object, Natoms, 1); // new updated state
*dL = Lnew - Lold; // I will update Object->Lhood
return 1;
}
/** @brief User defined functions required by BayeSys3
// Purpose: Build Lhood.
//-----------------------------------------------------------------------------
*/
static int
UserBuild ( // O positive values = OK, 0 = DO NOT USE, negative values = error
double* Lhood, // O loglikelihood
CommonStr* Common, // I General information
ObjectStr* Object, // I(O) Cubes in
int Natoms, // I # atoms
int output) // I positive values if Mock to be written out, negative values if not
{
double** Cubes = Object->Cubes; // I Cubes in [0,1) [Natoms][Ndim]
double tmpCube[Common->Ndim]; // temporary Cube
int valid, error;
double chi2;
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
valid = 0; // Initialize valid to 0 (=false)
*Lhood = 0.; // Initialize the loglikelihood to 0
error = 0; // Initialize error to 0 (= false)
// Copy by value of the cube (cube is modified in BayeSys3_rescaleCube_1Atom)
for (int i = 0; i < Common->Ndim; i++)
tmpCube[i] = Cubes[0][i];
if (Natoms == 1)
{
// Set the clumps parameters from the cube and possibly reject the whole object
valid = BayeSys3_rescaleCube_1Atom (tmpCube, Common);
// Compute the likelihood for this object
if (valid)
{
UserCommon->Nchi2++; // Increase the number of computed chi^2 by 1
*Lhood = BayeSys3_calculateLogLikelihood (chi2, error, Common);
}
}
return valid;
}
/** @brief Rescale the cube values to their corresponding parameter values
// in Lenstool. A single atom contains all the clumps
// Return 1 if the rescaled parameters have significant meaning,
// 0 otherwise.
// If a parameter is added, it has to be added here
*/
int
BayeSys3_rescaleCube_1Atom (double *cube, CommonStr* Common)
{
int valid = 1; // Set valid to 1 (=true)
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
std::string stringVariable; // Intitialize string
// Rescale the potentials: Loop over all lenses and their parameters.
// For each parameter that has to be optimized, we make
// valid *= BayeSys3_rescaleParam(haloID, parameterID, &cube[ipar]);
// The value of the cube[] element is rescaled in BayeSys3_rescaleParam
// and the corresponding potentials[] element is set to the same value (this way
// we can pass the potentials struct to the GPU function instead of the cube)
int parameterCounter = 0; // Initialize parameter counter
int haloCounter = 0; // Initialize halo counter
while (parameterCounter < Common->Ndim) // Loop over all halos and their parameters
{
if (UserCommon->potentialsToOptimize[haloCounter].position.x.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "position.x";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->position_x[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].position.y.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "position.y";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].position_y[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
/*
if (UserCommon->potentialsToOptimize[haloCounter].weight.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "weight";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].weight = cube[parameterCounter];
parameterCounter++;
};*/
if (UserCommon->potentialsToOptimize[haloCounter].vdisp.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "vdisp";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->vdisp[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity_angle.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity_angle";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->ellipticity_angle[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->ellipticity[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity_potential.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity_potential";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->ellipticity_potential[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].rcore.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "rcore";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->rcore[haloCounter] = cube[parameterCounter];
//printf ("rcore %f \n", UserCommon->potentials[haloCounter].rcore);
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].rcut.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "rcut";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->rcut[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
/*
if (UserCommon->potentialsToOptimize[haloCounter].rscale.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "rscale";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].rscale = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].exponent.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "exponent";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].exponent = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].alpha.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "alpha";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].alpha = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].einasto_kappacritic.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "einasto_kappacritic";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].einasto_kappacritic = cube[parameterCounter];
parameterCounter++;
};
*/
if (UserCommon->potentialsToOptimize[haloCounter].z.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "z";
valid *= BayeSys3_rescaleParam (haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials->z[haloCounter] = cube[parameterCounter];
parameterCounter++;
};
module_readParameters_calculatePotentialparameter_SOA(UserCommon->potentials, haloCounter);
haloCounter++; //We have finished a loop over 1 halo and go to the next in the next while loop
}
// If there is no need to go further:
if (valid == 0)
return valid;
return valid;
}
// BayeSys3_rescaleParam function
/** @brief Rescale the lens parameters for each lens (using *parametervalue) and test the
// physical validity of the rescaled values
// Return 1 if it's ok, 0 if it's not physically valid
*/
int
BayeSys3_rescaleParam (int haloID, std::string parameterID, double *parameterValue, CommonStr* Common)
{
double value = *parameterValue;
int valid = 1; // Intitialize to 1 (=true)
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
// Loop over all possible parameters, if the matching parameter is found
// we make
// value = BayeSys3_prior( Common, value, min, max );
// For certain parameters, e.g. r_cut, we check that they make sense, e.g.
// if r_cut < 0 -> valid = 0 (= false)
if (parameterID == "position.x")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].position.x.min,
UserCommon->potentialsToOptimize[haloID].position.x.max);
}
if (parameterID == "position.y")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].position.y.min,
UserCommon->potentialsToOptimize[haloID].position.y.max);
}
if (parameterID == "weight")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].weight.min, UserCommon->potentialsToOptimize[haloID].weight.max);
}
if (parameterID == "vdisp")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].vdisp.min, UserCommon->potentialsToOptimize[haloID].vdisp.max);
}
if (parameterID == "ellipticity_angle")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].ellipticity_angle.min,
UserCommon->potentialsToOptimize[haloID].ellipticity_angle.max);
}
if (parameterID == "ellipticity")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].ellipticity.min,
UserCommon->potentialsToOptimize[haloID].ellipticity.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "ellipticity_potential")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].ellipticity_potential.min,
UserCommon->potentialsToOptimize[haloID].ellipticity_potential.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "rcore")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].rcore.min, UserCommon->potentialsToOptimize[haloID].rcore.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "rcut")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].rcut.min, UserCommon->potentialsToOptimize[haloID].rcut.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "rscale")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].rscale.min, UserCommon->potentialsToOptimize[haloID].rscale.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "exponent")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].exponent.min,
UserCommon->potentialsToOptimize[haloID].exponent.max);
}
if (parameterID == "alpha")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].alpha.min, UserCommon->potentialsToOptimize[haloID].alpha.max);
}
if (parameterID == "einasto_kappacritic")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].einasto_kappacritic.min,
UserCommon->potentialsToOptimize[haloID].einasto_kappacritic.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "z")
{
value = BayeSys3_prior (Common, value, UserCommon->potentialsToOptimize[haloID].z.min, UserCommon->potentialsToOptimize[haloID].z.max);
if (value < 0)
{
valid = 0;
}
}
return valid;
}
/** @brief prior function
// value is a deviate from the unit Cube.
// Return the corresponding deviate drawn from the desired distribution
*/
double
BayeSys3_prior (CommonStr* Common, double parameterValue, double min, double max)
{
// To be done: get the switchIndex from an additional parameter in parametersoptimize: .distribution (has yet to be added to the struct)
// Now we set it always to a default of 1 (= uniform prior)
int switchIndex = 1;
switch (switchIndex)
{
case (1):
return BayeSys3_uniformPrior (parameterValue, min, max);
// To be done: implement other priors, e.g. Gaussian
default:
fprintf (stderr, "ERROR : Prior switchIndex %d is undefined.\n", switchIndex);
exit (-1);
}
return 0;
}
// BayeSys3_uniformPrior function
/** @brief Rescale Uniform [0:1] -> Uniform[x1:x2]
*/
double
BayeSys3_uniformPrior (double val, double min, double max)
{
return min + val * (max - min);
}
/** @brief function to calculate the log-likelihood
* */
double
BayeSys3_calculateLogLikelihood (double &chi2, int &error, CommonStr* Common)
{
double Lhood;
error = BayeSys3_calculateChi2 (chi2, Common);
Lhood = -0.5 * chi2;
return (Lhood);
}
//=============================================================================
// Dummy procedure to link to BayeSys3 when not running MassInf
//=============================================================================
int
UserFoot (double* Cube, CommonStr* Common, int* ibits, double* zbits, int* nbits)
{
return (Cube && Common && ibits && zbits && nbits) ? 1 : 1;
}
//=============================================================================
// Function: UserMonitor
//=============================================================================
// Purpose: There is a new ensemble of objects.
//
// 1. For probabilistic exploration, restrict Common->cool to <= 1.
// cool = 0 is prior, seen at the first call.
// 0 < cool < 1 is "burn-in"
// cool = 1 is evolve posterior,
// cool -> infinity is maximum likelihood
// Return with a POSITIVE return code to exit BayeSys,
// usually after the end of an evolution stage with cool = 1.
// You can use Common->Nsystem, which counts calls to UserMonitor.
//
// 2. You should reset any nuisance parameters x in UserCommon and/or
// each of your UserObject structures, preferably by sampling from
// their conditonal probabilities Pr(x|Objects) and/or Pr(x|Object).
//
// 3. Collect your statistics and display your diagnostics.
//-----------------------------------------------------------------------------
//
int
UserMonitor ( // O 0 = continue, positive values = finish, negative values = abort
CommonStr* Common, // I Full general information
ObjectStr* Objects) // I Full ensemble of sample objects [ENSEMBLE]
{
double tmpCube[Common->Ndim]; // temporary cube for calculations
int NumberRemaining; // Number of iterations to go until finish
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
ObjectStr* Object; // individual object
double chi2; // Chi^2
double logLikelihood;
int errorLogLikelihood;
clock_t currentTime;
clock_t timeDifference;
// Uncomment if we use the save lens parameters feature below
//long int totalNumberSamples;
//.............................................................................
// USER IMPOSES FINAL TEMPERATURE AND TERMINATION CONDITION !
// To be done: if we set in the parameter file inverse mode = 3, we want MCMC, not Bayesian max likelihood.
// In this case, we make:
// if ( Common->cool >= 1.0 && inverse == 3 )
// Common->cool = 1.0;
//.............................................................................
// Save the models in burnin.dat/bayes.dat
std::ofstream printFile;
std::string outputFile = UserCommon->outputPath;
if (Common->cool < 1.)
{
outputFile.append ("/Burnin.dat");
printFile.open (outputFile.c_str (), std::ios::app);
}
else
{
outputFile.append ("Bayes.dat");
printFile.open (outputFile.c_str (), std::ios::app);
}
UserCommon->sum = 0;
for (int j = 0; j < Common->ENSEMBLE; j++)
{
// Rescale the cube of object j
// The cube is changing during the run of BayeSys3 to
// reflect the new posterior, which is done automatically by BayeSys3
// We only rescale it to fit to our parameter range
for (int i = 0; i < Common->Ndim; i++)
{
tmpCube[i] = Objects[j].Cubes[0][i];
}
BayeSys3_rescaleCube_1Atom (tmpCube, Common);
// Take a test chi2 from our current sample distribution "Nsample"
// Save the chi2 and Lhood
UserCommon->Nchi2++;
logLikelihood = BayeSys3_calculateLogLikelihood (chi2, errorLogLikelihood, Common);
printFile << "Nsample = " << UserCommon->Nsample << ", logLikelihood = " << logLikelihood << "\n";
// TO BE DONE: Include statistics from addStat(val, UserCommon, param++, samp) and output it;
// append the evidence
if (Common->cool < 1.)
printFile << "Evidence = " << Common->Evidence << "\n";
//append Chi2 in last column in bayes.dat
if (Common->cool >= 1.)
printFile << "Chi^2 = " << chi2 << "\n";
}
printFile.close ();
//.............................................................................
// Print the status line to STDOUT
if (Common->cool < 1.)
{
NumberRemaining = Common->Nsystem; // NumberRemaining first increases with the number of iterations Nsystem
fprintf (stdout, "\rBurn-in: cool = %f, Nsystem = %u, chi2 = %f, Evidence = %f\n", Common->cool, NumberRemaining, chi2, Common->Evidence);
}
else
{
// To be done: If we manually set the max number of iterations in the parameter file
// use this value here = maxNumberIterations
int maxNumberIterations = Common->Nsystem;
NumberRemaining = maxNumberIterations - 2 * UserCommon->Nsample; // We sample until the number of samples == number of annealing iterations,
// as recommended in the BayeSys3 manual, we need the factor 2 because
// Nsystem continues to increase during the sampling phase, it is not
// fixed after the burn-in
UserCommon->sum /= Common->ENSEMBLE * Common->Ndim;
fprintf (stdout, "\rSampling: cool = %f, Finished (%) = %d/%d, chi2 = %f, sqrt(sum) = %f\n", Common->cool, UserCommon->Nsample + 1,
(maxNumberIterations - UserCommon->Nsample) + 1, chi2, sqrt (UserCommon->sum)); // We print maxNumberIterations - Nsample +1
// because this is the constant, total number of
// iterations we need to finish
}
// Print a performance measurement (Number chi2/sec)
currentTime = clock (); // get current time
timeDifference = currentTime - UserCommon->referenceTime;
fprintf (stdout, "Currently we have %lf chi2/sec\n", UserCommon->Nchi2 / (((float) timeDifference) / CLOCKS_PER_SEC)); // We have 1000 chi2 since the last time
// we called this function
UserCommon->referenceTime = clock (); // make current time reference time
UserCommon->Nchi2 = 0; // Reset number of chi2 counter
fflush (stdout); // Make sure the text is printed to stdout
//.............................................................................
// Set any NUISANCE PARAMETERS you may have here
// Could over-ride system's choices of FluxUnit (both Common and Objects) here
// Could CALIBRATE data here (corrupting Evidence and Information)
//.............................................................................
// BURN-IN
if (Common->cool < 1.0) // final temperature not reached..
return 0; // ...normal return, keep going
//.............................................................................
// EVOLUTION: Accumulate any desired statistics
for (int k = 0; k < Common->ENSEMBLE; k++)
{
Object = &Objects[k];
// # atoms
UserCommon->atoms += Object->Natoms;
}
// # samples of full ensemble
UserCommon->Nsample++;
if (NumberRemaining > 0) // not finished...
return 0; // ..normal return
//.............................................................................
// EXIT WITH POSITIVE FINISH CODE WHEN EVOLUTION IS DEEMED COMPLETE
fprintf (stdout, "\nSampling complete, exiting BayeSys3 \n");
fprintf (stdout, "Returning lens model with lowest chi2 = %lf\n\n", UserCommon->lowest_chi2);
return 1;
}
// BayeSys3_setReferencePotential function
//
/** @brief Set the lens parameters of the current lowest chi2 as a future reference
*/
int
BayeSys3_setReferencePotential (CommonStr* Common)
{
int valid = 1; // Set valid to 1 (=true)
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
for (int haloNumber = 0; haloNumber < UserCommon->nhalos; haloNumber++) // Loop over all halos and their parameters
{
UserCommon->reference_potential->position_x[haloNumber] = UserCommon->potentials->position_x[haloNumber];
UserCommon->reference_potential->position_y[haloNumber] = UserCommon->potentials->position_y[haloNumber];
//UserCommon->reference_potential->weight[haloNumber] = UserCommon->potentials->weight[haloNumber];
UserCommon->reference_potential->vdisp[haloNumber] = UserCommon->potentials->vdisp[haloNumber];
UserCommon->reference_potential->ellipticity_angle[haloNumber] = UserCommon->potentials->ellipticity_angle[haloNumber];
UserCommon->reference_potential->ellipticity[haloNumber] = UserCommon->potentials->ellipticity[haloNumber];
UserCommon->reference_potential->ellipticity_potential[haloNumber] = UserCommon->potentials->ellipticity_potential[haloNumber];
UserCommon->reference_potential->rcore[haloNumber] = UserCommon->potentials->rcore[haloNumber];
UserCommon->reference_potential->rcut[haloNumber] = UserCommon->potentials->rcut[haloNumber];
//UserCommon->reference_potential->rscale[haloNumber] = UserCommon->potentials->rscale[haloNumber];
//UserCommon->reference_potential->exponent[haloNumber] = UserCommon->potentials->exponent[haloNumber];
//UserCommon->reference_potential->alpha[haloNumber] = UserCommon->potentials->alpha[haloNumber];
//UserCommon->reference_potential->einasto_kappacritic[haloNumber] = UserCommon->potentials->einasto_kappacritic[haloNumber];
UserCommon->reference_potential->z[haloNumber] = UserCommon->potentials->z[haloNumber];
//never forget to calculate the parameter depended characteristic of the potential
//readParameters_calculatePotentialparameter(&UserCommon->reference_potential[haloNumber]);
module_readParameters_calculatePotentialparameter_SOA(UserCommon->reference_potential, haloNumber);
}
return valid;
}
/** @brief function to update the statistics
*/
int
BayeSys3_addStatistics (double value, UserCommonStr *UserCommon, int param, long int totalNumberSamples)
{
double var, mean;
UserCommon->err[param] += value * value;
UserCommon->avg[param] += value;
mean = UserCommon->avg[param] / totalNumberSamples;
var = UserCommon->err[param] / totalNumberSamples - mean * mean;
UserCommon->sum += var / mean / mean;
return 0;
}
// BayeSys3_returnBestPotential function
//
/** @brief Set the lens parameters of the lowest chi2 as output potential parameters
*/
int
BayeSys3_returnBestPotential (Potential_SOA* potentials, CommonStr* Common)
{
UserCommonStr* UserCommon = (UserCommonStr*) Common->UserCommon;
printf ("finished: returning best potential\n");
for (int haloNumber = 0; haloNumber < UserCommon->nhalos; haloNumber++) // Loop over all halos and their parameters
{
potentials->position_x[haloNumber] = UserCommon->reference_potential->position_x[haloNumber];
potentials->position_y[haloNumber] = UserCommon->reference_potential->position_y[haloNumber];
//potentials->weight[haloNumber] = UserCommon->reference_potential->weight[haloNumber];
potentials->vdisp[haloNumber] = UserCommon->reference_potential->vdisp[haloNumber];
potentials->ellipticity_angle[haloNumber] = UserCommon->reference_potential->ellipticity_angle[haloNumber];
potentials->ellipticity[haloNumber] = UserCommon->reference_potential->ellipticity[haloNumber];
potentials->ellipticity_potential[haloNumber] = UserCommon->reference_potential->ellipticity_potential[haloNumber];
potentials->rcore[haloNumber] = UserCommon->reference_potential->rcore[haloNumber];
potentials->rcut[haloNumber] = UserCommon->reference_potential->rcut[haloNumber];
//potentials->rscale[haloNumber] = UserCommon->reference_potential->rscale[haloNumber];
//potentials->exponent[haloNumber] = UserCommon->reference_potential->exponent[haloNumber];
//potentials->alpha[haloNumber] = UserCommon->reference_potential->alpha[haloNumber];
//potentials->einasto_kappacritic[haloNumber] = UserCommon->reference_potential->einasto_kappacritic[haloNumber];
potentials->z[haloNumber] = UserCommon->reference_potential->z[haloNumber];
//module_readParameters_calculatePotentialparameter_SOA (potentials, haloNumber);
}
printf ("returned best potential\n");
return 0; // We are no longer in the BayeSys3, so we return 0, not a "valid" code
}
c4science · Help