newton optimizer compiles, not yet tested. test is faulty

code/include/entropies/entropybase.h

@@ -19,10 +19,8 @@ class EntropyBase
      *  @returns: value of entropy functional at z */
     virtual double Entropy( double z ) = 0;
 
-    /*! @brief: computes the derivative of the entropy functional
-     *  @param: z = point where the derivative should be evaluated.
-     *          z must be in domain of the functional
-     *  @returns: value of derivative at z */
+    /*! @brief: computes eta'(z).z must be in domain of the functional.
+     *  @param: z = point where the derivative should be evaluated. */
     virtual double EntropyPrime( double z ) = 0;
 
     /*! @brief: computes the dual of the  entropy functional
@@ -31,16 +29,22 @@ class EntropyBase
      *  @returns: value of entropy functional at z */
     virtual double EntropyDual( double y ) = 0;
 
-    /*! @brief: computes the  dual of the derivative of the entropy functional
-     *  @param: z = point where the dual of the  derivative should be evaluated.
-     *          z must be in domain of the dual functional
-     *  @returns: value of dual of the derivative at z */
-    virtual double EntropyPrimeDual( double y ) = 0;
+    // /*! @brief: computes m*eta_*'(alpha*m). alpha*m must be in domain of the functional
+    //  *          alpha, m and grad must be of same size.
+    //  *  @param: alpha = point where the derivative should be evaluated.
+    //  *  @param: m = moment vector
+    //  *  @param: grad: vector in which the resulting gradient is stored */
+    // virtual void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) = 0;
+
+    /*! @brief: computes eta_*'(y).
+     *  @param: y = point where the derivative should be evaluated.
+     *  @returns: value of the derivative of the  entropy functional at y */
+    virtual double EntropyPrimeDual( double z ) = 0;
 
     /*! @brief: computes the hessian of the dual entropy functional
      *  @param: y = point where the hessian should be evaluated;
      *  @returns: value of the hessian at y */
-    // virtual double EntropyHessianDual( double y ) = 0;
+    virtual double EntropyHessianDual( double y ) = 0;
 
     /*! @brief: checks, if value is in domain of entropy */
     virtual bool CheckDomain( double z ) = 0;

code/include/entropies/maxwellboltzmannentropy.h

@@ -15,10 +15,16 @@ class MaxwellBoltzmannEntropy : public EntropyBase
 
     inline double EntropyPrime( double z ) override { return log( z ); }
 
+    // inline void EntropyPrime( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * log( dot( alpha, m ) ); }
+
     inline double EntropyDual( double y ) override { return exp( y ); }
 
+    // inline void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * exp( dot( alpha, m ) ); }
+
     inline double EntropyPrimeDual( double y ) override { return exp( y ); }
 
+    inline double EntropyHessianDual( double y ) override { return exp( y ); }
+
     inline bool CheckDomain( double z ) override { return ( z >= 0 && std::isnan( z ) ); }
 };
 

code/include/entropies/quadraticentropy.h

@@ -17,9 +17,11 @@ class QuadraticEntropy : public EntropyBase
 
     inline double EntropyDual( double y ) override { return 0.5 * y * y; }
 
+    // inline void EntropyPrimeDual( Vector& alpha, Vector& m, Vector& grad ) override { grad = m * dot( alpha, m ); }
+
     inline double EntropyPrimeDual( double y ) override { return y; }
 
-    // inline double EntropyHessianDual( double y ) override { return y; }
+    inline double EntropyHessianDual( double y ) override { return 1.0; }
 
     inline bool CheckDomain( double z ) override { return std::isnan( z ); }
 };

code/include/optimizers/newtonoptimizer.h

@@ -3,14 +3,34 @@
 
 #include "optimizerbase.h"
 
