datagenerator.cpp 13.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
/*!
 * \file datagenerator.cpp
 * \brief Class to generate data for the neural entropy closure
 * \author S. Schotthoefer
 */

#include "toolboxes/datagenerator.h"
#include "common/config.h"
#include "entropies/entropybase.h"
#include "optimizers/newtonoptimizer.h"
#include "quadratures/qlebedev.h"
#include "quadratures/quadraturebase.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalbase.h"

#include "spdlog/spdlog.h"

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

nnDataGenerator::nnDataGenerator( Config* settings ) {
    _settings = settings;
    _setSize  = settings->GetTrainingDataSetSize();
    _gridSize = _setSize;

    _LMaxDegree = settings->GetMaxMomentDegree();

    // Quadrature
    _quadrature       = QuadratureBase::Create( settings );
    _nq               = _quadrature->GetNq();
    _quadPoints       = _quadrature->GetPoints();
    _weights          = _quadrature->GetWeights();
    _quadPointsSphere = _quadrature->GetPointsSphere();

    // Spherical Harmonics
    if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS && _LMaxDegree > 0 ) {
        ErrorMessages::Error( "No sampling algorithm for spherical harmonics basis with degree higher than 0 implemented", CURRENT_FUNCTION );
    }
    _basis = SphericalBase::Create( _settings );

    _nTotalEntries = _basis->GetBasisSize();

    _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );

    ComputeMoments();

    // Optimizer
    _optimizer = new NewtonOptimizer( _settings );

    // Entropy
    _entropy = EntropyBase::Create( _settings );

    // Initialize Training Data
    if( _LMaxDegree == 0 ) {
    }
    else if( _LMaxDegree == 1 ) {
        // Sample points on unit sphere.
        QuadratureBase* quad = QuadratureBase::Create( _settings );
        unsigned long nq     = (unsigned long)quad->GetNq();

        // Allocate memory.
        _setSize = nq * _gridSize * ( _gridSize - 1 ) / 2;

        delete quad;
    }
    else if( _LMaxDegree == 2 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS && _settings->GetDim() == 1 ) {
        // Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1
        unsigned c = 1;
69
        double N1  = -1.0 + _settings->GetRealizableSetEpsilonU0();
70
71
        double N2;
        double dN = 2.0 / (double)_gridSize;
72
73
74
        while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
            N2 = N1 * N1 + _settings->GetRealizableSetEpsilonU0();
            while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
                c++;
                N2 += dN;
            }
            N1 += dN;
        }
        _setSize = c;
    }
    else {
        ErrorMessages::Error( "Sampling of training data of degree higher than 1 is not yet implemented.", CURRENT_FUNCTION );
    }

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

nnDataGenerator::~nnDataGenerator() {
    delete _quadrature;
    delete _entropy;
}

void nnDataGenerator::computeTrainingData() {
    // Prototype: Only for _LMaxDegree == 1
    // Prototype: u is sampled from [0,100]

    // --- sample u ---
    SampleSolutionU();

    PrintLoadScreen();

    // ---- Check realizability ---
    CheckRealizability();

    // --- compute alphas ---
    _optimizer->SolveMultiCell( _alpha, _uSol, _moments );

    // --- Postprocessing
    ComputeRealizableSolution();

    // --- compute entropy functional ---
    ComputeEntropyH_primal();

    // --- Print everything ----
    PrintTrainingData();
}

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

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
        my  = _quadPointsSphere[idx_quad][0];
        phi = _quadPointsSphere[idx_quad][1];

        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
    }
}

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

    // --- Determine stepsizes etc ---
136
    double du0 = _settings->GetMaxValFirstMoment() / (double)_gridSize;
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

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

        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
            _uSol[idx_set][0] = du0 * idx_set;
        }
    }
    else if( _LMaxDegree == 1 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
        // Sample points on unit sphere.
        QuadratureBase* quad = QuadratureBase::Create( _settings );
        VectorVector qpoints = quad->GetPoints();    // carthesian coordinates.
        unsigned long nq     = (unsigned long)quad->GetNq();

        // --- sample u in order 0 ---
        // u_0 = <1*psi>

156
        //#pragma omp parallel for schedule( guided )
157
        for( unsigned long idx_set = 0; idx_set < _gridSize; idx_set++ ) {
158
            Vector u1              = { 0, 0, 0 };
159
160
161
162
163
164
            unsigned long outerIdx = ( idx_set - 1 ) * ( idx_set ) / 2;    // sum over all radii up to current
            outerIdx *= nq;                                                // in each radius step, use all quad points

            //  --- sample u in order 1 ---
            /* order 1 has 3 elements. (omega_x, omega_y,  omega_z) = omega, let u_1 = (u_x, u_y, u_z) = <omega*psi>
             * Condition u_0 >= norm(u_1)   */
165
166
167
168
169
            double radiusU0 = du0 * ( idx_set + 1 );
            // Boundary correction
            if( radiusU0 < _settings->GetRealizableSetEpsilonU0() ) {
                radiusU0 = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
            }
170
171
172
173
174
175

            unsigned long localIdx = 0;
            double radius          = 0.0;
            unsigned long innerIdx = 0;
            // loop over all radii
            for( unsigned long idx_subset = 0; idx_subset < idx_set; idx_subset++ ) {
176
177
178
179
180
181
182
183
184
185
                radius = du0 * ( idx_subset + 1 );    // dont use radius 0 ==> shift by one
                // Boundary correction
                if( radius < _settings->GetRealizableSetEpsilonU0() ) {
                    radius = _settings->GetRealizableSetEpsilonU0();    // Boundary close to 0
                }
                if( radius / radiusU0 > 0.95 * _settings->GetRealizableSetEpsilonU1() ) {
                    // small offset to take care of rounding errors in quadrature
                    radius = radiusU0 * 0.95 * _settings->GetRealizableSetEpsilonU1();    // "outer" boundary
                }

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
                localIdx = outerIdx + idx_subset * nq;

                for( unsigned long quad_idx = 0; quad_idx < nq; quad_idx++ ) {
                    innerIdx = localIdx + quad_idx;    // gives the global index

                    _uSol[innerIdx][0] = radiusU0;    // Prevent 1st order moments to be too close to the boundary
                    // scale quadpoints with radius
                    _uSol[innerIdx][1] = radius * qpoints[quad_idx][0];
                    _uSol[innerIdx][2] = radius * qpoints[quad_idx][1];
                    _uSol[innerIdx][3] = radius * qpoints[quad_idx][2];
                }
            }
        }
    }
    else if( _LMaxDegree == 2 && _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS && _settings->GetDim() == 1 ) {
        // Carefull: This computes only normalized moments, i.e. sampling for u_0 = 1
        unsigned c = 0;
203
        double N1  = -1.0 + _settings->GetRealizableSetEpsilonU0();
204
205
        double N2;
        double dN = 2.0 / (double)_gridSize;
206
207
208
        while( N1 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
            N2 = N1 * N1 + _settings->GetRealizableSetEpsilonU0();
            while( N2 < 1.0 - _settings->GetRealizableSetEpsilonU0() ) {
209
210
211
212
213
214
215
216
217
218
                _uSol[c][0] = 1;     // u0 (normalized i.e. N0) by Monreals notation
                _uSol[c][1] = N1;    // u1 (normalized i.e. N1) by Monreals notation
                _uSol[c][2] = N2;    // u2 (normalized i.e. N2) by Monreals notation
                N2 += dN;
                c++;
            }
            N1 += dN;
        }
    }
    else {
219
        ErrorMessages::Error( "Sampling for this configuration is not yet supported", CURRENT_FUNCTION );
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
    }
}

