#100

Committed 01 Jul 2025 05:43PM UTC coverage: 46.655%. Remained the same

Build # #100

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

Merge 7e11939c4 into 1c5a8e2d5

Run Details

4422 of 9478 relevant lines covered (46.66%)

57163.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/datagenerator/datageneratorregression2D.cpp

/*!
 * \file datagenerator3D.cpp
 * \brief Class to generate data for the neural entropy closure
 * \author S. Schotthoefer
 */

#include "datagenerator/datageneratorregression2D.hpp"
#include "common/config.hpp"
#include "quadratures/quadraturebase.hpp"
#include "toolboxes/errormessages.hpp"
#include "velocitybasis/sphericalbase.hpp"

#include <iostream>
#include <omp.h>

DataGeneratorRegression2D::DataGeneratorRegression2D( Config* settings ) : DataGeneratorRegression( settings ) {
    ComputeMoments();

    // Initialize Training Data
    if( _settings->GetAlphaSampling() )
        ComputeSetSizeAlpha();
    else
        ComputeSetSizeU();

    _uSol     = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
    _alpha    = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );
    _hEntropy = std::vector<double>( _setSize, 0.0 );
}

DataGeneratorRegression2D::~DataGeneratorRegression2D() {}

void DataGeneratorRegression2D::ComputeMoments() {
    double my, phi;

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
        my                     = _quadPointsSphere[idx_quad][0];
        phi                    = _quadPointsSphere[idx_quad][1];
        _momentBasis[idx_quad] = _basisGenerator->ComputeSphericalBasis( my, phi );
    }

}

void DataGeneratorRegression2D::SampleSolutionU() {
    // Use necessary conditions from Monreal, Dissertation, Chapter 3.2.1, Page 26

    // different processes for different
    if( _maxPolyDegree == 0 ) {
        // --- sample u in order 0 ---
        // u_0 = <1*psi>

        // --- Determine stepsizes etc ---
        double du0 = _settings->GetMaxValFirstMoment() / (double)_gridSize;

        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
            _uSol[idx_set][0] = du0 * idx_set;
        }
    }
    else if( _settings->GetNormalizedSampling() ) {
        if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            // Sample points on unit sphere.
            // Use MC sampling: x randUniform ==> r = sqrt(x)*u_0Max
            //                  y radnUniform ==> a = 2*pi*y,

            // Create generator
            std::default_random_engine generator;
            std::uniform_real_distribution<double> distribution( 0.0, 1.0 );

#pragma omp parallel for schedule( guided )
            for( unsigned long idx_set = 0; idx_set < _setSize; idx_set++ ) {
                double mu  = std::sqrt( distribution( generator ) );
                double phi = 2 * M_PI * distribution( generator );

                if( std::sqrt( 1 - mu * mu ) > 1 - _settings->GetRealizableSetEpsilonU0() ) {
                    idx_set--;
                    continue;
                }
                else {
                    _uSol[idx_set][0] = 1.0;
                    _uSol[idx_set][1] = std::sqrt( 1 - mu * mu ) * std::cos( phi );
                    _uSol[idx_set][2] = std::sqrt( 1 - mu * mu ) * std::sin( phi );
                }
            }
        }
        else {
            ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
        }
    }
    else if( !_settings->GetNormalizedSampling() ) {
        if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            // Sample points on unit sphere.
            // Use MC sampling: x randUniform ==> r = sqrt(x)*u_0Max
            //                  y radnUniform ==> a = 2*pi*y,

            // Create generator
            std::default_random_engine generator;
            std::uniform_real_distribution<double> distribution( 0.0, 1.0 );

            double du     = 0.001;
            long gridSize = (long)( (double)_settings->GetMaxValFirstMoment() / du );
            long c        = 0;
            for( long idx_set = 0; idx_set < gridSize; idx_set++ ) {

                double radiusU0 = du * ( idx_set + 1 );
                // Boundary correction
                if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {
                    radiusU0 = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
                }
                // number of points for unit circle should scale with (u_0)^2,
                long currSetSize = _gridSize;    //(long)( ( radiusU0 * radiusU0 ) * (double)_gridSize );

                for( long idx_currSet = 0; idx_currSet < currSetSize; idx_currSet++ ) {
                    double mu  = std::sqrt( distribution( generator ) );
                    double phi = 2 * M_PI * distribution( generator );
                    if( std::sqrt( 1 - mu * mu ) < _settings->GetRealizableSetEpsilonU1() &&
                        std::sqrt( 1 - mu * mu ) > radiusU0 - _settings->GetRealizableSetEpsilonU1() ) {
                        ErrorMessages::Error( "sampling set is empty. Boundaries overlap", CURRENT_FUNCTION );    // empty set
                    }
                    if( std::sqrt( 1 - mu * mu ) > 1 - _settings->GetRealizableSetEpsilonU1() * radiusU0 ) {
                        idx_currSet--;
                        continue;
                    }
                    else if( std::sqrt( 1 - mu * mu ) < _settings->GetRealizableSetEpsilonU1() * radiusU0 ) {
                        idx_currSet--;
                        continue;
                    }
                    else {
                        _uSol[c][0] = radiusU0;
                        _uSol[c][1] = radiusU0 * std::sqrt( 1 - mu * mu ) * std::cos( phi );
                        _uSol[c][2] = radiusU0 * std::sqrt( 1 - mu * mu ) * std::sin( phi );
                        c++;
                    }
                }
            }
        }
        else {
            ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
        }
    }
}