+class QuadratureBase;
+
 class NewtonOptimizer : public OptimizerBase
 {
   public:
-    inline NewtonOptimizer( Config* settings ) : OptimizerBase( settings ) {}
+    NewtonOptimizer( Config* settings );
 
     inline ~NewtonOptimizer() {}
 
-    void Solve( Vector& lambda, Vector& u ) override;
+    void Solve( Vector& lambda, Vector& u, VectorVector& moments ) override;
+
+  private:
+    /*! @brief: Computes gradient of objective function and stores it in grad
+                grad = <m*eta*'(alpha*m)> - sol */
+    void ComputeGradient( Vector& alpha, Vector& sol, VectorVector& moments, Vector& grad );
+
+    /*! @brief: Computes hessian of objective function and stores it in hessian
+                grad = <mXm*eta*'(alpha*m)> */
+    void ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian );
+
+    QuadratureBase* _quadrature; /*! @brief: used quadrature */    // THis is memory doubling! Try to use a pointer.
+    unsigned _nq;                                                  /*! @brief: number of quadrature points */
+    Vector _weights;                                               /*!  @brief quadrature weights, dim(_weights) = (_nq) */
+    VectorVector _quadPointsSphere;                                /*!  @brief (my,phi), dim(_quadPoints) = (_nq,2) */
+
+    unsigned short _maxIterations;   /*! @brief: Max iterations of the newton solver */
+    double _alpha;                   /*! @brief: Newton Step Size */
+    unsigned short _maxLineSearches; /*! @brief: Max amount of line searches for Newton Algo */
 };
 
 #endif    // NEWTONOPTIMIZER_H

code/include/optimizers/optimizerbase.h

@@ -19,7 +19,7 @@ class OptimizerBase
     /*! @brief  : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem
      *  @param  : Vector  u = pointer to vector of given moments. // Maybe use pointer for performance?
      *  @return : Vector  alpha = optimal lagrange multipliers. Has the same length as Vector u. */
-    virtual void Solve( Vector& lambda, Vector& u ) = 0;
+    virtual void Solve( Vector& lambda, Vector& u, VectorVector& moments ) = 0;
 
   protected:
     EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */

code/include/settings/config.h

@@ -82,7 +82,9 @@ class Config
     KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
 
     // Optimizer
-    OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
+    OPTIMIZER_NAME _entropyOptimizerName;   /*!< @brief Choice of optimizer */
+    double _optimizerEpsilon;               /*!< @brief termination criterion epsilon for Newton Optmizer */
+    unsigned short _maxIterNewtonOptimizer; /*!< @brief Maximal Number of newton iterations */
 
     // --- Parsing Functionality and Initializing of Options ---
     /*!
@@ -228,7 +230,8 @@ class Config
 
     // Optimizer
     OPTIMIZER_NAME inline GetOptimizerName() const { return _entropyOptimizerName; }
-
+    double inline GetNewtonOptimizerEpsilon() const { return _optimizerEpsilon; }
+    unsigned inline GetMaxIterNewtonOptimizer() const { return _maxIterNewtonOptimizer; }
     // Boundary Conditions
     BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
 

code/src/optimizers/newtonoptimizer.cpp

 
 #include "optimizers/newtonoptimizer.h"
+#include "quadratures/quadraturebase.h"
 #include "settings/config.h"
 
-void NewtonOptimizer::Solve( Vector& lambda, Vector& u ) {
+NewtonOptimizer::NewtonOptimizer( Config* settings ) : OptimizerBase( settings ) {
+    _quadrature       = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
+    _nq               = _quadrature->GetNq();
+    _weights          = _quadrature->GetWeights();
+    _quadPointsSphere = _quadrature->GetPointsSphere();
+    _maxIterations    = settings->GetMaxIterNewtonOptimizer();
+    _alpha            = 0.1;     // settings->GetNewtonStepSize();
+    _maxLineSearches  = 1000;    // settings->GetMaxLineSearches();
+}
+
+void NewtonOptimizer::ComputeGradient( Vector& alpha, Vector& sol, VectorVector& moments, Vector& grad ) {
+
+    // Reset Vector
+    for( unsigned idx_sys; idx_sys < grad.size(); idx_sys++ ) {
+        grad[idx_sys] = 0.0;
+    }
+
+    // Integrate
+    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+        grad += moments[idx_quad] * ( _entropy->EntropyPrimeDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] );
+    }
+    grad -= sol;
+}
+
+void NewtonOptimizer::ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian ) {
+    // Reset Matrix
+    for( unsigned idx_Row; idx_Row < moments.size(); idx_Row++ ) {
+        for( unsigned idx_Col; idx_Col < moments.size(); idx_Col++ ) {
+            hessian( idx_Row, idx_Col ) = 0.0;
+        }
+    }
+
+    // Integrate
+    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+        hessian +=
+            outer( moments[idx_quad], moments[idx_quad] ) * ( _entropy->EntropyHessianDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] );
+    }
+}
+
+void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments ) {
+
+    /* solve the problem argmin ( <eta(alpha*m)>-alpha*u))
+     * where alpha = Lagrange multiplier
+     *           m = moment basis
+     *           u = current "moment solution"
+     */
 
     // if we have quadratic entropy, then alpha = u;
-    if( _settings->GetEntropyName() == QUADRATIC ) return;
+    if( _settings->GetEntropyName() == QUADRATIC ) {
+        lambda = sol;
+        return;
+    }
 
     // Start Newton Algorithm
     //
-    //    int maxRefinements = 1000;
-    //    unsigned nTotal    = _nTotalForRef[refLevel];
-    //
-    //    Vector g( _nStates * nTotal );
-    //
-    //    // check if initial guess is good enough
-    //    Gradient( g, lambda, u, refLevel );
-    //    if( CalcNorm( g, nTotal ) < _settings->GetEpsilon() ) {
-    //        return;
-    //    }
-    //    Matrix H( _nStates * nTotal, _nStates * nTotal );
-    //    Vector dlambdaNew( _nStates * nTotal );
-    //    // calculate initial Hessian and gradient
-    //    Vector dlambda = -g;
-    //    // std::cout << g << std::endl;
-    //    Hessian( H, lambda, refLevel );
-    //
-    //    // double eps     = 0.00001;
-    //    // lambda( 1, 0 ) = lambda( 1, 0 ) + eps;
-    //
-    //    // Vector gEps( _nStates * nTotal );
-    //
-    //    // Gradient( gEps, lambda, u, refLevel );
-    //
-    //    // std::cout << "H_FD = " << ( gEps - g ) / eps << std::endl;
-    //    // std::cout << "H = " << H << std::endl;
-    //    // exit( EXIT_FAILURE );
-    //
-    //    posv( H, g );
-    //    if( _maxIterations == 1 ) {
-    //        AddMatrixVectorToMatrix( lambda, -_alpha * g, lambda, nTotal );
-    //        return;
-    //    }
-    //    Matrix lambdaNew( _nStates, nTotal );
-    //    AddMatrixVectorToMatrix( lambda, -_alpha * g, lambdaNew, nTotal );
-    //    Gradient( dlambdaNew, lambdaNew, u, refLevel );
-    //    // perform Newton iterations
-    //    for( unsigned l = 0; l < _maxIterations; ++l ) {
-    //        double stepSize = 1.0;
-    //        if( l != 0 ) {
-    //            Gradient( g, lambda, u, refLevel );
-    //            dlambda = -g;
-    //            Hessian( H, lambda, refLevel );
-    //            posv( H, g );
-    //            AddMatrixVectorToMatrix( lambda, -stepSize * _alpha * g, lambdaNew, nTotal );
-    //            Gradient( dlambdaNew, lambdaNew, u, refLevel );
-    //        }
-    //        int refinementCounter = 0;
-    //        // std::cout << CalcNorm( dlambdaNew, nTotal ) << std::endl;
-    //        while( CalcNorm( dlambda, nTotal ) < CalcNorm( dlambdaNew, nTotal ) || !std::isfinite( CalcNorm( dlambdaNew, nTotal ) ) ) {
-    //            stepSize *= 0.5;
-    //            AddMatrixVectorToMatrix( lambda, -stepSize * _alpha * g, lambdaNew, nTotal );
-    //            Gradient( dlambdaNew, lambdaNew, u, refLevel );
-    //            if( CalcNorm( dlambdaNew, nTotal ) < _settings->GetEpsilon() ) {
-    //                lambda = lambdaNew;
-    //                return;
-    //            }
-    //            else if( ++refinementCounter > maxRefinements ) {
-    //                _log->error( "[closure] Newton needed too many refinement steps!" );
-    //                exit( EXIT_FAILURE );
-    //            }
-    //        }
-    //        lambda = lambdaNew;
-    //        if( CalcNorm( dlambdaNew, nTotal ) < _settings->GetEpsilon() ) {
-    //            lambda = lambdaNew;
-    //            return;
-    //        }
-    //    }
-    //    _log->error( "[closure] Newton did not converge!" );
-    //    exit( EXIT_FAILURE );
+    unsigned nSize = lambda.size();
 
