renamed sphericalbase class. code streaminlinging. Datagnerator is now able to use monomial basis

Former-commit-id: 22078d08

code/include/toolboxes/datagenerator.h

@@ -10,7 +10,7 @@
 #include "common/typedef.h"
 #include <vector>
 
-class SphericalHarmonics;
+class SphericalBase;
 class QuadratureBase;
 class Config;
 class NewtonOptimizer;
@@ -53,8 +53,8 @@ class nnDataGenerator
     Vector _weights;                /*!  @brief quadrature weights, dim(_weights) = (_nq) */
     VectorVector _quadPointsSphere; /*!  @brief (my,phi), dim(_quadPoints) = (_nq,2) */
 
-    SphericalHarmonics* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
-    VectorVector _moments;      /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
+    SphericalBase* _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
+    VectorVector _moments; /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
 
     NewtonOptimizer* _optimizer; /*! @brief: Class to solve minimal entropy problem */
     EntropyBase* _entropy;       /*! @brief: Class to handle entropy functional evaluations */

code/include/toolboxes/sphericalbasisbase.h → code/include/toolboxes/sphericalbase.h

 /*!
- * @file sphericalbasisbase.h
+ * @file sphericalbase.h
  * @brief Base Class to handle basis classes on the unit sphere
  * @author S. Schotthöfer
  *
  */
 
-#ifndef SPHERICALBASISBASE_H
-#define SPHERICALBASISBASE_H
+#ifndef SPHERICALBASE_H
+#define SPHERICALBASE_H
 
 #include "common/typedef.h"
 class Config;
 
-class SphericalBasisBase
+class SphericalBase
 {
   public:
-    SphericalBasisBase() {}
-    ~SphericalBasisBase() {}
+    SphericalBase() {}
+    ~SphericalBase() {}
 
     /*! @brief: Create a set of basis functions on the unit sphere defined in settings
      *  @param: Pointer to the config file
      *  @returns: Pointer to the createt basis class */
-    static SphericalBasisBase* Create( Config* settings );
+    static SphericalBase* Create( Config* settings );
 
     /*! @brief  : Computes all N basis functions at point (my, phi)
      *  @param  : my = cos(theta) - spherical coordinate, -1 <= x <= 1
@@ -34,6 +34,8 @@ class SphericalBasisBase
      *  @return : vector of basis functions at point (x,y,z) with size N
      */
     virtual Vector ComputeSphericalBasis( double x, double y, double z ) = 0;
+
+    virtual unsigned GetBasisSize() = 0;
 };
 
-#endif    // SPHERICALBASISBASE_H
+#endif    // SPHERICALBASE_H

code/include/toolboxes/sphericalharmonics.h

@@ -12,10 +12,10 @@
 #define SPHERICALHARMONICS_H
 
 #include "common/typedef.h"
-#include "toolboxes/sphericalbasisbase.h"
+#include "toolboxes/sphericalbase.h"
 #include <vector>
 
