diff --git a/include/optimizers/neuralnetworkoptimizer.hpp b/include/optimizers/neuralnetworkoptimizer.hpp index 4a66390b..366374c9 100644 --- a/include/optimizers/neuralnetworkoptimizer.hpp +++ b/include/optimizers/neuralnetworkoptimizer.hpp @@ -47,6 +47,8 @@ class NeuralNetworkOptimizer : public OptimizerBase Vector RotateM1( Vector& vec, Matrix& R ); /*!< @brief Rotates the M1 part of a 2D moment vector using a rotation matrix R */ /*!< @brief Rotates the tensorized M2 part of a 2D moment vector using a rotation matrix R */ Matrix RotateM2( Matrix& vec, Matrix& R, Matrix& Rt ); + /*!< @brief Rotates the tensorized M3 part of a 2D moment vector using a rotation matrix R */ + std::vector RotateM3( std::vector& vec, Matrix& R ); }; #else // Dummy class diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index b89d4ec6..9b0fb103 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -65,7 +65,7 @@ NeuralNetworkOptimizer::NeuralNetworkOptimizer( Config* settings ) : OptimizerBa _tfModel = new cppflow::model( tfModelPath ); // load model unsigned servingSize = _settings->GetNCells(); if( _settings->GetEnforceNeuralRotationalSymmetry() ) { - if( _settings->GetMaxMomentDegree() > 2 ) { + if( _settings->GetMaxMomentDegree() > 3 ) { ErrorMessages::Error( "This postprocessing step is currently only for M1 and M2 models available.", CURRENT_FUNCTION ); } servingSize *= 2; // Double number of vectors, since we mirror the rotated vector @@ -130,6 +130,46 @@ Vector NeuralNetworkOptimizer::RotateM1( Vector& vec, Matrix& R ) { return R * v Matrix NeuralNetworkOptimizer::RotateM2( Matrix& vec, Matrix& R, Matrix& Rt ) { return R * vec * Rt; } +std::vector NeuralNetworkOptimizer::RotateM3( std::vector& vec, Matrix& R ) { + std::vector res( vec ); + + //----- t3 tensor-rotation + for( unsigned j_1 = 0; j_1 < 2; j_1++ ) { + for( unsigned j_2 = 0; j_2 < 2; j_2++ ) { + for( unsigned i_3 = 0; i_3 < 2; i_3++ ) { + res[j_1][j_2][i_3] = 0; + for( unsigned j_3 = 0; j_3 < 2; j_3++ ) { + res[j_1][j_2][i_3] += R( i_3, j_3 ) * vec[j_1][j_2][j_3]; + } + } + } + } + std::vector res2( res ); + //----- t2 tensor-rotation + for( unsigned j_1 = 0; j_1 < 2; j_1++ ) { + for( unsigned i_2 = 0; i_2 < 2; i_2++ ) { + for( unsigned i_3 = 0; i_3 < 2; i_3++ ) { + res2[j_1][i_2][i_3] = 0; + for( unsigned j_2 = 0; j_2 < 2; j_2++ ) { + res2[j_1][i_2][i_3] += R( i_2, j_2 ) * res[j_1][j_2][i_3]; + } + } + } + } + //----- t1 tensor-rotation + for( unsigned i_1 = 0; i_1 < 2; i_1++ ) { + for( unsigned i_2 = 0; i_2 < 2; i_2++ ) { + for( unsigned i_3 = 0; i_3 < 2; i_3++ ) { + res[i_1][i_2][i_3] = 0; + for( unsigned j_1 = 0; j_1 < 2; j_1++ ) { + res[i_1][i_2][i_3] += R( i_1, j_1 ) * res2[j_1][i_2][i_3]; + } + } + } + } + return res; +} + void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& /*moments*/ ) { unsigned servingSize = _settings->GetNCells(); @@ -154,9 +194,8 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 1] = (float)( u1[1] ); // should be zero } } - else { - // Rotate everything to x-Axis of first moment tensor - //#pragma omp parallel for + else if( _settings->GetMaxMomentDegree() == 2 ) { +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { Vector u1{ u[idx_cell][1], u[idx_cell][2] }; Matrix u2{ { u[idx_cell][3], u[idx_cell][4] }, { u[idx_cell][4], u[idx_cell][5] } }; @@ -185,9 +224,56 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 4] = (float)( u2( 1, 1 ) ); } } + else { // M3 +#pragma omp parallel for + for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { + Vector u1{ u[idx_cell][1], u[idx_cell][2] }; + Matrix u2{ { u[idx_cell][3], u[idx_cell][4] }, { u[idx_cell][4], u[idx_cell][5] } }; + std::vector u3( 2, VectorVector( 2, Vector( 2, 0.0 ) ) ); + // fill u3 tensor + u3[0] = { { u[idx_cell][6], u[idx_cell][7] }, { u[idx_cell][7], u[idx_cell][8] } }; + u3[1] = { { u[idx_cell][7], u[idx_cell][8] }, { u[idx_cell][8], u[idx_cell][9] } }; + + _rotationMats[idx_cell] = CreateRotator( u1 ); + _rotationMatsT[idx_cell] = blaze::trans( _rotationMats[idx_cell] ); + + u1 = RotateM1( u1, _rotationMats[idx_cell] ); + u2 = RotateM2( u2, _rotationMats[idx_cell], _rotationMatsT[idx_cell] ); + u3 = RotateM3( u3, _rotationMats[idx_cell] ); + + _modelServingVectorU[idx_cell * ( _nSystem - 1 )] = (float)( u1[0] ); + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 1] = (float)( u1[1] ); // should be zero + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 2] = (float)( u2( 0, 0 ) ); // u2 part + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 3] = (float)( u2( 0, 1 ) ); + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 4] = (float)( u2( 1, 1 ) ); + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 5] = (float)( u3[0][0][0] ); // u3 part + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 6] = (float)( u3[1][0][0] ); + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 7] = (float)( u3[1][1][0] ); + _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 8] = (float)( u3[1][1][1] ); + + // Rotate Moment by 180 degrees and save mirrored moment + u1 = RotateM1( u1, rot180 ); + u2 = RotateM2( u2, rot180, rot180 ); + u3 = RotateM3( u3, rot180 ); + + // mirror matrix is symmetric + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 )] = (float)( u1[0] ); // Only first moment is mirrored + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 1] = (float)( u1[1] ); // should be zero + + // Even Moments cancel rotation + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 2] = (float)( u2( 0, 0 ) ); + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 3] = (float)( u2( 0, 1 ) ); + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 4] = (float)( u2( 1, 1 ) ); + // u3 part + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 5] = (float)( u3[0][0][0] ); + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 6] = (float)( u3[1][0][0] ); + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 7] = (float)( u3[1][1][0] ); + _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + 8] = (float)( u3[1][1][1] ); + } + } servingSize *= 2; } - else { // No Postprocessing + else { // No Preprocessing #pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { for( unsigned idx_sys = 0; idx_sys < _nSystem - 1; idx_sys++ ) { @@ -239,7 +325,7 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } } } - else { // Degree 2 + else if( _settings->GetMaxMomentDegree() == 2 ) { // Degree 2 #pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { Vector alphaRed = Vector( _nSystem - 1, 0.0 ); // local reduced alpha @@ -274,6 +360,56 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } alpha[idx_cell][0] = -log( integral ); // log trafo + // Store output + for( unsigned idx_sys = 0; idx_sys < _nSystem - 1; idx_sys++ ) { + alpha[idx_cell][idx_sys + 1] = alphaRed[idx_sys]; + } + } + } + else { // Degree 3 +#pragma omp parallel for + for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { + Vector alphaRed = Vector( _nSystem - 1, 0.0 ); // local reduced alpha + Vector alphaRedMirror = Vector( _nSystem - 1, 0.0 ); // local reduced mirrored alpha + for( unsigned idx_sys = 0; idx_sys < _nSystem - 1; idx_sys++ ) { + alphaRed[idx_sys] = (double)_modelServingVectorAlpha[idx_cell * ( _nSystem - 1 ) + idx_sys]; + alphaRedMirror[idx_sys] = (double)_modelServingVectorAlpha[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + idx_sys]; + if( idx_sys < 2 || idx_sys > 4 ) { // Miror order 1 moments back and Miror order 3 moments back + alphaRedMirror[idx_sys] *= -1; + } + alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) + } + + // Rotate Back + Vector alpha1{ alphaRed[0], alphaRed[1] }; + Matrix alpha2{ { alphaRed[2], alphaRed[3] * 0.5 }, { alphaRed[3] * 0.5, alphaRed[4] } }; + std::vector alpha3( 2, VectorVector( 2, Vector( 2, 0.0 ) ) ); + // fill u3 tensor + alpha3[0] = { { u[idx_cell][5], u[idx_cell][6] / 3.0 }, { u[idx_cell][6] / 3.0, u[idx_cell][7] / 3.0 } }; + alpha3[1] = { { u[idx_cell][6] / 3.0, u[idx_cell][7] / 3.0 }, { u[idx_cell][7] / 3.0, u[idx_cell][8] } }; + + alpha1 = RotateM1( alpha1, _rotationMatsT[idx_cell] ); // Rotate Back + alpha2 = RotateM2( alpha2, _rotationMatsT[idx_cell], _rotationMats[idx_cell] ); // Rotate Back + alpha3 = RotateM3( alpha3, _rotationMatsT[idx_cell] ); // Rotate Back + + // Store back-rotated alpha + alphaRed[0] = alpha1[0]; + alphaRed[1] = alpha1[1]; + alphaRed[2] = alpha2( 0, 0 ); + alphaRed[3] = 2 * alpha2( 1, 0 ); + alphaRed[4] = alpha2( 1, 1 ); + alphaRed[5] = alpha3[0][0][0]; + alphaRed[6] = 3 * alpha3[1][0][0]; + alphaRed[7] = 3 * alpha3[1][1][0]; + alphaRed[8] = alpha3[1][1][1]; + + // Restore alpha_0 + double integral = 0.0; + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + integral += _entropy->EntropyPrimeDual( dot( alphaRed, _reducedMomentBasis[idx_quad] ) ) * _weights[idx_quad]; + } + alpha[idx_cell][0] = -log( integral ); // log trafo + // Store output for( unsigned idx_sys = 0; idx_sys < _nSystem - 1; idx_sys++ ) { alpha[idx_cell][idx_sys + 1] = alphaRed[idx_sys];