-    // return Vector( 1, 0.0 );    // dummy
+    // int maxRefinements = 1000;
+    // unsigned nTotal    = _nTotalForRef[refLevel];
+
+    Vector grad( nSize, 0.0 );
+
+    // check if initial guess is good enough
+    ComputeGradient( lambda, sol, moments, grad );
+
+    // Gradient( g, lambda, u, refLevel );
+
+    if( norm( grad ) < _settings->GetNewtonOptimizerEpsilon() ) {
+        return;
+    }
+
+    // If not, compute the Hessian
+    Matrix H( nSize, nSize, 0.0 );
+
+    // calculate initial Hessian and gradient
+    Vector dlambda = -grad;
+
+    ComputeHessian( lambda, moments, H );
+
+    // invert Hessian
+    invert( H );
+    if( _maxIterations == 1 ) {
+        lambda = lambda - _alpha * H * grad;
+        return;
+    }
+
+    Vector lambdaNew( nSize, 0.0 );    //  New newton step
+    Vector dlambdaNew( nSize );        //  Gradient at New newton step
+
+    lambdaNew = lambda - _alpha * H * grad;
+
+    // Compute Gradient of new point;
+    ComputeGradient( lambdaNew, sol, moments, dlambdaNew );
+
+    // perform Newton iterations
+    for( unsigned l = 0; l < _maxIterations; ++l ) {
+        double stepSize = 1.0;
+        if( l != 0 ) {
+            ComputeGradient( lambda, sol, moments, grad );
+
+            dlambda = -grad;
+            ComputeHessian( lambda, moments, H );
+            invert( H );
+            lambdaNew = lambda - _alpha * H * grad;
+            ComputeGradient( lambdaNew, sol, moments, dlambdaNew );
+        }
+
+        int lineSearchCounter = 0;
+
+        // Line Search
+        while( norm( dlambda ) < norm( dlambdaNew ) || !std::isfinite( norm( dlambdaNew ) ) ) {
+            stepSize *= 0.5;
+
+            lambdaNew = lambda - stepSize * _alpha * H * grad;
+            ComputeGradient( lambdaNew, sol, moments, dlambdaNew );
+
+            // Check if FONC is fullfilled
+            if( norm( dlambdaNew ) < _settings->GetNewtonOptimizerEpsilon() ) {
+                lambda = lambdaNew;
+                return;
+            }
+            else if( ++lineSearchCounter > _maxLineSearches ) {
+                //_log->error( "[closure] Newton needed too many refinement steps!" );
+                std::cout << "[closure] Newton needed too many refinement steps!";
+                exit( EXIT_FAILURE );
+            }
+        }
+        lambda = lambdaNew;
+        if( norm( dlambdaNew ) < _settings->GetNewtonOptimizerEpsilon() ) {
+            lambda = lambdaNew;
+            return;
+        }
+    }
+    //    _log->error( "[closure] Newton did not converge!" );
+    std::cout << "[closure] Newton did not converge!";
+    exit( EXIT_FAILURE );
 }

code/src/settings/config.cpp

@@ -221,6 +221,10 @@ void Config::SetConfigOptions() {
     AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, QUADRATIC );
     /*! @brief Optimizer Name \n DESCRIPTION:  Optimizer used to determine the minimal Entropy reconstruction \n DEFAULT NEWTON \ingroup Config */
     AddEnumOption( "ENTROPY_OPTIMIZER", _entropyOptimizerName, Optimizer_Map, NEWTON );
+    /*! @brief Newton Optimizer Epsilon \n DESCRIPTION:  Convergencce Epsilon for Newton Optimizer \n DEFAULT 1e-3 \ingroup Config */
+    AddDoubleOption( "NEWTON_EPSILON", _optimizerEpsilon, 1e-3 );
+    /*! @brief Max Iter Newton Optmizers \n DESCRIPTION: Max number of newton iterations \n DEFAULT 10 \ingroup Config */
+    AddUnsignedShortOption( "MAX_ITER_NEWTON_OPTIMIZER", _maxIterNewtonOptimizer, 10 );
 
     // Mesh related options
     // Boundary Markers

code/src/solvers/mnsolver.cpp

