From 9ac9217b751bc9d48591ed302c3c89aa8c1d869b Mon Sep 17 00:00:00 2001 From: Steffen Date: Mon, 14 Nov 2022 13:11:26 +0100 Subject: [PATCH 1/5] added M3 rotation --- include/optimizers/neuralnetworkoptimizer.hpp | 2 + src/optimizers/neuralnetworkoptimizer.cpp | 147 +++++++++++++++++- 2 files changed, 144 insertions(+), 5 deletions(-) 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..1b8d4561 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -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,9 @@ 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 ) { +// Rotate everything to x-Axis of first moment tensor +#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 +225,55 @@ 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,57 @@ 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 * (double)_modelServingVectorAlpha[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + idx_sys]; + } + 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] }, { u[idx_cell][6], u[idx_cell][7] } }; + alpha3[1] = { { u[idx_cell][6], u[idx_cell][7] }, { u[idx_cell][7], 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]; From 027c343a0a8104d1a94e2137ad6601f6ab653a3b Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Mon, 14 Nov 2022 13:41:52 +0100 Subject: [PATCH 2/5] update version warning for neural optimizer --- src/optimizers/neuralnetworkoptimizer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index 1b8d4561..60081910 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 From 5023eba52b8bc3aac84d937deddf86311b0905c4 Mon Sep 17 00:00:00 2001 From: Steffen Date: Mon, 14 Nov 2022 17:20:10 +0100 Subject: [PATCH 3/5] some debugging output --- src/optimizers/neuralnetworkoptimizer.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index 1b8d4561..957b0a4e 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -378,6 +378,8 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& alphaRedMirror[idx_sys] = -1 * (double)_modelServingVectorAlpha[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + idx_sys]; } + + std::cout << alphaRed[idx_sys] << ":" << alphaRedMirror[idx_sys] << "\n ----------\n"; alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } From 796d631d75d55fd5e766a1c665dfdbb72611d058 Mon Sep 17 00:00:00 2001 From: Steffen Date: Mon, 14 Nov 2022 18:42:24 +0100 Subject: [PATCH 4/5] added scaling bug fix --- src/optimizers/neuralnetworkoptimizer.cpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index 66047d9f..6a96e262 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -252,10 +252,17 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 7] = (float)( u3[1][1][0] ); _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 8] = (float)( u3[1][1][1] ); + // std::cout << u1[0] << "," << u1[2] << "," << u2( 0, 0 ) << "," << u2( 0, 1 ) << "," << u2( 1, 1 ) << "," << u3[0][0][0] << "," + // << u3[1][0][0] << "," << u3[1][1][0] << "," << u3[1][1][1] << "\n"; + // Rotate Moment by 180 degrees and save mirrored moment u1 = RotateM1( u1, rot180 ); u2 = RotateM2( u2, rot180, rot180 ); u3 = RotateM3( u3, rot180 ); + // std::cout << u1[0] << "," << u1[2] << "," << u2( 0, 0 ) << "," << u2( 0, 1 ) << "," << u2( 1, 1 ) << "," << u3[0][0][0] << "," + //<< u3[1][0][0] << "," << u3[1][1][0] << "," << u3[1][1][1] << "\n"; + // std::cout << "-----\n"; + // 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 @@ -339,6 +346,7 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } + // std::cout << alphaRed - alphaRedMirror << "-----\n"; // Rotate Back Vector alpha1{ alphaRed[0], alphaRed[1] }; @@ -367,7 +375,7 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } } else { // Degree 3 - //#pragma omp parallel for +#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 @@ -375,21 +383,19 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& 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 * (double)_modelServingVectorAlpha[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 ) + idx_sys]; + alphaRedMirror[idx_sys] *= -1; } - - std::cout << alphaRed[idx_sys] << ":" << alphaRedMirror[idx_sys] << "\n ----------\n"; alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } + // std::cout << alphaRed - alphaRedMirror << "-----\n"; // 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] }, { u[idx_cell][6], u[idx_cell][7] } }; - alpha3[1] = { { u[idx_cell][6], u[idx_cell][7] }, { u[idx_cell][7], u[idx_cell][8] } }; + 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 From 9a1c1b2a6a453f348bb6c0080cebd7352ea18a50 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steffen=20Schotth=C3=B6fer?= Date: Wed, 16 Nov 2022 11:51:38 +0100 Subject: [PATCH 5/5] remove debugging output --- src/optimizers/neuralnetworkoptimizer.cpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index 6a96e262..9b0fb103 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -195,7 +195,6 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } } else if( _settings->GetMaxMomentDegree() == 2 ) { -// Rotate everything to x-Axis of first moment tensor #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] }; @@ -252,16 +251,10 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 7] = (float)( u3[1][1][0] ); _modelServingVectorU[idx_cell * ( _nSystem - 1 ) + 8] = (float)( u3[1][1][1] ); - // std::cout << u1[0] << "," << u1[2] << "," << u2( 0, 0 ) << "," << u2( 0, 1 ) << "," << u2( 1, 1 ) << "," << u3[0][0][0] << "," - // << u3[1][0][0] << "," << u3[1][1][0] << "," << u3[1][1][1] << "\n"; - // Rotate Moment by 180 degrees and save mirrored moment u1 = RotateM1( u1, rot180 ); u2 = RotateM2( u2, rot180, rot180 ); u3 = RotateM3( u3, rot180 ); - // std::cout << u1[0] << "," << u1[2] << "," << u2( 0, 0 ) << "," << u2( 0, 1 ) << "," << u2( 1, 1 ) << "," << u3[0][0][0] << "," - //<< u3[1][0][0] << "," << u3[1][1][0] << "," << u3[1][1][1] << "\n"; - // std::cout << "-----\n"; // mirror matrix is symmetric _modelServingVectorU[( _settings->GetNCells() + idx_cell ) * ( _nSystem - 1 )] = (float)( u1[0] ); // Only first moment is mirrored @@ -346,7 +339,6 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } - // std::cout << alphaRed - alphaRedMirror << "-----\n"; // Rotate Back Vector alpha1{ alphaRed[0], alphaRed[1] }; @@ -387,7 +379,6 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } - // std::cout << alphaRed - alphaRedMirror << "-----\n"; // Rotate Back Vector alpha1{ alphaRed[0], alphaRed[1] };