Diffusion Lenstool-HPC (develop_chris) Browse

wrapper_bayesis.cpp
8664617a3ed5develop_chris

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
			}

R1448/src/wrapper_bayesis.cpp8664617a3ed5develop_chris

wrapper_bayesis.cppView Options

R1448/src/wrapper_bayesis.cpp
8664617a3ed5develop_chris

wrapper_bayesis.cpp
View Options