@@ -81,21 +81,16 @@ void MNSolver::ComputeMoments() {
 Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
 
     // ---- Integration of Moment of flux ----
-    double w, entropyL, entropyR, entropyFlux;
+    double entropyL, entropyR, entropyFlux;
 
     Vector flux( _nTotalEntries, 0.0 );
-    Vector omega( 3, 0.0 );
 
     for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
 
-        w = _weights[idx_quad];
-
         entropyFlux = 0.0;    // Reset temorary flux
 
         entropyL = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
 
-        omega = _quadPoints[idx_quad];    // Use Pointer!
-
         for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
             // Store fluxes in psiNew, to save memory
             if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neigh] == _nCells )
@@ -103,9 +98,9 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
             else {
                 entropyR = _entropy->EntropyPrimeDual( blaze::dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] ) );
             }
-            entropyFlux += _g->Flux( omega, entropyL, entropyR, _normals[idx_cell][idx_neigh] );
+            entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
         }
-        flux += _moments[idx_quad] * ( w * entropyFlux );
+        flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
     }
     return flux;
 }
@@ -150,7 +145,8 @@ void MNSolver::Solve() {
 
             // ------- Reconstruction Step -------
 
-            _alpha[idx_cell] = _sol[idx_cell];    // _optimizer->Solve( _psi[idx_cell] );
+            // _alpha[idx_cell] = _sol[idx_cell];
+            _optimizer->Solve( _alpha[idx_cell], _sol[idx_cell], _moments );
 
             // ------- Flux Computation Step ---------
 

code/tests/input/unit_newtonOptimizer.cfg

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Newton Solver Unit Test Config File %
+%  Author <Steffen Schotthöfer>		   %
+%  Date   12.07.2020                   %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%
+% ---- File specifications ----
+OUTPUT_DIR = ../../result
+OUTPUT_FILE = rtsn_test_unit
+LOG_DIR = ../../result/logs
+MESH_FILE = unit.su2
+
+% ---- Solver specifications ----
+ENTROPY_OPTIMIZER = NEWTON
+NEWTON_EPSILON = 1e-3
+MAX_ITER_NEWTON_OPTIMIZER = 100
+ENTROPY_FUNCTIONAL = QUADRATIC
+MAX_MOMENT_SOLVER = 1
+
+
+% ---- Quadrature Specs ----
+
+QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
+QUAD_ORDER = 4

code/tests/test_optimizer.cpp

+#include <numeric>
+
+#include "catch.hpp"
+#include "optimizers/optimizerbase.h"
+#include "quadratures/quadraturebase.h"
+#include "settings/config.h"
+#include "solvers/sphericalharmonics.h"
+
+TEST_CASE( "Test the Newton Optimizer", "[optimizers]" ) {
+
+    char filename[] = "../tests/input/unit_newtonOptimizer.cfg";
+
+    // Load Settings from File
+    Config* config = new Config( filename );
+
+    // Get Basis
+    SphericalHarmonics basis( config->GetMaxMomentDegree() );
+
+    // Get Quadrature
+    QuadratureBase* quad = QuadratureBase::CreateQuadrature( config->GetQuadName(), config->GetQuadOrder() );
+
+    // Get Optimizer (Newton)
+    OptimizerBase* optimizer = OptimizerBase::Create( config );
+
+    // Get dummy Moment Vector
+    unsigned nTotalEntries = basis.GlobalIdxBasis( config->GetMaxMomentDegree(), config->GetMaxMomentDegree() ) + 1;    // = 4
+    Vector u( nTotalEntries, 1.0 );
+
+    // Get inital guess for solution
+    Vector alpha( nTotalEntries, 0.0 );
+
+    // Get Moments
+    VectorVector moments = VectorVector( config->GetNQuadPoints(), Vector( nTotalEntries, 0.0 ) );
+    double my, phi;
+    VectorVector quadPointsSphere = quad->GetPointsSphere();
+    for( unsigned idx_quad = 0; idx_quad < config->GetNQuadPoints(); idx_quad++ ) {
+        my  = quadPointsSphere[idx_quad][0];
+        phi = quadPointsSphere[idx_quad][1];
+
+        moments[idx_quad] = basis.ComputeSphericalBasis( my, phi );
+    }
+
+    optimizer->Solve( alpha, u, moments );
+
+    REQUIRE( std::fabs( norm( alpha - u ) ) < 1e2 * std::numeric_limits<double>::epsilon() );    // alpha = u for quadratic entropy
+}