Homec4science
Diffusion Lenstool-HPC (develop_chris) Browse

R1448/Benchmarks/BayesBenchmark/main.cpp
8664617a3ed5develop_chris

Last Change
Raw File

main.cpp
View Options

/**
Lenstool-HPC: HPC based massmodeling software and Lens-map generation
Copyright (C) 2017 Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch), Gilles Fourestey (gilles.fourestey@epfl.ch)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
@brief: Main Lenstool-HPC executable
*/
#include <iostream>
#include <iomanip>
#include <string.h>
#include <math.h>
#include <sys/time.h>
#include <fstream>
#include <sys/stat.h>
#include <unistd.h>
//
//#include <mm_malloc.h>
#include <omp.h>
//
//#include <cuda_runtime.h>
#include <structure_hpc.hpp>
//#include <cuda.h>
#include "timer.h"
#include "gradient.hpp"
#include "gradient2.hpp"
#include "chi_CPU.hpp"
#include "module_cosmodistances.hpp"
#include "module_readParameters.hpp"
#include "grid_gradient2_CPU.hpp"
#include "grid_amplif_CPU.hpp"
#include "module_writeFits.hpp"
#include "write_output.hpp"
#include "delense_CPU.hpp"
#include "allocation.hpp"
#include "wrapper_bayesis.hpp"
#ifdef __WITH_MPI
#include<mpi.h>
#endif
#ifdef _OPENMP
#include<omp.h>
#endif
#ifdef __WITH_GPU
#include "image_prediction_GPU.cuh"
#include "grid_gradient_GPU.cuh"
#include "grid_map_ampli_GPU.cuh"
#include "grid_map_pot_GPU.cuh"
#include "grid_map_shear_GPU.cuh"
#include "grid_map_mass_GPU.cuh"
#include "grid_map_dpl_GPU.cuh"
#include "grid_gradient2_GPU.cuh"
#include "allocation_GPU.cuh"
#include "chi_barycenter_GPU.cuh"
//#include "gradient_GPU.cuh"
#endif
void
gradient_grid_GPU_sorted(type_t *grid_grad_x, type_t *grid_grad_y, const struct grid_param *frame, const struct Potential_SOA *lens, int Nlens, int nbgridcells);
//
//
int module_readCheckInput_readInput(int argc, char *argv[], std::string *outdir)
{
/// check if there is a correct number of arguments, and store the name of the input file in infile
char* infile;
struct stat file_stat;
// If we do not have 3 arguments, stop
if ( argc != 3 )
{
fprintf(stderr, "\nUnexpected number of arguments\n");
fprintf(stderr, "\nUSAGE:\n");
fprintf(stderr, "lenstool input_file output_directorypath [-n]\n\n");
exit(-1);
}
else if ( argc == 3 )
infile=argv[1];
std::ifstream ifile(infile,std::ifstream::in); // Open the file
int ts = (int) time (NULL);
char buffer[10];
std::stringstream ss;
ss << ts;
std::string trimstamp = ss.str();
//
//std::string outdir = argv[2];
*outdir = argv[2];
//*outdir += "-";
//*outdir += trimstamp;
std::cout << *outdir << std::endl;
// check whether the output directory already exists
if (stat(outdir->c_str(), &file_stat) < 0){
mkdir(outdir->c_str(), S_IRUSR | S_IWUSR | S_IXUSR | S_IRGRP | S_IWGRP | S_IXGRP | S_IROTH );
}
else {
printf("Error : Directory %s already exists. Specify a non existing directory.\n",argv[2]);
exit(-1);
}
// check whether the input file exists. If it could not be opened (ifile = 0), it does not exist
if(ifile){
ifile.close();
}
else{
printf("The file %s does not exist, please specify a valid file name\n",infile);
exit(-1);
}
return 0;
}
//
//
//
int main(int argc, char *argv[])
{
int world_size = 1;
int world_rank = 0;
#ifdef __WITH_MPI
MPI_Init(NULL, NULL);
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
char processor_name[MPI_MAX_PROCESSOR_NAME];
int name_len;
MPI_Get_processor_name(processor_name, &name_len);
MPI_Barrier(MPI_COMM_WORLD);
// Print off a hello world message
#endif
int numthreads = 1;
#ifdef _OPENMP
#warning "using openmp"
#pragma omp parallel
numthreads = omp_get_num_threads();
#endif
//
int verbose = (world_rank == 0);
//
if (verbose) printf("\nLenstool-HPC\n\n"); fflush(stdout);
//
#ifdef __WITH_MPI
printf("Hello world from processor %s, rank %d out of %d processors and %d threads per rank\n", processor_name, world_rank, world_size, numthreads); fflush(stdout);
MPI_Barrier(MPI_COMM_WORLD);
#endif
//
double wallclock = myseconds();
if (world_rank == 0) printf("Reading parameter file at time %f s... ", myseconds() - wallclock);
// Setting Up the problem
//===========================================================================================================
// This module function reads the terminal input when calling LENSTOOL and checks that it is correct
// Otherwise it exits LENSTOOL
std::string path;
if (world_rank == 0) module_readCheckInput_readInput(argc, argv, &path);
// This module function reads the cosmology parameters from the parameter file
// Input: struct cosmologicalparameters cosmology, parameter file
// Output: Initialized cosmology struct
cosmo_param cosmology; // Cosmology struct to store the cosmology data from the file
std::string inputFile = argv[1]; // Input file
module_readParameters_readCosmology(inputFile, cosmology);
// This module function reads the runmode paragraph and the number of sources, arclets, etc. in the parameter file.
// The runmode_param stores the information of what exactly the user wants to do with lenstool.
struct runmode_param runmode;
module_readParameters_readRunmode(inputFile, &runmode);
runmode.debug = 1;
if (world_rank == 0)
{
module_readParameters_debug_cosmology(runmode.debug, cosmology);
module_readParameters_debug_runmode(runmode.debug, runmode);
}
//=== Declaring variables
struct grid_param frame;
struct galaxy images [runmode.nimagestot];
struct galaxy sources[runmode.nsets];
//struct Potential lenses[runmode.nhalos+runmode.npotfile-1];
struct Potential_SOA* lenses_SOA;
struct cline_param cline;
struct potfile_param potfile[runmode.Nb_potfile];
//struct Potential potfilepotentials[runmode.npotfile];
struct potentialoptimization host_potentialoptimization[runmode.nhalos];
int nImagesSet[runmode.nsets]; // Contains the number of images in each set of imagesnano
// This module function reads in the potential form and its parameters (e.g. NFW)
// Input: input file
// Output: Potentials and its parameters
#ifdef __WITH_GPU
PotentialSOAAllocation_GPU(&lenses_SOA, runmode.n_tot_halos);
#else
PotentialSOAAllocation(&lenses_SOA, runmode.n_tot_halos);
#endif
module_readParameters_PotentialSOA_noalloc(inputFile, lenses_SOA, runmode.nhalos, runmode.n_tot_halos, cosmology);
//
module_readParameters_debug_potential_SOA(1, lenses_SOA, runmode.nhalos);
module_readParameters_limit(inputFile, host_potentialoptimization, runmode.nhalos );
// This module function reads in the potfiles parameters
// Input: input file
// Output: Potentials from potfiles and its parameters
if (runmode.potfile == 1 )
{
module_readParameters_readpotfiles_param (inputFile, potfile, cosmology);
for(int ii = 0; ii < runmode.Nb_potfile; ii++){
module_readParameters_debug_potfileparam (1, &potfile[ii]);
}
module_readParameters_readpotfiles_SOA (&runmode , &cosmology, potfile, lenses_SOA);
module_readParameters_debug_potential_SOA(1, lenses_SOA, runmode.n_tot_halos);
}
module_readParameters_lens_dslds_calculation(&runmode,&cosmology,lenses_SOA);
// This module function reads in the grid form and its parameters
// Input: input file
// Output: grid and its parameters
module_readParameters_Grid(inputFile, &frame);
if (runmode.image == 1 or runmode.inverse == 1 or runmode.time > 0)
{
printf("Initialization of the images..."); fflush(stdout);
// This module function reads in the strong lensing images
module_readParameters_readImages(&runmode, images, nImagesSet);
//runmode.nsets = runmode.nimagestot;
for(int i = 0; i < runmode.nimagestot; ++i)
{
images[i].dls = module_cosmodistances_objectObject(lenses_SOA->z[0], images[i].redshift, cosmology);
images[i].dos = module_cosmodistances_observerObject(images[i].redshift, cosmology);
images[i].dr = module_cosmodistances_lensSourceToObserverSource(lenses_SOA->z[0], images[i].redshift, cosmology);
}
printf("end.\n"); fflush(stdout);
if (verbose) module_readParameters_debug_image(runmode.debug, images, nImagesSet, runmode.nsets);
}
printf("done.\n");fflush(stdout);
if (runmode.source == 1)
{
printf("Initialization of the sources..."); fflush(stdout);
//Initialisation to default values.(Setting sources to z = 1.5 default value)
for(int i = 0; i < runmode.nsets; ++i)
{
sources[i].shape.a = sources[i].shape.b = sources[i].shape.theta = (type_t) 0.;
sources[i].redshift = 1.5;
}
// This module function reads in the strong lensing sources
module_readParameters_readSources(&runmode, sources);
//Calculating cosmoratios
for(int i = 0; i < runmode.nsets; ++i)
{
sources[i].dls = module_cosmodistances_objectObject(lenses_SOA->z[0], sources[i].redshift, cosmology);
sources[i].dos = module_cosmodistances_observerObject(sources[i].redshift, cosmology);
sources[i].dr = module_cosmodistances_lensSourceToObserverSource(lenses_SOA->z[0], sources[i].redshift, cosmology);
}
module_readParameters_debug_source(runmode.debug, sources, runmode.nsets);
printf("end.\n"); fflush(stdout);
}
//
#ifdef __WITH_MPI
MPI_Barrier(MPI_COMM_WORLD);
#endif
//
//
//
std::cout << "--------------------------" << std::endl << std::endl; fflush(stdout);
#ifdef __WITH_GPU
runBayeSys3 (&mychi_bruteforce_SOA_CPU_grid_gradient_barycentersource,
argv[2], &runmode, &frame, lenses_SOA, host_potentialoptimization, images,
nImagesSet, 0.1, 1 );
std::cout << "Finished Bayesian Statistics " << std::endl;
module_readParameters_debug_potential_SOA(1, lenses_SOA, runmode.nhalos);
#endif
}
c4science · Help