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
69
70
71
72
73
74
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
// Copyright (C) 2016 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef LOCALSAMESIZELINEARREGRESSION_H
#define LOCALSAMESIZELINEARREGRESSION_H
#include <vector>
#include <memory>
#include <array>
#include <Eigen/Dense>
#include "StOpt/regression/LocalSameSizeRegression.h"
#include "StOpt/core/grids/InterpolatorSpectral.h"
/** \file LocalSameSizeLinearRegression.h
* \brief Compute conditional expectations by using linear local regressions with fixed size meshes.
* \author Xavier Warin
*/
namespace StOpt
{
/**
* \defgroup LocalSameSizeLinear Piecewise linear regressions with constant mesh size
* \brief It implements local linear functions for performing local regressions with meshes with same size
*@{
*/
/// \class LocalSameSizeLinearRegression LocalSameSizeLinearRegression.h
/// Linear regression on each cell of same size
class LocalSameSizeLinearRegression : public LocalSameSizeRegression
{
private:
Eigen::ArrayXXd m_matReg ; ///< Regression matrix (rank \f$ (n_d+1)^2 \f$ by the number of meshes ). Also store lower part of factorized matrix.
Eigen::ArrayXXd m_diagReg ; ///< Diagonal part of the factorized matrix during Choleski factorization
Eigen::Array<bool, Eigen::Dynamic, 1> m_bSingular ; ///< for each cell true if the matrix is singular
std::vector< std::shared_ptr< Eigen::MatrixXd > > m_pseudoInverse ; ///< when the Cholesky matrix is singular use Courrieu 2005 : Fast Computation of Moore-Penrose Inverse Matrices to calculate the pseudoInverse
/// \brief construct local matrices and factorize them
void constructAndFactorize();
/// \brief construct the second member for the regression
/// \param p_fToRegress function to regress
/// \return second member for the linear regression
Eigen::ArrayXd secondMember(const Eigen::ArrayXd &p_fToRegress) const;
/// \brief matrix inversion or pseudo inversion to get function basis
/// \param p_secMem second member associated to regression
Eigen::ArrayXd inversion(const Eigen::ArrayXd &p_secMem) const ;
/// \brief Reconstruction of the once the function basis coefficients are calculated
/// \param p_foncBasisCoef coefficients associated to the basis functions.
Eigen::ArrayXd reconstructionAllPoints(const Eigen::ArrayXd &p_foncBasisCoef) const;
/// \brief Reconstruction for one given
/// \param p_coord Point coordinate
/// \param p_cell Cell where it belong to
/// \param p_basisCoefficients Basis function coefficients
double reconstructionOnlyOnePoint(const Eigen::ArrayXd &p_coord, const int &p_cell, const Eigen::ArrayXd &p_basisCoefficients) const ;
public :
/// \brief Default constructor
LocalSameSizeLinearRegression() {}
/// \brief Constructor
/// \param p_lowValues in each dimension minimal value of the grid
/// \param p_step in each dimension the step size
/// \param p_nbStep in each dimension the number of steps
LocalSameSizeLinearRegression(const Eigen::ArrayXd &p_lowValues, const Eigen::ArrayXd &p_step, const Eigen::ArrayXi &p_nbStep) : LocalSameSizeRegression(p_lowValues, p_step, p_nbStep) {} ;
/// \brief Constructor for object constructed at each time step
/// \param p_bZeroDate first date is 0?
/// \param p_particles particles used for the meshes.
/// First dimension : dimension of the problem,
/// second dimension : the number of particles
/// \param p_lowValues in each dimension minimal value of the grid
/// \param p_step in each dimension the step size
/// \param p_nbStep in each dimension the number of steps
LocalSameSizeLinearRegression(const bool &p_bZeroDate,
const Eigen::ArrayXXd &p_particles,
const Eigen::ArrayXd &p_lowValues,
const Eigen::ArrayXd &p_step,
const Eigen::ArrayXi &p_nbStep);
/// \brief Constructor only used for serialization for simulation part
/// \param p_bZeroDate first date is 0?
/// \param p_lowValues in each dimension minimal value of the grid
/// \param p_step in each dimension the step size
/// \param p_nbStep in each dimension the number of steps
LocalSameSizeLinearRegression(const bool &p_bZeroDate,
const Eigen::ArrayXd &p_lowValues,
const Eigen::ArrayXd &p_step,
const Eigen::ArrayXi &p_nbStep):
LocalSameSizeRegression(p_bZeroDate, p_lowValues, p_step, p_nbStep) {}
/// \brief Copy constructor
/// \param p_objet object to copy
LocalSameSizeLinearRegression(const LocalSameSizeLinearRegression &p_objet);
/// \brief update the particles used in regression and construct the matrices
/// \param p_bZeroDate first date is 0?
/// \param p_particles particles used for the meshes.
/// First dimension : dimension of the problem,
/// second dimension : the number of particles
void updateSimulations(const bool &p_bZeroDate, const Eigen::ArrayXXd &p_particles);
/// \brief Get some local accessors
///@{
inline const Eigen::ArrayXXd &getMatReg() const
{
return m_matReg;
}
inline const Eigen::ArrayXXd &getDiagReg() const
{
return m_diagReg;
}
inline const Eigen::Array<bool, Eigen::Dynamic, 1> &getBSingular() const
{
return m_bSingular;
}
inline const std::vector< std::shared_ptr< Eigen::MatrixXd > > &getPseudoInverse() const
{
return m_pseudoInverse;
}
///@}
/// \brief conditional expectation basis function coefficient calculation
/// \param p_fToRegress function to regress associated to each simulation used in optimization
/// \return regression coordinates on the basis (size : number of meshes multiplied by the dimension plus one)
/// @{
Eigen::ArrayXd getCoordBasisFunction(const Eigen::ArrayXd &p_fToRegress) const;
Eigen::ArrayXXd getCoordBasisFunctionMultiple(const Eigen::ArrayXXd &p_fToRegress) const ;
///@}
/// \brief conditional expectation calculation
/// \param p_fToRegress simulations to regress used in optimization
/// \return regressed value function
/// @{
Eigen::ArrayXd getAllSimulations(const Eigen::ArrayXd &p_fToRegress) const ;
Eigen::ArrayXXd getAllSimulationsMultiple(const Eigen::ArrayXXd &p_fToRegress) const;
///@}
/// \brief Use basis functions to reconstruct the solution
/// \param p_basisCoefficients basis coefficients
///@{
Eigen::ArrayXd reconstruction(const Eigen::ArrayXd &p_basisCoefficients) const;
Eigen::ArrayXXd reconstructionMultiple(const Eigen::ArrayXXd &p_basisCoefficients) const;
/// @}
/// \brief use basis function to reconstruct a given simulation
/// \param p_isim simulation number
/// \param p_basisCoefficients basis coefficients to reconstruct a given conditional expectation
double reconstructionASim(const int &p_isim, const Eigen::ArrayXd &p_basisCoefficients) const ;
/// \brief conditional expectation reconstruction
/// \param p_coordinates coordinates to interpolate (uncertainty sample)
/// \param p_coordBasisFunction regression coordinates on the basis (size: number of meshes multiplied by the dimension plus one)
/// \return regressed value function reconstructed for each simulation
double getValue(const Eigen::ArrayXd &p_coordinates,
const Eigen::ArrayXd &p_coordBasisFunction) const;
/// \brief permits to reconstruct a function with basis functions coefficients values given on a grid
/// \param p_coordinates coordinates (uncertainty sample)
/// \param p_ptOfStock grid point
/// \param p_interpFuncBasis spectral interpolator to interpolate the basis functions coefficients used in regression on the grid (given for each basis function)
double getAValue(const Eigen::ArrayXd &p_coordinates, const Eigen::ArrayXd &p_ptOfStock,
const std::vector< std::shared_ptr<InterpolatorSpectral> > &p_interpFuncBasis) const;
/// \brief get the number of basis functions
inline int getNumberOfFunction() const
{
if (m_bZeroDate)
return 1;
else
return m_nbMeshTotal * (m_lowValues.size() + 1) ;
}
/// \brief Clone the regressor
std::shared_ptr<BaseRegression> clone() const
{
return std::static_pointer_cast<BaseRegression>(std::make_shared<LocalSameSizeLinearRegression>(*this));
}
};
/**@}*/
}
#endif /*LOCALSAMESIZELINEARREGRESSION_H*/
|
