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 DAKOTA_DATA_UTIL_H
10 #define DAKOTA_DATA_UTIL_H
11 #include <exception>
12 #include "dakota_system_defs.hpp"
13 #include "dakota_global_defs.hpp" // for Cerr
14 #include "dakota_data_types.hpp"
15 #include "pecos_data_types.hpp"
16 #include <boost/algorithm/string/case_conv.hpp>
17 #include <boost/algorithm/string/predicate.hpp>
18 #include <boost/functional/hash/hash.hpp>
19 #include <boost/regex.hpp>
20 #include <algorithm>
21 #include "Teuchos_SerialDenseHelpers.hpp"
22
23
24 // --------------
25 // hash functions
26 // --------------
27 namespace boost {
28
hash_value(const Dakota::StringMultiArray & sma)29 inline size_t hash_value(const Dakota::StringMultiArray& sma)
30 { return boost::hash_range( sma.begin(), sma.end() ); }
31
32 } // namespace boost
33
34 namespace Teuchos {
35
36 template<typename OrdinalType, typename ScalarType>
hash_value(const SerialDenseVector<OrdinalType,ScalarType> & sdv)37 size_t hash_value(const SerialDenseVector<OrdinalType, ScalarType>& sdv)
38 { return boost::hash_range( sdv.values(), sdv.values() + sdv.length() ); }
39
40 } // namespace Teuchos
41
42 namespace Dakota {
43
44
45 // --------------------------------------------
46 // Equality operator definitions for data types
47 // --------------------------------------------
48
49 /// tolerance-based equality operator for RealVector
50 bool nearby(const RealVector& rv1, const RealVector& rv2, Real rel_tol);
51 /// equality operator for IntArray
52 bool operator==(const IntArray& dia1, const IntArray& dia2);
53 /// equality operator for ShortArray
54 bool operator==(const ShortArray& dsa1, const ShortArray& dsa2);
55 /// equality operator for StringArray
56 bool operator==(const StringArray& dsa1, const StringArray& dsa2);
57
58 /// equality operator for std::vector and boost::multi_array::const_array_view
59 template <typename T>
operator ==(const std::vector<T> & vec,typename boost::multi_array<T,1>::template const_array_view<1>::type mav)60 bool operator==(const std::vector<T>& vec,
61 typename boost::multi_array<T, 1>::template
62 const_array_view<1>::type mav)
63 {
64 // Check for equality in array lengths
65 size_t len = vec.size();
66 if ( mav.size() != len )
67 return false;
68 // Check each size_t
69 for (size_t i=0; i<len; ++i)
70 if ( mav[i] != vec[i] )
71 return false;
72 return true;
73 }
74
75 /// equality operator for boost::multi_array::const_array_view and std::vector
76 template <typename T>
operator ==(typename boost::multi_array<T,1>::template const_array_view<1>::type mav,const std::vector<T> & vec)77 bool operator==(typename boost::multi_array<T, 1>::template
78 const_array_view<1>::type mav,
79 const std::vector<T>& vec)
80 { return (vec == mav); } // flip order
81
82 /// equality operator for boost::multi_array and
83 /// boost::multi_array::const_array_view
84 template <typename T>
operator ==(const boost::multi_array<T,1> & ma,typename boost::multi_array<T,1>::template const_array_view<1>::type mav)85 bool operator==(const boost::multi_array<T, 1>& ma,
86 typename boost::multi_array<T, 1>::template
87 const_array_view<1>::type mav)
88 {
89 if (ma.size() != mav.size())
90 return false;
91 // Note: template dependent names require disambiguation
92 typename boost::multi_array<T, 1>::template
93 const_array_view<1>::type::const_iterator cit1 = mav.begin(),
94 end1 = mav.end();
95 typename boost::multi_array<T, 1>::const_iterator cit2 = ma.begin();
96 for ( ; cit1 != end1; ++cit1, ++cit2)
97 if (*cit1 != *cit2)
98 return false;
99 return true;
100 }
101
102 /// equality operator for boost::multi_array::const_array_view and
103 /// boost::multi_array
104 template <typename T>
operator ==(typename boost::multi_array<T,1>::template const_array_view<1>::type mav,const boost::multi_array<T,1> & ma)105 bool operator==(typename boost::multi_array<T, 1>::template
106 const_array_view<1>::type mav,
107 const boost::multi_array<T, 1>& ma)
108 { return (ma == mav); } // flip order
109
110
111 // ----------------------------------------------
112 // Inequality operator definitions for data types
113 // ----------------------------------------------
114
115 /// inequality operator for IntArray
operator !=(const IntArray & dia1,const IntArray & dia2)116 inline bool operator!=(const IntArray& dia1, const IntArray& dia2)
117 { return !(dia1 == dia2); }
118
119 /// inequality operator for ShortArray
operator !=(const ShortArray & dsa1,const ShortArray & dsa2)120 inline bool operator!=(const ShortArray& dsa1, const ShortArray& dsa2)
121 { return !(dsa1 == dsa2); }
122
123 /// inequality operator for StringArray
operator !=(const StringArray & dsa1,const StringArray & dsa2)124 inline bool operator!=(const StringArray& dsa1, const StringArray& dsa2)
125 { return !(dsa1 == dsa2); }
126
127 /// inequality operator for std::vector and boost::multi_array::const_array_view
128 template <typename T>
operator !=(const std::vector<T> & vec,typename boost::multi_array<T,1>::template const_array_view<1>::type mav)129 bool operator!=(const std::vector<T>& vec,
130 typename boost::multi_array<T, 1>::template
131 const_array_view<1>::type mav)
132 { return !(vec == mav); }
133
134 /// inequality operator for boost::multi_array::const_array_view and std::vector
135 template <typename T>
operator !=(typename boost::multi_array<T,1>::template const_array_view<1>::type mav,const std::vector<T> & vec)136 bool operator!=(typename boost::multi_array<T, 1>::template
137 const_array_view<1>::type mav,
138 const std::vector<T>& vec)
139 { return !(vec == mav); } // flip order
140
141 /// inequality operator for boost::multi_array and
142 /// boost::multi_array::const_array_view
143 template <typename T>
operator !=(const boost::multi_array<T,1> & ma,typename boost::multi_array<T,1>::template const_array_view<1>::type mav)144 bool operator!=(const boost::multi_array<T, 1>& ma,
145 typename boost::multi_array<T, 1>::template
146 const_array_view<1>::type mav)
147 { return !(ma == mav); }
148
149 /// inequality operator for boost::multi_array::const_array_view and
150 /// boost::multi_array
151 template <typename T>
operator !=(typename boost::multi_array<T,1>::template const_array_view<1>::type mav,const boost::multi_array<T,1> & ma)152 bool operator!=(typename boost::multi_array<T, 1>::template
153 const_array_view<1>::type mav,
154 const boost::multi_array<T, 1>& ma)
155 { return !(ma == mav); } // flip order
156
157 // ---------------------------------
158 // miscellaneous numerical utilities
159 // ---------------------------------
160
161 /// checks for any non-zero value in std::vector(); useful for determining
162 /// whether an array of request codes (e.g., an ASV) has any actionable content
163 template <typename OrdinalType>
non_zero(const std::vector<OrdinalType> & vec)164 bool non_zero(const std::vector<OrdinalType>& vec)
165 {
166 size_t i, len = vec.size();
167 for (i=0; i<len; ++i)
168 if (vec[i] != 0)
169 return true;
170 return false; // includes case of empty vector (no actions)
171 }
172
173 /// Computes relative change between RealVectors using Euclidean L2 norm
174 Real rel_change_L2(const RealVector& curr_rv, const RealVector& prev_rv);
175
176 /// Computes relative change between Real/int/Real vector triples using
177 /// Euclidean L2 norm
178 Real rel_change_L2(const RealVector& curr_rv1, const RealVector& prev_rv1,
179 const IntVector& curr_iv, const IntVector& prev_iv,
180 const RealVector& curr_rv2, const RealVector& prev_rv2);
181
182 /// equality function for RealVector and a vector of arbitrary type
183 template <typename VectorType>
is_equal_vec(const RealVector & vec1,const VectorType & vec2)184 bool is_equal_vec( const RealVector & vec1,
185 const VectorType & vec2)
186 {
187 // Check for equality in array lengths
188 int len = vec1.length();
189 if ( (int)vec2.size() != len )
190 return false;
191 // Check each size_t
192 for (int i=0; i<len; ++i)
193 if ( vec1[i] != vec2[i] )
194 return false;
195 return true;
196 }
197
198 // ---------------------
199 // Misc matrix utilities
200 // ---------------------
201
202 /// Computes means of columns of matrix
203 void compute_col_means(RealMatrix& matrix, RealVector& avg_vals);
204 /// Computes standard deviations of columns of matrix
205 void compute_col_stdevs(RealMatrix& matrix, RealVector& avg_vals,
206 RealVector& std_devs);
207 /// Removes column from matrix
208 void remove_column(RealMatrix& matrix, int index);
209
210 /// Applies a RealMatrix to a vector (or subset of vector) v1
211 /** Optionally works with a subset of the passed vectors; applies the
212 matrix M to the first M.numCols() entries in v1, and populates the
213 first M.numRows entries in v2. */
214 template<typename MatrixType, typename VectorType>
apply_matrix_partial(const MatrixType & M,const VectorType & v1,VectorType & v2)215 void apply_matrix_partial(const MatrixType& M, const VectorType & v1, VectorType & v2)
216 {
217 if( M.numCols() > v1.size() ) {
218 Cerr << "apply_matrix Error: incoming vector size is inconsistent with matrix column dimension."
219 << std::endl;
220 abort_handler(-1);
221 }
222
223 // Resize target vector if needed
224 if( M.numRows() > v2.size() )
225 v2.resize(M.numRows());
226
227 // Apply the matrix
228 for(size_t i=0; i<M.numRows(); ++i) {
229 v2[i] = 0.0;
230 for (size_t j=0; j<M.numCols(); ++j)
231 v2[i] += M(i,j) * v1[j];
232 }
233 }
234
235 /// Applies transpose of a RealMatrix to a vector (or subset of vector) v1
236 /** Optionally works with a subset of the passed vectors; applies the
237 matrix M^T to the first M.numRows() entries in v1, and populates the
238 first M.numCols() entries in v2. */
239 template<typename VectorType>
apply_matrix_transpose_partial(const RealMatrix & M,const VectorType & v1,VectorType & v2)240 void apply_matrix_transpose_partial(const RealMatrix& M, const VectorType & v1, VectorType & v2)
241 {
242 if( M.numRows() > v1.size() ) {
243 Cerr << "apply_matrix_transpose Error: incoming vector size is inconsistent with matrix row dimension."
244 << std::endl;
245 abort_handler(-1);
246 }
247
248 // Resize target vector if needed
249 if( M.numCols() > v2.size() )
250 v2.resize(M.numCols());
251
252 // Apply the matrix
253 for(size_t j=0; j<M.numCols(); ++j) {
254 v2[j] = 0.0;
255 for (size_t i=0; i<M.numRows(); ++i)
256 v2[j] += M(i,j) * v1[i];
257 }
258 }
259
260 // -----
261 // Utility functions for manipulating or searching strings
262 // -----
263
264 /// Return lowercase copy of string s
strtolower(const std::string & s)265 inline std::string strtolower(const std::string& s)
266 { return boost::to_lower_copy(s); }
267
268 /// Return true if input string begins with string test
strbegins(const std::string & input,const std::string & test)269 inline bool strbegins(const std::string& input, const std::string& test)
270 { return(boost::starts_with(input, test)); }
271
272 /// Return true if input string ends with string test
strends(const std::string & input,const std::string & test)273 inline bool strends(const std::string& input, const std::string& test)
274 { return(boost::ends_with(input, test)); }
275
276 /// Return true if input string contains string test
strcontains(const std::string & input,const std::string & test)277 inline bool strcontains(const std::string& input, const std::string& test)
278 { return(boost::contains(input, test)); }
279
280 /// Trim then split a string on {space, tab} and return as vector of strings
281 std::vector<std::string> strsplit(const std::string& input);
282
283 /// Return the length of the longest string in the passed vector
284 std::string::size_type longest_strlen(const std::vector<std::string>& vecstr);
285
286
287 // --------------------------------------------
288 // Utility functions for creating string arrays
289 // --------------------------------------------
290
291 /// create a label by appending a numerical tag to the root_label, o
build_label(String & label,const String & root_label,size_t tag,const String & separator="")292 inline void build_label(String& label, const String& root_label, size_t tag,
293 const String& separator = "")
294 {
295 label = root_label + separator + std::to_string(tag);
296 }
297
298 /// create an array of labels by tagging root_label for each entry in
299 /// label_array. Uses build_label().
build_labels(StringArray & label_array,const String & root_label)300 inline void build_labels(StringArray& label_array, const String& root_label)
301 {
302 size_t len = label_array.size();
303 for (size_t i=0; i<len; ++i)
304 build_label(label_array[i], root_label, i+1);
305 }
306
307 /// create an array of labels by tagging root_label for each entry in
308 /// label_array. Uses build_label().
build_labels(StringMultiArray & label_array,const String & root_label)309 inline void build_labels(StringMultiArray& label_array,
310 const String& root_label)
311 {
312 size_t len = label_array.size();
313 for (size_t i=0; i<len; ++i)
314 build_label(label_array[i], root_label, i+1);
315 }
316
317 /// create a partial array of labels by tagging root_label for a subset
318 /// of entries in label_array. Uses build_label().
build_labels_partial(StringArray & label_array,const String & root_label,size_t start_index,size_t num_items)319 inline void build_labels_partial(StringArray& label_array,
320 const String& root_label, size_t start_index,
321 size_t num_items)
322 {
323 for (size_t i=0; i<num_items; ++i)
324 build_label(label_array[start_index+i], root_label, i+1);
325 }
326
327
328 // --------------------------
329 // templated assign functions
330 // --------------------------
331
332 /// assign a value to an arbitrary vector
333 template <typename vecType, typename valueType>
assign_value(vecType & target,valueType val)334 void assign_value(vecType& target, valueType val)
335 {
336 size_t i, len = target.size();
337 for (i=0; i<len; ++i)
338 target[i] = val;
339 }
340
341 /// assign a value to a portion of an arbitrary vector
342 template <typename vecType, typename valueType>
assign_value(vecType & target,valueType val,size_t start,size_t len)343 void assign_value(vecType& target, valueType val, size_t start, size_t len)
344 {
345 size_t i, end = start + len;
346 for (i=start; i<end; ++i)
347 target[i] = val;
348 }
349
350 // ----------------------------
351 // non-templated copy functions
352 // ----------------------------
353
354 // ------------------------
355 // templated copy functions
356 // ------------------------
357
358 // Templated conversion functions between/among DAKOTA data types and
359 // built in/pointer data types
360
361
362 ///// copy Array<T> to T*
363 //template <typename T>
364 //void copy_data(const std::vector<T>& vec, T* ptr, const size_t ptr_len)
365 //{
366 // if (ptr_len != vec.size()) { // could use <, but enforce exact match
367 // Cerr << "Error: bad ptr_len in copy_data(Dakota::Array<T>, T* ptr)."
368 // << std::endl;
369 // abort_handler(-1);
370 // }
371 // for (size_t i=0; i<ptr_len; ++i)
372 // ptr[i] = vec[i];
373 //}
374 //
375 //template <typename ScalarType1> //, typename ScalarType1, typename ScalarType2>
376 //void copy_data(const RealVector& vec, ScalarType1* Sptr, const size_t ptr_len, size_t dummy)
377 //{
378 // if (ptr_len != vec.length()) { // could use <, but enforce exact match
379 // Cerr << "Error: bad ptr_len in copy_data(Dakota::Array<T>, T* ptr)."
380 // << std::endl;
381 // abort_handler(-1);
382 // }
383 // for (size_t i=0; i<ptr_len; ++i)
384 // Sptr[i] = vec[i];
385 //}
386 //
387 ///// copy T* to Array<T>
388 //template <typename T>
389 //void copy_data(const T* ptr, const size_t ptr_len, std::vector<T>& vec)
390 //{
391 // if (ptr_len != vec.size())
392 // vec.resize(ptr_len);
393 // for (size_t i=0; i<ptr_len; ++i)
394 // vec[i] = ptr[i];
395 //}
396 //
397 //
398 ///// copy Array<Teuchos::SerialDenseVector<OT,ST> > to ST*
399 //template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
400 //void copy_data(
401 // const std::vector<Teuchos::SerialDenseVector<OrdinalType1, ScalarType> >& va,
402 // ScalarType* ptr, const OrdinalType2 ptr_len, const String& ptr_type)
403 //{
404 // bool c_type;
405 // if (strtolower(ptr_type) == "c") // case insensitive
406 // c_type = true;
407 // else if (strtolower(ptr_type) == "fortran") // case insensitive
408 // c_type = false;
409 // else {
410 // Cerr << "Error: invalid ptr_type in copy_data(Dakota::Array<SDV<OT,ST> >, "
411 // << "ST* ptr)" << std::endl;
412 // abort_handler(-1);
413 // }
414 // OrdinalType2 i, j, num_vec = va.size(), total_len = 0, max_vec_len = 0;
415 // for (i=0; i<num_vec; ++i) { // loop over vectors in array
416 // OrdinalType1 vec_len = va[i].length();
417 // total_len += vec_len;
418 // if (!c_type && vec_len > max_vec_len)
419 // max_vec_len = vec_len;
420 // }
421 // if (ptr_len != total_len) {
422 // Cerr << "Error: bad ptr_len in copy_data(Dakota::Array<Vector<T> >, T* "
423 // << "ptr)." << std::endl;
424 // abort_handler(-1);
425 // }
426 // int cntr = 0;
427 // if (c_type) {
428 // for (i=0; i<num_vec; ++i) { // loop over rows
429 // OrdinalType1 vec_len = va[i].length(); // allowed to vary
430 // for (j=0; j<vec_len; ++j) // loop over columns
431 // ptr[cntr++] = va[i][j];
432 // }
433 // }
434 // else {
435 // for (j=0; j<max_vec_len; ++j) // loop over longest column
436 // for (i=0; i<num_vec; ++i) // loop over rows
437 // if (j < va[i].length())
438 // ptr[cntr++] = va[i][j];
439 // }
440 //}
441
442
443 /// copy Array<Teuchos::SerialDenseVector<OT,ST> > to
444 /// Teuchos::SerialDenseMatrix<OT,ST> - used by read_data_tabular - RWH
445 template <typename OrdinalType, typename ScalarType>
copy_data(const std::vector<Teuchos::SerialDenseVector<OrdinalType,ScalarType>> & sdva,Teuchos::SerialDenseMatrix<OrdinalType,ScalarType> & sdm)446 void copy_data(
447 const std::vector<Teuchos::SerialDenseVector<OrdinalType, ScalarType> >& sdva,
448 Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm)
449 {
450 OrdinalType i, j, num_vec = sdva.size(), max_vec_len = 0;
451 for (i=0; i<num_vec; ++i) { // loop over vectors in array
452 OrdinalType vec_len = sdva[i].length();
453 if (vec_len > max_vec_len)
454 max_vec_len = vec_len;
455 }
456
457 // each vector in array becomes a row in the matrix, with shorter
458 // rows padded with trailing zeros
459 sdm.shape(num_vec, max_vec_len);
460 for (i=0; i<num_vec; ++i) {
461 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec_i = sdva[i];
462 OrdinalType vec_len = vec_i.length(); // allowed to vary
463 for (j=0; j<vec_len; ++j)
464 sdm(i,j) = vec_i[j];
465 }
466 }
467
468
469 /// copy Array<Teuchos::SerialDenseVector<OT,ST> > to transposed
470 /// Teuchos::SerialDenseMatrix<OT,ST> - used by read_data_tabular - RWH
471 template <typename OrdinalType, typename ScalarType>
copy_data_transpose(const std::vector<Teuchos::SerialDenseVector<OrdinalType,ScalarType>> & sdva,Teuchos::SerialDenseMatrix<OrdinalType,ScalarType> & sdm)472 void copy_data_transpose(
473 const std::vector<Teuchos::SerialDenseVector<OrdinalType, ScalarType> >& sdva,
474 Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm)
475 {
476 OrdinalType i, j, num_vec = sdva.size(), max_vec_len = 0;
477 for (i=0; i<num_vec; ++i) { // loop over vectors in array
478 OrdinalType vec_len = sdva[i].length();
479 if (vec_len > max_vec_len)
480 max_vec_len = vec_len;
481 }
482
483 // each vector in array becomes a column in the matrix, with shorter
484 // columns padded with trailing zeros
485 sdm.shape(max_vec_len, num_vec);
486 for (i=0; i<num_vec; ++i) {
487 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec_i = sdva[i];
488 OrdinalType vec_len = vec_i.length(); // allowed to vary
489 ScalarType* sdm_i = sdm[i]; // ith column
490 for (j=0; j<vec_len; ++j)
491 sdm_i[j] = vec_i[j];
492 }
493 }
494
495
496 ///// copy Teuchos::SerialDenseMatrix<OT,ST> to
497 ///// Array<Teuchos::SerialDenseVector<OT,ST> >
498 //template <typename OrdinalType, typename ScalarType>
499 //void copy_data(
500 // const Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm,
501 // std::vector<Teuchos::SerialDenseVector<OrdinalType, ScalarType> >& sdva)
502 //{
503 // OrdinalType i, j, num_vec = sdm.numRows(), vec_len = sdm.numCols();
504 // sdva.resize(num_vec);
505 // for (i=0; i<num_vec; ++i) {
506 // Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec_i = sdva[i];
507 // vec_i.sizeUninitialized(vec_len);
508 // for (j=0; j<vec_len; ++j)
509 // vec_i[j] = sdm(i,j);
510 // }
511 //}
512 //
513 //
514 ///// copy Teuchos::SerialDenseMatrix<OT,ST> to transposed
515 ///// Array<Teuchos::SerialDenseVector<OT,ST> >
516 //template <typename OrdinalType, typename ScalarType>
517 //void copy_data_transpose(
518 // const Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm,
519 // std::vector<Teuchos::SerialDenseVector<OrdinalType, ScalarType> >& sdva)
520 //{
521 // OrdinalType i, j, num_vec = sdm.numCols(), vec_len = sdm.numRows();
522 // sdva.resize(num_vec);
523 // for (i=0; i<num_vec; ++i) {
524 // Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec_i = sdva[i];
525 // const ScalarType* sdm_i = sdm[i]; // ith column
526 // vec_i.sizeUninitialized(vec_len);
527 // for (j=0; j<vec_len; ++j)
528 // vec_i[j] = sdm_i[j];
529 // }
530 //}
531
532
533 /// copy Teuchos::SerialDenseVector<OT,ST> to Teuchos::SerialDenseMatrix<OT,ST> - used by NestedModel::update_sub_iterator - RWH
534 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv,Teuchos::SerialDenseMatrix<OrdinalType1,ScalarType> & sdm,OrdinalType2 nr,OrdinalType2 nc)535 void copy_data(const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv,
536 Teuchos::SerialDenseMatrix<OrdinalType1, ScalarType>& sdm,
537 OrdinalType2 nr, OrdinalType2 nc)
538 {
539 OrdinalType1 size_sdv = sdv.length();
540
541 // This function is set up to do the transformation with either nr or nc or
542 // both specified. To omit nr or nc specification, a 0 is passed.
543 if (nr && nc) { // both specified
544 if (size_sdv != nr*nc) {
545 Cerr << "Error: sdv length (" << size_sdv << ") does not equal nr*nc ("
546 << nr << '*' << nc << ") in copy_data(Teuchos_SerialDenseVector<>, "
547 << "Teuchos_SerialDenseMatrix<>)." << std::endl;
548 abort_handler(-1);
549 }
550 }
551 else if (nr) { // only nr is non-zero
552 if (size_sdv%nr) {
553 Cerr << "Error: sdv length (" << size_sdv << ") not evenly divisible by "
554 << "number of rows (" << nr << ") in copy_data(Teuchos_"
555 << "SerialDenseVector<>, Teuchos_SerialDenseMatrix<>)." << std::endl;
556 abort_handler(-1);
557 }
558 nc = size_sdv/nr;
559 }
560 else if (nc) { // only nc is non-zero
561 if (size_sdv%nc) {
562 Cerr << "Error: sdv length (" << size_sdv << ") not evenly divisible by "
563 << "number of columns (" << nc << ") in copy_data(Teuchos_"
564 << "SerialDenseVector<>, Teuchos_SerialDenseMatrix<>)." << std::endl;
565 abort_handler(-1);
566 }
567 nr = size_sdv/nc;
568 }
569 else { // neither specified
570 Cerr << "Error: either nr or nc must be specified in copy_data(Teuchos_"
571 << "SerialDenseVector<>, Teuchos_SerialDenseMatrix<>)." << std::endl;
572 abort_handler(-1);
573 }
574
575 if (sdm.numRows() != nr || sdm.numCols() != nc)
576 sdm.shapeUninitialized(nr, nc);
577 // sdv is head to tail by rows, which matches the visual matrix format that
578 // a user would employ in specifying a matrix as a <LISTof><REAL>
579 OrdinalType1 counter = 0;
580 for (OrdinalType1 i=0; i<nr; ++i)
581 for (OrdinalType1 j=0; j<nc; ++j, ++counter)
582 sdm(i,j) = sdv[counter];
583 }
584
585 ///// copy std::list<T> to std::vector<T>
586 //template <typename T>
587 //void copy_data(const std::list<T>& dl, std::vector<T>& da)
588 //{
589 // size_t size_dl = dl.size();
590 // if (size_dl != da.size())
591 // da.resize(size_dl);
592 // std::copy(dl.begin(), dl.end(), da.begin());
593 //}
594 //
595 ///// copy std::list<T> to std::vector<std::vector<T> >
596 //template <typename T>
597 //void copy_data(const std::list<T>& dl, std::vector<std::vector<T> >& d2a,
598 // size_t num_a, size_t a_len)
599 //{
600 // size_t i, j, size_dl = dl.entries();
601 //
602 // // This function is set up to do the copy with either num_a or a_len or both
603 // // specified. To omit num_a or a_len specification, a 0 is passed.
604 // if (num_a && a_len) { // both specified
605 // if (size_dl != num_a*a_len) {
606 // Cerr << "Error: dl length (" << size_dl <<") does not equal num_a*a_len ("
607 // << num_a << '*' << a_len << ") in copy_data(std::list<T>, "
608 // << "std::vector<vector<T> >)." << std::endl;
609 // abort_handler(-1);
610 // }
611 // }
612 // else if (num_a) { // only num_a is non-zero
613 // if (size_dl%num_a) {
614 // Cerr << "Error: dl length (" << size_dl << ") not evenly divisible by "
615 // << "number of arrays (" << num_a << ") in copy_data(std::list<T>"
616 // << ", std::vector<vector<T> >)." << std::endl;
617 // abort_handler(-1);
618 // }
619 // a_len = size_dl/num_a;
620 // }
621 // else if (a_len) { // only a_len is non-zero
622 // if (size_dl%a_len) {
623 // Cerr << "Error: dl length (" << size_dl << ") not evenly divisible by "
624 // << "array length (" << a_len << ") in copy_data(std::list<T>, "
625 // << "std::vector<vector<T> >)." << std::endl;
626 // abort_handler(-1);
627 // }
628 // num_a = size_dl/a_len;
629 // }
630 // else { // neither specified
631 // Cerr << "Error: either num_a or a_len must be specified in "
632 // << "copy_data(std::list<T>,std::vector<vector<T> >)." << std::endl;
633 // abort_handler(-1);
634 // }
635 //
636 // if (d2a.size() != num_a)
637 // d2a.reshape(num_a);
638 // typename std::list<T>::const_iterator dl_cit = dl.begin();
639 // for (i=0; i<num_a; ++i) {
640 // if (d2a[i].size() != a_len)
641 // d2a[i].reshape(a_len);
642 // for (j=0; j<a_len; ++j, ++dl_cit)
643 // d2a[i][j] = *dl_cit;
644 // }
645 //}
646
647 /// copy std::vector<vector<T> > to std::vector<T>(unroll vecOfvecs into vector) - used by ProcessApplicInterface::write_parameters_files - RWH
648 template <typename T>
copy_data(const std::vector<std::vector<T>> & d2a,std::vector<T> & da)649 void copy_data(const std::vector<std::vector<T> >& d2a, std::vector<T>& da)
650 {
651 typename std::vector<T>::size_type i, j, size_d2a = 0,
652 num_arrays = d2a.size(), cntr = 0;
653 for (i=0; i<num_arrays; ++i)
654 size_d2a += d2a[i].size();
655 if (size_d2a != da.size())
656 da.resize(size_d2a);
657 for (i=0; i<num_arrays; ++i) {
658 typename std::vector<T>::size_type array_len = d2a[i].size();
659 for (j=0; j<array_len; ++j)
660 da[cntr++] = d2a[i][j];
661 }
662 }
663
664 /// copy map<int, T> to std::vector<T> (discard integer keys) - used by SurrBasedGlobalMinimizer::core_run - RWH
665 template <typename T>
copy_data(const std::map<int,T> & im,std::vector<T> & da)666 void copy_data(const std::map<int, T>& im, std::vector<T>& da)
667 {
668 size_t i, size_im = im.size();
669 if (size_im != da.size())
670 da.resize(size_im);
671 typename std::map<int, T>::const_iterator im_cit = im.begin();
672 for (i=0; i<size_im; ++i, ++im_cit)
673 da[i] = im_cit->second;
674 }
675
676 /// copy Teuchos::SerialDenseVector<OrdinalType, ScalarType> to same
677 /// (used in place of operator= when a deep copy is required) -
678 /// used by Response - MSE
679 template <typename OrdinalType, typename ScalarType>
copy_data(const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & sdv1,Teuchos::SerialDenseVector<OrdinalType,ScalarType> & sdv2)680 void copy_data(const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& sdv1,
681 Teuchos::SerialDenseVector<OrdinalType, ScalarType>& sdv2)
682 {
683 OrdinalType size_sdv1 = sdv1.length();
684 if (sdv2.length() != size_sdv1)
685 sdv2.sizeUninitialized(size_sdv1);
686 sdv2.assign(sdv1);
687 }
688
689 /// copy Teuchos::SerialDenseMatrix<OrdinalType, ScalarType> to same
690 /// (used in place of operator= when a deep copy is required) -
691 /// used by Response - MSE
692 template <typename OrdinalType, typename ScalarType>
copy_data(const Teuchos::SerialDenseMatrix<OrdinalType,ScalarType> & sdm1,Teuchos::SerialDenseMatrix<OrdinalType,ScalarType> & sdm2)693 void copy_data(const Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm1,
694 Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& sdm2)
695 {
696 OrdinalType nr1 = sdm1.numRows(), nc1 = sdm1.numCols();
697 if (sdm2.numRows() != nr1 || sdm2.numCols() != nc1)
698 sdm2.shapeUninitialized(nr1, nc1);
699 sdm2.assign(sdm1);
700 }
701
702 /// copy Teuchos::SerialSymDenseMatrix<OrdinalType, ScalarType> to same
703 /// (used in place of operator= when a deep copy is required) -
704 /// used by Response - MSE
705 template <typename OrdinalType, typename ScalarType>
copy_data(const Teuchos::SerialSymDenseMatrix<OrdinalType,ScalarType> & ssdm1,Teuchos::SerialSymDenseMatrix<OrdinalType,ScalarType> & ssdm2)706 void copy_data(const Teuchos::SerialSymDenseMatrix<OrdinalType, ScalarType>& ssdm1,
707 Teuchos::SerialSymDenseMatrix<OrdinalType, ScalarType>& ssdm2)
708 {
709 OrdinalType nr1 = ssdm1.numRows();
710 if (ssdm2.numRows() != nr1)
711 ssdm2.shapeUninitialized(nr1);
712 ssdm2.assign(ssdm1);
713 }
714
715 /// Taken from pecos/src/MathTools.hpp, BUT
716 /// not templated because the implementation is specific to RealMatrix
copy_data(const RealMatrix & source,RealMatrix & dest,int num_rows,int num_cols,int start_row=0,int start_col=0)717 inline void copy_data( const RealMatrix &source, RealMatrix &dest,
718 int num_rows, int num_cols, int start_row=0, int start_col=0 )
719 {
720 RealMatrix source_subset( Teuchos::View, source, num_rows, num_cols,
721 start_row, start_col );
722 dest.reshape( num_rows, num_cols );
723 dest.assign( source_subset );
724 }
725
726 /// copy Teuchos::SerialDenseVector<OrdinalType, ScalarType> to
727 /// VecType - used by APPS for HOPS vector types
728 template <typename OrdinalType, typename ScalarType, typename VecType>
copy_data(const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & sdv,VecType & vec)729 void copy_data(const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& sdv,
730 VecType& vec)
731 {
732 OrdinalType size_sdv = sdv.length();
733 if (size_sdv != vec.size())
734 vec.resize(size_sdv);
735 for (OrdinalType i=0; i<size_sdv; ++i)
736 vec[i] = sdv[i];
737 }
738
739 /// copy Teuchos::SerialDenseVector<OrdinalType, ScalarType> to
740 /// std::vector<ScalarType> - used by DakotaModel
741 template <typename OrdinalType, typename ScalarType1, typename ScalarType2>
copy_data(const Teuchos::SerialDenseVector<OrdinalType,ScalarType1> & sdv,std::vector<ScalarType2> & vec)742 void copy_data(const Teuchos::SerialDenseVector<OrdinalType, ScalarType1>& sdv,
743 std::vector<ScalarType2>& vec)
744 {
745 OrdinalType size_sdv = sdv.length();
746 if (size_sdv != vec.size())
747 vec.resize(size_sdv);
748 for (OrdinalType i=0; i<size_sdv; ++i)
749 vec[i] = (ScalarType2)sdv[i];
750 }
751
752 /// copy Array<ScalarType> to
753 /// Teuchos::SerialDenseVector<OrdinalType, ScalarType> - used by NOWPACOptimizer - MSE
754 template <typename OrdinalType, typename ScalarType>
copy_data(const std::vector<ScalarType> & da,Teuchos::SerialDenseVector<OrdinalType,ScalarType> & sdv)755 void copy_data(const std::vector<ScalarType>& da,
756 Teuchos::SerialDenseVector<OrdinalType, ScalarType>& sdv)
757 {
758 size_t size_da = da.size();
759 if (sdv.length() != size_da)
760 sdv.sizeUninitialized(size_da);
761 for (OrdinalType i=0; i<size_da; ++i)
762 sdv[i] = da[i];
763 }
764
765 /// copy ScalarType* to Teuchos::SerialDenseVector<OrdinalType, ScalarType> - used by ScalingModel::response_modify_n2s - RWH
766 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data(const ScalarType * ptr,const OrdinalType2 ptr_len,Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv)767 void copy_data(const ScalarType* ptr, const OrdinalType2 ptr_len,
768 Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv)
769 {
770 if (sdv.length() != ptr_len)
771 sdv.sizeUninitialized(ptr_len);
772 for (OrdinalType2 i=0; i<ptr_len; ++i)
773 sdv[i] = ptr[i];
774
775 //sdv = RealVector(Copy, ptr, ptr_len);
776 // not advisable since relying on default operator=. Would be a good
777 // approach if a reference counting scheme was used for operator=.
778 }
779
780 /// copy ScalarType* to Teuchos::SerialDenseVector<OrdinalType, ScalarType> - used by NL2SOLLeastSq::core_run - RWH
781 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv,ScalarType * ptr,const OrdinalType2 ptr_len)782 void copy_data(const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv,
783 ScalarType* ptr, const OrdinalType2 ptr_len)
784 {
785 if (ptr_len != sdv.length()) { // could use <, but enforce exact match
786 Cerr << "Error: bad ptr_len in copy_data(Teuchos::SerialDenseVector<>, "
787 << "T* ptr)." << std::endl;
788 abort_handler(-1);
789 }
790 for (OrdinalType2 i=0; i<ptr_len; ++i)
791 ptr[i] = sdv[i];
792 }
793
794
795 /// copy SerialDenseVector<> to Array<SerialDenseVector<> > - used by ConcurrentMetaIterator constructor - RWH
796 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv,std::vector<Teuchos::SerialDenseVector<OrdinalType1,ScalarType>> & sdva,OrdinalType2 num_vec,OrdinalType2 vec_len)797 void copy_data(const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv,
798 std::vector<Teuchos::SerialDenseVector<OrdinalType1, ScalarType> >& sdva,
799 OrdinalType2 num_vec, OrdinalType2 vec_len)
800 {
801 // This function is set up to do the transformation with either num_vec or
802 // vec_len or both specified. To omit num_vec or vec_len specification, a 0
803 // is passed.
804 OrdinalType1 size_sdv = sdv.length();
805 if (num_vec && vec_len) { // both specified
806 if (size_sdv != num_vec*vec_len) {
807 Cerr << "Error: sdv length (" << size_sdv << ") does not equal num_vec*"
808 << "vec_len (" << num_vec << '*' << vec_len << ") in copy_data("
809 << "Teuchos::SerialDenseVector<>, Dakota::Array<Teuchos::"
810 << "SerialDenseVector<> >)." << std::endl;
811 abort_handler(-1);
812 }
813 }
814 else if (num_vec) { // only num_vec is non-zero
815 if (size_sdv%num_vec) {
816 Cerr << "Error: sdv length (" << size_sdv << ") not evenly divisible by "
817 << "number of vectors (" << num_vec << ") in copy_data("
818 << "Teuchos::SerialDenseVector<>, Dakota::Array<Teuchos::"
819 << "SerialDenseVector<> >)." << std::endl;
820 abort_handler(-1);
821 }
822 vec_len = size_sdv/num_vec;
823 }
824 else if (vec_len) { // only vec_len is non-zero
825 if (size_sdv%vec_len) {
826 Cerr << "Error: sdv length (" << size_sdv << ") not evenly divisible by "
827 << "vector length (" << vec_len << ") in copy_data(Teuchos::"
828 << "SerialDenseVector<>, Dakota::Array<Teuchos::"
829 << "SerialDenseVector<> >)." << std::endl;
830 abort_handler(-1);
831 }
832 num_vec = size_sdv/vec_len;
833 }
834 else { // neither specified
835 Cerr << "Error: either num_vec or vec_len must be specified in "
836 << "copy_data(Teuchos::SerialDenseVector<>, Dakota::Array<Teuchos::"
837 << "SerialDenseVector<> >)." << std::endl;
838 abort_handler(-1);
839 }
840
841 if (sdva.size() != num_vec)
842 sdva.resize(num_vec);
843 // sdv is head to tail by rows, which matches the visual matrix format that
844 // a user would employ in specifying a matrix as a <LISTof><REAL>
845 OrdinalType2 counter = 0;
846 for (OrdinalType2 i=0; i<num_vec; ++i) {
847 if (sdva[i].length() != vec_len)
848 sdva[i].sizeUninitialized(vec_len);
849 for (OrdinalType2 j=0; j<vec_len; ++j)
850 sdva[i][j] = sdv[counter++];
851 }
852 }
853
854
855 // Partial copy functions
856
857 /// copy a portion arbitrary vector to all of another arbitrary vector
858 template <typename vecType1, typename vecType2>
copy_data_partial(const vecType1 & source,size_t source_start_idx,vecType2 & target,size_t target_start_idx,size_t len)859 void copy_data_partial(
860 const vecType1& source,
861 size_t source_start_idx,
862 vecType2& target,
863 size_t target_start_idx,
864 size_t len)
865 {
866 // This requires that the target type supports the size() method, which RealVectors don't
867 //if( len != target.size()) { // could use <, but enforce exact match
868 // Cerr << "Error: bad target vector length copy_data_partial." << std::endl;
869 // abort_handler(-1);
870 //}
871 for( size_t i=0; i<len; ++i)
872 target[i+target_start_idx] = source[i+source_start_idx];
873 }
874
875 /// copy portion of first SerialDenseVector to all of second SerialDenseVector - used by DataTransformModel::vars_mapping - RWH
876 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data_partial(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv1,OrdinalType2 start_index1,OrdinalType2 num_items,Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv2)877 void copy_data_partial(
878 const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv1,
879 OrdinalType2 start_index1, OrdinalType2 num_items,
880 Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv2)
881 {
882 // sdv1 will be indexed from start_index1 to start_index1+num_items-1
883 if (start_index1 + num_items > sdv1.length()) {
884 Cerr << "Error: indexing out of bounds in copy_data_partial("
885 << "Teuchos::SerialDenseVector<OrdinalType, ScalarType>, size_t, "
886 << "size_t, Teuchos::SerialDenseVector<OrdinalType, ScalarType>)."
887 << std::endl;
888 abort_handler(-1);
889 }
890 if (num_items != sdv2.length())
891 sdv2.sizeUninitialized(num_items);
892 for (OrdinalType2 i=0; i<num_items; ++i)
893 sdv2[i] = sdv1[start_index1+i];
894 }
895
896 /// copy all of first SerialDenseVector to portion of second SerialDenseVector - used by MixedVariables - RWH, NLSSOLLeastSq - BMA
897 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data_partial(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv1,Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv2,OrdinalType2 start_index2)898 void copy_data_partial(
899 const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv1,
900 Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv2,
901 OrdinalType2 start_index2)
902 {
903 OrdinalType1 num_items = sdv1.length();
904 // In this case, incoming sdv2 must already be sized and will be
905 // indexed from start_index2 to start_index2+num_items-1
906 if (start_index2 + num_items > sdv2.length()) {
907 Cerr << "Error: indexing out of bounds in copy_data_partial("
908 << "Teuchos::SerialDenseVector<OrdinalType, ScalarType>, "
909 << "Teuchos::SerialDenseVector<OrdinalType, ScalarType>, OrdinalType)."
910 << std::endl;
911 abort_handler(-1);
912 }
913 for (OrdinalType1 i=0; i<num_items; ++i)
914 sdv2[start_index2+i] = sdv1[i];
915 }
916
917 /// copy portion of first SerialDenseVector to portion of second
918 /// SerialDenseVector - used by ScalingModel::secondary_resp_scaled2native - RWH
919 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data_partial(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv1,OrdinalType2 start_index1,OrdinalType2 num_items,Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv2,OrdinalType2 start_index2)920 void copy_data_partial(
921 const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv1,
922 OrdinalType2 start_index1, OrdinalType2 num_items,
923 Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv2,
924 OrdinalType2 start_index2)
925 {
926 // In this case, incoming sdv2 must already be sized and will be
927 // indexed from start_index2 to start_index2+num_items-1. sdv1 will
928 // be indexed from start_index1 to start_index1+num_items-1
929 if (start_index1 + num_items > sdv1.length() ||
930 start_index2 + num_items > sdv2.length()) {
931 Cerr << "Error: indexing out of bounds in copy_data_partial("
932 << "Teuchos::SerialDenseVector<OrdinalType, ScalarType>, OrdinalType, "
933 << "OrdinalType, Teuchos::SerialDenseVector<OrdinalType, ScalarType>, "
934 << "OrdinalType)." << std::endl;
935 abort_handler(-1);
936 }
937 for (OrdinalType2 i=0; i<num_items; ++i)
938 sdv2[start_index2+i] = sdv1[start_index1+i];
939 }
940
941 /// copy all of first SerialDenseVector to portion of second SerialDenseVector - used by SharedSurfpackApproxData::merge_variable_arrays - RWH
942 template <typename OrdinalType1, typename OrdinalType2, typename ScalarType>
copy_data_partial(const Teuchos::SerialDenseVector<OrdinalType1,ScalarType> & sdv1,std::vector<ScalarType> & da2,OrdinalType2 start_index2)943 void copy_data_partial(
944 const Teuchos::SerialDenseVector<OrdinalType1, ScalarType>& sdv1,
945 std::vector<ScalarType>& da2, OrdinalType2 start_index2)
946 {
947 OrdinalType1 num_items = sdv1.length();
948 // In this case, incoming da2 must already be sized and will be
949 // indexed from start_index2 to start_index2+num_items-1
950 if (start_index2 + num_items > da2.size()) {
951 Cerr << "Error: indexing out of bounds in copy_data_partial(Teuchos::"
952 << "SerialDenseVector<OrdinalType, ScalarType>, "
953 << "std::vector<ScalarType>, OrdinalType)." << std::endl;
954 abort_handler(-1);
955 }
956 for (OrdinalType1 i=0; i<num_items; ++i)
957 da2[start_index2+i] = sdv1[i];
958 }
959
960 /// copy portion of first Array<T> to all of second Array<T> - used by SharedResponseDataRep constructor - RWH
961 template <typename T>
copy_data_partial(const std::vector<T> & da1,size_t start_index1,size_t num_items,std::vector<T> & da2)962 void copy_data_partial(const std::vector<T>& da1, size_t start_index1,
963 size_t num_items, std::vector<T>& da2)
964 {
965 // da1 will be indexed from start_index1 to start_index1+num_items-1
966 if (start_index1 + num_items > da1.size()) {
967 Cerr << "Error: indexing out of bounds in copy_data_partial("
968 << "Dakota::Array<T>, size_t, size_t, Dakota::Array<T>)." << std::endl;
969 abort_handler(-1);
970 }
971 if (num_items != da2.size())
972 da2.resize(num_items);
973 for (size_t i=0; i<num_items; ++i)
974 da2[i] = da1[start_index1+i];
975 }
976
977 /// copy all of first Array<T> to portion of second Array<T> - used by ParamStudy::multidim_loop - RWH
978 template <typename T>
copy_data_partial(const std::vector<T> & da1,std::vector<T> & da2,size_t start_index2)979 void copy_data_partial(const std::vector<T>& da1, std::vector<T>& da2,
980 size_t start_index2)
981 {
982 size_t i, num_items = da1.size();
983 // In this case, incoming da2 must already be sized and will be
984 // indexed from start_index2 to start_index2+num_items-1
985 if (start_index2 + num_items > da2.size()) {
986 Cerr << "Error: indexing out of bounds in copy_data_partial("
987 << "Dakota::Array<T>, Dakota::Array<T>, size_t)." << std::endl;
988 abort_handler(-1);
989 }
990 for (i=0; i<num_items; ++i)
991 da2[start_index2+i] = da1[i];
992 }
993
994 /// copy all of first Array<T> to portion of boost::multi_array<T, 1> - used by RelaxedVariables - RWH
995 template <typename T>
copy_data_partial(const std::vector<T> & da,boost::multi_array<T,1> & bma,size_t start_index_bma)996 void copy_data_partial(const std::vector<T>& da, boost::multi_array<T, 1>& bma,
997 size_t start_index_bma)
998 {
999 size_t i, num_items = da.size();
1000 // In this case, incoming bma must already be sized and will be
1001 // indexed from start_index_bma to start_index_bma+num_items-1
1002 if (start_index_bma + num_items > bma.size()) {
1003 Cerr << "Error: indexing out of bounds in copy_data_partial("
1004 << "Dakota::Array<T>, boost::multi_array<T, 1>, size_t)." << std::endl;
1005 abort_handler(-1);
1006 }
1007 for (i=0; i<num_items; ++i)
1008 bma[start_index_bma+i] = da[i];
1009 }
1010
1011 ///// copy portion of first Array<T> to portion of second Array<T>
1012 //template <typename T>
1013 //void copy_data_partial(const std::vector<T>& da1, size_t start_index1,
1014 // size_t num_items,std::vector<T>& da2,size_t start_index2)
1015 //{
1016 // // In this case, incoming da2 must already be sized and will be
1017 // // indexed from start_index2 to start_index2+num_items-1. da1 will
1018 // // be indexed from start_index1 to start_index1+num_items-1
1019 // if (start_index1 + num_items > da1.size() ||
1020 // start_index2 + num_items > da2.size()) {
1021 // Cerr << "Error: indexing out of bounds in copy_data_partial(Dakota::"
1022 // << "Array<T>, size_t, size_t, Dakota::Array<T>, size_t)." << std::endl;
1023 // abort_handler(-1);
1024 // }
1025 // for (size_t i=0; i<num_items; ++i)
1026 // da2[start_index2+i] = da1[start_index1+i];
1027 //}
1028
1029 /// Copies a column of a Teuchos_SerialDenseMatrix<int,Real> to std::vector<Real>
1030 template<typename VectorType>
copy_column_vector(const RealMatrix & m,RealMatrix::ordinalType j,VectorType & col)1031 void copy_column_vector(const RealMatrix& m,
1032 RealMatrix::ordinalType j,
1033 VectorType& col)
1034 {
1035 RealMatrix::ordinalType i, num_items = m.numRows();
1036 if (col.size() != num_items)
1037 col.resize(num_items);
1038 for(i=0; i<num_items; ++i)
1039 col[i] = m(i,j);
1040 }
1041
1042 /// Copies a row of a Teuchos_SerialDenseMatrix<int,Real> to std::vector<Real>
1043 template<typename VectorType>
copy_row_vector(const RealMatrix & m,RealMatrix::ordinalType i,VectorType & row)1044 void copy_row_vector(const RealMatrix& m, RealMatrix::ordinalType i,
1045 VectorType& row)
1046 {
1047 RealMatrix::ordinalType j, num_items = m.numCols();
1048 if (row.size() != num_items)
1049 row.resize(num_items);
1050 for(j=0; j<num_items; ++j)
1051 row[j] = m(i,j);
1052 }
1053
1054
1055 /// Inserts a std::vector<Real> into a row of a Teuchos_SerialDenseMatrix<int,Real>
1056 template<typename ScalarType>
insert_row_vector(const std::vector<ScalarType> & row,RealMatrix::ordinalType i,RealMatrix & m)1057 void insert_row_vector(const std::vector<ScalarType>& row,
1058 RealMatrix::ordinalType i,
1059 RealMatrix& m)
1060 {
1061 if( (m.numRows() < i+1) || (m.numCols() != row.size()) )
1062 m.reshape(i+1, row.size());
1063 for( size_t j=0; j<row.size(); ++j)
1064 m(i,j) = row[j];
1065 }
1066
1067
1068 // ----------------------------------------------------
1069 // Non-templated functions for creating relaxed vectors
1070 // ----------------------------------------------------
1071
1072 /// merge a discrete integer vector into a single continuous vector
merge_data_partial(const IntVector & d_vec,RealVector & m_vec,size_t start_index_ma)1073 inline void merge_data_partial(const IntVector& d_vec,
1074 RealVector& m_vec, size_t start_index_ma)
1075 {
1076 size_t i, num_items = d_vec.length();
1077 // In this case, incoming m_vec must already be sized and will be
1078 // indexed from start_index_ma to start_index_ma+num_items-1
1079 size_t m_vec_len = m_vec.length();
1080 if (start_index_ma + num_items > m_vec_len) {
1081 Cerr << "Error: indexing out of bounds in merge_data_partial(IntVector, "
1082 << "RealVector, size_t)." << std::endl;
1083 abort_handler(-1);
1084 }
1085 for (i=0; i<num_items; ++i)
1086 m_vec[start_index_ma+i] = (Real)d_vec[i];
1087 }
1088
1089 /// merge a discrete integer vector into a single continuous array
merge_data_partial(const IntVector & d_vec,RealArray & m_array,size_t start_index_ma)1090 inline void merge_data_partial(const IntVector& d_vec,
1091 RealArray& m_array, size_t start_index_ma)
1092 {
1093 size_t i, num_items = d_vec.length();
1094 // In this case, incoming m_array must already be sized and will be
1095 // indexed from start_index_ma to start_index_ma+num_items-1
1096 if (start_index_ma + num_items > m_array.size()) {
1097 Cerr << "Error: indexing out of bounds in merge_data_partial(IntVector, "
1098 << "RealArray, size_t)." << std::endl;
1099 abort_handler(-1);
1100 }
1101 for (i=0; i<num_items; ++i)
1102 m_array[start_index_ma+i] = (Real)d_vec[i];
1103 }
1104
1105
1106 /// round entries of a RealVector yielding an IntVector
1107 void iround(const RealVector& input_vec, IntVector& rounded_vec);
1108
1109
1110 // -------------------------------
1111 // templated set utility functions
1112 // -------------------------------
1113
1114 /// retrieve the set value corresponding to the passed index
1115 template <typename OrdinalType, typename ScalarType>
set_index_to_value(OrdinalType index,const std::set<ScalarType> & values)1116 const ScalarType& set_index_to_value(OrdinalType index,
1117 const std::set<ScalarType>& values)
1118 {
1119 // TO DO: conditional activation for automatic bounds checking
1120 if (index < 0 || index >= values.size())
1121 throw std::out_of_range(String("Error: index ") + std::to_string(index) +
1122 " must be between 0 and " + std::to_string(values.size() - 1) +
1123 " in set_index_to_value()");
1124
1125 typename std::set<ScalarType>::const_iterator cit = values.begin();
1126 std::advance(cit, index);
1127 return *cit;
1128 }
1129
1130
1131 /// calculate the set index corresponding to the passed value
1132 template <typename ScalarType>
set_value_to_index(const ScalarType & value,const std::set<ScalarType> & values)1133 size_t set_value_to_index(const ScalarType& value,
1134 const std::set<ScalarType>& values)
1135 {
1136 typename std::set<ScalarType>::const_iterator cit = values.find(value);
1137 return (cit == values.end()) ? _NPOS : std::distance(values.begin(), cit);
1138
1139 // linear search provides index in one pass, but find() plus distance()
1140 // should be faster for sorted associative containers
1141 //size_t index = 0;
1142 //typename std::set<ScalarType>::const_iterator cit;
1143 //for (cit=values.begin(); cit!=values.end(); ++cit, ++index)
1144 // if (*cit == value)
1145 // return index;
1146 //return _NPOS;
1147 }
1148
1149
1150 /// retrieve the set value corresponding to the passed index
1151 template <typename OrdinalType, typename KeyType, typename ValueType>
map_index_to_key(OrdinalType index,const std::map<KeyType,ValueType> & pairs)1152 const KeyType& map_index_to_key(OrdinalType index,
1153 const std::map<KeyType, ValueType>& pairs)
1154 {
1155 // TO DO: conditional activation for automatic bounds checking
1156 if (index < 0 || index >= pairs.size()) {
1157 Cerr << "\nError: index out of range in map_index_to_key()" << std::endl;
1158 abort_handler(-1);
1159 }
1160 typename std::map<KeyType, ValueType>::const_iterator cit = pairs.begin();
1161 std::advance(cit, index);
1162 return cit->first;
1163 }
1164
1165
1166 /// retrieve the set value corresponding to the passed index
1167 template <typename OrdinalType, typename KeyType, typename ValueType>
map_index_to_value(OrdinalType index,const std::map<KeyType,ValueType> & pairs)1168 const ValueType& map_index_to_value(OrdinalType index,
1169 const std::map<KeyType, ValueType>& pairs)
1170 {
1171 // TO DO: conditional activation for automatic bounds checking
1172 if (index < 0 || index >= pairs.size()) {
1173 Cerr << "\nError: index out of range in map_index_to_value()" << std::endl;
1174 abort_handler(-1);
1175 }
1176 typename std::map<KeyType, ValueType>::const_iterator cit = pairs.begin();
1177 std::advance(cit, index);
1178 return cit->second;
1179 }
1180
1181
1182 /// calculate the map index corresponding to the passed key
1183 template <typename KeyType, typename ValueType>
map_keys_to_set(const std::map<KeyType,ValueType> & source_map,std::set<KeyType> & target_set)1184 void map_keys_to_set(const std::map<KeyType, ValueType>& source_map,
1185 std::set<KeyType>& target_set)
1186 {
1187 target_set.clear();
1188 typename std::map<KeyType, ValueType>::const_iterator cit;
1189 for (cit=source_map.begin(); cit!=source_map.end(); ++cit)
1190 target_set.insert(cit->first);
1191 }
1192
1193
1194 /// calculate the map index corresponding to the passed key
1195 template <typename KeyType, typename ValueType>
map_key_to_index(const KeyType & key,const std::map<KeyType,ValueType> & pairs)1196 size_t map_key_to_index(const KeyType& key,
1197 const std::map<KeyType, ValueType>& pairs)
1198 {
1199 typename std::map<KeyType, ValueType>::const_iterator cit = pairs.find(key);
1200 return (cit == pairs.end()) ? _NPOS : std::distance(pairs.begin(), cit);
1201
1202 // linear search provides index in one pass, but find() plus distance()
1203 // should be faster for sorted associative containers
1204 //size_t index = 0;
1205 //typename std::map<KeyType, ValueType>::const_iterator cit;
1206 //for (cit=pairs.begin(); cit!=pairs.end(); ++cit, ++index)
1207 // if (cit->first == value)
1208 // return index;
1209 //return _NPOS;
1210 }
1211
1212
1213 /// calculate the map index corresponding to the passed value
1214 template <typename KeyType, typename ValueType>
map_value_to_index(const ValueType & value,const std::map<KeyType,ValueType> & pairs)1215 size_t map_value_to_index(const ValueType& value,
1216 const std::map<KeyType, ValueType>& pairs)
1217 {
1218 size_t index = 0;
1219 typename std::map<KeyType, ValueType>::const_iterator cit;
1220 for (cit=pairs.begin(); cit!=pairs.end(); ++cit, ++index)
1221 if (cit->second == value)
1222 return index;
1223 return _NPOS;
1224 }
1225
1226
1227 /// convert a SerialDenseVector of head-to-tail (x,y) pairs into a
1228 /// std::set of (x), discarding the y values
1229 template <typename OrdinalType, typename ScalarType>
x_y_pairs_to_x_set(const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & xy_pairs,std::set<ScalarType> & x_set)1230 void x_y_pairs_to_x_set(
1231 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& xy_pairs,
1232 std::set<ScalarType>& x_set)
1233 {
1234 x_set.clear();
1235 OrdinalType i, num_pairs = xy_pairs.length()/2;
1236 for (i=0; i<num_pairs; ++i)
1237 x_set.insert(xy_pairs[2*i]);
1238 }
1239
1240
1241 // ---------------------------------
1242 // templated array utility functions
1243 // ---------------------------------
1244
1245
1246 template <typename ScalarType>
find_min(const std::vector<ScalarType> & vec)1247 inline ScalarType find_min(const std::vector<ScalarType>& vec)
1248 {
1249 size_t i, len = vec.size();
1250 ScalarType min = (len) ? vec[0] : std::numeric_limits<ScalarType>::max();
1251 for (i=1; i<len; ++i)
1252 if (vec[i] < min)
1253 min = vec[i];
1254 return min;
1255 }
1256
1257
1258 template <typename ScalarType>
find_max(const std::vector<ScalarType> & vec)1259 inline ScalarType find_max(const std::vector<ScalarType>& vec)
1260 {
1261 size_t i, len = vec.size();
1262 ScalarType max = (len) ? vec[0] : std::numeric_limits<ScalarType>::min();
1263 for (i=1; i<len; ++i)
1264 if (vec[i] > max)
1265 max = vec[i];
1266 return max;
1267 }
1268
1269
1270 #if defined(_MSC_VER)
1271 // MSE: this may be too generic and could hide special cases:
1272 // can we rely on partial template specialization?
1273
1274 /// generic find_index (inactive)
1275 template <typename ContainerType>
find_index(const ContainerType & c,const typename ContainerType::value_type & search_data)1276 size_t find_index(const ContainerType& c,
1277 const typename ContainerType::value_type& search_data)
1278 {
1279 // should be more efficient than find() + distance()
1280 size_t cntr = 0;
1281 for(const typename ContainerType::value_type& entry : c) {
1282 if (entry == search_data)
1283 return cntr;
1284 else
1285 ++cntr;
1286 }
1287 return _NPOS;
1288 }
1289
1290
1291 ///// generic copy (inactive)
1292 //template <typename MultiArrayType, typename DakArrayType>
1293 //void copy_data(const MultiArrayType& ma, DakArrayType& da)
1294 //{
1295 // size_t size_ma = ma.size();
1296 // if (size_ma != da.size())
1297 // da.resize(size_ma);
1298 // for (size_t i=0; i<size_ma; ++i)
1299 // da[i] = ma[i];
1300 //}
1301
1302
1303 #else
1304 /// compute the index of an entry within a boost::multi_array
1305 template <typename T>
find_index(const boost::multi_array<T,1> & bma,const T & search_data)1306 size_t find_index(const boost::multi_array<T, 1>& bma, const T& search_data)
1307 {
1308 // should be more efficient than find() + distance()
1309 size_t i, len = bma.size();
1310 for (i=0; i<len; i++)
1311 if (bma[i] == search_data)
1312 return i;
1313 return _NPOS;
1314 }
1315
1316
1317 //template <typename T>
1318 //size_t find_index(boost::multi_array<T, 1>::const_array_view<1>::type bmav,
1319 // const T& search_data)
1320 // TO DO: difficulty with compilation of ::type --> work-around by enumeration
1321
1322
1323 /// compute the index of an entry within a boost::multi_array view
find_index(SizetMultiArrayConstView bmacv,size_t search_data)1324 inline size_t find_index(SizetMultiArrayConstView bmacv, size_t search_data)
1325 {
1326 // should be more efficient than find() + distance()
1327 size_t len = bmacv.size();
1328 for (size_t i=0; i<len; i++)
1329 if (bmacv[i] == search_data)
1330 return i;
1331 return _NPOS;
1332 }
1333
1334
1335 /// compute the index of an entry within a boost::multi_array view
find_index(StringMultiArrayConstView bmacv,const String & search_data)1336 inline size_t find_index(StringMultiArrayConstView bmacv,
1337 const String& search_data)
1338 {
1339 // should be more efficient than find() + distance()
1340 size_t len = bmacv.size();
1341 for (size_t i=0; i<len; i++)
1342 if (bmacv[i] == search_data)
1343 return i;
1344 return _NPOS;
1345 }
1346
1347
1348 /// compute the index of an entry within a std::list
1349 template <typename ListT>
find_index(const ListT & l,const typename ListT::value_type & val)1350 size_t find_index(const ListT& l, const typename ListT::value_type& val)
1351 {
1352 // should be more efficient than find() + distance()
1353 size_t cntr = 0;
1354 for(const typename ListT::value_type& entry : l) {
1355 if (entry == val)
1356 return cntr;
1357 else
1358 ++cntr;
1359 }
1360 return _NPOS;
1361 }
1362
1363
1364 // copy MultiArrayView<T> to Array<T>
1365 //template <typename T>
1366 //void copy_data(boost::multi_array<T, 1>::const_array_view<1>::type ma,
1367 // Array<T>& da)
1368 // TO DO: difficulty with compilation of ::type --> work-around by enumeration
1369
1370
1371 /// return an iterator to the first list element satisfying the
1372 /// predicate test_fn w.r.t. the passed test_fn_data; end if not found
1373 template <typename ListT>
1374 typename ListT::const_iterator
find_if(const ListT & c,bool (* test_fn)(const typename ListT::value_type &,const std::string &),const std::string & test_fn_data)1375 find_if(const ListT& c,
1376 bool (*test_fn)(const typename ListT::value_type&, const std::string&),
1377 const std::string& test_fn_data)
1378 {
1379 // a hand-coded list traversal is faster than the Functor approach
1380 for (typename ListT::const_iterator it = c.begin(); it != c.end(); ++it)
1381 if (test_fn(*it, test_fn_data))
1382 return it;
1383 return c.end();
1384
1385 // while only a single list traversal is needed with find_if() in this case,
1386 // it is still a fair amount slower, presumably due to Functor overhead.
1387 //return find_if(begin(), end(), FunctionCompare<T>(test_fn, test_fn_data));
1388 }
1389
1390 #endif
1391
1392 // copy std::vector<VecType> to Real*
1393 // VectorType must support the length() method.
1394 template<typename VectorType, typename ScalarType>
copy_data(const std::vector<VectorType> & va,ScalarType * ptr,int ptr_len)1395 void copy_data(const std::vector<VectorType>& va, ScalarType * ptr, int ptr_len)
1396 {
1397 size_t total_len=0, cntr=0, num_vec = va.size();
1398 for( size_t i=0; i<num_vec; ++i)
1399 total_len += va[i].length();
1400 if (total_len != ptr_len) {
1401 Cerr << "copy_data Error: pointer allocation (" << ptr_len << ") does not equal "
1402 << "total std::vector<VecType> length (" << total_len << ")." << std::endl;
1403 abort_handler(-1);
1404 }
1405 for( size_t i=0; i<num_vec; ++i)
1406 {
1407 int vec_len = va[i].length();
1408 for(int j=0; j<vec_len; ++j)
1409 ptr[cntr++] = va[i][j];
1410 }
1411 }
1412
1413 /// copy boost::multi_array view to Array - used by ActiveSet::derivative_vector - RWH
copy_data(SizetMultiArrayConstView ma,SizetArray & da)1414 inline void copy_data(SizetMultiArrayConstView ma, SizetArray& da)
1415 {
1416 size_t size_ma = ma.size();
1417 if (size_ma != da.size())
1418 da.resize(size_ma);
1419 for (size_t i=0; i<size_ma; ++i)
1420 da[i] = ma[i];
1421 }
1422
1423
1424 /// copy boost::multi_array view to Array - used by Pecos::copy_data - RWH
copy_data(StringMultiArrayConstView ma,StringArray & da)1425 inline void copy_data(StringMultiArrayConstView ma, StringArray& da)
1426 {
1427 size_t size_ma = ma.size();
1428 if (size_ma != da.size())
1429 da.resize(size_ma);
1430 for (size_t i=0; i<size_ma; ++i)
1431 da[i] = ma[i];
1432 }
1433
1434
1435 /// return true if the item val appears in container v
1436 template <typename DakContainerType>
contains(const DakContainerType & v,const typename DakContainerType::value_type & val)1437 inline bool contains(const DakContainerType& v,
1438 const typename DakContainerType::value_type& val)
1439 {
1440 return ( std::find(v.begin(), v.end(), val) != v.end() ) ? true : false;
1441 }
1442
1443
1444 } // namespace Dakota
1445
1446 /// return true if the string contains a floating point value
isfloat(const Dakota::String token)1447 inline bool isfloat(const Dakota::String token) {
1448 static boost::regex float_regex("[\\+-]?[0-9]*\\.?[0-9]+\\.?[0-9]*[eEdD]?[\\+-]?[0-9]*|[Nn][Aa][Nn]|[\\+-]?[Ii][Nn][Ff](?:[Ii][Nn][Ii][Tt][Yy])?");
1449 return boost::regex_match(token, float_regex);
1450 }
1451
1452
1453 #endif // DAKOTA_DATA_UTIL_H
1454