Commit c0c71a9b authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

now able to sample multipliers alpha and able to sample without realizability reconstruction

parent dceec721
Pipeline #146098 failed with stage
in 1 second
......@@ -128,6 +128,8 @@ class Config
double _RealizableSetEpsilonU0; /*!< @brief Distance to 0 of the sampled moments to the boundary of the realizable set */
double _RealizableSetEpsilonU1; /*!< @brief norm(u_1)/u_0 !< _RealizableSetEpsilonU1 */
bool _normalizedSampling; /*!< @brief Flag for sampling of normalized moments, i.e. u_0 =1 */
bool _alphaSampling; /*!< @brief Flag for sampling alpha instead of u */
bool _realizabilityRecons; /*!< @brief Turns realizability reconstruction on/off for u sampling */
// --- Parsing Functionality and Initializing of Options ---
/*!
......@@ -326,6 +328,8 @@ class Config
double GetRealizableSetEpsilonU0() { return _RealizableSetEpsilonU0; }
double GetRealizableSetEpsilonU1() { return _RealizableSetEpsilonU1; }
bool inline GetNormalizedSampling() { return _normalizedSampling; }
bool inline GetAlphaSampling() { return _alphaSampling; }
bool inline GetRelizabilityReconsU() { return _realizabilityRecons; }
// ---- Setters for option structure
// This section is dangerous
......
......@@ -68,6 +68,7 @@ class DataGeneratorBase
// Main methods
virtual void SampleSolutionU() = 0; /*!< @brief Samples solution vectors u */
void SampleMultiplierAlpha(); /*!< @brief Sample Lagrange multipliers alpha */
void ComputeEntropyH_dual(); /*!< @brief Compute the entropy functional at (u,alpha) in dual formulation */
void ComputeEntropyH_primal(); /*!< @brief Compute the entropy functional at (u,alpha) in primal formulation */
void ComputeRealizableSolution(); /*!< @brief make u the realizable moment to alpha, since Newton has roundoff errors. */
......
......@@ -7,13 +7,16 @@
%
% ---- Datagenerator settings ----
DATA_GENERATOR_MODE = YES
TRAINING_SET_SIZE = 100
TRAINING_SET_SIZE = 100000
MAX_VALUE_FIRST_MOMENT = 1
SPATIAL_DIM = 1
REALIZABLE_SET_EPSILON_U1 = 0.003
REALIZABLE_SET_EPSILON_U0 = 0.003
MAX_MOMENT_SOLVER = 1
REALIZABLE_SET_EPSILON_U1 = 0.01
REALIZABLE_SET_EPSILON_U0 = 0.01
NORMALIZED_SAMPLING = YES
MAX_MOMENT_SOLVER = 2
ALPHA_SAMPLING = YES
REALIZABILITY_RECONS_U = NO
%
% ---- File specifications ----
%
......
......@@ -334,6 +334,12 @@ void Config::SetConfigOptions() {
/*! @brief Flag for sampling of normalized moments, i.e. u_0 =1 \n DESCRIPTION: Flag for sampling of normalized moments, i.e. u_0 =1 \n
* DEFAULT False \ingroup Config */
AddBoolOption( "NORMALIZED_SAMPLING", _normalizedSampling, false );
/*! @brief Flag for sampling the space of Legendre multipliers instead of the moments \n DESCRIPTION: Sample alpha instead of u \n DEFAULT False
* \ingroup Config */
AddBoolOption( "ALPHA_SAMPLING", _alphaSampling, false );
/*! @brief Flag for sampling the space of Legendre multipliers instead of the moments \n DESCRIPTION: Sample alpha instead of u \n DEFAULT False
* \ingroup Config */
AddBoolOption( "REALIZABILITY_RECONS_U", _realizabilityRecons, true );
}
void Config::SetConfigParsing( string case_filename ) {
......@@ -694,7 +700,7 @@ bool Config::TokenizeString( string& str, string& option_name, vector<string>& o
pos = str.find( "=" );
if( pos == string::npos ) {
string errmsg = "Error in Config::TokenizeString(): line in the configuration file with no \"=\" sign. ";
errmsg += "\nLook for: \n str.length() = " + std::to_string(str.length());
errmsg += "\nLook for: \n str.length() = " + std::to_string( str.length() );
spdlog::error( errmsg );
throw( -1 );
}
......
......@@ -81,24 +81,38 @@ DataGeneratorBase* DataGeneratorBase::Create( Config* settings ) {
void DataGeneratorBase::ComputeTrainingData() {
PrintLoadScreen();
// --- sample u ---
SampleSolutionU();
auto log = spdlog::get( "event" );
log->info( "| Moments sampled." );
log->info( "| Start solving the optimization problems. This may take some minutes." );
// ---- Check realizability ---
CheckRealizability();
if( _settings->GetAlphaSampling() ) {
// --- sample alpha ---
SampleMultiplierAlpha();
log->info( "| Multipliers sampled." );
log->info( "| Making moments realizable problems." );
// --- compute alphas ---
// --- Postprocessing
ComputeRealizableSolution();
}
else {
// --- sample u ---
SampleSolutionU();
log->info( "| Moments sampled." );
log->info( "| Start solving the optimization problems. This may take some minutes." );
_optimizer->SolveMultiCell( _alpha, _uSol, _moments );
// ---- Check realizability ---
CheckRealizability();
log->info( "| Making moments realizable problems." );
// --- compute alphas ---
// --- Postprocessing
ComputeRealizableSolution();
_optimizer->SolveMultiCell( _alpha, _uSol, _moments );
log->info( "| Making moments realizable problems." );
// --- Postprocessing
if( _settings->GetRelizabilityReconsU() ) {
ComputeRealizableSolution();
}
}
log->info( "| Compute entropies." );
......@@ -111,6 +125,58 @@ void DataGeneratorBase::ComputeTrainingData() {
PrintTrainingData();
}
void DataGeneratorBase::SampleMultiplierAlpha() {
double minAlphaValue = -50;
double maxAlphaValue = 50;
if( _settings->GetNormalizedSampling() ) {
// compute reduced version of alpha and m
if( _maxPolyDegree == 0 ) {
ErrorMessages::Error( "Normalized sampling not meaningful for M0 closure", CURRENT_FUNCTION );
}
VectorVector alphaRed = VectorVector( _setSize, Vector( _nTotalEntries - 1, 0.0 ) );
VectorVector momentsRed = VectorVector( _nq, Vector( _nTotalEntries - 1, 0.0 ) );
for( unsigned idx_nq = 0; idx_nq < _nq; idx_nq++ ) { // copy (reduced) moments
for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) {
momentsRed[idx_nq][idx_sys - 1] = _moments[idx_nq][idx_sys];
}
}
double dalpha = 0;
switch( _maxPolyDegree ) {
case 1:
// Sample alpha1 from [minAlphaValue, maxAlphaValue], then compute alpha_0 s.t. u_0 = 1
dalpha = ( maxAlphaValue - minAlphaValue ) / (double)_setSize;
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
alphaRed[idx_set][0] = minAlphaValue + idx_set * dalpha;
}
break;
default: ErrorMessages::Error( "Not yet implemented!", CURRENT_FUNCTION );
}
// Compute alpha_0 = log(<exp(alpha m )>) // for maxwell boltzmann! only
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
double integral = 0.0;
// Integrate (eta(eta'_*(alpha*m))
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
integral += _entropy->EntropyPrimeDual( dot( alphaRed[idx_set], momentsRed[idx_quad] ) ) * _weights[idx_quad];
}
_alpha[idx_set][0] = -log( integral ); // log trafo
// copy all other alphas to the member
for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) {
_alpha[idx_set][idx_sys] = alphaRed[idx_set][idx_sys - 1];
}
}
}
else {
// non normalized sampling
// TODO
ErrorMessages::Error( "Not yet implemented!", CURRENT_FUNCTION );
}
}
void DataGeneratorBase::ComputeEntropyH_dual() {
#pragma omp parallel for schedule( guided )
for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment