diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index 07fd013d..374afa38 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -1683,179 +1683,6 @@ void MmUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm UpgradeThread::ensemble_solution(iter,verbose_level,maxsing,thread_id,t_count, use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid,par_diff,Am,obs_resid, obs_diff,upgrade_1,obs_err,weights,parcov_inv,empty_obs_names,empty_par_names); -// Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, d_dash; -// -// //---------------------------------- -// //es-mda solution -// //---------------------------------- -// -// if (!use_glm_form) -// { -// if (true) -// { -// // Low rank Cee. Section 14.3.2 Evenson Book -// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); -// obs_err = sqrt(cur_lam) * obs_err; -// Eigen::MatrixXd s0, V0, U0, s0_i; -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); -// s0_i = s0.asDiagonal().inverse(); -// Eigen::MatrixXd X0 = U0.transpose() * obs_err; -// X0 = s0_i * X0; -// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; -// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); -// -// s1_2 = s1.cwiseProduct(s1); -// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); -// Eigen::MatrixXd X1 = s0_i * U1; -// X1 = U0 * X1; -// -// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); -// Eigen::MatrixXd X2 = X4 * obs_resid; -// Eigen::MatrixXd X3 = X1 * X2; -// -// X3 = obs_diff.transpose() * X3; -// upgrade_1 = -1 * par_diff * X3; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// upgrade_1.transposeInPlace(); -// -// } -// else -// { -// // Use when ensemble size is larger than number of observation. -// // This is a good option when number of observation is small -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); -// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); -// -// Eigen::MatrixXd s2_, s, V, U, cum_sum; -// SVD_REDSVD rsvd; -// Eigen::MatrixXd C; -// C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); -// Eigen::VectorXd s2; -// -// -// rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); -// s2 = s.cwiseProduct(s); -// s2_ = s2.asDiagonal().inverse(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); -// -// -// Eigen::MatrixXd X1 = s2_ * U.transpose(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); -// -// X1 = X1 * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); -// -// X1 = U * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); -// -// X1 = obs_diff.transpose() * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); -// -// upgrade_1 = -1 * par_diff * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// upgrade_1.transposeInPlace(); -// } -// -// -// } -// -// -//// ---------------------------------- -//// glm solution -//// ---------------------------------- -// else -// { -// -// obs_resid = weights * obs_resid; -// obs_diff = scale * (weights * obs_diff); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); -// if (use_prior_scaling) -// par_diff = scale * parcov_inv * par_diff; -// else -// par_diff = scale * par_diff; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); -// -// Ut.transposeInPlace(); -// obs_diff.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); -// -// Eigen::MatrixXd s2 = s.cwiseProduct(s); -// -// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); -// -// Eigen::MatrixXd X1 = Ut * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); -// -// Eigen::MatrixXd X2 = ivec * X1; -// X1.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); -// -// Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; -// X2.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); -// upgrade_1 = -1.0 * par_diff * X3; -// -// if (use_prior_scaling) -// { -// //upgrade_1 = parcov_inv * upgrade_1; -// } -// -// upgrade_1.transposeInPlace(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// X3.resize(0, 0); -// -// Eigen::MatrixXd upgrade_2; -// if ((!use_approx) && (iter > 1)) -// { -// if (use_prior_scaling) -// { -// par_resid = parcov_inv * par_resid; -// } -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); -// Eigen::MatrixXd x4 = Am.transpose() * par_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); -// -// par_resid.resize(0, 0); -// -// Eigen::MatrixXd x5 = Am * x4; -// x4.resize(0, 0); -// //Am.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); -// Eigen::MatrixXd x6 = par_diff.transpose() * x5; -// x5.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); -// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; -// x6.resize(0, 0); -// -// if (use_prior_scaling) -// { -// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; -// } -// else -// { -// upgrade_2 = -1.0 * (par_diff * x7); -// } -// x7.resize(0, 0); -// -// upgrade_1 = upgrade_1 + upgrade_2.transpose(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); -// upgrade_2.resize(0, 0); -// -// } -// } //assuming that the fist row is the realization we are after... row_vec = upgrade_1.row(0); @@ -1900,493 +1727,6 @@ MmUpgradeThread::MmUpgradeThread(PerformanceLog* _performance_log, unordered_map } -//void CovLocalizationUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector par_names, -// vector obs_names) -//{ -// -// //declare these helpers in here so they are thread safe... -// class local_utils -// { -// public: -// -// static void save_names(int verbose_level, int tid, int iter, int t_count, string prefix, vector& names) -// { -// if (verbose_level < 2) -// return; -// -// if (verbose_level < 3) -// return; -// stringstream ss; -// -// ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; -// string fname = ss.str(); -// ofstream f(fname); -// if (!f.good()) -// cout << "error getting ofstream " << fname << endl; -// else -// { -// for (const auto& name : names) -// { -// f << name << endl; -// } -// -// } -// f.close(); -// return; -// } -// static Eigen::DiagonalMatrix get_matrix_from_map(vector& names, unordered_map& dmap) -// { -// Eigen::VectorXd vec(names.size()); -// int i = 0; -// for (auto name : names) -// { -// vec[i] = dmap.at(name); -// ++i; -// } -// Eigen::DiagonalMatrix m = vec.asDiagonal(); -// return m; -// } -// static Eigen::MatrixXd get_matrix_from_map(int num_reals, vector& names, unordered_map& emap) -// { -// Eigen::MatrixXd mat(num_reals, names.size()); -// mat.setZero(); -// -// for (int j = 0; j < names.size(); j++) -// { -// mat.col(j) = emap[names[j]]; -// } -// -// return mat; -// } -// static void save_mat(int verbose_level, int tid, int iter, int t_count, string prefix, Eigen::MatrixXd& mat) -// { -// if (verbose_level < 2) -// return; -// -// if (verbose_level < 3) -// return; -// //cout << "thread: " << tid << ", " << t_count << ", " << prefix << " rows:cols" << mat.rows() << ":" << mat.cols() << endl; -// stringstream ss; -// -// ss << "thread_" << tid << ".count_ " << t_count << ".iter_" << iter << "." << prefix << ".dat"; -// string fname = ss.str(); -// ofstream f(fname); -// if (!f.good()) -// cout << "error getting ofstream " << fname << endl; -// else -// { -// -// try -// { -// f << mat << endl; -// f.close(); -// } -// catch (...) -// { -// cout << "error saving matrix " << fname << endl; -// } -// } -// } -// }; -// -// stringstream ss; -// -// unique_lock ctrl_guard(ctrl_lock, defer_lock); -// int maxsing, num_reals, verbose_level, pcount = 0, t_count=0; -// double eigthresh; -// bool use_approx; -// bool use_prior_scaling; -// bool use_localizer = false; -// bool loc_by_obs = true; -// -// while (true) -// { -// if (ctrl_guard.try_lock()) -// { -// -// maxsing = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().maxsing; -// eigthresh = pe_upgrade.get_pest_scenario_ptr()->get_svd_info().eigthresh; -// use_approx = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_approx(); -// use_prior_scaling = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_use_prior_scaling(); -// num_reals = pe_upgrade.shape().first; -// verbose_level = pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_verbose_level(); -// ctrl_guard.unlock(); -// //if (pe_upgrade.get_pest_scenario_ptr()->get_pestpp_options().get_ies_localize_how()[0] == 'P') -// if (how == Localizer::How::PARAMETERS) -// loc_by_obs = false; -// else -// { -// throw runtime_error("Covariance localization only supporte for localization by parameters..."); -// } -// break; -// } -// } -// ofstream f_thread; -// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_obs_names",obs_names); -// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_par_names",par_names); -// if (verbose_level > 2) -// { -// ss.str(""); -// ss << "thread_" << thread_id << ".part_map.csv"; -// f_thread.open(ss.str()); -// ss.str(""); -// } -// Eigen::MatrixXd par_resid, par_diff, Am; -// Eigen::MatrixXd obs_resid, obs_diff, obs_err, loc; -// Eigen::DiagonalMatrix weights, parcov_inv; -// vector case_par_names, case_obs_names; -// string key; -// -// -// //these locks are used to control (thread-safe) access to the fast look up containers -// unique_lock next_guard(next_lock, defer_lock); -// unique_lock obs_diff_guard(obs_diff_lock, defer_lock); -// unique_lock obs_resid_guard(obs_resid_lock, defer_lock); -// unique_lock obs_err_guard(obs_err_lock, defer_lock); -// unique_lock par_diff_guard(par_diff_lock, defer_lock); -// unique_lock par_resid_guard(par_resid_lock, defer_lock); -// unique_lock loc_guard(loc_lock, defer_lock); -// unique_lock weight_guard(weight_lock, defer_lock); -// unique_lock parcov_guard(parcov_lock, defer_lock); -// unique_lock am_guard(am_lock, defer_lock); -// -// //reset all the solution parts -// par_resid.resize(0, 0); -// par_diff.resize(0, 0); -// obs_resid.resize(0, 0); -// obs_diff.resize(0, 0); -// obs_err.resize(0, 0); -// loc.resize(0, 0); -// Am.resize(0, 0); -// weights.resize(0); -// parcov_inv.resize(0); -// Am.resize(0, 0); -// -// -// -// obs_diff = local_utils::get_matrix_from_map(num_reals, obs_names, obs_diff_map); -// obs_resid = local_utils::get_matrix_from_map(num_reals, obs_names, obs_resid_map); -// obs_err = local_utils::get_matrix_from_map(num_reals, obs_names, obs_err_map); -// par_diff = local_utils::get_matrix_from_map(num_reals, par_names, par_diff_map); -// par_resid = local_utils::get_matrix_from_map(num_reals, par_names, par_resid_map); -// weights = local_utils::get_matrix_from_map(obs_names, weight_map); -// parcov_inv = local_utils::get_matrix_from_map(par_names, parcov_inv_map); -// -// if ((!use_approx) && (Am.rows() == 0)) -// { -// //Am = local_utils::get_matrix_from_map(num_reals, par_names, Am_map).transpose(); -// int am_cols = Am_map[par_names[0]].size(); -// Am.resize(par_names.size(), am_cols); -// Am.setZero(); -// -// for (int j = 0; j < par_names.size(); j++) -// { -// Am.row(j) = Am_map[par_names[j]]; -// } -// } -// -// par_diff.transposeInPlace(); -// obs_diff.transposeInPlace(); -// obs_resid.transposeInPlace(); -// par_resid.transposeInPlace(); -// obs_err.transposeInPlace(); -// -// -// //container to quickly look indices -// map par2col_map; -// for (int i = 0; i < par_names.size(); i++) -// par2col_map[par_names[i]] = i; -// -// map obs2row_map; -// for (int i = 0; i < obs_names.size(); i++) -// obs2row_map[obs_names[i]] = i; -// -// -// //form the scaled obs resid matrix -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); -// //Eigen::MatrixXd scaled_residual = weights * obs_resid; -// -// //form the (optionally) scaled par resid matrix -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); -// -// ss.str(""); -// double scale = (1.0 / (sqrt(double(num_reals - 1)))); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); -// -// //calculate some full solution components first... -// -// Eigen::MatrixXd ivec, upgrade_1, s, s2, V, Ut, t; -// Eigen::MatrixXd upgrade_2; -// Eigen::VectorXd loc_vec; -// -// upgrade_2.resize(num_reals, pe_upgrade.shape().second); -// upgrade_2.setZero(); -// -// if (!use_glm_form) -// { -// if (true) -// { -// // Low rank Cee. Section 14.3.2 Evenson Book -// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); -// obs_err = sqrt(cur_lam) * obs_err; -// Eigen::MatrixXd s0, V0, U0, s0_i; -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); -// s0_i = s0.asDiagonal().inverse(); -// Eigen::MatrixXd X0 = U0.transpose() * obs_err; -// X0 = s0_i * X0; -// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; -// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); -// -// s1_2 = s1.cwiseProduct(s1); -// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); -// Eigen::MatrixXd X1 = s0_i * U1; -// X1 = U0 * X1; -// -// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); -// Eigen::MatrixXd X2 = X4 * obs_resid; -// Eigen::MatrixXd X3 = X1 * X2; -// -// X3 = obs_diff.transpose() * X3; -// //upgrade_1 = -1 * par_diff * X3; -// //local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// //upgrade_1.transposeInPlace(); -// t = obs_diff * X3; -// Ut.resize(0, 0); -// obs_diff.resize(0, 0); -// X3.resize(0,0); -// X2.resize(0,0); -// X4.resize(0,0); -// X1.resize(0,0); -// -// } -// else { -// obs_diff = scale * obs_diff;// (H-Hm)/sqrt(N-1) -// par_diff = scale * par_diff;// (K-Km)/sqrt(N-1) -// obs_err = scale * obs_err; // (E-Em)/sqrt(N-1) -// -// SVD_REDSVD rsvd; -// Eigen::MatrixXd C = obs_diff + (cur_lam * obs_err); // curr_lam is the inflation factor -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); -// rsvd.solve_ip(C, s, Ut, V, eigthresh, maxsing); -// Ut.transposeInPlace(); -// V.resize(0, 0); -// C.resize(0, 0); -// obs_err.resize(0, 0); -// -// // s2 = s.asDiagonal().inverse(); -// s2 = s.cwiseProduct(s).asDiagonal().inverse(); -// for (int i = 0; i < s.size(); i++) { -// if (s(i) < 1e-50) { -// s2(i, i) = 0; -// } -// } -// -// t = obs_diff.transpose() * Ut.transpose() * s2 * Ut; -// Ut.resize(0, 0); -// obs_diff.resize(0, 0); -// } -// } -// -// -// //---------------------------------- -// //glm solution -// //---------------------------------- -// else -// { -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); -// obs_diff = scale * (weights * obs_diff); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_diff", obs_diff); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); -// if (verbose_level > 1) { -// if (parcov_inv.size() < 10000) { -// Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); -// } -// ss.str(""); -// ss << "solution scaling factor: " << scale; -// cout << ss.str() << endl; -// } -// if (use_prior_scaling) -// par_diff = scale * parcov_inv * par_diff; -// else -// par_diff = scale * par_diff; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); -// -// Ut.transposeInPlace(); -// obs_diff.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); -// -// Eigen::MatrixXd s2 = s.cwiseProduct(s); -// -// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); -// -// obs_resid = weights * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_resid", obs_resid); -// t = V * s.asDiagonal() * ivec * Ut; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "t", t); -// -// if ((!use_approx) && (iter > 1)) -// { -// if (use_prior_scaling) -// { -// par_resid = parcov_inv * par_resid; -// } -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); -// Eigen::MatrixXd x4 = Am.transpose() * par_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); -// -// par_resid.resize(0, 0); -// -// Eigen::MatrixXd x5 = Am * x4; -// x4.resize(0, 0); -// Am.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); -// Eigen::MatrixXd x6 = par_diff.transpose() * x5; -// x5.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); -// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; -// x6.resize(0, 0); -// -// if (use_prior_scaling) -// { -// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; -// } -// else -// { -// upgrade_2 = -1.0 * (par_diff * x7); -// } -// x7.resize(0, 0); -// -// //add upgrade_2 piece -// //upgrade_1 = upgrade_1 + upgrade_2.transpose(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); -// //upgrade_2.resize(0, 0); -// } -// } -// -// //This is the main thread loop - it continues until all upgrade pieces have been completed -// map loc_map; -// //solve for each case par_name -// Eigen::VectorXd pt = parcov_inv.diagonal(); -// string name; -// while (true) -// { -// //clear the pieces used last time... -// key = ""; -// loc_map.clear(); -// -// //get new pieces... -// while (true) -// { -// if ((key != "") && (loc_map.size() > 0)) -// break; -// if (next_guard.try_lock()) -// { -// //if all the pieces have been completed, return -// if (count == keys.size()) -// { -// if (verbose_level > 3) -// { -// cout << "upgrade thread: " << thread_id << " processed " << pcount << " upgrade parts" << endl; -// } -// if (f_thread.good()) -// f_thread.close(); -// return; -// } -// key = keys[count]; -// pair, vector> p = cases.at(key); -// //we can potentially optimize the speed of this loc type by changing how many pars are -// //passed in each "case" so herein, we support a generic number of pars per case... -// case_par_names = p.second; -// //In this solution, we ignore case obs names since we are using the full set of obs for the solution... -// case_obs_names = p.first; -// -// if (count % 1000 == 0) -// { -// ss.str(""); -// ss << "upgrade thread progress: " << count << " of " << total << " parts done"; -// if (verbose_level > 3) -// cout << ss.str() << endl; -// performance_log->log_event(ss.str()); -// } -// count++; -// t_count = count; -// pcount++; -// next_guard.unlock(); -// -// } -// //get access to the localizer -// if ((key != "") && (loc_map.size() == 0) && (loc_guard.try_lock())) -// { -// //get the nobs-length localizing vector for each case par name -// for (auto& par_name : case_par_names) -// { -// loc_map[par_name] = localizer.get_obs_hadamard_vector(par_name, obs_names); -// } -// loc_guard.unlock(); -// } -// } -// -// if (verbose_level > 2) -// { -// f_thread << t_count << "," << iter; -// for (auto name : case_par_names) -// f_thread << "," << name; -// for (auto name : case_obs_names) -// f_thread << "," << name; -// f_thread << endl; -// } -// upgrade_1.resize(num_reals,case_par_names.size()); -// upgrade_1.setZero(); -// -// for (int i = 0; i < case_par_names.size(); i++) -// { -// name = case_par_names[i]; -// loc_vec = loc_map[name]; -// Eigen::VectorXd par_vec = par_diff.row(par2col_map[name]) * t; -// //apply the localizer -// par_vec = par_vec.cwiseProduct(loc_vec); -// if (!use_glm_form) -// par_vec = -1.0 * par_vec.transpose() * obs_resid; -// else -// { -// par_vec = -1.0 * par_vec.transpose() * obs_resid; -// //if (use_prior_scaling) -// // par_vec *= pt[i]; -// } -// upgrade_1.col(i) += par_vec.transpose(); -// //add the par change part for the full glm solution -// if ((use_glm_form) && (!use_approx) && (iter > 1)) -// { -// //update_2 is transposed relative to upgrade_1 -// upgrade_1.col(i) += upgrade_2.row(par2col_map[name]); -// } -// -// } -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// -// //put this piece of the upgrade vector in -// unique_lock put_guard(put_lock, defer_lock); -// while (true) -// { -// if (put_guard.try_lock()) -// { -// pe_upgrade.add_2_cols_ip(case_par_names, upgrade_1); -// put_guard.unlock(); -// break; -// } -// } -// } -//} void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, bool use_glm_form, @@ -2590,19 +1930,6 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b } } -// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_obs_names",obs_names); -// local_utils::save_names(verbose_level,thread_id, iter, t_count,"act_par_names",par_names); -// -// if (verbose_level > 2) -// { -// -// f_thread << t_count << "," << iter; -// for (auto name : par_names) -// f_thread << "," << name; -// for (auto name : obs_names) -// f_thread << "," << name; -// f_thread << endl; -// } //reset all the solution parts par_resid.resize(0, 0); @@ -2673,227 +2000,6 @@ void LocalAnalysisUpgradeThread::work(int thread_id, int iter, double cur_lam, b parcov_inv, obs_names,par_names); - -// par_diff.transposeInPlace(); -// obs_diff.transposeInPlace(); -// obs_resid.transposeInPlace(); -// par_resid.transposeInPlace(); -// obs_err.transposeInPlace(); -// -// //form the scaled obs resid matrix -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_resid", obs_resid); -// //Eigen::MatrixXd scaled_residual = weights * obs_resid; -// -// //form the (optionally) scaled par resid matrix -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_resid", par_resid); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); -// -// stringstream ss; -// -// double scale = (1.0 / (sqrt(double(num_reals - 1)))); -// -// if (verbose_level > 1) { -// ss.str(""); -// ss << "solution scaling factor: " << scale; -// cout << ss.str() << endl; -// if (parcov_inv.size() < 10000) -// { -// Eigen::MatrixXd temp = parcov_inv.toDenseMatrix(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); -// } -// } -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff", obs_diff); -// //apply the localizer here... -// if (use_localizer) -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "loc", loc); -// if (use_localizer) -// { -// if (loc_by_obs) -// par_diff = par_diff.cwiseProduct(loc); -// else -// obs_diff = obs_diff.cwiseProduct(loc); -// } -// -// //Eigen::MatrixXd upgrade_1; -//// ensemble_solution(iter,verbose_level,maxsing,thread_id,t_count, use_prior_scaling,use_approx,use_glm_form,cur_lam,eigthresh,par_resid,par_diff,Am,obs_resid, -//// obs_diff,upgrade_1,obs_err,weights,parcov_inv); -// -// -// -// //---------------------------------- -// //es-mda solution -// //---------------------------------- -// -// if (!use_glm_form) -// { -// if (true) -// { -// // Low rank Cee. Section 14.3.2 Evenson Book -// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); -// obs_err = sqrt(cur_lam) * obs_err; -// Eigen::MatrixXd s0, V0, U0, s0_i; -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s0, U0, V0, eigthresh, maxsing); -// s0_i = s0.asDiagonal().inverse(); -// Eigen::MatrixXd X0 = U0.transpose() * obs_err; -// X0 = s0_i * X0; -// Eigen::MatrixXd s1, V1, U1, s1_2, s1_2i; -// rsvd.solve_ip(X0, s1, U1, V1, 0, maxsing); -// -// s1_2 = s1.cwiseProduct(s1); -// s1_2i = (Eigen::VectorXd::Ones(s1_2.size()) + s1_2).asDiagonal().inverse(); -// Eigen::MatrixXd X1 = s0_i * U1; -// X1 = U0 * X1; -// -// Eigen::MatrixXd X4 = s1_2i * X1.transpose(); -// Eigen::MatrixXd X2 = X4 * obs_resid; -// Eigen::MatrixXd X3 = X1 * X2; -// -// X3 = obs_diff.transpose() * X3; -// upgrade_1 = -1 * par_diff * X3; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// upgrade_1.transposeInPlace(); -// -// } -// else -// { -// // Use when ensemble size is larger than number of observation. -// // This is a good option when number of observation is small -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_err", obs_err); -// obs_err = obs_err.colwise() - obs_err.rowwise().mean(); -// -// Eigen::MatrixXd s2_, s, V, U, cum_sum; -// SVD_REDSVD rsvd; -// Eigen::MatrixXd C; -// C = obs_diff + (sqrt(cur_lam) * obs_err); // curr_lam is the inflation factor -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "C", C); -// Eigen::VectorXd s2; -// -// -// rsvd.solve_ip(C, s, U, V, eigthresh, maxsing); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U", U); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); -// s2 = s.cwiseProduct(s); -// s2_ = s2.asDiagonal().inverse(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "inv_s2_", s2_); -// -// -// Eigen::MatrixXd X1 = s2_ * U.transpose(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); -// -// X1 = X1 * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1_obs_resid", X1); -// -// X1 = U * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "U_X1", X1); -// -// X1 = obs_diff.transpose() * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "obs_diff_X1", X1); -// -// upgrade_1 = -1 * par_diff * X1; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// upgrade_1.transposeInPlace(); -// } -// -// -// } -// -// -//// ---------------------------------- -//// glm solution -//// ---------------------------------- -// else -// { -// -// obs_resid = weights * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_resid", par_diff); -// obs_diff = scale * (weights * obs_diff); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_obs_diff", par_diff); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "par_diff", par_diff); -// if (use_prior_scaling) -// par_diff = scale * parcov_inv * par_diff; -// else -// par_diff = scale * par_diff; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "scaled_par_diff", par_diff); -// SVD_REDSVD rsvd; -// rsvd.solve_ip(obs_diff, s, Ut, V, eigthresh, maxsing); -// -// Ut.transposeInPlace(); -// obs_diff.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Ut", Ut); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "s", s); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "V", V); -// -// Eigen::MatrixXd s2 = s.cwiseProduct(s); -// -// ivec = ((Eigen::VectorXd::Ones(s2.size()) * (cur_lam + 1.0)) + s2).asDiagonal().inverse(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "ivec", ivec); -// -// Eigen::MatrixXd X1 = Ut * obs_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X1", X1); -// -// Eigen::MatrixXd X2 = ivec * X1; -// X1.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X2", X2); -// -// Eigen::MatrixXd X3 = V * s.asDiagonal() * X2; -// X2.resize(0, 0); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X3", X3); -// upgrade_1 = -1.0 * par_diff * X3; -// -// if (use_prior_scaling) -// { -// //upgrade_1 = parcov_inv * upgrade_1; -// } -// -// upgrade_1.transposeInPlace(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_1", upgrade_1); -// X3.resize(0, 0); -// -// Eigen::MatrixXd upgrade_2; -// if ((!use_approx) && (iter > 1)) -// { -// if (use_prior_scaling) -// { -// par_resid = parcov_inv * par_resid; -// } -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "Am", Am); -// Eigen::MatrixXd x4 = Am.transpose() * par_resid; -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X4", x4); -// -// par_resid.resize(0, 0); -// -// Eigen::MatrixXd x5 = Am * x4; -// x4.resize(0, 0); -// Am.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X5", x5); -// Eigen::MatrixXd x6 = par_diff.transpose() * x5; -// x5.resize(0, 0); -// -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "X6", x6); -// Eigen::MatrixXd x7 = V * ivec * V.transpose() * x6; -// x6.resize(0, 0); -// -// if (use_prior_scaling) -// { -// upgrade_2 = -1.0 * parcov_inv * par_diff * x7; -// } -// else -// { -// upgrade_2 = -1.0 * (par_diff * x7); -// } -// x7.resize(0, 0); -// -// upgrade_1 = upgrade_1 + upgrade_2.transpose(); -// local_utils::save_mat(verbose_level, thread_id, iter, t_count, "upgrade_2", upgrade_2); -// upgrade_2.resize(0, 0); -// -// } -// } while (true) { @@ -3710,62 +2816,6 @@ void L2PhiHandler::report(bool echo, bool group_report) -//void PhiHandler::report(bool echo) -//{ -// ofstream &f = file_manager->rec_ofstream(); -// f << get_summary_header(); -// if (echo) -// cout << get_summary_header(); -// string s = get_summary_string(PhiHandler::phiType::COMPOSITE); -// f << s; -// if (echo) -// cout << s; -// s = get_summary_string(PhiHandler::phiType::MEAS); -// f << s; -// if (echo) -// cout << s; -// if (org_reg_factor != 0.0) -// { -// s = get_summary_string(PhiHandler::phiType::REGUL); -// f << s; -// if (echo) -// cout << s; -// } -// s = get_summary_string(PhiHandler::phiType::ACTUAL); -// f << s; -// if (echo) -// cout << s; -// if (org_reg_factor != 0.0) -// { -// if (*reg_factor == 0.0) -// { -// f << " (note: reg_factor is zero; regularization phi reported but not used)" << endl; -// if (echo) -// cout << " (note: reg_factor is zero; regularization phi reported but not used)" << endl; -// } -// else -// { -// f << " current reg_factor: " << *reg_factor << endl; -// if (echo) -// cout << " current reg_factor: " << *reg_factor << endl; -// } -// if (*reg_factor != 0.0) -// { -// -// f << " note: regularization phi reported above does not " << endl; -// f << " include the effects of reg_factor, " << endl; -// f << " but composite phi does." << endl; -// if (echo) -// { -// cout << " note: regularization phi reported above does not " << endl; -// cout << " include the effects of reg_factor, " << endl; -// cout << " but composite phi does." << endl; -// } -// } -// } -// f << endl << endl; -// f.flush(); -//} void L2PhiHandler::write(int iter_num, int total_runs, bool write_group) { diff --git a/src/libs/pestpp_common/MOEA.cpp b/src/libs/pestpp_common/MOEA.cpp index 9db7a5b2..61847e5b 100644 --- a/src/libs/pestpp_common/MOEA.cpp +++ b/src/libs/pestpp_common/MOEA.cpp @@ -3389,17 +3389,20 @@ vector MOEA::get_pso_gbest_solutions(int num_reals, ParameterEnsemble& _ DomPair dompair = objectives.get_nsga2_pareto_dominance(-999, _op, _dp, &constraints, false); vector nondom_solutions = dompair.first; vector gbest_solutions; + int num_objs = pi_obj_names.size()+obs_obj_names.size(); + //if no non dom solutions, then use the dominated ones... if (nondom_solutions.size() == 0) { ss.str(""); - ss << "WARNING: no nondom solutions for pst gbest calculation, using dominated solutions" << endl; + ss << "WARNING: no nondom solutions for pso gbest calculation, using dominated solutions" << endl; nondom_solutions = dompair.second; } - else if (nondom_solutions.size() == 1) + //todo: should we warn for nondom > 1 and objs == 1? + else if ((nondom_solutions.size() == 1) && (num_objs > 1)) { ss.str(""); - ss << "WARNING: only one nondom solution for pst gbest calculation" << endl; + ss << "WARNING: only one nondom solution for pso gbest calculation" << endl; file_manager.rec_ofstream() << ss.str(); cout << ss.str(); for (int i = 0; i < num_reals; i++) diff --git a/src/libs/pestpp_common/Pest.cpp b/src/libs/pestpp_common/Pest.cpp index a6aaf3e9..6c94eefe 100644 --- a/src/libs/pestpp_common/Pest.cpp +++ b/src/libs/pestpp_common/Pest.cpp @@ -1141,7 +1141,7 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec) //try to use this as a control data arg if (stat == PestppOptions::ARG_STATUS::ARG_NOTFOUND) { - stat = control_info.assign_value_by_key(kv.first,kv.second); + stat = control_info.assign_value_by_key(kv.first,kv.second,f_rec); check_report_assignment(f_rec, stat, kv.first, kv.second); } diff --git a/src/libs/pestpp_common/pest_data_structs.cpp b/src/libs/pestpp_common/pest_data_structs.cpp index ef28e356..11de6594 100644 --- a/src/libs/pestpp_common/pest_data_structs.cpp +++ b/src/libs/pestpp_common/pest_data_structs.cpp @@ -2251,7 +2251,7 @@ ostream& operator<< (ostream &os, const SVDInfo& val) return os; } -PestppOptions::ARG_STATUS ControlInfo::assign_value_by_key(const string key, const string org_value) +PestppOptions::ARG_STATUS ControlInfo::assign_value_by_key(const string key, const string org_value, ofstream& f_rec) { /*enum PestMode { ESTIMATION, REGUL, PARETO, UNKNOWN }; double phiredswh; @@ -2266,6 +2266,7 @@ PestppOptions::ARG_STATUS ControlInfo::assign_value_by_key(const string key, con int noptswitch; double splitswh; PestMode pestmode;*/ + stringstream ss; set valid_args{"RSTFLE","PESTMODE","NPAR","NOBS","NPARGP","NPRIOR","NOBSGP","MAXCOMPDIM","NTPLFLE","NINSFLE", "PRECIS","DPOINT","NUMCOM","JACFILE","MESSFILE","OBSREREF","RLAMBDA1","RLAMFAC","PHIREDLAM", "PHIRATSUF","NUMLAM","JACUPDATE","LAMFORGIVE","DERFORGIVE","RELPARMAX","FACPARMAX", @@ -2298,6 +2299,14 @@ PestppOptions::ARG_STATUS ControlInfo::assign_value_by_key(const string key, con convert_ip(value, noptswitch); else if (key == "SPLITSWH") convert_ip(value, splitswh); + else if (key == "PHIREDSWH") + { + convert_ip(value,phiredswh); + } + else if (key == "NRELPAR") + { + convert_ip(value, nrelpar); + } else if (key == "PESTMODE") { if (value == "ESTIMATION") @@ -2311,7 +2320,10 @@ PestppOptions::ARG_STATUS ControlInfo::assign_value_by_key(const string key, con } else if (valid_args.find(key) != valid_args.end()) { - //a valid but unused ctl data arg + ss.str(""); + ss << "WARNING: unused 'control data keyword' option,value:" << key << "," << value; + f_rec << ss.str(); + cout << ss.str(); } else return PestppOptions::ARG_STATUS::ARG_NOTFOUND; diff --git a/src/libs/pestpp_common/pest_data_structs.h b/src/libs/pestpp_common/pest_data_structs.h index a569a323..0236ef2e 100644 --- a/src/libs/pestpp_common/pest_data_structs.h +++ b/src/libs/pestpp_common/pest_data_structs.h @@ -879,7 +879,7 @@ class ControlInfo { int noptswitch; double splitswh; PestMode pestmode; - PestppOptions::ARG_STATUS assign_value_by_key(const string key, const string org_value); + PestppOptions::ARG_STATUS assign_value_by_key(const string key, const string org_value, ofstream& f_rec); ControlInfo() { ; } void set_defaults();