Homec4science
Diffusion Lenstool-HPC (BayesImplementation) Browse

R1448/Benchmarks/BayesBenchmark/bayes_app.cpp
c34e73e12be4BayesImplementation

Last Change
Raw File

bayes_app.cpp
View Options

/**
* @file module_BayeSys3.cu
* @Author Markus Rexroth, EPFL; Christoph Schaefer "christophernstrerne.schaefer@epfl.ch" EPFL
* @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"bayesys3.h"
#include"userstr.h"
#include"bayes_app.hpp"
#include<typeinfo>
/// 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 module_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 module_BayeSys3_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
#if 0
BayeSys3_gpuStruct gpuStruct[1]; // Struct with variables to be put on the gpu
#endif
int numberDim = 0; // number of dimensions of the object, initialize to 0
// Initialize data structs
module_BayeSys3_getNumberDim(numberDim, 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_param *)malloc((int) runmode->nhalos * sizeof(potential_param));
UserCommon->reference_potential = (potential_param *)malloc((int) runmode->nhalos * sizeof(potential_param));
// 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;
#if 0
// Add our gpuStruct to UserCommon
UserCommon->gpuStruct = gpuStruct;
#endif
// Objects points to the object structs of the ensemble
Objects = (ObjectStr *)malloc((int) Common->ENSEMBLE * sizeof(ObjectStr));
#if 0
// 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
// Allocate the memory on the GPU only if we use a GPU
if (UserCommon->use_gpu == 1)
{
std::cout << "\nAllocating memory on the GPU and CPU..." << std::endl;
// Allocate memory on the CPU memory
gpuStruct->nImagesSet_gpu = (int *)malloc((int) runmode->nsets * sizeof(int));
gpuStruct->chi2_gpu = (double *)malloc(sizeof(double));
gpuStruct->images_gpu = (galaxy *)malloc( runmode->nimagestot* sizeof(galaxy));
gpuStruct->runmode_gpu = (runmode_param *)malloc((int) sizeof(runmode_param));
gpuStruct->frame_gpu = (grid_param *)malloc((int) sizeof(grid_param));
gpuStruct->lens_gpu = (potential_param *)malloc((int) runmode->nhalos * sizeof(potential_param));
gpuStruct->error_gpu = (int *)malloc(sizeof(int));
// Allocate variables on the GPU
cudasafe(cudaMalloc( (void**)&(gpuStruct), sizeof(BayeSys3_gpuStruct)),"Alloc nImagesSet_gpu: " );
cudasafe(cudaMalloc( (void**)&(gpuStruct->nImagesSet_gpu), runmode->nsets*sizeof(int)),"Alloc nImagesSet_gpu: " );
cudasafe(cudaMalloc( (void**)&(gpuStruct->images_gpu), runmode->nimagestot*sizeof(galaxy)),"Alloc images_gpu: " );
cudasafe(cudaMalloc( (void**)&(gpuStruct->chi2_gpu), sizeof(double) ),"Alloc chi2_gpu: ");
cudasafe(cudaMalloc( (void**)&(gpuStruct->runmode_gpu), sizeof(runmode_param) ),"Alloc runmode_gpu: ");
cudasafe(cudaMalloc( (void**)&(gpuStruct->frame_gpu), sizeof(grid_param) ),"Alloc frame_gpu: ");
cudasafe(cudaMalloc( (void**)&(gpuStruct->lens_gpu), runmode->nhalos * sizeof(potential_param) ),"Alloc lens_gpu: ");
cudasafe(cudaMalloc( (void**)&(gpuStruct->error_gpu), sizeof(int) ),"Alloc error_gpu: ");
// Copy fixed values to the GPU (only once copied to save computation time in the loop)
cudaMemcpy( gpuStruct->nImagesSet_gpu, nImagesSet, runmode->nsets * sizeof(int), cudaMemcpyHostToDevice );
cudaMemcpy( gpuStruct->images_gpu, images, runmode->nimagestot * sizeof(galaxy), cudaMemcpyHostToDevice );
cudaMemcpy(gpuStruct->runmode_gpu, runmode, sizeof(runmode_param), cudaMemcpyHostToDevice);
cudaMemcpy(gpuStruct->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;
}
#endif
//Run BayeSys3
BayeSys3(Common, Objects);
// Return potential parameters for lowest chi2
module_BayeSys3_returnBestPotential(potentials, Common);
// Free the memory
free( UserCommon->err );
free( UserCommon->avg );
free( UserCommon->potentials );
free( UserCommon->potentialsToOptimize );
free( Objects );
#if 0
// Free GPU memory
if (UserCommon->use_gpu == 1)
{
cudaFree(gpuStruct->nImagesSet_gpu);
cudaFree(gpuStruct->chi2_gpu);
cudaFree(gpuStruct->runmode_gpu);
cudaFree(gpuStruct->frame_gpu);
cudaFree(gpuStruct->lens_gpu);
cudaFree(gpuStruct->error_gpu);
cudaFree(gpuStruct->images_gpu);
cudaFree(gpuStruct);
}
#endif
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 module_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 module_BayeSys3_testchi(void (*chi2_function)(double*, int*, struct runmode_param *, const struct potential_param*, const struct grid_param *,const int*, galaxy*), char* FOLDER, runmode_param *runmode, grid_param *frame, potential_param* 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
module_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_param *)malloc((int) runmode->nhalos * sizeof(potential_param));
UserCommon->reference_potential = (potential_param *)malloc((int) runmode->nhalos * sizeof(potential_param));
// 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;
#if 0
// Add our gpuStruct to UserCommon
UserCommon->gpuStruct = gpuStruct;
#endif
// 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);
#if 0
// 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_param *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_param) ),"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;
}*/
#endif
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;
#if 0
// Copy the new data from BayeSys3 into the GPU memory
cudaMemcpy(lens_gpu, UserCommon->potentials, (UserCommon->nhalos) * sizeof(potential_param), 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 module_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);
#endif
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 module_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 update non independent parameters
*
* These include b0 (impact parameter) and ellipticity_potential for the moment.
*
*/
void update_non_independent_parameter(Potential_SOA *lens, int nhalos){
int N_type[100];
int Indice_type[100];
int ind;
for(int i=0;i < 100; ++i){
N_type[i] = 0;
Indice_type[i] = 0;
}
for (int i = 0; i < nhalos; ++i){
N_type[lens.type[i]] += 1;
}
for(int i=1;i < 100; ++i){
Indice_type[i] = N_type[i]+Indice_type[i-1];
//printf("%d %d \n ",N_type[i], Indice_type[i]);
}
for (int i = 0; i < nhalos; ++i){
if(Indice_type[lens.type[i]-1] <nhalos){
ind = Indice_type[lens.type[i]-1];
switch (lens->type[ind])
{
case(5): /*Elliptical Isothermal Sphere*/
//impact parameter b0
lens->b0[ind] = 4* pi_c2 * lens->vdisp[i] * lens->vdisp[i] ;
//ellipticity_potential
lens->ellipticity_potential[ind] = lens->ellipticity[i]/3 ;
break;
case(8): /* PIEMD */
//impact parameter b0
lens->b0[ind] = 6.*pi_c2 * lens->vdisp[i] * lens->vdisp[i];
//ellipticity_parameter
if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] != 0. ){
// emass is (a2-b2)/(a2+b2)
lens->ellipticity[ind] = 2.*lens->ellipticity_potential[ind] / (1. + lens->ellipticity_potential[ind] * lens->ellipticity_potential[ind]);
printf("1 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]);
}
else if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] == 0. ){
lens->ellipticity_potential[ind] = 0.00001;
printf("2 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]);
}
else{
// epot is (a-b)/(a+b)
lens->ellipticity_potential[ind] = (1. - sqrt(1 - lens->ellipticity[ind] * lens->ellipticity[ind])) / lens->ellipticity[ind];
printf("3 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]);
}
break;
default:
printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name[ind], lens->type[ind]);
break;
};
lens->b0[ind] = lens[i].b0;
lens->ellipticity_potential[ind] = lens[i].ellipticity_potential;
Indice_type[lens[i].type-1] += 1;
}
}
}
/** @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 module_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;
#if 0
BayeSys3_gpuStruct* gpuStruct = (BayeSys3_gpuStruct*)UserCommon->gpuStruct;
#endif
//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;
update_non_independent_parameter(&UserCommon->potentials);
for (int i = 0; i < UserCommon->nhalos; ++i){
printf("Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d, type_name = %s, name = %s,\n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n \t rscale = %f, exponent = %f, alpha = %f, einasto kappa critical = %f, z = %f\n", i,lenses.position_x[i], lenses.position_y[i], lenses.vdisp[i], lenses.type[i], lenses.type_name[i], lenses.name[i], lenses.ellipticity[i], lenses.ellipticity_potential[i], lenses.ellipticity_angle[i], lenses.rcore[i], lenses.rcut[i], lenses.rscale[i], lenses.exponent[i], lenses.alpha[i], lenses.einasto_kappacritic[i], lenses.z[i]);
}
#if 0
// Copy the new data from BayeSys3 into the GPU memory
cudasafe(cudaMemcpy(gpuStruct->lens_gpu, UserCommon->potentials, (UserCommon->nhalos) * sizeof(potential_param), cudaMemcpyHostToDevice), "Baysis3, lens-gpu");
cudasafe(cudaMemcpy(gpuStruct->chi2_gpu, &chi2, sizeof(double), cudaMemcpyHostToDevice),"Baysis3, chi2_gpu");
cudasafe(cudaMemcpy(gpuStruct->error_gpu, &error, sizeof(int), cudaMemcpyHostToDevice),"Baysis3, error_gpu");
printf("allocated\n");
#endif
// Run the chi2 function on the GPU
// Call the chi2_function specified in the module_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>>>(gpuStruct->chi2_gpu, gpuStruct->error_gpu, gpuStruct->runmode_gpu, gpuStruct->lens_gpu, gpuStruct->frame_gpu, gpuStruct->nImagesSet_gpu, gpuStruct->images_gpu);
UserCommon->chi2_function(chi2, error, UserCommon->runmode, UserCommon->potentials, UserCommon->grid, UserCommon->nImagesSet_gpu, UserCommon->constraint);
printf("finished running\n");
cudaError_t errAsync = cudaDeviceSynchronize();
printf("finished synchronising\n");
/*
//For debugging use remove it afterwards because it is relly slow
cudaError_t errSync = cudaGetLastError();
cudaError_t errAsync = cudaDeviceSynchronize();
if (errSync != cudaSuccess)
printf("Sync kernel error: %s\n", cudaGetErrorString(errSync));
if (errAsync != cudaSuccess)
printf("Async kernel error: %s\n", cudaGetErrorString(errAsync));*/
// 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(&error, gpuStruct->error_gpu, sizeof(int), cudaMemcpyDeviceToHost),"Baysis3, error_gpu");
printf("copied 1\n");
cudasafe(cudaMemcpy(&chi2, gpuStruct->chi2_gpu, sizeof(double), cudaMemcpyDeviceToHost),"Baysis3, chi2 result");
printf("copied 2\n");
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
module_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 module_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 = module_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 = module_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 module_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 *= module_BayeSys3_rescaleParam(haloID, parameterID, &cube[ipar]);
// The value of the cube[] element is rescaled in module_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 *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].position.x = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].position.y.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "position.y";
valid *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].position.y = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].weight.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "weight";
valid *= module_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 *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].vdisp = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity_angle.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity_angle";
valid *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].ellipticity_angle = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity";
valid *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].ellipticity = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].ellipticity_potential.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "ellipticity_potential";
valid *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].ellipticity_potential = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].rcore.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "rcore";
valid *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].rcore = 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 *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].rcut = cube[parameterCounter];
parameterCounter++;
};
if (UserCommon->potentialsToOptimize[haloCounter].rscale.block == 1) // 1=true, the user wants to optimize this parameter
{
stringVariable = "rscale";
valid *= module_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 *= module_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 *= module_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 *= module_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 *= module_BayeSys3_rescaleParam(haloCounter, stringVariable, &cube[parameterCounter], Common);
UserCommon->potentials[haloCounter].z = cube[parameterCounter];
parameterCounter++;
};
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;
}
// module_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 module_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 = module_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 = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].position.x.min, UserCommon->potentialsToOptimize[haloID].position.x.max);
}
if (parameterID == "position.y")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].position.y.min, UserCommon->potentialsToOptimize[haloID].position.y.max);
}
if (parameterID == "weight")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].weight.min, UserCommon->potentialsToOptimize[haloID].weight.max);
}
if (parameterID == "vdisp")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].vdisp.min, UserCommon->potentialsToOptimize[haloID].vdisp.max);
}
if (parameterID == "ellipticity_angle")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].ellipticity_angle.min, UserCommon->potentialsToOptimize[haloID].ellipticity_angle.max);
}
if (parameterID == "ellipticity")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].ellipticity.min, UserCommon->potentialsToOptimize[haloID].ellipticity.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "ellipticity_potential")
{
value = module_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 = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].rcore.min, UserCommon->potentialsToOptimize[haloID].rcore.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "rcut")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].rcut.min, UserCommon->potentialsToOptimize[haloID].rcut.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "rscale")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].rscale.min, UserCommon->potentialsToOptimize[haloID].rscale.max);
if (value < 0)
{
valid = 0;
}
}
if (parameterID == "exponent")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].exponent.min, UserCommon->potentialsToOptimize[haloID].exponent.max);
}
if (parameterID == "alpha")
{
value = module_BayeSys3_prior(Common, value, UserCommon->potentialsToOptimize[haloID].alpha.min, UserCommon->potentialsToOptimize[haloID].alpha.max);
}
if (parameterID == "einasto_kappacritic")
{
value = module_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 = module_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 module_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 module_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;
}
// module_BayeSys3_uniformPrior function
/** @brief Rescale Uniform [0:1] -> Uniform[x1:x2]
*/
double module_BayeSys3_uniformPrior(double val, double min, double max)
{
return min + val * (max - min);
}
/** @brief function to calculate the log-likelihood
* */
double module_BayeSys3_calculateLogLikelihood(double &chi2, int &error, CommonStr* Common)
{
double Lhood;
error = module_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];
}
module_BayeSys3_rescaleCube_1Atom(tmpCube, Common);
// Take a test chi2 from our current sample distribution "Nsample"
// Save the chi2 and Lhood
UserCommon->Nchi2++;
logLikelihood = module_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;
}
// module_BayeSys3_setReferencePotential function
//
/** @brief Set the lens parameters of the current lowest chi2 as a future reference
*/
int module_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[haloNumber].position.x = UserCommon->potentials[haloNumber].position.x;
UserCommon->reference_potential[haloNumber].position.y = UserCommon->potentials[haloNumber].position.y;
UserCommon->reference_potential[haloNumber].weight = UserCommon->potentials[haloNumber].weight;
UserCommon->reference_potential[haloNumber].vdisp = UserCommon->potentials[haloNumber].vdisp;
UserCommon->reference_potential[haloNumber].ellipticity_angle = UserCommon->potentials[haloNumber].ellipticity_angle;
UserCommon->reference_potential[haloNumber].ellipticity = UserCommon->potentials[haloNumber].ellipticity;
UserCommon->reference_potential[haloNumber].ellipticity_potential = UserCommon->potentials[haloNumber].ellipticity_potential;
UserCommon->reference_potential[haloNumber].rcore = UserCommon->potentials[haloNumber].rcore;
UserCommon->reference_potential[haloNumber].rcut = UserCommon->potentials[haloNumber].rcut;
UserCommon->reference_potential[haloNumber].rscale = UserCommon->potentials[haloNumber].rscale;
UserCommon->reference_potential[haloNumber].exponent = UserCommon->potentials[haloNumber].exponent;
UserCommon->reference_potential[haloNumber].alpha = UserCommon->potentials[haloNumber].alpha;
UserCommon->reference_potential[haloNumber].einasto_kappacritic = UserCommon->potentials[haloNumber].einasto_kappacritic;
UserCommon->reference_potential[haloNumber].z = UserCommon->potentials[haloNumber].z;
//never forget to calculate the parameter depended characteristic of the potential
//module_readParameters_calculatePotentialparameter(&UserCommon->reference_potential[haloNumber]);
}
return valid;
}
/** @brief function to update the statistics
*/
int module_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;
}
// module_BayeSys3_returnBestPotential function
//
/** @brief Set the lens parameters of the lowest chi2 as output potential parameters
*/
int module_BayeSys3_returnBestPotential(potential_param potentials[], CommonStr* Common)
{
UserCommonStr* UserCommon = (UserCommonStr*)Common->UserCommon;
printf( "finished: return best potential\n");
for (int haloNumber = 0; haloNumber < UserCommon->nhalos; haloNumber++) // Loop over all halos and their parameters
{
potentials[haloNumber].position.x = UserCommon->reference_potential[haloNumber].position.x;
potentials[haloNumber].position.y = UserCommon->reference_potential[haloNumber].position.y;
potentials[haloNumber].weight = UserCommon->reference_potential[haloNumber].weight;
potentials[haloNumber].vdisp = UserCommon->reference_potential[haloNumber].vdisp;
potentials[haloNumber].ellipticity_angle = UserCommon->reference_potential[haloNumber].ellipticity_angle;
potentials[haloNumber].ellipticity = UserCommon->reference_potential[haloNumber].ellipticity;
potentials[haloNumber].ellipticity_potential = UserCommon->reference_potential[haloNumber].ellipticity_potential;
potentials[haloNumber].rcore = UserCommon->reference_potential[haloNumber].rcore;
potentials[haloNumber].rcut = UserCommon->reference_potential[haloNumber].rcut;
potentials[haloNumber].rscale = UserCommon->reference_potential[haloNumber].rscale;
potentials[haloNumber].exponent = UserCommon->reference_potential[haloNumber].exponent;
potentials[haloNumber].alpha = UserCommon->reference_potential[haloNumber].alpha;
potentials[haloNumber].einasto_kappacritic = UserCommon->reference_potential[haloNumber].einasto_kappacritic;
potentials[haloNumber].z = UserCommon->reference_potential[haloNumber].z;
module_readParameters_calculatePotentialparameter(&UserCommon->potentials[haloNumber]);
}
return 0; // We are no longer in the BayeSys3, so we return 0, not a "valid" code
}
c4science · Help