void nnDataGenerator::ComputeEntropyH_dual() {
    for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
        _hEntropy[idx_set] = _optimizer->ComputeObjFunc( _alpha[idx_set], _uSol[idx_set], _moments );
    }
}

void nnDataGenerator::ComputeEntropyH_primal() {
    double result = 0.0;

    for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
        result = 0.0;
        // Integrate (eta(eta'_*(alpha*m))
        for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
            result += _entropy->Entropy( _entropy->EntropyPrimeDual( dot( _alpha[idx_set], _moments[idx_quad] ) ) ) * _weights[idx_quad];
        }
        _hEntropy[idx_set] = result;
    }
}

void nnDataGenerator::PrintTrainingData() {
    auto log    = spdlog::get( "event" );
    auto logCSV = spdlog::get( "tabular" );
    log->info( "---------------------- Data Generation Successful ------------------------" );

    std::string uSolString  = "";
    std::string alphaString = "";
    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
        uSolString += "u_" + std::to_string( idx_sys ) + ",";
        alphaString += "alpha_" + std::to_string( idx_sys ) + ",";
    }
    // log->info( uSolString + alphaString + "h" );
    logCSV->info( uSolString + alphaString + "h" );

    for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
        std::string uSolString  = "";
        std::string alphaString = "";
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
            uSolString += std::to_string( _uSol[idx_set][idx_sys] ) + ",";
            alphaString += std::to_string( _alpha[idx_set][idx_sys] ) + ",";
        }
        // log->info( uSolString + alphaString + "{}", _hEntropy[idx_set] );
        logCSV->info( uSolString + alphaString + "{}", _hEntropy[idx_set] );
    }
}

void nnDataGenerator::CheckRealizability() {
269
    double epsilon = _settings->GetRealizableSetEpsilonU0();
270
271
272
273
274
275
    if( _LMaxDegree == 1 ) {
        double normU1 = 0.0;
        Vector u1( 3, 0.0 );
        for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) {
            if( _uSol[idx_set][0] < epsilon ) {
                if( std::abs( _uSol[idx_set][1] ) > 0 || std::abs( _uSol[idx_set][2] ) > 0 || std::abs( _uSol[idx_set][3] ) > 0 ) {
276
277
278
279
                    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] ) + "|" +
                                              std::to_string( _uSol[idx_set][3] ) + ")",
                                          CURRENT_FUNCTION );
280
281
282
283
284
                }
            }
            else {
                u1     = { _uSol[idx_set][1], _uSol[idx_set][2], _uSol[idx_set][3] };
                normU1 = norm( u1 );
285
286
287
                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";
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
                    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*/ ) {
                    // std::cout << "_uSol" << _uSol[idx_set][1] << " | " << _uSol[idx_set][2] << " | " << _uSol[idx_set][3] << " \n";
                    // 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 2].\nBoundary ratio: " +
                                              std::to_string( normU1 / _uSol[idx_set][0] ),
                                          CURRENT_FUNCTION );
                }
            }
        }
    }
}

void nnDataGenerator::ComputeRealizableSolution() {
    for( unsigned idx_sol = 0; idx_sol < _setSize; idx_sol++ ) {
        double entropyReconstruction = 0.0;
        _uSol[idx_sol]               = 0;
        for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
            // Make entropyReconstruction a member vector, s.t. it does not have to be re-evaluated in ConstructFlux
            entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_sol], _moments[idx_quad] ) );
            _uSol[idx_sol] += _moments[idx_quad] * ( _weights[idx_quad] * entropyReconstruction );
        }
    }
}
void nnDataGenerator::PrintLoadScreen() {
    auto log = spdlog::get( "event" );
    log->info( "------------------------ Data Generation Starts --------------------------" );
    log->info( "| Generating {} datapoints.", _setSize );
}