Skip to content

Commit

Permalink
Random classes refactoring - alternative with minimal changes (#2364)
Browse files Browse the repository at this point in the history
* Random class refactoring

* Fixes

* Minimal modifications

* Reset modifications

* Address Mayeul's comments

* pseudo instead of psedo in Random*.h

* trigger CI: cdash server was restarted

---------

Co-authored-by: Anna Shlyaeva <[email protected]>
  • Loading branch information
benjaminmenetrier and shlyaeva authored Sep 19, 2023
1 parent f09288b commit a86234e
Show file tree
Hide file tree
Showing 14 changed files with 136 additions and 38 deletions.
8 changes: 3 additions & 5 deletions l95/src/lorenz95/LocalizationMatrixL95.cc
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,9 @@ namespace lorenz95 {
// -----------------------------------------------------------------------------
LocalizationMatrixL95::LocalizationMatrixL95(const Resolution & resol,
const oops::Variables & vars,
const eckit::Configuration & config,
const size_t & timeRank)
const eckit::Configuration & config)
: resol_(resol.npoints()),
rscale_(1.0/config.getDouble("length_scale")),
timeRank_(timeRank)
rscale_(1.0/config.getDouble("length_scale"))
{
// Gaussian structure function
unsigned int size = resol_/2+1;
Expand All @@ -51,7 +49,7 @@ LocalizationMatrixL95::LocalizationMatrixL95(const Resolution & resol,
}
// -----------------------------------------------------------------------------
void LocalizationMatrixL95::randomize(IncrementL95 & dx) const {
dx.random(1+timeRank_);
dx.random();
unsigned int size = resol_/2+1;
Eigen::FFT<double> fft;
std::vector<std::complex<double> > four(size);
Expand Down
3 changes: 1 addition & 2 deletions l95/src/lorenz95/LocalizationMatrixL95.h
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,7 @@ class LocalizationMatrixL95: public oops::interface::LocalizationBase<lorenz95::

LocalizationMatrixL95(const Resolution &,
const oops::Variables &,
const eckit::Configuration &,
const size_t &);
const eckit::Configuration &);
void randomize(IncrementL95 &) const override;
void multiply(IncrementL95 &) const override;

Expand Down
2 changes: 1 addition & 1 deletion qg/model/ErrorCovarianceQG.cc
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ namespace qg {
ErrorCovarianceQG::ErrorCovarianceQG(const GeometryQG & resol, const oops::Variables & vars,
const eckit::Configuration & conf,
const StateQG &, const StateQG &) {
qg_error_covariance_setup_f90(keyConfig_, conf, resol.toFortran(), 0);
qg_error_covariance_setup_f90(keyConfig_, conf, resol.toFortran());
oops::Log::trace() << "ErrorCovarianceQG created" << std::endl;
}
// -----------------------------------------------------------------------------
Expand Down
5 changes: 2 additions & 3 deletions qg/model/LocalizationMatrixQG.cc
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,8 @@ namespace qg {
// -----------------------------------------------------------------------------
LocalizationMatrixQG::LocalizationMatrixQG(const GeometryQG & resol,
const oops::Variables & vars,
const eckit::Configuration & config,
const size_t & timeRank) {
qg_error_covariance_setup_f90(keyLocal_, config, resol.toFortran(), timeRank);
const eckit::Configuration & config) {
qg_error_covariance_setup_f90(keyLocal_, config, resol.toFortran());
}
// -----------------------------------------------------------------------------
LocalizationMatrixQG::~LocalizationMatrixQG() {
Expand Down
3 changes: 1 addition & 2 deletions qg/model/LocalizationMatrixQG.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,7 @@ class LocalizationMatrixQG: public oops::interface::LocalizationBase<qg::QgTrait

LocalizationMatrixQG(const GeometryQG &,
const oops::Variables &,
const eckit::Configuration &,
const size_t &);
const eckit::Configuration &);
~LocalizationMatrixQG();

void randomize(IncrementQG &) const override;
Expand Down
2 changes: 1 addition & 1 deletion qg/model/QgFortran.h
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ extern "C" {
// Error covariance
// -----------------------------------------------------------------------------
void qg_error_covariance_setup_f90(F90error_covariance &, const eckit::Configuration &,
const F90geom &, const int &);
const F90geom &);
void qg_error_covariance_delete_f90(F90error_covariance &);
void qg_error_covariance_mult_f90(const F90error_covariance &, const F90flds &, const F90flds &);
void qg_error_covariance_inv_mult_f90(const F90error_covariance &, const F90flds &,
Expand Down
5 changes: 2 additions & 3 deletions qg/model/qg_error_covariance_interface.F90
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,14 @@ module qg_error_covariance_interface
contains
! ------------------------------------------------------------------------------
!> Setup error covariance matrix
subroutine qg_error_covariance_setup_c(c_key_self,c_conf,c_key_geom,c_time_rank) bind (c,name='qg_error_covariance_setup_f90')
subroutine qg_error_covariance_setup_c(c_key_self,c_conf,c_key_geom) bind (c,name='qg_error_covariance_setup_f90')

implicit none

! Passed variables
integer(c_int),intent(inout) :: c_key_self !< Error covariance configuration
type(c_ptr),value,intent(in) :: c_conf !< Configuration
integer(c_int),intent(in) :: c_key_geom !< Geometry
integer(c_int),intent(in) :: c_time_rank !< Time rank

! Local variables
type(fckit_configuration) :: f_conf
Expand All @@ -47,7 +46,7 @@ subroutine qg_error_covariance_setup_c(c_key_self,c_conf,c_key_geom,c_time_rank)
call qg_error_covariance_registry%get(c_key_self,self)

! Call Fortran
call qg_error_covariance_setup(self,f_conf,geom,c_time_rank)
call qg_error_covariance_setup(self,f_conf,geom)

end subroutine qg_error_covariance_setup_c
! ------------------------------------------------------------------------------
Expand Down
6 changes: 1 addition & 5 deletions qg/model/qg_error_covariance_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -58,15 +58,14 @@ module qg_error_covariance_mod
#include "oops/util/linkedList_c.f"
! ------------------------------------------------------------------------------
!> Setup error covariance matrix
subroutine qg_error_covariance_setup(self,f_conf,geom,time_rank)
subroutine qg_error_covariance_setup(self,f_conf,geom)

implicit none

! Passed variables
type(qg_error_covariance_config),intent(inout) :: self !< Error covariance configuration
type(fckit_configuration),intent(in) :: f_conf !< FCKIT configuration
type(qg_geom),intent(in) :: geom !< Geometry
integer,intent(in) :: time_rank !< Time rank

! Local variables
integer :: ix,iy,jy,ky,iz,jz,kz,info
Expand All @@ -91,9 +90,6 @@ subroutine qg_error_covariance_setup(self,f_conf,geom,time_rank)
self%seed = rseed
end if

! Modify the random seed based on the rank
self%seed = self%seed+time_rank

! Check nx
if (mod(geom%nx,2)/=0) then
write(record,*) 'qg_error_covariance_setup: number of zonal gridpoints nx=',geom%nx
Expand Down
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -397,6 +397,7 @@ oops/util/random_f.h
oops/util/random_mod.F90
oops/util/Random.h
oops/util/random.intfb.h
oops/util/RandomField.h
oops/util/Range.h
oops/util/ScalarOrMap.h
oops/util/Serializable.h
Expand Down
2 changes: 1 addition & 1 deletion src/oops/interface/LocalizationBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ class LocalizationMaker : public oops::LocalizationFactory<MODEL> {
std::unique_ptr<oops::LocalizationBase<MODEL>> make(const Geometry_ & geometry,
const oops::Variables & vars,
const eckit::Configuration & conf) override
{ return std::make_unique<T>(geometry.geometry(), vars, conf, geometry.timeComm().rank()); }
{ return std::make_unique<T>(geometry.geometry(), vars, conf); }
public:
explicit LocalizationMaker(const std::string & name) : oops::LocalizationFactory<MODEL>(name) {}
};
Expand Down
12 changes: 5 additions & 7 deletions src/oops/util/FieldSetHelpers.cc
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
#include "oops/util/FloatCompare.h"
#include "oops/util/Logger.h"
#include "oops/util/missingValues.h"
#include "oops/util/Random.h"
#include "oops/util/RandomField.h"

#define ERR(e) {ABORT(nc_strerror(e));}

Expand All @@ -33,8 +33,7 @@ namespace util {
atlas::FieldSet createRandomFieldSet(const eckit::mpi::Comm & comm,
const atlas::FunctionSpace & fspace,
const std::vector<size_t> & variableSizes,
const std::vector<std::string> & vars,
const size_t & timeRank) {
const std::vector<std::string> & vars) {
oops::Log::trace() << "createRandomFieldSet starting" << std::endl;

// Local ghost points
Expand Down Expand Up @@ -123,7 +122,7 @@ atlas::FieldSet createRandomFieldSet(const eckit::mpi::Comm & comm,
std::vector<double> rand_vec_glb(nglb);
if (comm.rank() == 0) {
// Generate global random vector
util::NormalDistribution<double> dist(nglb, 0.0, 1.0, 1+timeRank);
util::NormalDistributionField dist(nglb, 0.0, 1.0);
for (size_t i = 0; i < nglb; ++i) {
rand_vec_glb[i] = dist[i];
}
Expand Down Expand Up @@ -1288,13 +1287,12 @@ void writeFieldSet(const eckit::mpi::Comm & comm,

atlas::FieldSet createRandomFieldSet(const eckit::mpi::Comm & comm,
const atlas::FunctionSpace & fspace,
const oops::Variables & vars,
const size_t & timeRank) {
const oops::Variables & vars) {
std::vector<size_t> variableSizes;
for (const std::string & var : vars.variables()) {
variableSizes.push_back(vars.getLevels(var));
}
return createRandomFieldSet(comm, fspace, variableSizes, vars.variables(), timeRank);
return createRandomFieldSet(comm, fspace, variableSizes, vars.variables());
}

// -----------------------------------------------------------------------------
Expand Down
6 changes: 2 additions & 4 deletions src/oops/util/FieldSetHelpers.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,7 @@ namespace util {
atlas::FieldSet createRandomFieldSet(const eckit::mpi::Comm &,
const atlas::FunctionSpace &,
const std::vector<size_t> &,
const std::vector<std::string> &,
const size_t & timeRank = 0);
const std::vector<std::string> &);
/// Returns a fieldset with the same smooth field for all variables.
/// Useful for testing interpolation.
atlas::FieldSet createSmoothFieldSet(const eckit::mpi::Comm &,
Expand Down Expand Up @@ -63,8 +62,7 @@ void writeFieldSet(const eckit::mpi::Comm &,

atlas::FieldSet createRandomFieldSet(const eckit::mpi::Comm &,
const atlas::FunctionSpace &,
const oops::Variables &,
const size_t & timeRank = 0);
const oops::Variables &);
/// Returns a fieldset with the same smooth field for all variables.
/// Useful for testing interpolation.
atlas::FieldSet createSmoothFieldSet(const eckit::mpi::Comm &,
Expand Down
8 changes: 4 additions & 4 deletions src/oops/util/Random.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
namespace util {

// -----------------------------------------------------------------------------
/*! Base Class for generating compiler-indpependent random numbers
/*! Base Class for generating compiler-independent random numbers
*
* \details The *Random* class and its derived classes are intended to provide
* a compiler independent pseudo random number generator. So, if you use the
Expand Down Expand Up @@ -66,7 +66,7 @@ class Random : public util::Printable {
// -----------------------------------------------------------------------------
/*! Class for generating uniformly-distributed random numbers
*
* \details *util::UniformDistribution* creates a vector of psedo-random real
* \details *util::UniformDistribution* creates a vector of pseudo-random real
* numbers that are uniformly distributed across a specified interval of values.
*
* \param[in] N The size of the desired array (default 1)
Expand Down Expand Up @@ -116,7 +116,7 @@ class UniformDistribution : public Random<datatype> {
// -----------------------------------------------------------------------------
/*! Class for generating uniformly-distributed random integers
*
* \details *util::UniformIntDistribution* creates a vector of psedo-random integers
* \details *util::UniformIntDistribution* creates a vector of pseudo-random integers
* that are uniformly distributed across a specified interval of values.
*
* \param[in] N The size of the desired array (default 1)
Expand Down Expand Up @@ -164,7 +164,7 @@ class UniformIntDistribution : public Random<datatype> {
// ------------------------------------------------------------------------------
/*! Class for generating Gaussian-distributed random numbers
*
* \details *util::NormalDistribution* creates a vector of psedo-random numbers
* \details *util::NormalDistribution* creates a vector of pseudo-random numbers
* with a normal (Gaussian) distribution.
*
* \param[in] N The size of the desired array (default 1)
Expand Down
111 changes: 111 additions & 0 deletions src/oops/util/RandomField.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
/*
* (C) Copyright 2023 Meteorologisk Institutt
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
*/

#ifndef OOPS_UTIL_RANDOMFIELD_H_
#define OOPS_UTIL_RANDOMFIELD_H_

#include <algorithm>
#include <iostream>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
#include <boost/random.hpp>
#include "eckit/exception/Exceptions.h"
#include "oops/mpi/mpi.h"
#include "oops/util/abor1_cpp.h"
#include "oops/util/formats.h"
#include "oops/util/Logger.h"
#include "oops/util/Printable.h"

namespace util {

// -----------------------------------------------------------------------------
/*! Base Class for generating compiler-independent random numbers for fields
*
* \details The *RandomField* class and its derived classes are intended to provide
* a compiler independent pseudo random number generator. So, if you use the
* same random seed, you should get the same results, regardless of compiler
* or platform. Its sub-classes implement Gaussian (Normal) deviates.
* For usage tips, see each individual sub-class.
*
* \date Aug, 2023 (copy from Random.h, remove reset option)
*
* \sa util::NormalDistribution
*/

class RandomField : public util::Printable {
public:
const double & operator[](const size_t ii) const {return data_[ii];}
std::vector<double> data() {return data_;}

protected:
explicit RandomField(size_t N): N_(N) {data_.reserve(N);}
virtual ~RandomField() {}

size_t N_;
std::vector<double> data_;

private:
/*! This prints in a format that can be easily inserted into a yaml file for testing */
void print(std::ostream & os) const {
for (size_t jj=0; jj < N_; ++jj) {
os << " - " << util::full_precision(data_[jj]) << std::endl;
}
}
};

// ------------------------------------------------------------------------------
/*! Class for generating Gaussian-distributed random numbers
*
* \details *util::NormalDistributionField* creates a vector of pseudo-random numbers
* with a normal (Gaussian) distribution. Specific seed handling for fields (different
* on each MPI task).
*
* \param[in] N The size of the desired array
* \param[in] mean The mean of the distribution
* \param[in] sdev The standard deviation of the destribution
* \param[in] seed The seed to use for the random number generator (optional).
*
* \example Example usage:
* util::NormalDistributionField x(N,0.0,20.0)
* std::cout << x[i] << std::endl; // access one element
* std::cout << x << std::endl; // print full array
*/

class NormalDistributionField : public RandomField {
public:
NormalDistributionField(size_t N, double mean, double sdev):
RandomField(N), mean_(mean), sdev_(sdev) {
static boost::random::mt19937 generator(oops::mpi::world().rank()+1);
boost::random::normal_distribution<double> distribution(mean_, sdev_);
for (size_t jj=0; jj < this->N_; ++jj) this->data_.push_back(distribution(generator));
}
NormalDistributionField(size_t N, double mean, double sdev, unsigned int seed):
RandomField(N), mean_(mean), sdev_(sdev) {
static boost::random::mt19937 generator(oops::mpi::world().rank()+1);
std::stringstream ss;
ss << generator;
generator.seed(seed);
boost::random::normal_distribution<double> distribution(mean_, sdev_);
for (size_t jj=0; jj < this->N_; ++jj) this->data_.push_back(distribution(generator));
ss >> generator;
}

virtual ~NormalDistributionField() {}

private:
double mean_;
double sdev_;
};

// ------------------------------------------------------------------------------

} // namespace util

#endif // OOPS_UTIL_RANDOMFIELD_H_

0 comments on commit a86234e

Please sign in to comment.