wrapper_bayesis.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Dec 1, 07:26

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
	}

wrapper_bayesis.cppNo OneTemporaryActions

File Metadata

wrapper_bayesis.cppView Options

Event Timeline

wrapper_bayesis.cpp
No OneTemporary
Actions

wrapper_bayesis.cpp
View Options