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

Former-commit-id: c0c71a9b

code/include/common/config.h

@@ -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

code/include/toolboxes/datageneratorbase.h

@@ -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. */

code/input/DataGenerator1D.cfg

@@ -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 ----
 %

code/src/common/config.cpp

@@ -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 );
     }

code/src/toolboxes/datageneratorbase.cpp

@@ -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++ ) {