-class SphericalHarmonics : public SphericalBasisBase
+class SphericalHarmonics : public SphericalBase
 {
   public:
     /*! @brief : Sets up class for spherical harmonics basis based on legendre
@@ -44,7 +44,7 @@ class SphericalHarmonics : public SphericalBasisBase
      *  @param : l_degree - current degree of basis function, 0 <= l <= L
      *  @param : k_order  - current order of basis function,  -l <= k <= l
      */
-    unsigned inline GlobalIdxBasis( unsigned l_degree, unsigned k_order ) { return k_order + l_degree + l_degree * l_degree; }
+    unsigned GlobalIdxBasis( unsigned l_degree, unsigned k_order );
 
     /*! @brief : Computes an entire set of (komplex congjugate) P_l^k and stores
      *           it in the vector _assLegendreP
@@ -53,6 +53,9 @@ class SphericalHarmonics : public SphericalBasisBase
      */
     std::vector<double> GetAssLegendrePoly( const double my );
 
+    /*! @brief: Returns length of the basis, i.e. number of elements of the basis */
+    unsigned GetBasisSize() override;
+
   private:
     /*! @brief: maximal degree of the spherical harmonics basis (this is "L" in the comments)*/
     unsigned _LMaxDegree;

code/include/toolboxes/sphericalmonomials.h

@@ -8,10 +8,10 @@
 #define SPHERICALMONOMIALS_H
 
 #include "common/typedef.h"
-#include "toolboxes/sphericalbasisbase.h"
+#include "toolboxes/sphericalbase.h"
 #include <vector>
 
-class SphericalMonomials : public SphericalBasisBase
+class SphericalMonomials : public SphericalBase
 {
   public:
     /*! @brief : Sets up class for monomial basis on sphere up to degree L.
@@ -34,6 +34,11 @@ class SphericalMonomials : public SphericalBasisBase
      */
     Vector ComputeSphericalBasis( double x, double y, double z ) override;
 
+    /*! @brief: Computes the length of the basis vector for a given max degree and
+     *          spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
+     *  @return: lenght of whole basis */
+    unsigned GetBasisSize() override;
+
   private:
     /*! @brief: maximal degree of the spherical monomial basis (this is "L" in the comments)*/
     unsigned _LMaxDegree;
@@ -52,11 +57,6 @@ class SphericalMonomials : public SphericalBasisBase
      *  @return: lenght of a single dimension */
     unsigned ComputeDimensionSize( unsigned degree );
 
-    /*! @brief: Computes the length of the basis vector for a given max degree and
-     *          spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
-     *  @return: lenght of whole basis */
-    unsigned ComputeBasisSize( unsigned degree );
-
     /*! @brief: Function to compute factorial of n (n!) in recursive manner */
     unsigned Factorial( unsigned n );
 

code/src/main.cpp

@@ -12,6 +12,8 @@
 #include "common/io.h"
 #include "solvers/solverbase.h"
 
+#include "toolboxes/datagenerator.h"
+
 int main( int argc, char** argv ) {
     MPI_Init( &argc, &argv );
     wchar_t* program = Py_DecodeLocale( argv[0], NULL );
@@ -25,12 +27,20 @@ int main( int argc, char** argv ) {
     // Print input file and run info to file
     PrintLogHeader( filename );
 
-    // Build solver
-    Solver* solver = Solver::Create( config );
-
-    // Run solver and export
-    solver->Solve();
-    solver->PrintVolumeOutput();
+    if( config->GetDataGeneratorMode() ) {
+        // Build Data generator
+        nnDataGenerator* datagen = new nnDataGenerator( config );
+        // Generate Data and export
+        datagen->computeTrainingData();
+    }
+    else {
+        // Build solver
+        Solver* solver = Solver::Create( config );
+
+        // Run solver and export
+        solver->Solve();
+        solver->PrintVolumeOutput();
+    }
 
     MPI_Finalize();
     return EXIT_SUCCESS;

code/src/toolboxes/datagenerator.cpp

@@ -9,7 +9,8 @@
 #include "entropies/entropybase.h"
 #include "optimizers/newtonoptimizer.h"
 #include "quadratures/quadraturebase.h"
-#include "toolboxes/sphericalharmonics.h"
+#include "toolboxes/errormessages.h"
+#include "toolboxes/sphericalbase.h"
 
 #include "spdlog/spdlog.h"
 
@@ -19,8 +20,7 @@ nnDataGenerator::nnDataGenerator( Config* settings ) {
     _settings = settings;
     _setSize  = settings->GetTrainingDataSetSize();
 
-    _LMaxDegree    = settings->GetMaxMomentDegree();
-    _nTotalEntries = (unsigned)GlobalIndex( _LMaxDegree, _LMaxDegree ) + 1;
+    _LMaxDegree = settings->GetMaxMomentDegree();
 
     // Quadrature
     _quadrature       = QuadratureBase::Create( settings );
@@ -31,8 +31,15 @@ nnDataGenerator::nnDataGenerator( Config* settings ) {
     _settings->SetNQuadPoints( _nq );
 
     // Spherical Harmonics
-    _basis   = new SphericalHarmonics( _LMaxDegree );
+    if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS && _LMaxDegree > 0 ) {
+        ErrorMessages::Error( "No sampling algorithm for spherical harmonics basis with degree higher than 0 implemented", CURRENT_FUNCTION );
+    }
+    _basis = SphericalBase::Create( _settings );
+
+    _nTotalEntries = _basis->GetBasisSize();
+
     _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
+
     ComputeMoments();
 
     // Optimizer

code/src/toolboxes/sphericalbasisbase.cpp → code/src/toolboxes/sphericalbase.cpp

 /*!
- * @file sphericalbasisbase.cpp
+ * @file sphericalbase.cpp
  * @brief Base Class to handle basis classes on the unit sphere
  * @author S. Schotthöfer
- *
  */
 
-#include "toolboxes/sphericalbasisbase.h"
 #include "common/config.h"
 #include "toolboxes/errormessages.h"
+#include "toolboxes/sphericalbase.h"
 #include "toolboxes/sphericalharmonics.h"
 #include "toolboxes/sphericalmonomials.h"
 
-SphericalBasisBase* SphericalBasisBase::Create( Config* settings ) {
+SphericalBase* SphericalBase::Create( Config* settings ) {
     SPHERICAL_BASIS_NAME name = settings->GetSphericalBasisName();
     unsigned maxMomentDegree  = settings->GetMaxMomentDegree();
 

code/src/toolboxes/sphericalharmonics.cpp

@@ -12,10 +12,11 @@ SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
 
     ComputeCoefficients();
 
-    unsigned basisSize = GlobalIdxBasis( _LMaxDegree, _LMaxDegree ) + 1;
-    _YBasis            = Vector( basisSize, 0.0 );
+    _YBasis = Vector( GetBasisSize(), 0.0 );
 }
 
+unsigned SphericalHarmonics::GetBasisSize() { return GlobalIdxBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */ }
+
 Vector SphericalHarmonics::ComputeSphericalBasis( double my, double phi ) {
     ComputeAssLegendrePoly( my );
     ComputeYBasis( phi );
@@ -38,6 +39,8 @@ Vector SphericalHarmonics::ComputeSphericalBasis( double x, double y, double z )
     return _YBasis;
 }
 
+unsigned SphericalHarmonics::GlobalIdxBasis( unsigned l_degree, unsigned k_order ) { return k_order + l_degree + l_degree * l_degree; }
+
 std::vector<double> SphericalHarmonics::GetAssLegendrePoly( const double my ) {
     ComputeAssLegendrePoly( my );
     return _assLegendreP;
@@ -117,12 +120,4 @@ void SphericalHarmonics::ComputeYBasis( const double phi ) {
             _YBasis[GlobalIdxBasis( l_idx, k_idx )]  = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * c;
         }
     }
-
-    // slower version
-    // for( int l_idx = 1; l_idx <= int( _LMaxDegree ); l_idx++ ) {
-    //     for( int k_idx = 1; k_idx <= l_idx; k_idx++ ) {
-    //         _YBasis[GlobalIdxBasis( l_idx, -k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * sin( k_idx * phi );
-    //         _YBasis[GlobalIdxBasis( l_idx, k_idx )]  = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * cos( k_idx * phi );
-    //     }
-    // }
 }

code/src/toolboxes/sphericalmonomials.cpp

@@ -10,8 +10,7 @@ SphericalMonomials::SphericalMonomials( unsigned L_degree ) {
     _LMaxDegree = L_degree;
     _spatialDim = 3;    // Spatial dimension 2 is hardcoded
 
-    unsigned basisSize = ComputeBasisSize( L_degree );
-    _YBasis            = Vector( basisSize );
+    _YBasis = Vector( GetBasisSize() );
 }
 
 Vector SphericalMonomials::ComputeSphericalBasis( double my, double phi ) {
@@ -22,11 +21,9 @@ Vector SphericalMonomials::ComputeSphericalBasis( double my, double phi ) {
         double omega_X = Omega_x( my, phi );
         double omega_Y = Omega_y( my, phi );
         double omega_Z = Omega_z( my );
-        double tmp     = .0;
         for( unsigned a = 0; a <= idx_degree; a++ ) {
             for( unsigned b = 0; b <= idx_degree - a; b++ ) {
                 unsigned c          = idx_degree - a - b;    // c uniquely defined
-                tmp                 = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
                 _YBasis[idx_vector] = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
                 idx_vector++;
             }
@@ -48,21 +45,18 @@ Vector SphericalMonomials::ComputeSphericalBasis( double x, double y, double z )
     return ComputeSphericalBasis( my, phi );
 }
 
-/*! @brief: Computes the length of the basis vector for a given max degree and
- *          spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
- *  @return: lenght of whole basis */
-unsigned SphericalMonomials::ComputeDimensionSize( unsigned degree ) {
-    return Factorial( degree + _spatialDim - 1 ) / ( Factorial( degree ) * Factorial( _spatialDim - 1 ) );
-}
-
-unsigned SphericalMonomials::ComputeBasisSize( unsigned degree ) {
+unsigned SphericalMonomials::GetBasisSize() {
     unsigned basisLen = 0;
-    for( unsigned idx_degree = 0; idx_degree <= degree; idx_degree++ ) {
+    for( unsigned idx_degree = 0; idx_degree <= _LMaxDegree; idx_degree++ ) {
         basisLen += ComputeDimensionSize( idx_degree );
     }
     return basisLen;
 }
 
+unsigned SphericalMonomials::ComputeDimensionSize( unsigned degree ) {
+    return Factorial( degree + _spatialDim - 1 ) / ( Factorial( degree ) * Factorial( _spatialDim - 1 ) );
+}
+
 unsigned SphericalMonomials::Factorial( unsigned n ) { return ( n == 1 || n == 0 ) ? 1 : Factorial( n - 1 ) * n; }
 
 double SphericalMonomials::Omega_x( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }

code/tests/test_sphericalBasis.cpp

@@ -6,7 +6,6 @@
 #include "toolboxes/sphericalmonomials.h"
 
 #include <fstream>
-#include <iostream>
 #include <sstream>
 
 double Y0_0( double, double ) { return sqrt( 1 / ( 4 * M_PI ) ); }