mirror/music-panphasia
1
0
Fork 0
mirror of https://github.com/glatterf42/music-panphasia.git synced 2024-09-13 09:13:46 +02:00
Code Releases Activity
music-panphasia/random.hh
glatterf42 1e32e12349 First commit of original files.
2022-04-29 14:37:23 +02:00

495 lines
13 KiB
C++
Raw Permalink Blame History

/*
random.hh - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
*/
//... for testing purposes.............
//#define DEGRADE_RAND1
//#define DEGRADE_RAND2
//.....................................
#ifndef __RANDOM_HH
#define __RANDOM_HH
#define DEF_RAN_CUBE_SIZE 32
#include <fstream>
#include <algorithm>
#include <map>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "general.hh"
#include "mesh.hh"
#include "mg_operators.hh"
#include "constraints.hh"
#include "convolution_kernel.hh"
class RNG_plugin{
protected:
config_file *pcf_; //!< pointer to config_file from which to read parameters
const refinement_hierarchy *prefh_; //!< pointer to refinement hierarchy structure containing the grid sizes
public:
explicit RNG_plugin( config_file& cf )//, const refinement_hierarchy& refh )
: pcf_( &cf )//, prefh_( & refh )
{ }
virtual ~RNG_plugin() { }
virtual bool is_multiscale() const = 0;
virtual void fill_grid( int level, DensityGrid<real_t>& R ) = 0;
virtual void initialize_for_grid_structure( const refinement_hierarchy& refh ) = 0;
};
struct RNG_plugin_creator
{
virtual RNG_plugin * create( config_file& cf ) const = 0;
virtual ~RNG_plugin_creator() { }
};
std::map< std::string, RNG_plugin_creator *>&
get_RNG_plugin_map();
void print_RNG_plugins( void );
template< class Derived >
struct RNG_plugin_creator_concrete : public RNG_plugin_creator
{
//! register the plugin by its name
RNG_plugin_creator_concrete( const std::string& plugin_name )
{
get_RNG_plugin_map()[ plugin_name ] = this;
}
//! create an instance of the plugin
RNG_plugin* create( config_file& cf ) const //, const refinement_hierarchy& refh ) const
{
return new Derived( cf ); //, refh );
}
};
typedef RNG_plugin RNG_instance;
RNG_plugin *select_RNG_plugin( config_file& cf );//, const refinement_hierarchy& refh );
/*!
* @brief encapsulates all things for multi-scale white noise generation
*/
template< typename T >
class random_number_generator
{
protected:
config_file * pcf_;
//const refinement_hierarchy * prefh_;
RNG_plugin * generator_;
int levelmin_, levelmax_;
public:
//! constructor
random_number_generator( config_file& cf, transfer_function *ptf = NULL )
: pcf_( &cf )//, prefh_( &refh )
{
levelmin_ = pcf_->getValue<int>("setup","levelmin");
levelmax_ = pcf_->getValue<int>("setup","levelmax");
generator_ = select_RNG_plugin( cf );
}
//! destructor
~random_number_generator()
{
}
//! initialize_for_grid_structure
void initialize_for_grid_structure( const refinement_hierarchy& refh )
{
generator_->initialize_for_grid_structure( refh );
}
//! load random numbers to a new array
template< typename array >
void load( array& A, int ilevel )
{
generator_->fill_grid( ilevel, A );
}
};
//typedef random_numbers<real_t> rand_nums;
//typedef random_number_generator< rand_nums,real_t> rand_gen;
typedef random_number_generator<real_t> noise_generator;
#if 0
constraint_set constraints;
int levelmin_,
levelmax_,
levelmin_seed_;
std::vector<long> rngseeds_;
std::vector<std::string> rngfnames_;
bool disk_cached_;
bool restart_;
std::vector< std::vector<T>* > mem_cache_;
unsigned ran_cube_size_;
protected:
//! checks if the specified string is numeric
bool is_number(const std::string& s);
//! parses the random number parameters in the conf file
void parse_rand_parameters( void );
//! correct coarse grid averages for the change in small scale when using Fourier interpolation
void correct_avg( int icoarse, int ifine );
//! the main driver routine for multi-scale white noise generation
void compute_random_numbers( void );
//! store the white noise fields in memory or on disk
void store_rnd( int ilevel, rng* prng );
public:
//! constructor
random_number_generator( config_file& cf, refinement_hierarchy& refh, transfer_function *ptf = NULL );
//! destructor
~random_number_generator();
//! load random numbers to a new array
template< typename array >
void load( array& A, int ilevel )
{
if( restart_ )
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.",ilevel,ilevel);
if( disk_cached_ )
{
char fname[128];
sprintf(fname,"wnoise_%04d.bin",ilevel);
LOGUSER("Loading white noise from file \'%s\'...",fname);
std::ifstream ifs( fname, std::ios::binary );
if( !ifs.good() )
{
LOGERR("White noise file \'%s\'was not found.",fname);
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
}
int nx,ny,nz;
ifs.read( reinterpret_cast<char*> (&nx), sizeof(int) );
ifs.read( reinterpret_cast<char*> (&ny), sizeof(int) );
ifs.read( reinterpret_cast<char*> (&nz), sizeof(int) );
if( nx!=(int)A.size(0) || ny!=(int)A.size(1) || nz!=(int)A.size(2) )
{
if( nx==(int)A.size(0)*2 && ny==(int)A.size(1)*2 && nz==(int)A.size(2)*2 )
{
std::cerr << "CHECKPOINT" << std::endl;
int ox = nx/4, oy = ny/4, oz = nz/4;
std::vector<T> slice( ny*nz, 0.0 );
for( int i=0; i<nx; ++i )
{
ifs.read( reinterpret_cast<char*> ( &slice[0] ), ny*nz*sizeof(T) );
if( i<ox ) continue;
if( i>=3*ox ) break;
#pragma omp parallel for
for( int j=oy; j<3*oy; ++j )
for( int k=oz; k<3*oz; ++k )
A(i-ox,j-oy,k-oz) = slice[j*nz+k];
}
ifs.close();
}
else
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].",
nx,ny,nz,A.size(0),A.size(1),A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad.");
}
}else{
for( int i=0; i<nx; ++i )
{
std::vector<T> slice( ny*nz, 0.0 );
ifs.read( reinterpret_cast<char*> ( &slice[0] ), ny*nz*sizeof(T) );
#pragma omp parallel for
for( int j=0; j<ny; ++j )
for( int k=0; k<nz; ++k )
A(i,j,k) = slice[j*nz+k];
}
ifs.close();
}
}
else
{
LOGUSER("Copying white noise from memory cache...");
if( mem_cache_[ilevel-levelmin_] == NULL )
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n",ilevel);
int nx( A.size(0) ), ny( A.size(1) ), nz( A.size(2) );
if ( (size_t)nx*(size_t)ny*(size_t)nz != mem_cache_[ilevel-levelmin_]->size() )
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].",nx,ny,nz,A.size(0),A.size(1),A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
}
#pragma omp parallel for
for( int i=0; i<nx; ++i )
for( int j=0; j<ny; ++j )
for( int k=0; k<nz; ++k )
A(i,j,k) = (*mem_cache_[ilevel-levelmin_])[((size_t)i*ny+(size_t)j)*nz+(size_t)k];
std::vector<T>().swap( *mem_cache_[ilevel-levelmin_] );
delete mem_cache_[ilevel-levelmin_];
mem_cache_[ilevel-levelmin_] = NULL;
}
}
};
/*!
* @brief encapsulates all things random number generator related
*/
template< typename T >
class random_numbers
{
public:
unsigned
res_, //!< resolution of the full mesh
cubesize_, //!< size of one independent random number cube
ncubes_; //!< number of random number cubes to cover the full mesh
long baseseed_; //!< base seed from which cube seeds are computed
protected:
//! vector of 3D meshes (the random number cubes) with random numbers
std::vector< Meshvar<T>* > rnums_;
//! map of 3D indices to cube index
std::map<size_t,size_t> cubemap_;
typedef std::map<size_t,size_t>::iterator cubemap_iterator;
protected:
//! register a cube with the hash map
void register_cube( int i, int j, int k);
//! fills a subcube with random numbers
double fill_cube( int i, int j, int k);
//! subtract a constant from an entire cube
void subtract_from_cube( int i, int j, int k, double val );
//! copy random numbers from a cube to a full grid array
template< class C >
void copy_cube( int i, int j, int k, C& dat )
{
int offi, offj, offk;
offi = i*cubesize_;
offj = j*cubesize_;
offk = k*cubesize_;
i = (i+ncubes_)%ncubes_;
j = (j+ncubes_)%ncubes_;
k = (k+ncubes_)%ncubes_;
size_t icube = (i*ncubes_+j)*ncubes_+k;
cubemap_iterator it = cubemap_.find( icube );
if( it == cubemap_.end() )
{
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d",i,j,k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
for( int ii=0; ii<(int)cubesize_; ++ii )
for( int jj=0; jj<(int)cubesize_; ++jj )
for( int kk=0; kk<(int)cubesize_; ++kk )
dat(offi+ii,offj+jj,offk+kk) = (*rnums_[cubeidx])(ii,jj,kk);
}
//! free the memory associated with a subcube
void free_cube( int i, int j, int k );
//! initialize member variables and allocate memory
void initialize( void );
//! fill a cubic subvolume of the full grid with random numbers
double fill_subvolume( int *i0, int *n );
//! fill an entire grid with random numbers
double fill_all( void );
//! fill an external array instead of the internal field
template< class C >
double fill_all( C& dat )
{
double sum = 0.0;
for( int i=0; i<(int)ncubes_; ++i )
for( int j=0; j<(int)ncubes_; ++j )
for( int k=0; k<(int)ncubes_; ++k )
{
int ii(i),jj(j),kk(k);
register_cube(ii,jj,kk);
}
#pragma omp parallel for reduction(+:sum)
for( int i=0; i<(int)ncubes_; ++i )
for( int j=0; j<(int)ncubes_; ++j )
for( int k=0; k<(int)ncubes_; ++k )
{
int ii(i),jj(j),kk(k);
ii = (ii+ncubes_)%ncubes_;
jj = (jj+ncubes_)%ncubes_;
kk = (kk+ncubes_)%ncubes_;
sum+=fill_cube(ii, jj, kk);
copy_cube(ii,jj,kk,dat);
free_cube(ii, jj, kk);
}
return sum/(ncubes_*ncubes_*ncubes_);
}
//! write the number of allocated random number cubes to stdout
void print_allocated( void );
public:
//! constructor
random_numbers( unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx );
//! constructor for constrained fine field
random_numbers( random_numbers<T>& rc, unsigned cubesize, long baseseed,
bool kspace=false, bool isolated=false, int *x0_=NULL, int *lx_=NULL, bool zeromean=true );
//! constructor
random_numbers( unsigned res, unsigned cubesize, long baseseed, bool zeromean=true );
//! constructor to read white noise from file
random_numbers( unsigned res, std::string randfname, bool rndsign );
//! copy constructor for averaged field (not copying) hence explicit!
explicit random_numbers( /*const*/ random_numbers <T>& rc, bool kdegrade = true );
//! destructor
~random_numbers()
{
for( unsigned i=0; i<rnums_.size(); ++i )
if( rnums_[i] != NULL )
delete rnums_[i];
rnums_.clear();
}
//! access a random number, this allocates a cube and fills it with consistent random numbers
inline T& operator()( int i, int j, int k, bool fillrand=true )
{
int ic, jc, kc, is, js, ks;
if( ncubes_ == 0 )
throw std::runtime_error("random_numbers: internal error, not properly initialized");
//... determine cube
ic = (int)((double)i/cubesize_ + ncubes_) % ncubes_;
jc = (int)((double)j/cubesize_ + ncubes_) % ncubes_;
kc = (int)((double)k/cubesize_ + ncubes_) % ncubes_;
size_t icube = ((size_t)ic*ncubes_+(size_t)jc)*ncubes_+(size_t)kc;
cubemap_iterator it = cubemap_.find( icube );
if( it == cubemap_.end() )
{
LOGERR("Attempting to copy data from non-existing RND cube %d,%d,%d @ %d,%d,%d",ic,jc,kc,i,j,k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
if( rnums_[ cubeidx ] == NULL )
{
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d",ic,jc,kc);
throw std::runtime_error("attempting to access data from non-allocated RND cube");
}
//... determine cell in cube
is = (i - ic * cubesize_ + cubesize_) % cubesize_;
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
return (*rnums_[ cubeidx ])(is,js,ks);
}
//! free all cubes
void free_all_mem( void )
{
for( unsigned i=0; i<rnums_.size(); ++i )
if( rnums_[i] != NULL )
{
delete rnums_[i];
rnums_[i] = NULL;
}
}
};
#endif
#endif //__RANDOM_HH
View git blame Copy permalink