implemented dimension reduction strategies

Former-commit-id: 300652ec

code/include/toolboxes/datageneratorbase.h

@@ -52,7 +52,7 @@ class DataGeneratorBase
         _gridSize; /*!< @brief Size of the grid discretizing moment U0 for higher order sampling (has different uses for different samplers)*/
 
     unsigned short _maxPolyDegree; /*!< @brief Max Order of Spherical Harmonics */
-    unsigned _nTotalEntries;    /*!< @brief Total number of equations in the system */
+    unsigned _nTotalEntries;       /*!< @brief Total number of equations in the system */
 
     QuadratureBase* _quadrature;    /*!< @brief quadrature to create members below */
     unsigned _nq;                   /*!< @brief number of quadrature points */
@@ -80,5 +80,6 @@ class DataGeneratorBase
     virtual void ComputeMoments()     = 0; /*!< @brief Pre-Compute Moments at all quadrature points. */
     virtual void CheckRealizability() = 0; /*!< @brief Debugging helper. Will be removed */
     virtual void ComputeSetSize()     = 0; /*!< @brief Computes the size of the training set, depending on the chosen settings.*/
+    void AdaptBasisSize();                 /*!< @brief In case of normal sampling, deletes zero order basis for dimension reduction. */
 };
 #endif    // DATAGENERATOR_H

code/input/DataGenerator1D.cfg

@@ -7,20 +7,20 @@
 %
 % ---- Datagenerator settings ----
 DATA_GENERATOR_MODE = YES
-TRAINING_SET_SIZE = 100000
+TRAINING_SET_SIZE = 100
 MAX_VALUE_FIRST_MOMENT = 1
 SPATIAL_DIM = 1
 REALIZABLE_SET_EPSILON_U1 = 0.003
 REALIZABLE_SET_EPSILON_U0 = 0.003
 NORMALIZED_SAMPLING = YES
-MAX_MOMENT_SOLVER = 1
+MAX_MOMENT_SOLVER = 2
 %
 % ---- File specifications ----
 %
 % Output directory
 OUTPUT_DIR = ../result
 % Output file
-OUTPUT_FILE = trainM01D
+OUTPUT_FILE = trainM11D
 % Log directory
 LOG_DIR = ../result/logs
 LOG_FILE = dataGen1D

code/src/optimizers/newtonoptimizer.cpp

@@ -9,6 +9,7 @@
 #include "entropies/entropybase.h"
 #include "quadratures/quadraturebase.h"
 #include "toolboxes/errormessages.h"
+#include "toolboxes/textprocessingtoolbox.h"
 
 #include <omp.h>
 
