1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 #ifndef EXPERIMENT_DATA_UTILS
10 #define EXPERIMENT_DATA_UTILS
11
12 #include "dakota_data_types.hpp"
13 #include "Teuchos_SerialSpdDenseSolver.hpp"
14
15 namespace Dakota {
16
17 class Response;
18
19
20 /**
21 * \brief find the interval containing a target value.
22 * This function assumes the data is in ascending order.
23 */
24 template<typename O, typename T>
binary_search(T target,Teuchos::SerialDenseVector<O,T> & data)25 int binary_search( T target, Teuchos::SerialDenseVector<O,T> &data )
26 {
27 O low = 0, high = data.length()-1, mid;
28 while ( low <= high )
29 {
30 mid = low + ( high - low ) / 2;
31 if ( target < data[mid] ) high = mid - 1;
32 else if ( target > data[mid] ) low = mid + 1;
33 else return mid;
34 }
35 if ( high < 0 ) return 0;
36 else if ( high < low ) return high;
37 else return low;
38
39 }
40
41 void copy_field_data(const RealVector& fn_vals, RealMatrix& fn_grad,
42 const RealSymMatrixArray& fn_hess, size_t offset,
43 size_t num_fns, Response& response);
44
45 void interpolate_simulation_field_data( const Response &sim_resp,
46 const RealMatrix &exp_coords,
47 size_t field_num, short total_asv,
48 size_t interp_resp_offset,
49 Response &interp_resp );
50
51 /**
52 * \brief Returns the value of at 1D function f and its gradient and hessians
53 * (if available) at the points of vector pred_pts using linear interpolation.
54 * The vector build_pts specifies the coordinates of the underlying interval at
55 * which the values (build_vals) of the function f are known. The length of
56 * output pred_vals is equal to the length of pred_pts.
57 * This function assumes the build_pts is in ascending order */
58 void linear_interpolate_1d( const RealMatrix &build_pts,
59 const RealVector &build_vals,
60 const RealMatrix &build_grads,
61 const RealSymMatrixArray &build_hessians,
62 const RealMatrix &pred_pts,
63 RealVector &pred_vals,
64 RealMatrix &pred_grads,
65 RealSymMatrixArray &pred_hessians );
66
67 class CovarianceMatrix
68 {
69 private:
70
71 /// The number of rows in the covariance matrix
72 int numDOF_;
73
74 /// The covariance matrix where each element (i,j) is the covariance
75 /// between points Xi and Xj in the initial set of samples
76 RealSymMatrix covMatrix_;
77
78 /// The diagonal entries of a diagonal covariance matrix.
79 RealVector covDiagonal_;
80
81 /// The inverse of the covariance matrix
82 RealSymMatrix covCholFactor_;
83
84 /// The inverse of the Cholesky factor of the covariance matrix
85 RealMatrix cholFactorInv_;
86
87 /// Flag specifying if the covariance matrix is diagonal
88 bool covIsDiagonal_;
89
90 /// The global solver for all computations involving the inverse of
91 /// the covariance matrix
92 Teuchos::SerialSpdDenseSolver<int, Real> covSlvr_;
93
94 /// Compute the Cholesky factorization of the covariance matrix
95 void factor_covariance_matrix();
96
97 /// Compute the inverse of the Cholesky factor of the covariance matrix
98 void invert_cholesky_factor();
99
100 /// Copy the values from one existing CovarianceMatrix to another.
101 void copy( const CovarianceMatrix &source );
102
103 public:
104
105 /// Specification for covariance type
106 enum FORMAT { CONSTANT,
107 VECTOR,
108 MATRIX };
109
110 /// Default Constructor
111 CovarianceMatrix();
112
113 /// Copy Constructor. Needed for creating std::vector<CovarianceMatrix>
114 CovarianceMatrix( const CovarianceMatrix &cov_mat );
115
116 /// Deconstructor
117 ~CovarianceMatrix();
118
119 /// The operator= copies the values from one existing CovarianceMatrix to
120 /// another.
121 CovarianceMatrix& operator=( const CovarianceMatrix &source );
122
123 /// Return the full dense matrix specifying the covariance among the
124 /// field data
125 void dense_covariance(RealSymMatrix &cov) const;
126
127 /// Retrieve the covariance as a correlation matrix
128 void as_correlation(RealSymMatrix& corr_mat) const;
129
130 /// Set the matrix specifying the correlation between the field data
131 void set_covariance( const RealMatrix &cov );
132
133 /// Set the matrix specifying the correlation between the field data.
134 // A scalar covariance is a degenerate case of a diagonal covariance matrix
135 void set_covariance( Real cov );
136
137 /// Set the diagonal of the matrix specifying the correlation between
138 // the field data
139 void set_covariance( const RealVector & cov );
140
141 /// Compute the triple product of r'*inv(C)*r where r is a vector r and C
142 /// is the covariance matrix
143 Real apply_covariance_inverse( const RealVector &vector ) const;
144
145 /// Multiply a vector r by the sqrt of the inverse covariance matrix C, i.e.
146 /// compute L'*r where L is the cholesky factor of the positive definite
147 /// covariance matrix
148 void apply_covariance_inverse_sqrt( const RealVector &vector,
149 RealVector &result ) const;
150
151 /// Multiply a matrix of gradients g (each column is a gradient vector)
152 /// by the sqrt of the inverse covariance matrix C, i.e.
153 /// compute L'*g where L is the cholesky factor of the positive definite
154 /// covariance matrix
155 void apply_covariance_inverse_sqrt_to_gradients( const RealMatrix &gradients,
156 RealMatrix &result ) const;
157
158 /// Apply the sqrt of the inverse covariance matrix to a list of Hessians.
159 /// the argument hessians is a numDOF_ list of num_grads x num_grads Hessian
160 /// matrices. \param start is an index pointing to the first Hessian to
161 /// consider in the list. This helps avoid copying large Hessian matrix
162 /// when applying block covariances using ExperimentCovariance class
163 void apply_covariance_inverse_sqrt_to_hessian( RealSymMatrixArray &hessians,
164 int start) const;
165
166 /// Return the number of rows in the covariance matrix
167 int num_dof() const;
168
169 /// Print a covariance matrix
170 void print() const;
171
172 /// Return a (copy) vector containing the entries of the main diagonal of the
173 /// covariance matrix
174 void get_main_diagonal( RealVector &diagonal ) const;
175
176 /// The determinant of the covariance matrix
177 Real determinant() const;
178
179 /// The log(determinant) of the covariance matrix
180 Real log_determinant() const;
181
182 };
183
184
185 class ExperimentCovariance
186 {
187 private:
188 /// A list of block covariance matrices. ExperimentCovariance consists of
189 /// multiple block covariance matrices. These are stored in a list
190 /// rather than a dense matrix to save memory
191 std::vector<CovarianceMatrix> covMatrices_;
192
193 /// The number of block covariance matrices
194 int numBlocks_;
195
196 /// The number of entries along the diagonal;
197 int numDOF_;
198
199 public:
200
201 /// Default Constructor
ExperimentCovariance()202 ExperimentCovariance() : numBlocks_(0), numDOF_(0) {};
203
204 /// Assignment operator
205 ExperimentCovariance & operator=(const ExperimentCovariance& source);
206
207 /// Deconstructor
~ExperimentCovariance()208 ~ExperimentCovariance(){numBlocks_=0;numDOF_=0;};
209
210 /// Set the experiment covariance matrix blocks
211 void set_covariance_matrices( std::vector<RealMatrix> &matrices,
212 std::vector<RealVector> &diagonals,
213 RealVector &scalars,
214 IntVector matrix_map_indices,
215 IntVector diagonal_map_indices,
216 IntVector scalar_map_indices );
217
218 /// Compute the triple product v'*inv(C)*v
219 Real apply_experiment_covariance( const RealVector &vector ) const;
220
221 /// Compute the product inv(L)*v where L is the Cholesky factor of the
222 /// covariance matrix C
223 void apply_experiment_covariance_inverse_sqrt( const RealVector &vector,
224 RealVector &result ) const;
225
226 /// Compute the product inv(L)*G where L is the Cholesky factor of the
227 /// covariance matrix C and G is a matrix whose columns are gradient vectors
228 /// for each degree of freedom
229 void apply_experiment_covariance_inverse_sqrt_to_gradients(
230 const RealMatrix &grads,
231 RealMatrix &result) const;
232
233 /// Compute the products inv(L)*H where L is the Cholesky factor of the
234 /// covariance matrix C and H is a Hessian matrix. The product is computed
235 /// for each Hessian of every degree of freedom.
236 void apply_experiment_covariance_inverse_sqrt_to_hessians(
237 const RealSymMatrixArray &hesians, RealSymMatrixArray &result ) const;
238
239 /// Print each block of the covariance matrix
240 void print_covariance_blocks() const;
241
242 /// Return a (copy) vector containing the main diagonal entries of the
243 /// experimental covariance matrix
244 void get_main_diagonal( RealVector &diagonal ) const;
245
246 /// Retrieve a full/dense representation of the covariance
247 void dense_covariance(RealSymMatrix& cov_mat) const;
248
249 /// Retrieve the covariance as a (dense) correlation matrix
250 void as_correlation(RealSymMatrix& corr_mat) const;
251
252 /// Determinant of the total ExperimentCovariance structure
253 Real determinant() const;
254
255 /// Log of determinant of the total ExperimentCovariance structure
256 Real log_determinant() const;
257
258 /// Return the number of diagonal blocks in the entire matrix
num_blocks() const259 int num_blocks() const{
260 return numBlocks_;
261 };
262
num_dof() const263 int num_dof() const{
264 return numDOF_;
265 }
266
267 };
268
269 /**
270 * \brief Computes the eigenvalues and, optionally, eigenvectors of a
271 * real symmetric matrix A.
272 *
273 * Eigenvalues are returned in ascending order.
274 */
275 void symmetric_eigenvalue_decomposition( const RealSymMatrix &matrix,
276 RealVector &eigenvalues,
277 RealMatrix &eigenvectors );
278
279 /**
280 * \brief Compute the means of each column of an arbitrary matrix
281 */
282 void compute_column_means( RealMatrix & matrix, RealVector & avg_vals );
283
284 /**
285 * \brief Sort incoming vector with result and corresponding indices returned in passed arguments
286 */
287 void sort_vector( const RealVector & vec, RealVector & sort_vec, IntVector & indices );
288
289 /**
290 * \brief Sort incoming matrix columns with result and corresponding indices returned in passed arguments
291 */
292 void sort_matrix_columns( const RealMatrix & mat, RealMatrix & sort_mat, IntMatrix & indices );
293
294 /**
295 * \brief Test if incoming matrix is symmetric
296 */
297 bool is_matrix_symmetric( const RealMatrix & matrix );
298
299 } // namespace dakota
300
301 #endif //EXPERIMENT_DATA_UTILS
302