void DataGeneratorRegression2D::CheckRealizability() {
    double epsilon = _settings->GetRealizableSetEpsilonU0();
    if( _maxPolyDegree == 1 ) {
        //#pragma omp parallel for schedule( guided )
        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
            double normU1 = 0.0;
            Vector u1( 3, 0.0 );
            if( _uSol[idx_set][0] < epsilon ) {
                if( std::abs( _uSol[idx_set][1] ) > 0 || std::abs( _uSol[idx_set][2] ) > 0 ) {
                    ErrorMessages::Error( "Moment not realizable [code 0]. Values: (" + std::to_string( _uSol[idx_set][0] ) + "|" +
                                              std::to_string( _uSol[idx_set][1] ) + "|" + std::to_string( _uSol[idx_set][2] ) + ")",
                                          CURRENT_FUNCTION );
                }
            }
            else {
                u1     = { _uSol[idx_set][1], _uSol[idx_set][2] };
                normU1 = norm( u1 );
                if( normU1 / _uSol[idx_set][0] > _settings->GetRealizableSetEpsilonU1() ) {    // Danger Hardcoded
                    std::cout << "normU1 / _uSol[" << idx_set << "][0]: " << normU1 / _uSol[idx_set][0] << "\n";
                    std::cout << "normU1: " << normU1 << " | _uSol[idx_set][0] " << _uSol[idx_set][0] << "\n";
                    ErrorMessages::Error( "Moment to close to boundary of realizable set [code 1].\nBoundary ratio: " +
                                              std::to_string( normU1 / _uSol[idx_set][0] ),
                                          CURRENT_FUNCTION );
                }
                if( normU1 / _uSol[idx_set][0] <= 0 /*+ 0.5 * epsilon*/ ) {
                    ErrorMessages::Error( "Moment to close to boundary of realizable set [code 2].\nBoundary ratio: " +
                                              std::to_string( normU1 / _uSol[idx_set][0] ),
                                          CURRENT_FUNCTION );
                }
            }
        }
    }
}

void DataGeneratorRegression2D::ComputeSetSizeU() {
    if( _maxPolyDegree == 0 ) {
    }
    else if( _settings->GetNormalizedSampling() ) {
        if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            // Do nothing. Setsize is already given.
        }
        else {
            ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
        }
    }
    else if( !_settings->GetNormalizedSampling() ) {
        if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {

            double du     = 0.001;
            long gridSize = (long)( (double)_settings->GetMaxValFirstMoment() / du );    // Hardcoded
            long c        = 0;
            for( long idx_set = 0; idx_set < gridSize; idx_set++ ) {

                double radiusU0 = du * ( idx_set + 1 );
                // Boundary correction
                if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {
                    radiusU0 = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
                }
                // number of points for unit circle should scale with (u_0)^2,
                long currSetSize = _gridSize;    // (long)( ( radiusU0 * radiusU0 ) * (double)_gridSize );

                for( long idx_currSet = 0; idx_currSet < currSetSize; idx_currSet++ ) {
                    c++;
                }
            }
            _setSize = c;
        }
        else {
            ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
        }
    }
}

