diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f38c62d4..b0719d46 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ jobs: fail-fast: false matrix: os: [windows-latest, ubuntu-latest] - python-version: [3.9] # , 3.7, 3.6] + python-version: [3.10] # , 3.7, 3.6] run-type: [std] # test_repo: [""] # test_dir: [""] @@ -46,7 +46,7 @@ jobs: steps: - - uses: actions/checkout@v2.3.4 + - uses: actions/checkout@v4 # - name: Setup Ninja # if: ${{ runner.os == 'Windows' }} - uses: seanmiddleditch/gha-setup-ninja@master @@ -179,9 +179,15 @@ jobs: ls -l cd ${{ env.test_dir }} - nosetests --nocapture -v ${{ env.test_script }} - cd - + nosetests -v ${{ env.test_script }} + status=$? + echo $status + + if [ $status -eq 0 ]; then + exit 0 + else + exit 127 + fi diff --git a/benchmarks/basic_tests.py b/benchmarks/basic_tests.py index a9eec28e..0ae729ea 100644 --- a/benchmarks/basic_tests.py +++ b/benchmarks/basic_tests.py @@ -821,6 +821,7 @@ def ext_stdcol_test(): par.loc[pst.adj_par_names,"standard_deviation"] = (par.loc[pst.adj_par_names,"parubnd_trans"] - par.loc[pst.adj_par_names,"parlbnd_trans"]) / 4.0 #par.loc[pst.adj_par_names[0],"mean"] = par.loc[pst.adj_par_names[0],"parubnd"] pst.pestpp_options["ies_num_reals"] = 10 + pst.pestpp_options["ies_no_noise"] = False pst.control_data.noptmax = -1 pst.write(os.path.join(m_d,"pest_base.pst")) pyemu.os_utils.run("{0} pest_base.pst".format(exe_path),cwd=m_d) @@ -1481,14 +1482,16 @@ def sweep_bin_test(): print(np.abs(diff).max()) assert np.abs(diff).max() < 1e-7 - +#def fail_test(): +# raise Exception("fail please") if __name__ == "__main__": #run() #mf6_v5_ies_test() #prep_ends() - sweep_bin_test() + #sweep_bin_test() #mf6_v5_sen_test() + #ext_stdcol_test() #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-glm.exe"),os.path.join("..","bin","win","pestpp-glm.exe")) #shutil.copy2(os.path.join("..", "exe", "windows", "x64", "Debug", "pestpp-ies.exe"), # os.path.join("..", "bin", "win", "pestpp-ies.exe")) @@ -1526,7 +1529,7 @@ def sweep_bin_test(): #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-ies.exe"),os.path.join("..","bin","win","pestpp-ies.exe")) #tplins1_test() - #fr_timeout_test() + fr_timeout_test() #mf6_v5_ies_test() #mf6_v5_sen_test() diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index 15d9c7d4..eb876fd0 100644 --- a/src/libs/common/config_os.h +++ b/src/libs/common/config_os.h @@ -2,7 +2,7 @@ #define CONFIG_OS_H_ -#define PESTPP_VERSION "5.2.10"; +#define PESTPP_VERSION "5.2.11"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/common/system_variables.cpp b/src/libs/common/system_variables.cpp index 00eac76b..7caa60b6 100644 --- a/src/libs/common/system_variables.cpp +++ b/src/libs/common/system_variables.cpp @@ -165,7 +165,7 @@ int start(string &cmd_string) // argv[icmd] = cmds[icmd].data(); //} //argv[cmds.size() + 1] = NULL; //last arg must be NULL - + arg_v.push_back(NULL); pid_t pid = fork(); if (pid == 0) diff --git a/src/libs/common/utilities.cpp b/src/libs/common/utilities.cpp index 222878f7..3f716824 100644 --- a/src/libs/common/utilities.cpp +++ b/src/libs/common/utilities.cpp @@ -1127,12 +1127,14 @@ void read_dense_binary(const string& filename, vector& row_names, vector cout << "reading 'dense' format matrix with " << n_obs_and_pi << " columns" << endl; //first read the names of the columns and the rows col_names = read_dense_binary_col_names(in,n_obs_and_pi); + streampos first_record = in.tellg(); row_names = read_dense_binary_remaining_row_names(in,col_names); in.close(); - in.open(filename.c_str(), ifstream::binary); + in.open(filename.c_str(), ifstream::binary); + in.seekg(first_record); //resize the matrix now that we know big it should be matrix.resize(row_names.size(), col_names.size()); diff --git a/src/libs/pestpp_common/Ensemble.cpp b/src/libs/pestpp_common/Ensemble.cpp index 6f9b97f4..0701f1e0 100644 --- a/src/libs/pestpp_common/Ensemble.cpp +++ b/src/libs/pestpp_common/Ensemble.cpp @@ -2166,10 +2166,18 @@ ParameterEnsemble ParameterEnsemble::zeros_like(int nrows) if (nrows < 0) nrows = real_names.size(); Eigen::MatrixXd new_reals = Eigen::MatrixXd::Zero(nrows, var_names.size()); - + stringstream ss; vector new_real_names; for (int i = 0; i < nrows; i++) - new_real_names.push_back(real_names[i]); + if (i < real_names.size()) + new_real_names.push_back(real_names[i]); + else + { + ss.str(""); + ss << "zeros_like_real_" << i; + new_real_names.push_back(ss.str()); + } + ParameterEnsemble new_en(pest_scenario_ptr, rand_gen_ptr); new_en.from_eigen_mat(new_reals, new_real_names, var_names); diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index fe1181c2..aa260072 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -1374,13 +1374,13 @@ void UpgradeThread::ensemble_solution(const int iter, const int verbose_level,co Eigen::MatrixXd temp = parcov_inv.diagonal().matrix(); local_utils::save_mat(verbose_level, thread_id, iter, t_count, "parcov_inv", temp); } - if (act_obs_names.size() > 0) { //this works bc the mm solve doesnt pass these names... - ss.str(""); - ss << "solution scaling factor: " << scale << endl; - ss << "eigthresh: " << eigthresh << endl; - ss << "maxsing: " << maxsing << endl; - cout << ss.str() << endl; - } +// if (act_obs_names.size() > 0) { //this works bc the mm solve doesnt pass these names... +// ss.str(""); +// ss << "solution scaling factor: " << scale << endl; +// ss << "eigthresh: " << eigthresh << endl; +// ss << "maxsing: " << maxsing << endl; +// cout << ss.str() << endl; +// } } if (use_prior_scaling) @@ -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) { @@ -5676,7 +4726,8 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) //check to see if any explicit obs noise options are set bool reset_to_nonoise = true; - if (ppo->get_passed_args().find("IES_NO_NOISE") != ppo->get_passed_args().end()) + set passed = ppo->get_passed_args(); + if (passed.find("IES_NO_NOISE") != passed.end()) { reset_to_nonoise = false; } @@ -5684,12 +4735,15 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) reset_to_nonoise = false; else if (!ppo->get_obscov_filename().empty()) reset_to_nonoise = false; + else if (ppo->get_ies_use_mda()) + reset_to_nonoise = false; else { map obs_std = pest_scenario.get_ext_file_double_map("observation data external", "standard_deviation"); if (obs_std.size() > 0) reset_to_nonoise = false; } + if (reset_to_nonoise) { ss.str(""); @@ -6171,13 +5225,16 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) ss.str(""); ss << " WARNING: Prior-data conflict detected. Continuing with IES parameter" << endl; ss << " adjustment will likely result in parameter and forecast bias." << endl; - ss << " Consider using 'ies_drop_conflicts' as a quick fix."; + //ss << " Consider using 'ies_drop_conflicts' as a quick fix."; message(1, ss.str()); } else { - - //check that all obs are in conflict + ss.str(""); + ss << " WARNING: 'ies_drop_conflicts' is not an ideal approach. A better" << endl; + ss << " approach would be to dig in and figure out why these observations " << endl; + ss << " are in conflict." << endl; + //check that all obs are in conflict message(1, "dropping conflicted observations"); if (in_conflict.size() == act_obs_names.size()) { @@ -6456,7 +5513,14 @@ void EnsembleMethod::check_and_fill_phi_factors(map>& grou { for (auto& entry : row.second) { - if (entry.second <= 0.0) + if (entry.second == -999.) + { + ss.str(""); + ss << "phi factor tag '" << row.first << "' has -999 value, indicating that weights for this tag should not be adjusted"; + message(1,ss.str()); + + } + else if (entry.second <= 0.0) { problems.push_back(make_pair(row.first,entry.first)); } @@ -6585,7 +5649,13 @@ void EnsembleMethod::check_and_fill_phi_factors(map>& grou for (auto &g : group_map.at(pf.first)) ss << g << ","; message(2, ss.str()); - if (pf.second <= 0.0) { + if (pf.second == -999.0) + { + ss.str(""); + ss << "phi factor tag '" << pf.first << "' has -999 value, indicating that weights for this tag should not be adjusted"; + message(1,ss.str()); + } + else if (pf.second <= 0.0) { ss.str(""); ss << "adjust_weights(): phi factor '" << pf.first << "' less or equal 0.0 - this is not allowed"; throw_em_error(ss.str()); @@ -6691,6 +5761,13 @@ void EnsembleMethod::adjust_weights_by_real(map>& group_to sub_total=0; scale_fac = 0; for (auto& pf: phi_fracs) { + if (pf.second == -999) + { + ss.str(""); + ss << "NOTE: adjust_weights(): tag " << pf.first << " for realization " << swr_map.first << " has a factor value of -999, skipping this tag"; + message(1,ss.str()); + continue; + } total = 0; for (auto &g : group_map.at(pf.first)) { sub_total = 0; @@ -6822,6 +5899,13 @@ void EnsembleMethod::adjust_weights_single(map>& group_to_ double sub_total = 0; for (auto& pf: phi_fracs) { + if (pf.second == -999) + { + ss.str(""); + ss << "NOTE: adjust_weights(): tag " << pf.first << " has a factor value of -999, skipping this tag"; + message(1,ss.str()); + continue; + } total = 0; for (auto& g : group_map.at(pf.first)) { diff --git a/src/libs/pestpp_common/MOEA.cpp b/src/libs/pestpp_common/MOEA.cpp index 9f825f16..61847e5b 100644 --- a/src/libs/pestpp_common/MOEA.cpp +++ b/src/libs/pestpp_common/MOEA.cpp @@ -291,9 +291,11 @@ void ParetoObjectives::update(ObservationEnsemble& op, ParameterEnsemble& dp, Co vector onames = op.get_var_names(), pnames = dp.get_var_names(); set obs_obj_set(obs_obj_names_ptr->begin(), obs_obj_names_ptr->end()); set pi_obj_set(pi_obj_names_ptr->begin(), pi_obj_names_ptr->end()); + set::iterator end; ObservationInfo *oi = pest_scenario.get_observation_info_ptr(); PriorInformation *pi = pest_scenario.get_prior_info_ptr(); feas_member_struct.clear(); + for (auto real_name : real_names) { vsum = 0.0; obs.update_without_clear(onames, op.get_real_vector(real_name)); @@ -301,14 +303,16 @@ void ParetoObjectives::update(ObservationEnsemble& op, ParameterEnsemble& dp, Co // the 'false' arg is to not apply risk shifting to the satisfaction calcs since // 'op' has already been shifted violations = constraints_ptr->get_unsatified_obs_constraints(obs, 0.0, false); + end = obs_obj_set.end(); for (auto v : violations) { - if (obs_obj_set.find(v.first) == obs_obj_set.end()) + if (obs_obj_set.find(v.first) == end) vsum += pow(v.second * oi->get_weight(v.first), 2); } pars.update_without_clear(pnames, dp.get_real_vector(real_name)); violations = constraints_ptr->get_unsatified_pi_constraints(pars, 0.0); + end = pi_obj_set.end(); for (auto v : violations) { - if (pi_obj_set.find(v.first) == pi_obj_set.end()) + if (pi_obj_set.find(v.first) == end) vsum += pow(v.second * pi->get_pi_rec(v.first).get_weight(), 2); } if (vsum > 0.0) { @@ -2993,6 +2997,9 @@ void MOEA::iterate_to_solution() if ((q == 1) || (q == 2)) { message(0,"'pest.stp' found, quitting"); + message(1,"force-saving current populations"); + save_populations(dp,op,"",true); + save_populations(dp_archive, op_archive, "archive",true); break; } else if (q == 4) { @@ -3382,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++) @@ -4126,7 +4136,7 @@ ParameterEnsemble MOEA::generate_sbx_population(int num_members, ParameterEnsemb return tmp_dp; } -void MOEA::save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag) +void MOEA::save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag, bool force_save) { stringstream ss; @@ -4153,7 +4163,7 @@ void MOEA::save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, st ss << "saved decision variable population of size " << _dp.shape().first << " X " << _dp.shape().second << " to '" << name << "'"; message(1, ss.str()); ss.str(""); - if ((save_every > 0) && (iter % save_every == 0)) + if (((save_every > 0) && (iter % save_every == 0)) || (iter == pest_scenario.get_control_info().noptmax) || (force_save)) { ss << file_manager.get_base_filename() << "." << iter; if (tag.size() > 0) diff --git a/src/libs/pestpp_common/MOEA.h b/src/libs/pestpp_common/MOEA.h index cbecfdc3..265f1bb4 100644 --- a/src/libs/pestpp_common/MOEA.h +++ b/src/libs/pestpp_common/MOEA.h @@ -224,7 +224,7 @@ class MOEA string get_new_member_name(string tag = string()); - void save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag = string()); + void save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag = string(), bool force_save=false); void gauss_mutation_ip(ParameterEnsemble& _dp); pair sbx(double probability, double eta_m, int idx1, int idx2); pair sbx_new(double crossover_probability, double di, Eigen::VectorXd& parent1, diff --git a/src/libs/pestpp_common/Pest.cpp b/src/libs/pestpp_common/Pest.cpp index a6aaf3e9..738dc8eb 100644 --- a/src/libs/pestpp_common/Pest.cpp +++ b/src/libs/pestpp_common/Pest.cpp @@ -463,12 +463,6 @@ vector Pest::get_nonregul_obs() const return ret_val; } -//int Pest::process_ctl_file_old(ifstream &fin, string pst_filename) -//{ -// ofstream f_out("ctl_process.out"); -// return process_ctl_file_old(fin, pst_filename, f_out); -//} - int Pest::process_ctl_file(ifstream& fin, string _pst_filename) { ofstream f_out("ctl_process.out"); @@ -476,521 +470,6 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename) } -//int Pest::process_ctl_file_old(ifstream &fin, string _pst_filename, ofstream &f_rec) -//{ -// string line; -// string line_upper; -// string section(""); -// vector tokens; -// int sec_begin_lnum, sec_lnum; -// double value; -// string name; -// string *trans_type; -// string prior_info_string; -// pair pi_name_group; -// int lnum; -// int num_par; -// int num_tpl_file; -// int dercom; -// int i_tpl_ins = 0; -// double phimlim; -// double phimaccept; -// double fracphim; -// double wfinit = 1.0; -// double wfmin = numeric_limits::min(); -// double wfmax = numeric_limits::max(); -// double wffac; -// double wftol; -// bool use_dynamic_reg = false; -// bool reg_adj_grp_weights = false; -// vector pestpp_input; -// regul_scheme_ptr = new DynamicRegularization(use_dynamic_reg); -// -// //dont change these text names - they are used in ParamTransformSeq -// TranTied *t_tied = new TranTied("PEST to model tied transformation"); -// TranOffset *t_offset = new TranOffset("PEST to model offset transformation"); -// TranScale *t_scale = new TranScale("PEST to model scale transformation"); -// TranLog10 *t_log = new TranLog10("PEST to model log transformation"); -// TranFixed *t_fixed = new TranFixed("PEST to model fixed transformation"); -// //TranNormalize *t_auto_norm = new TranNormalize("PEST auto-normalization transformation"); -// -// base_par_transform.push_back_ctl2model(t_scale); -// base_par_transform.push_back_ctl2model(t_offset); -// base_par_transform.push_back_ctl2active_ctl(t_tied); -// base_par_transform.push_back_ctl2active_ctl(t_fixed); -// base_par_transform.push_back_active_ctl2numeric(t_log); -// -// -// set tied_names; -// pst_filename = _pst_filename; -// -//#ifndef _DEBUG -// try { -//#endif -// prior_info_string = ""; -// for(lnum=1, sec_begin_lnum=1; getline(fin, line); ++ lnum) -// { -// strip_ip(line); -// line_upper = upper_cp(line); -// tokens.clear(); -// tokenize(line_upper, tokens); -// sec_lnum = lnum - sec_begin_lnum; -// -// if (lnum == 1) -// { -// if (tokens[0] != "PCF") -// { -// cout << "WARNING: fist line of control file should be 'PCF' not " << tokens[0] << endl; -// } -// } -// -// else if (tokens.empty()) -// { -// //skip blank line -// } -// else if (line[0] == '#') -// { -// -// } -// else if (line_upper.substr(0,2) == "++") -// { -// pestpp_input.push_back(line); -// } -// -// else if (line_upper[0] == '*') -// { -// section = upper_cp(strip_cp(line_upper, "both", " *\t\n")); -// sec_begin_lnum = lnum; -// } -// else if (section == "CONTROL DATA") -// { -// if (sec_lnum == 1) -// { -// if (tokens[1] == "REGULARIZATION" || tokens[1] == "REGULARISATION") -// { -// use_dynamic_reg = true; -// control_info.pestmode = ControlInfo::PestMode::REGUL; -// -// } -// else if (tokens[1] == "ESTIMATION") -// control_info.pestmode = ControlInfo::PestMode::ESTIMATION; -// else if (tokens[1] == "PARETO") -// control_info.pestmode = ControlInfo::PestMode::PARETO; -// } -// -// -// else if (sec_lnum == 2) -// { -// convert_ip(tokens[0], num_par); -// } -// else if (sec_lnum == 3) -// { -// convert_ip(tokens[0], num_tpl_file); -// if(tokens.size() >= 5) { -// convert_ip(tokens[4], control_info.numcom); -// } -// else { -// control_info.numcom = 0; -// } -// if(tokens.size() >= 6) { -// convert_ip(tokens[5], control_info.jacfile); -// } -// else { -// control_info.jacfile = 0; -// } -// } -// else if (sec_lnum == 5) -// { -// convert_ip(tokens[0], control_info.relparmax); -// convert_ip(tokens[1], control_info.facparmax); -// convert_ip(tokens[2], control_info.facorig); -// } -// else if (sec_lnum == 6) -// { -// // remove text arguements from the line as these can be specified out of order -// // and PEST++ does not use them -// set remove_tags = { "aui", "auid", "noaui", "senreuse", "nsenreuse", "boundscale", "noboundscale" }; -// auto end_iter = std::remove_if(tokens.begin(), tokens.end(), -// [&remove_tags](string &str)->bool{return (remove_tags.find(upper_cp(str)) != remove_tags.end() -// || remove_tags.find(lower_cp(str)) != remove_tags.end()); }); -// tokens.resize(std::distance(tokens.begin(), end_iter)); -// -// convert_ip(tokens[0], control_info.phiredswh); -// if (tokens.size() >= 2) { -// convert_ip(tokens[1], control_info.noptswitch); -// } -// if (tokens.size() >= 3) { -// convert_ip(tokens[2], control_info.splitswh); -// } -// -// } -// else if (sec_lnum == 7) -// { -// convert_ip(tokens[0], control_info.noptmax); -// convert_ip(tokens[1], control_info.phiredstp); -// convert_ip(tokens[2], control_info.nphistp); -// convert_ip(tokens[3], control_info.nphinored); -// convert_ip(tokens[4], control_info.relparstp); -// convert_ip(tokens[5], control_info.nrelpar); -// } -// else -// { -// other_lines[lnum] = line; -// } -// } -// else if (section == "SINGULAR VALUE DECOMPOSITION") -// { -// if (sec_lnum == 2) { -// convert_ip(tokens[0], svd_info.maxsing); -// convert_ip(tokens[1], svd_info.eigthresh); -// } -// else if (sec_lnum == 3) { -// convert_ip(tokens[0], svd_info.eigwrite); -// } -// else -// { -// other_lines[lnum] = line; -// } -// } -// else if (section == "PARAMETER GROUPS") -// { -// ParameterGroupRec pgi; -// name = tokens[0]; -// size_t n_tokens = tokens.size(); -// pgi.name = name; -// ctl_ordered_par_group_names.push_back(name); -// convert_ip(tokens[1], pgi.inctyp); -// convert_ip(tokens[2], pgi.derinc); -// convert_ip(tokens[3], pgi.derinclb); -// convert_ip(tokens[4], pgi.forcen); -// convert_ip(tokens[5], pgi.derincmul); -// convert_ip(tokens[6], pgi.dermthd); -// if (n_tokens >= 8) convert_ip(tokens[7], pgi.splitthresh); -// if (n_tokens >= 9) convert_ip(tokens[8], pgi.splitreldiff); -// base_group_info.insert_group(name, pgi); -// } -// else if (section == "PARAMETER DATA") -// { -// if (sec_lnum <= num_par) { -// double scale; -// double offset; -// ParameterRec pi; -// name = tokens[0]; -// trans_type = &tokens[1]; -// convert_ip(tokens[2], pi.chglim); -// convert_ip(tokens[3], pi.init_value); -// convert_ip(tokens[4], pi.lbnd); -// convert_ip(tokens[5], pi.ubnd); -// convert_ip(tokens[6], pi.group); -// convert_ip(tokens[7], scale); -// convert_ip(tokens[8], offset); -// if (control_info.numcom > 1) -// convert_ip(tokens[9], pi.dercom); -// else -// pi.dercom = 1; -// pi.scale = scale; -// pi.offset = offset; -// // add parameters to model parameter and paramter_info datasets -// ctl_ordered_par_names.push_back(name); -// if (*trans_type == "FIXED") -// { -// pi.tranform_type = ParameterRec::TRAN_TYPE::FIXED; -// } -// else if (*trans_type == "LOG") -// { -// pi.tranform_type = ParameterRec::TRAN_TYPE::LOG; -// n_adj_par++; -// } -// else if (*trans_type == "TIED") -// { -// pi.tranform_type = ParameterRec::TRAN_TYPE::TIED; -// } -// else if (*trans_type == "NONE") -// { -// pi.tranform_type = ParameterRec::TRAN_TYPE::NONE; -// n_adj_par++; -// } -// else -// { -// //pi.tranform_type = ParameterRec::TRAN_TYPE::NONE; -// //assert(true); -// //n_adj_par++; -// throw PestError("unrecognized partrans for par " + name + ": " + *trans_type); -// } -// ctl_parameter_info.insert(name, pi); -// ctl_parameters.insert(name, pi.init_value); -// base_group_info.insert_parameter_link(name, pi.group); -// -// // build appropriate transformations -// if (*trans_type == "FIXED") { -// t_fixed->insert(name, pi.init_value);} -// else if (*trans_type == "LOG") { -// t_log->insert(name); -// } -// if (offset!=0) { -// t_offset->insert(name, offset); -// } -// if (scale !=1) { -// t_scale->insert(name, scale); -// } -// } -// // Get rest of information for tied paramters -// else { -// name = tokens[0]; -// string name_tied = tokens[1]; -// double ratio = ctl_parameters[name] / ctl_parameters[name_tied]; -// t_tied->insert(name, pair(name_tied, ratio)); -// tied_names.insert(name_tied); -// } -// } -// else if (section == "OBSERVATION GROUPS") -// { -// string name = tokens[0]; -// if (tokens.size() > 1) -// { -// stringstream ss; -// ss << "observation covariance matrix detected for group '" << tokens[0] << "' - these are not supported...yet!"; -// string s = ss.str(); -// throw PestError(s); -// } -// ObservationGroupRec group_rec; -// observation_info.groups[name] = group_rec; -// vector::iterator is = find(ctl_ordered_obs_group_names.begin(), ctl_ordered_obs_group_names.end(), name); -// if (is == ctl_ordered_obs_group_names.end()) -// { -// ctl_ordered_obs_group_names.push_back(name); -// } -// } -// else if (section == "OBSERVATION DATA") -// { -// ObservationRec obs_i; -// name = tokens[0]; -// convert_ip(tokens[1], value); -// convert_ip(tokens[2], obs_i.weight); -// obs_i.group = tokens[3]; -// ctl_ordered_obs_names.push_back(name); -// observation_info.observations[name] = obs_i; -// observation_values.insert(name, value); -// } -// -// else if (section == "PRIOR INFORMATION") -// { -// //This section processes the prior information. It does not write out the -// //last prior infomration. THis is because it must check for line continuations -// if (!prior_info_string.empty() && tokens[0] != "&"){ -// pi_name_group = prior_info.AddRecord(prior_info_string); -// ctl_ordered_pi_names.push_back(pi_name_group.first); -// vector::iterator is = find(ctl_ordered_obs_group_names.begin(), ctl_ordered_obs_group_names.end(), pi_name_group.second); -// if (is == ctl_ordered_obs_group_names.end()) -// { -// ctl_ordered_obs_group_names.push_back(pi_name_group.second); -// } -// prior_info_string.clear(); -// } -// else if (tokens[0] == "&") { -// prior_info_string.append(" "); -// } -// prior_info_string.append(line_upper); -// } -// -// else if (section == "PRIOR INFORMATION" ) -// { -// //This section processes the prior information. It does not write out the -// //last prior infomration -// if (!prior_info_string.empty() && tokens[0] != "&") { -// pi_name_group = prior_info.AddRecord(prior_info_string); -// ctl_ordered_pi_names.push_back(pi_name_group.first); -// vector::iterator is = find(ctl_ordered_obs_group_names.begin(), ctl_ordered_obs_group_names.end(), pi_name_group.second); -// if (is == ctl_ordered_obs_group_names.end()) -// { -// ctl_ordered_obs_group_names.push_back(pi_name_group.second); -// } -// prior_info_string.clear(); -// } -// else if (tokens[0] != "&") { -// prior_info_string.append(" "); -// } -// prior_info_string.append(line); -// } -// else if (section == "MODEL COMMAND LINE" ) -// { -// model_exec_info.comline_vec.push_back(line); -// } -// else if (section == "MODEL INPUT/OUTPUT" ) -// { -// vector tokens_case_sen; -// tokenize(line, tokens_case_sen); -// if(i_tpl_ins < num_tpl_file) -// { -// model_exec_info.tplfile_vec.push_back(tokens_case_sen[0]); -// model_exec_info.inpfile_vec.push_back(tokens_case_sen[1]); -// } -// else -// { -// model_exec_info.insfile_vec.push_back(tokens_case_sen[0]); -// model_exec_info.outfile_vec.push_back(tokens_case_sen[1]); -// } -// ++i_tpl_ins; -// } -// else if (section == "REGULARISATION" || section=="REGULARIZATION" ) -// { -// if (sec_lnum == 1) { -// convert_ip(tokens[0], phimlim); -// convert_ip(tokens[1], phimaccept); -// fracphim = 0.0; -// if(tokens.size() >=3) convert_ip(tokens[2], fracphim); -// } -// else if (sec_lnum == 2) { -// convert_ip(tokens[0], wfinit); -// convert_ip(tokens[1], wfmin); -// convert_ip(tokens[2], wfmax); -// } -// else if (sec_lnum == 3) { -// int iregadj; -// convert_ip(tokens[0], wffac); -// convert_ip(tokens[1], wftol); -// if (tokens.size() > 2) -// { -// convert_ip(tokens[2], iregadj); -// if (iregadj == 1) reg_adj_grp_weights = true; -// } -// delete regul_scheme_ptr; -// regul_scheme_ptr = new DynamicRegularization(use_dynamic_reg, reg_adj_grp_weights, phimlim, -// phimaccept, fracphim, wfmin, wfmax, wffac, wftol, wfinit); -// } -// else -// { -// other_lines[lnum] = line; -// } -// } -// else if (section == "PARETO") -// { -// if (sec_lnum == 1) -// { -// convert_ip(tokens[0], pareto_info.obsgroup); -// } -// else if (sec_lnum == 2) -// { -// convert_ip(tokens[0], pareto_info.wf_start); -// convert_ip(tokens[1], pareto_info.wf_fin); -// convert_ip(tokens[2], pareto_info.wf_inc); -// } -// else if (sec_lnum == 3) -// { -// convert_ip(tokens[0], pareto_info.niter_start); -// convert_ip(tokens[1], pareto_info.niter_gen); -// convert_ip(tokens[2], pareto_info.niter_fin); -// } -// -// } -// else -// { -// other_lines[lnum] = line; -// } -// } -// -// // write out last prior information record -// if (!prior_info_string.empty()) -// { -// pi_name_group = prior_info.AddRecord(prior_info_string); -// ctl_ordered_pi_names.push_back(pi_name_group.first); -// vector::iterator is = find(ctl_ordered_obs_group_names.begin(), ctl_ordered_obs_group_names.end(), pi_name_group.second); -// if (is == ctl_ordered_obs_group_names.end()) -// { -// ctl_ordered_obs_group_names.push_back(pi_name_group.second); -// } -// prior_info_string.clear(); -// } -//#ifndef _DEBUG -// } -// catch (PestConversionError &e) { -// std::stringstream out; -// out << "Error parsing \"" << pst_filename << "\" on line number " << lnum << endl; -// out << e.what() << endl; -// e.add_front(out.str()); -// e.raise(); -// } -//#endif -// fin.close(); -// -// map arg_map, line_arg_map; -// for(vector::const_iterator b=pestpp_input.begin(),e=pestpp_input.end(); -// b!=e; ++b) { -// -// try -// { -// line_arg_map = pestpp_options.parse_plusplus_line(*b); -// } -// catch (...) -// { -// throw runtime_error("error parsing ++ line '" + line + "'"); -// } -// arg_map.insert(line_arg_map.begin(),line_arg_map.end()); -// } -// //check for not "accepted" args here... -// -// -// regul_scheme_ptr->set_max_reg_iter(pestpp_options.get_max_reg_iter()); -// -// if (pestpp_options.get_tie_by_group()) -// { -// cout << "Note: ++tie_by_group(true) - tying adjustable parameters by groups" << endl; -// f_rec << "Note: ++tie_by_group(true) - tying adjustable parameters by groups" << endl; -// map> group_map; -// string gname; -// ParameterRec::TRAN_TYPE tlog = ParameterRec::TRAN_TYPE::LOG, -// tnone = ParameterRec::TRAN_TYPE::NONE, -// ttied = ParameterRec::TRAN_TYPE::TIED; -// int new_n_adj_par = 0; -// for (auto pname : ctl_ordered_par_names) -// { -// if ((ctl_parameter_info.get_parameter_rec_ptr(pname)->tranform_type == tlog) || -// (ctl_parameter_info.get_parameter_rec_ptr(pname)->tranform_type == tnone)) -// { -// -// if (tied_names.find(pname) != tied_names.end()) -// { -// new_n_adj_par++; -// continue; -// } -// gname = ctl_parameter_info.get_parameter_rec_ptr(pname)->group; -// if (group_map.find(gname) == group_map.end()) -// group_map[gname] = vector(); -// group_map[gname].push_back(pname); -// } -// -// } -// string tie_to_name; -// vector to_tie_names; -// vector::const_iterator first, last; -// for (auto gm : group_map) -// { -// tie_to_name = gm.second[0]; -// new_n_adj_par++; -// first = gm.second.begin() + 1; -// last = gm.second.end(); -// to_tie_names = vector(first, last); -// for (auto pname : to_tie_names) -// { -// double ratio = ctl_parameters[pname] / ctl_parameters[tie_to_name]; -// t_tied->insert(pname, pair(tie_to_name, ratio)); -// ctl_parameter_info.get_parameter_rec_ptr_4_mod(pname)->tranform_type = ttied; -// } -// -// } -// f_rec << "-->number of adjustable parameters reduced from " << n_adj_par << " to " << new_n_adj_par << endl; -// cout << "-->number of adjustable parameters reduced from " << n_adj_par << " to " << new_n_adj_par << endl; -// n_adj_par = new_n_adj_par; -// if (tied_names.size() > 0) -// { -// f_rec << "-->existing adjustable parameters that others tie to have been maintained:" << endl; -// for (auto tname : tied_names) -// f_rec << tname << endl; -// } -// } -// -// return 0; -//} int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec) { @@ -1141,7 +620,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/constraints.cpp b/src/libs/pestpp_common/constraints.cpp index f1f571c8..f90bc1f0 100644 --- a/src/libs/pestpp_common/constraints.cpp +++ b/src/libs/pestpp_common/constraints.cpp @@ -229,7 +229,7 @@ void Constraints::initialize(vector& ctl_ord_dec_var_names, double _dbl_ //current_constraints_sim_ptr = _current_constraints_sim_ptr; //same for dec var values //current_pars_and_dec_vars_ptr = _current_pars_and_dec_vars_ptr; - //check for a stack size arg - if postive, then use stacks for chances + //check for a stack size arg - if positive, then use stacks for chances stack_size = pest_scenario.get_pestpp_options().get_opt_stack_size(); //an existing parameter stack for chances string par_stack_name = pest_scenario.get_pestpp_options().get_opt_par_stack(); @@ -450,7 +450,12 @@ void Constraints::initialize(vector& ctl_ord_dec_var_names, double _dbl_ } set dec_set(ctl_ord_dec_var_names.begin(), ctl_ord_dec_var_names.end()); - for (auto& name : pest_scenario.get_ctl_ordered_par_names()) + vector prob_pars; + //disabled the exception for adj pars tied to dec vars + //string partied; + //map> tied_info = pest_scenario.get_base_par_tran_seq().get_tied_ptr()->get_items(); + + for (auto& name : pest_scenario.get_ctl_ordered_par_names()) { //if this parameter is not a decision var //if (find(start, end, name) == end) @@ -458,13 +463,26 @@ void Constraints::initialize(vector& ctl_ord_dec_var_names, double _dbl_ { ParameterRec::TRAN_TYPE tt = pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(name)->tranform_type; if ((tt == ParameterRec::TRAN_TYPE::LOG) || (tt == ParameterRec::TRAN_TYPE::NONE)) - adj_par_names.push_back(name); + { + adj_par_names.push_back(name); + } +// else if (tt == ParameterRec::TRAN_TYPE::TIED) +// { +// partied = tied_info[name].first; +// if (dec_set.find(partied) != dec_set.end()) +// prob_pars.push_back(name); +// } } + } +// if (!prob_pars.empty()) +// { +// throw_constraints_error("the following adjustable pars are 'tied' to decision variables",prob_pars); +// } //------------------------------------------ - // --- chance constratints --- + // --- chance constraints --- //------------------------------------------ risk = pest_scenario.get_pestpp_options().get_opt_risk(); if (risk != 0.5) 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(); diff --git a/src/libs/pestpp_common/sequential_lp.cpp b/src/libs/pestpp_common/sequential_lp.cpp index a2f2e788..cf30b41f 100644 --- a/src/libs/pestpp_common/sequential_lp.cpp +++ b/src/libs/pestpp_common/sequential_lp.cpp @@ -82,7 +82,7 @@ void sequentialLP::initial_report() f_rec << "-->objective function sense (direction): " << optobjfunc.get_obj_sense() << endl; - f_rec << "-->number of decision variable: " << num_dec_vars() << endl; + f_rec << "-->number of decision variables: " << num_dec_vars() << endl; f_rec << "-->number of observation constraints: " << constraints.num_obs_constraints() << endl; f_rec << "-->number of prior information constraints: " << constraints.num_pi_constraints() << endl; if (iter_derinc_fac != 1.0) @@ -326,12 +326,41 @@ void sequentialLP::initialize_and_check() //if any decision vars have a transformation that is not allowed vector problem_trans; - for (auto &name : dv_names) - if (pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(name)->tranform_type != ParameterRec::TRAN_TYPE::NONE) + vector temp_dv_names; + vector tied_dv_names; + //this should just skip over fixed dvs + for (auto &name : dv_names) + { + if (pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(name)->tranform_type == ParameterRec::TRAN_TYPE::NONE) + temp_dv_names.push_back(name); + else if (pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(name)->tranform_type == ParameterRec::TRAN_TYPE::TIED) + tied_dv_names.push_back(name); + else if (pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(name)->tranform_type == ParameterRec::TRAN_TYPE::LOG) problem_trans.push_back(name); - if (problem_trans.size() > 0) - throw_sequentialLP_error("the following decision variables don't have 'none' type parameter transformation: ", problem_trans); + } + if (!problem_trans.empty()) + throw_sequentialLP_error("the following decision variables have 'log' type parameter transformation: ", problem_trans); + if (!tied_dv_names.empty()) + { + set sdv_names(temp_dv_names.begin(),temp_dv_names.end()); + auto send = sdv_names.end(); + problem_trans.clear(); + string tied_name; + map> tied_info = pest_scenario.get_base_par_tran_seq().get_tied_ptr()->get_items(); + for (auto& name : tied_dv_names) + { + if (sdv_names.find(tied_info[name].first) == send) + problem_trans.push_back(name); + } + if (!problem_trans.empty()) + { + throw_sequentialLP_error("the following 'tied' decision variables are not tied to a decision variable: ", problem_trans); + } + + } + dv_names = temp_dv_names; + temp_dv_names.clear(); current_constraints_sim = pest_scenario.get_ctl_observations(); constraints.initialize(dv_names, COIN_DBL_MAX); @@ -958,11 +987,11 @@ bool sequentialLP::make_upgrade_run(Parameters &upgrade_pars, Observations &upgr { cout << " --- running the model once with optimal decision variables --- " << endl; - int run_id = run_mgr_ptr->add_run(par_trans.ctl2model_cp(upgrade_pars)); + int run_id = run_mgr_ptr->add_run(par_trans.active_ctl2model_cp(upgrade_pars)); run_mgr_ptr->run(); bool success = run_mgr_ptr->get_run(run_id, upgrade_pars, upgrade_obs); if (success) - par_trans.model2ctl_ip(upgrade_pars); + par_trans.model2active_ctl_ip(upgrade_pars); return success; } @@ -1044,7 +1073,7 @@ void sequentialLP::iter_presolve() } else { - //make the intial base run + //make the initial base run cout << " --- running the model once with initial decision variables --- " << endl; int run_id = run_mgr_ptr->add_run(par_trans.ctl2model_cp(current_pars)); //this would be only for stack runs since the fosm runs should have been in the jco diff --git a/src/libs/run_managers/abstract_base/model_interface.cpp b/src/libs/run_managers/abstract_base/model_interface.cpp index 7902b0b1..dfa314b4 100644 --- a/src/libs/run_managers/abstract_base/model_interface.cpp +++ b/src/libs/run_managers/abstract_base/model_interface.cpp @@ -274,7 +274,7 @@ void ThreadedTemplateProcess::work(int tid, vector& tpl_idx, Parameters par { if (tpl_idx.size() == 0) { - cout << "thread " << tid << " processed " << count << " template files" << endl; + //cout << "thread " << tid << " processed " << count << " template files" << endl; return; } i = tpl_idx[tpl_idx.size() - 1]; @@ -346,7 +346,8 @@ void ModelInterface::write_input_files(Parameters *pars_ptr) if (nnum_threads > tplfile_vec.size()) nnum_threads = tplfile_vec.size(); std::chrono::system_clock::time_point start_time = chrono::system_clock::now(); - cout << pest_utils::get_time_string() << " processing template files with " << nnum_threads << " threads..." << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " processing template files with " << nnum_threads << " threads..." << endl; vector threads; vector exception_ptrs; Parameters pro_pars = *pars_ptr; //copy @@ -454,7 +455,8 @@ void ModelInterface::write_input_files(Parameters *pars_ptr) } pars_ptr->update_without_clear(pro_pars.get_keys(), pro_pars.get_data_vec(pro_pars.get_keys())); - cout << pest_utils::get_time_string() << " done, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " done, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; } void ModelInterface::read_output_files(Observations *obs) @@ -463,7 +465,8 @@ void ModelInterface::read_output_files(Observations *obs) if (nnum_threads > insfile_vec.size()) nnum_threads = insfile_vec.size(); std::chrono::system_clock::time_point start_time = chrono::system_clock::now(); - cout << pest_utils::get_time_string() << " processing instruction files with " << nnum_threads << " threads..." << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " processing instruction files with " << nnum_threads << " threads..." << endl; vector threads; vector exception_ptrs; Observations temp_obs; @@ -600,7 +603,8 @@ void ModelInterface::read_output_files(Observations *obs) } t = temp_obs.get_keys(); obs->update(t, temp_obs.get_data_vec(t)); - cout << pest_utils::get_time_string() << " done, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " done, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; } @@ -661,7 +665,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ write_input_files(pars_ptr); std::chrono::system_clock::time_point start_time = chrono::system_clock::now(); - cout << pest_utils::get_time_string() << " calling forward run command(s)" << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " calling forward run command(s)" << endl; #ifdef OS_WIN //a flag to track if the run was terminated @@ -677,7 +682,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ } for (auto &cmd_string : comline_vec) { - cout << pest_utils::get_time_string() << " calling forward run command: '" << cmd_string << "' " << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " calling forward run command: '" << cmd_string << "' " << endl; //start the command PROCESS_INFORMATION pi; try @@ -694,7 +700,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ throw PestError("could not add process to job object: " + cmd_string); } DWORD pid = pi.dwProcessId; - cout << "...pid: " << pid << endl; + if (should_echo) + cout << "...pid: " << pid << endl; DWORD exitcode; while (true) { @@ -717,7 +724,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ //check for termination flag if (terminate->get()) { - std::cout << "received terminate signal" << std::endl; + if (should_echo) + std::cout << "received terminate signal" << std::endl; //try to kill the process bool success = (CloseHandle(job) != 0); @@ -745,10 +753,12 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ bool term_break = false; for (auto &cmd_string : comline_vec) { - cout << pest_utils::get_time_string() << " calling forward run command: '" << cmd_string << "' " << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " calling forward run command: '" << cmd_string << "' " << endl; //start the command int command_pid = start(cmd_string); - cout << "...pid: " << command_pid << endl; + if (should_echo) + cout << "...pid: " << command_pid << endl; while (true) { std::this_thread::sleep_for(std::chrono::milliseconds(sleep_ms)); @@ -770,7 +780,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ //check for termination flag if (terminate->get()) { - std::cout << "received terminate signal" << std::endl; + if (should_echo) + std::cout << "received terminate signal" << std::endl; //try to kill the process errno = 0; int success = kill(-command_pid, SIGKILL); @@ -789,8 +800,8 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ if (term_break) break; } #endif - - cout << pest_utils::get_time_string() << " foward run command(s) finished, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; + if (should_echo) + cout << pest_utils::get_time_string() << " forward run command(s) finished, took " << pest_utils::get_duration_sec(start_time) << " seconds" << endl; if (term_break) return; diff --git a/src/libs/run_managers/abstract_base/model_interface.h b/src/libs/run_managers/abstract_base/model_interface.h index 93d64edf..7206272a 100644 --- a/src/libs/run_managers/abstract_base/model_interface.h +++ b/src/libs/run_managers/abstract_base/model_interface.h @@ -113,6 +113,7 @@ class ModelInterface{ void set_tpl_force_decimal(bool _flag) {tpl_force_decimal = _flag;} void set_num_threads(int _num_threads) { num_threads = _num_threads; } void set_sleep_ms(int _sleep_ms){sleep_ms = _sleep_ms;} + void set_should_echo(bool _should_echo){should_echo=_should_echo;} private: int num_threads; @@ -128,6 +129,7 @@ class ModelInterface{ bool fill_tpl_zeros; bool tpl_force_decimal; string additional_ins_delimiters; + bool should_echo; void write_input_files(Parameters *pars_ptr); void read_output_files(Observations *obs_ptr); diff --git a/src/libs/run_managers/serial/RunManagerSerial.cpp b/src/libs/run_managers/serial/RunManagerSerial.cpp index d5d7a2ce..1d9c5e0c 100644 --- a/src/libs/run_managers/serial/RunManagerSerial.cpp +++ b/src/libs/run_managers/serial/RunManagerSerial.cpp @@ -40,7 +40,7 @@ RunManagerSerial::RunManagerSerial(const vector _comline_vec, const vector _insfile_vec, const vector _outfile_vec, const string &stor_filename, const string &_run_dir, int _max_run_fail, bool fill_tpl_zeros, string additional_ins_delimiters, int _num_threads, - bool tpl_force_decimal) + bool tpl_force_decimal, bool should_echo) : RunManagerAbstract(_comline_vec, _tplfile_vec, _inpfile_vec, _insfile_vec, _outfile_vec, stor_filename, _max_run_fail), run_dir(_run_dir), mi(_tplfile_vec,_inpfile_vec,_insfile_vec,_outfile_vec, _comline_vec) @@ -50,6 +50,7 @@ RunManagerSerial::RunManagerSerial(const vector _comline_vec, mi.set_num_threads(_num_threads); mi.set_tpl_force_decimal(tpl_force_decimal); mi.set_sleep_ms(5); + mi.set_should_echo(should_echo); cout << " starting serial run manager ..." << endl << endl; mgr_type = RUN_MGR_TYPE::SERIAL; } diff --git a/src/libs/run_managers/serial/RunManagerSerial.h b/src/libs/run_managers/serial/RunManagerSerial.h index aca8189b..7d7e1554 100644 --- a/src/libs/run_managers/serial/RunManagerSerial.h +++ b/src/libs/run_managers/serial/RunManagerSerial.h @@ -31,13 +31,14 @@ class RunManagerSerial : public RunManagerAbstract const std::vector _insfile_vec, const std::vector _outfile_vec, const std::string &stor_filename, const std::string &run_dir, int _max_run_fail=1, bool fill_tpl_zeros=false, string additional_ins_delimiters="", int _num_threads=1, - bool tpl_force_decimal=false); + bool tpl_force_decimal=false, bool should_echo=true); virtual void run(); ~RunManagerSerial(void); private: ModelInterface mi; std::string run_dir; + void run_async(pest_utils::thread_flag* terminate, pest_utils::thread_flag* finished, exception_ptr& run_exception, Parameters* pars, Observations* obs); diff --git a/src/programs/pestpp-ies/pestpp-ies.cpp b/src/programs/pestpp-ies/pestpp-ies.cpp index 41257816..cce7cdb8 100644 --- a/src/programs/pestpp-ies/pestpp-ies.cpp +++ b/src/programs/pestpp-ies/pestpp-ies.cpp @@ -239,7 +239,8 @@ int main(int argc, char* argv[]) pest_scenario.get_pestpp_options().get_fill_tpl_zeros(), pest_scenario.get_pestpp_options().get_additional_ins_delimiters(), pest_scenario.get_pestpp_options().get_num_tpl_ins_threads(), - pest_scenario.get_pestpp_options().get_tpl_force_decimal()); + pest_scenario.get_pestpp_options().get_tpl_force_decimal(), + pest_scenario.get_pestpp_options().get_panther_echo()); } const ParamTransformSeq &base_trans_seq = pest_scenario.get_base_par_tran_seq(); diff --git a/src/programs/sweep/pestpp-swp.vcxproj.user b/src/programs/sweep/pestpp-swp.vcxproj.user index 249a734e..a0d40ce8 100644 --- a/src/programs/sweep/pestpp-swp.vcxproj.user +++ b/src/programs/sweep/pestpp-swp.vcxproj.user @@ -1,8 +1,8 @@  - 10par_xsec.pst - C:\Dev\pyemu\autotest\moouu\10par_xsec\template + pest_forgive.pst + C:\Dev\pestpp\benchmarks\ies_10par_xsec\master_sweep_bin WindowsLocalDebugger diff --git a/src/programs/sweep/sweep.cpp b/src/programs/sweep/sweep.cpp index 4740ed38..c96b9333 100644 --- a/src/programs/sweep/sweep.cpp +++ b/src/programs/sweep/sweep.cpp @@ -716,6 +716,8 @@ int main(int argc, char* argv[]) { cout << " --- dense binary file detected for par_csv" << endl; fout_rec << " --- dense binary file detected for par_csv" << endl; + par_stream.close(); + par_stream.open(par_csv_file, ifstream::binary); vector col_names; header_info = prepare_parameter_dense_binary(pest_scenario.get_ctl_parameters(),par_stream, pest_scenario.get_pestpp_options().get_sweep_forgive(),col_names);