Merge branch 'develop' of https://github.com/usgs/pestpp

README.md

@@ -14,7 +14,7 @@ PEST++ is a software suite aimed at supporting complex numerical models in the d
 
 ## Documentation
 
-The lastest PEST++ users manual is available [here](https://github.com/jtwhite79/pestpp/tree/develop/documentation). Direct zip download [here](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/jtwhite79/pestpp/tree/develop/documentation)
+The lastest PEST++ users manual is available [here](https://github.com/usgs/pestpp/tree/develop/documentation). Direct zip download [here](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/usgs/pestpp/tree/develop/documentation)
 
 ## Links to latest binaries
 
@@ -44,11 +44,11 @@ All members of the software suite can be compiled for PC, MAC, or Linux and have
 
 ## Funding
 
-Funding for PEST++ has been provided by the U.S. Geologial Survey. The New Zealand Strategic Science Investment Fund as part of GNS Science’s (https://www.gns.cri.nz/) Groundwater Research Programme has also funded contributions 2018-present.
+Funding for PEST++ has been provided by the U.S. Geologial Survey. The New Zealand Strategic Science Investment Fund as part of GNS Science’s (https://www.gns.cri.nz/) Groundwater Research Programme has also funded contributions 2018-present.  Intera, Inc. also provides ongoing support for PEST++.
 
 ## Recent developements
 
-PEST++ version 5 will be released soon!
+PEST++ version 5 has been released!  Please see the users manual for current input instructions and options for all PEST++ tools.  Also, several new tools are in development, include PESTPP-DA (generalized data assimilation including iterative ensemble Kalman filter), PESTPP-MOU (single and multiple objective constrained optimization under uncertainty) and PESTPP-SQP (ensemble-based constrainted sequential quadratic programming under uncertainty).  If you would like to be an early adopter/beta tester, please let us know!
 
 The PEST++ suite has been refactored to remove the fortran dependancy.  This includes the template and instruction file processing routines.
 

benchmarks/basic_tests.py

@@ -752,7 +752,9 @@ def tplins1_test():
     pst = pyemu.Pst(os.path.join(t_d, "pest.pst"))
     dum_obs = ['h01_03', 'h01_07']
     pst.observation_data.drop(index=dum_obs, inplace=True)
-    pst.instruction_files = ['out1dum.dat.ins']
+    #pst.instruction_files = ['out1dum.dat.ins']
+    pst.model_output_data = pd.DataFrame({"pest_file":"out1dum.dat.ins","model_file":"out1.dat"},
+                                         index=["out1dum.dat.ins"])
     pst.write(os.path.join(t_d, "pest_dum.pst"))
     pyemu.os_utils.run("{0} pest_dum.pst".format(exe_path.replace("-ies", "-glm")), cwd=t_d)
     obf_df = pd.read_csv(os.path.join(t_d, "out1.dat.obf"), delim_whitespace=True, header=None,
@@ -1012,7 +1014,7 @@ def cmdline_test():
     #tie_by_group_test()
     #sen_basic_test()
     #salib_verf()
-    #tplins1_test()
+    tplins1_test()
     #ext_stdcol_test()
     mf6_v5_ies_test()
     #mf6_v5_sen_test()

benchmarks/ies_10par_xsec/template/10par_xsec.hds

@@ -1,2 +1,2 @@
-     1.000     1.200     1.400     1.600     1.800     2.000     2.200     2.400     2.600     2.800
-     1.000     1.400     1.800     2.200     2.600     3.000     3.400     3.800     4.200     4.600
+     0.100     0.120     0.140     0.160     0.180     0.200     0.220     0.240     0.260     0.280
+     0.100     0.140     0.180     0.220     0.260     0.300     0.340     0.380     0.420     0.460

benchmarks/ies_10par_xsec/template/hk_Layer_1.ref

@@ -1 +1 @@
- 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000
+ 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000 25.00000000000

benchmarks/ies_10par_xsec/template/strt_Layer_1.ref

@@ -1 +1 @@
-   1.00000000000000   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00
+   0.10000000000000   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00   1.000000E+00

src/libs/common/config_os.h

@@ -2,7 +2,7 @@
 #define CONFIG_OS_H_
 
 
-#define PESTPP_VERSION "5.0.4";
+#define PESTPP_VERSION "5.0.6";
 
 #if defined(_WIN32) || defined(_WIN64)
 #define OS_WIN

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -664,9 +664,17 @@ vector<int> L2PhiHandler::get_idxs_greater_than(double bad_phi, double bad_phi_s
 	double std = calc_std(&_meas);
 	vector<int> idxs;
 	vector<string> names = oe.get_real_names();
-	for (int i=0;i<names.size();i++)
+	for (int i = 0; i < names.size(); i++)
+	{
+
 		if ((_meas[names[i]] > bad_phi) || (_meas[names[i]] > mean + (std * bad_phi_sigma)))
-			idxs.push_back(i);
+		{
+			if (names[i] == BASE_REAL_NAME)
+				cout << "...not dropping 'base' real even though phi is 'bad'" << endl;
+			else
+				idxs.push_back(i);
+		}	
+	}
 	return idxs;
 }
 

src/libs/pestpp_common/EnsembleSmoother.cpp

@@ -718,12 +718,12 @@ void IterEnsembleSmoother::initialize_restart()
 				if (pest_scenario.get_pestpp_options().get_ies_save_binary())
 				{
 					ss << file_manager.get_base_filename() << ".obs+noise.jcb";
-					oe.to_binary(ss.str());
+					oe_base.to_binary(ss.str());
 				}
 				else
 				{
 					ss << file_manager.get_base_filename() << ".obs+noise.csv";
-					oe.to_csv(ss.str());
+					oe_base.to_csv(ss.str());
 				}
 				message(1, "re-saved obs+noise observation ensemble (obsval+noise) to ", ss.str());
 			}
@@ -2945,7 +2945,7 @@ bool IterEnsembleSmoother::solve_new()
 		}
 		else
 		{
-			message(0, "not updating lambda");
+			message(0, "not updating lambda (standard deviation reduction criteria not met)");
 		}
 		last_best_std = best_std;
 	}
@@ -2956,10 +2956,34 @@ bool IterEnsembleSmoother::solve_new()
 		message(0, "only updating realizations with reduced phi");
 		update_reals_by_phi(pe_lams[best_idx], oe_lam_best);
 		ph.update(oe, pe);
-		double new_lam = last_best_lam * lam_inc;
-		new_lam = (new_lam > lambda_max) ? lambda_max : new_lam;
-		message(0, "incresing lambda to: ", new_lam);
-		last_best_lam = new_lam;
+		//re-check phi
+		best_mean = ph.get_mean(L2PhiHandler::phiType::COMPOSITE);
+		best_std = ph.get_std(L2PhiHandler::phiType::COMPOSITE);
+		//replace the last entry in the best mean phi tracker
+		best_mean_phis[best_mean_phis.size() - 1] = best_mean;
+		message(1, "current best mean phi (after updating reduced-phi reals): ", best_mean);
+		if (best_mean < last_best_mean * acc_fac)
+		{
+			if (best_std < last_best_std * acc_fac)
+			{
+				double new_lam = lam_vals[best_idx] * lam_dec;
+				new_lam = (new_lam < lambda_min) ? lambda_min : new_lam;
+				message(0, "updating lambda to ", new_lam);
+				last_best_lam = new_lam;
+			}
+			else
+			{
+				message(0, "not updating lambda (standard deviation reduction criteria not met)");
+			}
+			last_best_std = best_std;
+		}
+		else
+		{
+			double new_lam = last_best_lam * lam_inc;
+			new_lam = (new_lam > lambda_max) ? lambda_max : new_lam;
+			message(0, "incresing lambda to: ", new_lam);
+			last_best_lam = new_lam;
+		}
 	}
 	//report_and_save();
 	return true;