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 //- Class: SharedC3ApproxData
10 //- Description: Implementation code for SharedC3ApproxData class
11 //-
12 //- Owner: Mike Eldred
13
14 #include "SharedC3ApproxData.hpp"
15 #include "ProblemDescDB.hpp"
16 #include "NonDIntegration.hpp"
17
18 #include "pecos_stat_util.hpp"
19 #include "pecos_global_defs.hpp"
20
21 #include <assert.h>
22 //#define DEBUG
23
24 namespace Dakota {
25
26
27 SharedC3ApproxData::
SharedC3ApproxData(ProblemDescDB & problem_db,size_t num_vars)28 SharedC3ApproxData(ProblemDescDB& problem_db, size_t num_vars):
29 SharedApproxData(BaseConstructor(), problem_db, num_vars),
30 kickOrder(problem_db.get_ushort("model.c3function_train.kick_order")),
31 maxOrder(problem_db.get_ushort("model.c3function_train.max_order")),
32 adaptOrder(problem_db.get_bool("model.c3function_train.adapt_order")),
33 startRank(problem_db.get_sizet("model.c3function_train.start_rank")),
34 kickRank(problem_db.get_sizet("model.c3function_train.kick_rank")),
35 maxRank(problem_db.get_sizet("model.c3function_train.max_rank")),
36 adaptRank(problem_db.get_bool("model.c3function_train.adapt_rank")),
37 regressType(problem_db.get_short("model.surrogate.regression_type")),
38 regressRegParam(problem_db.get_real("model.surrogate.regression_penalty")),
39 solverTol(problem_db.get_real("model.c3function_train.solver_tolerance")),
40 solverRoundingTol(
41 problem_db.get_real("model.c3function_train.solver_rounding_tolerance")),
42 statsRoundingTol(
43 problem_db.get_real("model.c3function_train.stats_rounding_tolerance")),
44 maxSolverIterations(problem_db.get_int("model.max_solver_iterations")),
45 crossMaxIter(
46 problem_db.get_int("model.c3function_train.max_cross_iterations")),
47 //adaptConstruct(false),
48 c3AdvancementType(NO_C3_ADVANCEMENT)
49 {
50 // This ctor used for user-spec of DataFitSurrModel (surrogate global FT
51 // used by generic surrogate-based UQ in NonDSurrogateExpansion)
52
53 RealVector dim_pref_spec; // isotropic for now, prior to XML support
54 NonDIntegration::dimension_preference_to_anisotropic_order(
55 problem_db.get_ushort("model.c3function_train.start_order"),
56 dim_pref_spec,//problem_db.get_rv("model.dimension_preference"),
57 numVars, startOrders);
58
59 multiApproxOpts = multi_approx_opts_alloc(num_vars);
60 oneApproxOpts.assign(num_vars, NULL);
61 }
62
63
64 SharedC3ApproxData::
SharedC3ApproxData(const String & approx_type,const UShortArray & approx_order,size_t num_vars,short data_order,short output_level)65 SharedC3ApproxData(const String& approx_type, const UShortArray& approx_order,
66 size_t num_vars, short data_order, short output_level):
67 SharedApproxData(NoDBBaseConstructor(), approx_type, num_vars, data_order,
68 output_level),
69 // default values overridden by set_parameter
70 startOrders(approx_order), kickOrder(1), maxOrder(USHRT_MAX),
71 adaptOrder(false), startRank(2), kickRank(1),
72 maxRank(std::numeric_limits<size_t>::max()), adaptRank(false),
73 regressType(FT_LS), // non-regularized least sq
74 solverTol(1.e-10), solverRoundingTol(1.e-10), statsRoundingTol(1.e-10),
75 maxSolverIterations(-1), crossMaxIter(5), //adaptConstruct(false),
76 c3AdvancementType(NO_C3_ADVANCEMENT)
77 {
78 // This ctor used by lightweight/on-the-fly DataFitSurrModel ctor
79 // (used to build an FT on top of a user model in NonDC3FuntionTrain)
80
81 // short basis_type; approx_type_to_basis_type(approxType, basis_type);
82
83 multiApproxOpts = multi_approx_opts_alloc(num_vars);
84 oneApproxOpts.assign(num_vars, NULL);
85 }
86
87
~SharedC3ApproxData()88 SharedC3ApproxData::~SharedC3ApproxData()
89 {
90 multi_approx_opts_free(multiApproxOpts); //multiApproxOpts = NULL;
91 for (size_t i=0; i<numVars; ++i) {
92 struct OneApproxOpts*& a_opts = oneApproxOpts[i]; // ref to ptr
93 if (a_opts) one_approx_opts_free_deep(&a_opts); //a_opts = NULL;
94 }
95 }
96
97
98 void SharedC3ApproxData::
construct_basis(const Pecos::MultivariateDistribution & mv_dist)99 construct_basis(const Pecos::MultivariateDistribution& mv_dist)
100 {
101 const ShortArray& rv_types = mv_dist.random_variable_types();
102 const BitArray& active_vars = mv_dist.active_variables();
103 bool no_mask = active_vars.empty();
104 size_t i, av_cntr = 0, num_rv = rv_types.size(),
105 num_active_rv = (no_mask) ? num_rv : active_vars.count();
106 assert (num_active_rv == numVars);
107
108 struct OpeOpts * o_opts;
109 // uses startOrders at construct time (initialize_u_space_model()) since
110 // startOrdersMap[key] not meaningful until run-time active key assignment
111 const UShortArray& so = start_orders();
112 size_t np, max_np = max_order(); ++max_np;//no overflow (default is USHRT_MAX)
113 for (i=0; i<num_rv; ++i)
114 if (no_mask || active_vars[i]) {
115 switch (rv_types[i]) {
116 case Pecos::STD_NORMAL:
117 o_opts = ope_opts_alloc(HERMITE); break;
118 case Pecos::STD_UNIFORM:
119 o_opts = ope_opts_alloc(LEGENDRE); break;
120 default:
121 o_opts = NULL;
122 PCerr << "Error: unsupported RV type (" << rv_types[i] << ") in "
123 << "SharedC3ApproxData::distribution_parameters()" << std::endl;
124 abort_handler(-1); break;
125 }
126
127 np = std::min((size_t)so[i] + 1, max_np);
128 ope_opts_set_nparams(o_opts, np); // startnum = startord + 1
129 // Note: maxOrder not used for C3 regression (limits increment_order()
130 // on Dakota side); to be used for adaptation by cross-approximation
131 ope_opts_set_maxnum(o_opts, max_np); // maxnum = maxord + 1
132
133 struct OneApproxOpts*& a_opts = oneApproxOpts[av_cntr]; // ref to ptr
134 if (a_opts) one_approx_opts_free_deep(&a_opts); // a_opts frees o_opts
135 a_opts = one_approx_opts_alloc(POLYNOMIAL, o_opts);
136 multi_approx_opts_set_dim(multiApproxOpts, av_cntr, a_opts);
137
138 ++av_cntr;
139 }
140
141 // don't bother to assign formUpdated prior to an active key definition
142 if (!activeKey.empty())
143 formUpdated[activeKey] = true;
144 }
145
146
147 void SharedC3ApproxData::
update_basis(const UShortArray & start_orders,unsigned short max_order)148 update_basis(const UShortArray& start_orders, unsigned short max_order)
149 {
150 size_t np, max_np = max_order; ++max_np; // no overflow (default is USHRT_MAX)
151 for (size_t v=0; v<numVars; ++v) {
152 struct OneApproxOpts*& a_opts = oneApproxOpts[v];
153 np = (size_t)std::min(start_orders[v], max_order) + 1;
154 one_approx_opts_set_nparams(a_opts, np);
155 one_approx_opts_set_maxnum( a_opts, max_np);
156 }
157
158 formUpdated[activeKey] = true;
159 }
160
161
162 void SharedC3ApproxData::
update_basis(size_t v,unsigned short start_order,unsigned short max_order)163 update_basis(size_t v, unsigned short start_order, unsigned short max_order)
164 {
165 size_t max_np = max_order; ++max_np; // no overflow (default is USHRT_MAX)
166 size_t np = (size_t)std::min(start_order, max_order) + 1;
167
168 struct OneApproxOpts*& a_opts = oneApproxOpts[v];
169 one_approx_opts_set_nparams(a_opts, np);
170 one_approx_opts_set_maxnum( a_opts, max_np);
171
172 //formUpdated[activeKey] = true; // elevate to clients to avoid redundancy
173 }
174
175
advancement_available()176 bool SharedC3ApproxData::advancement_available()
177 {
178 // START_*_ADVANCEMENT cases are tested first. If false, then each C3Approx
179 // is tested for MAX_*_ADVANCEMENT cases (see ApproximationInterface::
180 // advancement_available()). This distributes the c3AdvancementType cases.
181
182 switch (c3AdvancementType) {
183
184 // these two options only require the shared data config:
185 case START_ORDER_ADVANCEMENT: {
186 const UShortArray& s_ord = start_orders();// adapted orders
187 unsigned short m_ord = max_order();
188 size_t i, num_ord = s_ord.size();
189 for (size_t i=0; i<num_ord; ++i)
190 if (s_ord[i] < m_ord)
191 return true;
192 return false; break;
193 }
194 case START_RANK_ADVANCEMENT:
195 return (start_rank() < max_rank()); break;
196
197 // MAX_*_ADVANCEMENT refine types are QoI-dependent: shared false induces
198 // per-QoI checks (see ApproximationInterface::advancement_available()
199 default:
200 // clear state prior to accumulation across C3Approximations
201 c3MaxRankAdvance[activeKey] = c3MaxOrderAdvance[activeKey] = false;
202 return false; break;
203 }
204 }
205
206
increment_order()207 void SharedC3ApproxData::increment_order()
208 {
209 switch (c3AdvancementType) {
210 case START_ORDER_ADVANCEMENT: {
211 UShortArray& start_ord = start_orders();
212 unsigned short max_ord = max_order(); // default is USHRT_MAX
213 bool incremented = false;
214 for (size_t v=0; v<numVars; ++v) {
215 unsigned short &s_ord = start_ord[v], prev_ord = s_ord;
216 // unconditional increment (preserve reproducibility w/decrement)
217 // Note: dim_pref ratios are not preserved
218 s_ord += kickOrder;
219 if (prev_ord < max_ord) // increment occurs, but kick may be truncated
220 incremented = true;
221 }
222 if (incremented)
223 formUpdated[activeKey] = true;
224 else
225 Cerr << "Warning: SharedC3ApproxData::increment_order() cannot advance "
226 << "order beyond maximum allowable (" << max_ord << ") for key:\n"
227 << activeKey << std::endl;
228 break;
229 }
230 case START_RANK_ADVANCEMENT: {
231 // To ensure symmetry with decrement, don't saturate at maxRank
232 // > Must bound start_ranks vector in C3Approximation::build()
233 // > build() truncates kick to max: start_ranks = std::min(start_r, max_r)
234 size_t &start_r = start_rank();
235 if (start_r < max_rank()) // increment occurs but kick might be truncated
236 formUpdated[activeKey] = true;
237 start_r += kickRank; // invertible in decrement_order()
238 break;
239 }
240 case MAX_RANK_ADVANCEMENT: {
241 size_t& max_r = max_rank(); max_r += kickRank;
242 formUpdated[activeKey] = true; break;
243 }
244 case MAX_ORDER_ADVANCEMENT: {
245 unsigned short& max_o = max_order(); max_o += kickOrder;
246 formUpdated[activeKey] = true; break;
247 }
248 case MAX_RANK_ORDER_ADVANCEMENT: {
249 // Prior to implementing a multi-index approach, we use heuristics.
250 // > unconditionally advancing both is wasteful; only advance non-saturated
251 // > could also consider only advancing one when both bounds are active:
252 // least saturated first with tie break to max rank (recovered ranks are
253 // heterogeneous anyway)
254 bool r_advance = c3MaxRankAdvance[activeKey],
255 o_advance = c3MaxOrderAdvance[activeKey];
256 if (r_advance)
257 { size_t& max_r = max_rank(); max_r += kickRank; }
258 if (o_advance)
259 { unsigned short& max_o = max_order(); max_o += kickOrder; }
260 if (r_advance || o_advance) formUpdated[activeKey] = true;
261 break;
262 }
263 }
264 }
265
266
decrement_order()267 void SharedC3ApproxData::decrement_order()
268 {
269 bool bad_range = false;
270 switch (c3AdvancementType) {
271 case START_ORDER_ADVANCEMENT: {
272 // always decrement and only update_basis() if within bounds
273 // (preserve symmetry/reproducibility w/ increment)
274 UShortArray& start_ord = start_orders();
275 unsigned short max_ord = max_order(); // default is USHRT_MAX
276 bool decremented = false;
277 for (size_t v=0; v<numVars; ++v) {
278 unsigned short &s_ord = start_ord[v];
279 if (s_ord >= kickOrder) { // for completeness (should not happen)
280 s_ord -= kickOrder; // preserve symmetry/reproducibility w/increment
281 if (s_ord < max_ord) // only communicate if in bounds
282 decremented = true;
283 }
284 }
285 if (decremented) formUpdated[activeKey] = true;
286 else bad_range = true;
287 break;
288 }
289 case START_RANK_ADVANCEMENT: {
290 size_t& start_r = start_rank();
291 if (start_r < kickRank) bad_range = true;
292 else { start_r -= kickRank; formUpdated[activeKey] = true; }
293 break;
294 }
295 case MAX_RANK_ADVANCEMENT: {
296 size_t& max_r = max_rank();
297 if (max_r < kickRank) bad_range = true;
298 else { max_r -= kickRank; formUpdated[activeKey] = true; }
299 break;
300 }
301 case MAX_ORDER_ADVANCEMENT: {
302 unsigned short& max_o = max_order();
303 if (max_o < kickOrder) bad_range = true;
304 else { max_o -= kickOrder; formUpdated[activeKey] = true; }
305 break;
306 }
307 case MAX_RANK_ORDER_ADVANCEMENT: {
308 bool r_advance = c3MaxRankAdvance[activeKey],
309 o_advance = c3MaxOrderAdvance[activeKey];
310 if (r_advance) {
311 size_t& max_r = max_rank();
312 if (max_r < kickRank) bad_range = true;
313 else max_r -= kickRank;
314 }
315 if (o_advance) {
316 unsigned short& max_o = max_order();
317 if (max_o < kickOrder) bad_range = true;
318 else max_o -= kickOrder;
319 }
320 if (r_advance || o_advance) formUpdated[activeKey] = true;
321 break;
322 }
323 }
324
325 if (bad_range)
326 Cerr << "Warning: SharedC3ApproxData::decrement_order() outside of valid "
327 << "range for key:\n" << activeKey << std::endl;
328 }
329
330
pre_combine()331 void SharedC3ApproxData::pre_combine()
332 {
333 combinedOrders.assign(numVars, 0);
334 std::map<UShortArray, UShortArray>::const_iterator cit; size_t i;
335 for (cit=startOrdersMap.begin(); cit!=startOrdersMap.end(); ++cit) {
336 const UShortArray& so = cit->second;
337 for (i=0; i<numVars; ++i)
338 if (so[i] > combinedOrders[i])
339 combinedOrders[i] = so[i];
340 }
341
342 update_basis(combinedOrders, max_order());
343 }
344
345
346 //void SharedC3ApproxData::post_combine()
347 //{ update_basis(); } // restore to active
348
349 } // namespace Dakota
350