@@ -82,6 +83,9 @@ void NewtonOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& sol, co
 
 #pragma omp parallel for schedule( guided )
     for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {
+        // std::cout << "Sol Vector"
+        //          << "|" << sol[idx_cell] << "\n";
+
         Solve( alpha[idx_cell], sol[idx_cell], moments, idx_cell );
         // if( idx_cell % 10000 == 0 ) {
         //    printf( "%d\n", idx_cell );

code/src/toolboxes/datagenerator1D.cpp

@@ -16,6 +16,8 @@
 DataGenerator1D::DataGenerator1D( Config* settings ) : DataGeneratorBase( settings ) {
     ComputeMoments();
 
+    AdaptBasisSize();
+
     // Initialize Training Data
     ComputeSetSize();
 
@@ -140,8 +142,8 @@ void DataGenerator1D::SampleSolutionU() {
             double dN  = 2.0 / (double)_gridSize;
             double N1  = -1.0 + _settings->GetRealizableSetEpsilonU0();
             while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
-                _uSol[c][0] = 1;     // u0 (normalized i.e. N0) by Monreals notation
-                _uSol[c][1] = N1;    // u1 (normalized i.e. N1) by Monreals notation
+
+                _uSol[c][0] = N1;    // u1 (normalized i.e. N1) by Monreals notation
                 N1 += dN;
                 c++;
             }
@@ -155,9 +157,8 @@ void DataGenerator1D::SampleSolutionU() {
             while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
                 N2 = N1 * N1 + _settings->GetRealizableSetEpsilonU0();
                 while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
-                    _uSol[c][0] = 1;     // u0 (normalized i.e. N0) by Monreals notation
-                    _uSol[c][1] = N1;    // u1 (normalized i.e. N1) by Monreals notation
-                    _uSol[c][2] = N2;    // u2 (normalized i.e. N2) by Monreals notation
+                    _uSol[c][0] = N1;    // u1 (normalized i.e. N1) by Monreals notation
+                    _uSol[c][1] = N2;    // u2 (normalized i.e. N2) by Monreals notation
                     N2 += dN;
                     c++;
                 }
@@ -174,10 +175,9 @@ void DataGenerator1D::SampleSolutionU() {
             while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
                 N3 = -N2 + ( N1 + N2 ) * ( N1 + N2 ) / ( 1 + N1 ) + _settings->GetRealizableSetEpsilonU1();
                 while( N3 < N2 - ( N1 - N2 ) * ( N1 - N2 ) / ( 1 - N1 ) - _settings->GetRealizableSetEpsilonU1() ) {
-                    _uSol[c][0] = 1;     // u0  by Monreals notation
-                    _uSol[c][1] = N1;    // u1  by Monreals notation
-                    _uSol[c][2] = N2;    // u2  by Monreals notation
-                    _uSol[c][3] = N3;    // u3  by Monreals notation
+                    _uSol[c][0] = N1;    // u1  by Monreals notation
+                    _uSol[c][1] = N2;    // u2  by Monreals notation
+                    _uSol[c][2] = N3;    // u3  by Monreals notation
                     c++;
                     N3 += dN;
                 }
@@ -294,7 +294,7 @@ void DataGenerator1D::ComputeSetSize() {
         }
         else if( _maxPolyDegree == 2 ) {
             // Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1
-            unsigned c = 1;
+            unsigned c = 0;
             double N1  = -1.0 + _settings->GetRealizableSetEpsilonU0();
             double N2;
             double dN = 2.0 / (double)_gridSize;
@@ -310,7 +310,7 @@ void DataGenerator1D::ComputeSetSize() {
         }
         else if( _maxPolyDegree == 3 ) {
             // Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1, N1=0
-            unsigned c = 1;
+            unsigned c = 0;
             double N1  = 0 + _settings->GetRealizableSetEpsilonU0();
             double N2, N3;
             double dN = 1.0 / (double)_gridSize;

code/src/toolboxes/datagenerator2D.cpp

@@ -16,6 +16,8 @@
 DataGenerator2D::DataGenerator2D( Config* settings ) : DataGeneratorBase( settings ) {
     ComputeMoments();
 
+    AdaptBasisSize();
+
     // Initialize Training Data
     ComputeSetSize();
 
@@ -71,9 +73,8 @@ void DataGenerator2D::SampleSolutionU() {
                     continue;
                 }
                 else {
-                    _uSol[idx_set][0] = 1.0;
-                    _uSol[idx_set][1] = std::sqrt( 1 - mu * mu ) * std::cos( phi );
-                    _uSol[idx_set][2] = std::sqrt( 1 - mu * mu ) * std::sin( phi );
+                    _uSol[idx_set][0] = std::sqrt( 1 - mu * mu ) * std::cos( phi );
+                    _uSol[idx_set][1] = std::sqrt( 1 - mu * mu ) * std::sin( phi );
                 }
             }
         }

code/src/toolboxes/datagenerator3D.cpp

@@ -16,6 +16,8 @@
 DataGenerator3D::DataGenerator3D( Config* settings ) : DataGeneratorBase( settings ) {
     ComputeMoments();
 
+    AdaptBasisSize();
+
     // Initialize Training Data
     ComputeSetSize();
 

code/src/toolboxes/datageneratorbase.cpp

@@ -92,6 +92,7 @@ void DataGeneratorBase::ComputeTrainingData() {
     CheckRealizability();
 
     // --- compute alphas ---
+
     _optimizer->SolveMultiCell( _alpha, _uSol, _moments );
 
     log->info( "| Making moments realizable problems." );
@@ -199,3 +200,19 @@ void DataGeneratorBase::PrintLoadScreen() {
     log->info( "------------------------ Data Generation Starts --------------------------" );
     log->info( "| Generating {} datapoints.", _setSize );
 }
+
+void DataGeneratorBase::AdaptBasisSize() {
+    // Remove zero order Moment for dimension reduction
+
+    if( _settings->GetNormalizedSampling() ) {
+        VectorVector momentTemp = _moments;
+        _nTotalEntries          = _nTotalEntries - 1;
+        _moments                = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
+
+        for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+            for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+                _moments[idx_quad][idx_sys] = momentTemp[idx_quad][idx_sys + 1];
+            }
+        }
+    }
+}