1	/*!
2	* \file datagenerator3D.cpp
3	* \brief Class to generate data for the neural entropy closure
4	* \author S. Schotthoefer
5	*/
6
7	#include "datagenerator/datageneratorregression2D.hpp"
8	#include "common/config.hpp"
9	#include "quadratures/quadraturebase.hpp"
10	#include "toolboxes/errormessages.hpp"
11	#include "velocitybasis/sphericalbase.hpp"
12
13	#include <iostream>
14	#include <omp.h>
15
16	DataGeneratorRegression2D::DataGeneratorRegression2D( Config* settings ) : DataGeneratorRegression( settings ) {	×
17	ComputeMoments();	×
18
19	// Initialize Training Data
20	if( _settings->GetAlphaSampling() )	×
21	ComputeSetSizeAlpha();	×
22	else
23	ComputeSetSizeU();	×
24
25	_uSol = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );	×
26	_alpha = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) );	×
27	_hEntropy = std::vector<double>( _setSize, 0.0 );	×
28	}
29
30	DataGeneratorRegression2D::~DataGeneratorRegression2D() {}	×
31		×
32	void DataGeneratorRegression2D::ComputeMoments() {	×
33	double my, phi;
34		×
35	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
36	my = _quadPointsSphere[idx_quad][0];
37	phi = _quadPointsSphere[idx_quad][1];	×
38	_momentBasis[idx_quad] = _basisGenerator->ComputeSphericalBasis( my, phi );	×
39	}	×
40		×
41	}
42
43	void DataGeneratorRegression2D::SampleSolutionU() {
44	// Use necessary conditions from Monreal, Dissertation, Chapter 3.2.1, Page 26
45		×
46	// different processes for different
47	if( _maxPolyDegree == 0 ) {
48	// --- sample u in order 0 ---
49	// u_0 = <1*psi>	×
50
51	// --- Determine stepsizes etc ---
52	double du0 = _settings->GetMaxValFirstMoment() / (double)_gridSize;
53
54	for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {	×
55	_uSol[idx_set][0] = du0 * idx_set;
56	}	×
57	}	×
58	else if( _settings->GetNormalizedSampling() ) {
59	if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
60	// Sample points on unit sphere.	×
61	// Use MC sampling: x randUniform ==> r = sqrt(x)*u_0Max	×
62	// y radnUniform ==> a = 2piy,
63
64	// Create generator
65	std::default_random_engine generator;
66	std::uniform_real_distribution<double> distribution( 0.0, 1.0 );
67		×
68	#pragma omp parallel for schedule( guided )	×
69	for( unsigned long idx_set = 0; idx_set < _setSize; idx_set++ ) {
70	double mu = std::sqrt( distribution( generator ) );	×
71	double phi = 2 * M_PI * distribution( generator );
72
73	if( std::sqrt( 1 - mu * mu ) > 1 - _settings->GetRealizableSetEpsilonU0() ) {
74	idx_set--;
75	continue;
76	}
77	else {
78	_uSol[idx_set][0] = 1.0;
79	_uSol[idx_set][1] = std::sqrt( 1 - mu * mu ) * std::cos( phi );
80	_uSol[idx_set][2] = std::sqrt( 1 - mu * mu ) * std::sin( phi );
81	}
82	}
83	}
84	else {
85	ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
86	}
87	}	×
88	else if( !_settings->GetNormalizedSampling() ) {
89	if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
90	// Sample points on unit sphere.	×
91	// Use MC sampling: x randUniform ==> r = sqrt(x)*u_0Max	×
92	// y radnUniform ==> a = 2piy,
93
94	// Create generator
95	std::default_random_engine generator;
96	std::uniform_real_distribution<double> distribution( 0.0, 1.0 );
97		×
98	double du = 0.001;	×
99	long gridSize = (long)( (double)_settings->GetMaxValFirstMoment() / du );
100	long c = 0;	×
101	for( long idx_set = 0; idx_set < gridSize; idx_set++ ) {	×
102		×
103	double radiusU0 = du * ( idx_set + 1 );	×
104	// Boundary correction
105	if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {	×
106	radiusU0 = _settings->GetRealizableSetEpsilonU0(); // Boundary close to 0
107	}	×
108	// number of points for unit circle should scale with (u_0)^2,	×
109	long currSetSize = _gridSize; //(long)( ( radiusU0 * radiusU0 ) * (double)_gridSize );
110
111	for( long idx_currSet = 0; idx_currSet < currSetSize; idx_currSet++ ) {	×
112	double mu = std::sqrt( distribution( generator ) );
113	double phi = 2 * M_PI * distribution( generator );	×
114	if( std::sqrt( 1 - mu * mu ) < _settings->GetRealizableSetEpsilonU1() &&	×
115	std::sqrt( 1 - mu * mu ) > radiusU0 - _settings->GetRealizableSetEpsilonU1() ) {	×
116	ErrorMessages::Error( "sampling set is empty. Boundaries overlap", CURRENT_FUNCTION ); // empty set	×
117	}	×
118	if( std::sqrt( 1 - mu * mu ) > 1 - _settings->GetRealizableSetEpsilonU1() * radiusU0 ) {	×
119	idx_currSet--;
120	continue;	×
121	}	×
122	else if( std::sqrt( 1 - mu * mu ) < _settings->GetRealizableSetEpsilonU1() * radiusU0 ) {	×
123	idx_currSet--;
124	continue;	×
125	}	×
126	else {	×
127	_uSol[c][0] = radiusU0;
128	_uSol[c][1] = radiusU0 * std::sqrt( 1 - mu * mu ) * std::cos( phi );
129	_uSol[c][2] = radiusU0 * std::sqrt( 1 - mu * mu ) * std::sin( phi );	×
130	c++;	×
131	}	×
132	}	×
133	}
134	}
135	else {
136	ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
137	}
138	}	×
139	}
140
141	void DataGeneratorRegression2D::CheckRealizability() {
142	double epsilon = _settings->GetRealizableSetEpsilonU0();
143	if( _maxPolyDegree == 1 ) {	×
144	//#pragma omp parallel for schedule( guided )	×
145	for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {	×
146	double normU1 = 0.0;
147	Vector u1( 3, 0.0 );	×
148	if( _uSol[idx_set][0] < epsilon ) {	×
149	if( std::abs( _uSol[idx_set][1] ) > 0 \|\| std::abs( _uSol[idx_set][2] ) > 0 ) {	×
150	ErrorMessages::Error( "Moment not realizable [code 0]. Values: (" + std::to_string( _uSol[idx_set][0] ) + "\|" +	×
151	std::to_string( _uSol[idx_set][1] ) + "\|" + std::to_string( _uSol[idx_set][2] ) + ")",	×
152	CURRENT_FUNCTION );	×
153	}	×
154	}
155	else {
156	u1 = { _uSol[idx_set][1], _uSol[idx_set][2] };
157	normU1 = norm( u1 );
158	if( normU1 / _uSol[idx_set][0] > _settings->GetRealizableSetEpsilonU1() ) { // Danger Hardcoded	×
159	std::cout << "normU1 / _uSol[" << idx_set << "][0]: " << normU1 / _uSol[idx_set][0] << "\n";	×
160	std::cout << "normU1: " << normU1 << " \| _uSol[idx_set][0] " << _uSol[idx_set][0] << "\n";	×
161	ErrorMessages::Error( "Moment to close to boundary of realizable set [code 1].\nBoundary ratio: " +	×
162	std::to_string( normU1 / _uSol[idx_set][0] ),	×
163	CURRENT_FUNCTION );	×
164	}	×
165	if( normU1 / _uSol[idx_set][0] <= 0 /+ 0.5 epsilon*/ ) {
166	ErrorMessages::Error( "Moment to close to boundary of realizable set [code 2].\nBoundary ratio: " +
167	std::to_string( normU1 / _uSol[idx_set][0] ),	×
168	CURRENT_FUNCTION );	×
169	}	×
170	}
171	}
172	}
173	}
174
175	void DataGeneratorRegression2D::ComputeSetSizeU() {
176	if( _maxPolyDegree == 0 ) {
177	}	×
178	else if( _settings->GetNormalizedSampling() ) {	×
179	if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
180	// Do nothing. Setsize is already given.	×
181	}	×
182	else {
183	ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
184	}
185	}	×
186	else if( !_settings->GetNormalizedSampling() ) {
187	if( _maxPolyDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
188		×
189	double du = 0.001;	×
190	long gridSize = (long)( (double)_settings->GetMaxValFirstMoment() / du ); // Hardcoded
191	long c = 0;	×
192	for( long idx_set = 0; idx_set < gridSize; idx_set++ ) {	×
193		×
194	double radiusU0 = du * ( idx_set + 1 );	×
195	// Boundary correction
196	if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {	×
197	radiusU0 = _settings->GetRealizableSetEpsilonU0(); // Boundary close to 0
198	}	×
199	// number of points for unit circle should scale with (u_0)^2,	×
200	long currSetSize = _gridSize; // (long)( ( radiusU0 * radiusU0 ) * (double)_gridSize );
201
202	for( long idx_currSet = 0; idx_currSet < currSetSize; idx_currSet++ ) {	×
203	c++;
204	}	×
205	}	×
206	_setSize = c;
207	}
208	else {	×
209	ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
210	}
211	}	×
212	}

CSMMLab / KiT-RT / #100

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous