bayes_app.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Aug 9, 20:41

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->nsetssizeof(int)),"Alloc nImagesSet_gpu: " );
	cudasafe(cudaMalloc( (void*)&(gpuStruct->images_gpu), runmode->nimagestotsizeof(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->nsetssizeof(int)),"Alloc nImagesSet_gpu: " );
	cudasafe(cudaMalloc( (void*)&(images_gpu), runmode->nimagestotsizeof(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
	}

bayes_app.cppNo OneTemporaryActions

File Metadata

bayes_app.cppView Options

Event Timeline

bayes_app.cpp
No OneTemporary
Actions

bayes_app.cpp
View Options