Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

M3 rotations #35

Open
wants to merge 7 commits into
base: develop
Choose a base branch
from
Open

M3 rotations #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9ac9217
added M3 rotation
Nov 14, 2022
027c343
update version warning for neural optimizer
ScSteffen Nov 14, 2022
5023eba
some debugging output
Nov 14, 2022
3b1e55c
Merge branch 'M3_rotations' of github.com:CSMMLab/KiT-RT into M3_rota…
Nov 14, 2022
796d631
added scaling bug fix
Nov 14, 2022
9a1c1b2
remove debugging output
Nov 16, 2022
d94d933
Merge branch 'develop' into M3_rotations
ScSteffen Jan 21, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions include/optimizers/neuralnetworkoptimizer.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<VectorVector> RotateM3( std::vector<VectorVector>& vec, Matrix& R );
};
#else
// Dummy class
Expand Down
148 changes: 142 additions & 6 deletions src/optimizers/neuralnetworkoptimizer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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<VectorVector> NeuralNetworkOptimizer::RotateM3( std::vector<VectorVector>& vec, Matrix& R ) {
std::vector<VectorVector> 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<VectorVector> 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();
Expand All @@ -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] } };
Expand Down Expand Up @@ -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<VectorVector> 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++ ) {
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<VectorVector> 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];
Expand Down