Refs/heads/feat adjbnds

autotest/en_tests.py

@@ -769,7 +769,70 @@ def mixed_par_draw_2_test():
     assert pst.npar == npar
 
 
+def draw_new_test():
+
+    import os
+    import numpy as np
+    import pyemu
+    pst = pyemu.Pst(os.path.join("en", "pest.pst"))
+    cov = pyemu.Cov.from_binary(os.path.join("en", "cov.jcb"))
+    print(pst.npar, cov.shape)
+    num_reals = 10000
+
+    pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst, cov=cov, num_reals=num_reals, factor="cholesky")
+
+    sub_pe = pe.iloc[:10000,:]
+
+    new_pe_nonoise = sub_pe.draw_new_ensemble(num_reals=num_reals)
+    new_pe_stdnoise = sub_pe.draw_new_ensemble(num_reals=num_reals,include_noise=True)
+    new_pe_usernoise = sub_pe.draw_new_ensemble(num_reals=num_reals,include_noise=1./np.sqrt(sub_pe.shape[0]))
+    new_pe_ennoise = sub_pe.draw_new_ensemble(num_reals=num_reals,include_noise=True,noise_reals=pe)
+
+    #pes = [pe,sub_pe,new_pe_nonoise,new_pe_stdnoise,new_pe_usernoise,new_pe_ennoise]
+    colors = ["r","y","b","g","m","c"]
+
+    pes = [pe,new_pe_stdnoise,new_pe_usernoise,new_pe_ennoise]
+    for ppe in pes:
+        ppe.transform()
+        #ppe.enforce()
+
+    for pname in pst.adj_par_names:
+        std = [ppe.loc[:,pname].std() for ppe in pes]
+        mean = [ppe.loc[:,pname].mean() for ppe in pes]
+        sdiff = [np.abs(std[0] - s) for s in std[1:]]
+
+        print(pname,max(sdiff))
+        assert max(sdiff) < 0.1
+
+    # pst.add_transform_columns()
+    # ubnd = pst.parameter_data.parubnd_trans.to_dict()
+    # lbnd = pst.parameter_data.parlbnd_trans.to_dict()
+
+
+    # import matplotlib.pyplot as plt
+    # from matplotlib.backends.backend_pdf import PdfPages
+    # with PdfPages("check.pdf") as pdf:
+    #     for pname in pst.adj_par_names:
+    #         fig,ax = plt.subplots(1,1,figsize=(10,10))
+    #         for ppe,c in zip(pes,colors):
+
+    #             ax.hist(ppe.loc[:,pname],fc=c,alpha=0.5,bins=10,density=True)
+    #         ylim = ax.get_ylim()
+    #         ax.plot([ubnd[pname],ubnd[pname]],ylim,"k--",lw=3)
+    #         ax.plot([lbnd[pname],lbnd[pname]],ylim,"k--",lw=3)
+
+    #         ax.grid()
+    #         ax.set_title(pname)
+    #         plt.tight_layout()
+    #         pdf.savefig()
+    #         plt.close(fig)
+    #         print(pname)
+
+
+
+
 if __name__ == "__main__":
+    draw_new_test()
     #par_gauss_draw_consistency_test()
     #obs_gauss_draw_consistency_test()
     # phi_vector_test()
@@ -784,7 +847,7 @@ def mixed_par_draw_2_test():
     # uniform_draw_test()
     #fill_test()
     #factor_draw_test()
-    emp_cov_test()
+    #emp_cov_test()
     #emp_cov_draw_test()
     #mixed_par_draw_2_test()
     #binary_test()

autotest/pst_from_tests.py

@@ -6272,8 +6272,8 @@ def xsec_pars_as_obs_test(tmp_path):
     sr = m.modelgrid
     t_d = "template_xsec"
     pf = pyemu.utils.PstFrom(tmp_model_ws,t_d,remove_existing=True,spatial_reference=sr)
-    pf.add_parameters("hk_Layer_1.ref",par_type="grid",par_style="direct",upper_bound=25,
-                      lower_bound=0.25)
+    pf.add_parameters("hk_Layer_1.ref",par_type="grid",par_style="direct",upper_bound=3,
+                      lower_bound=2,transform="none")
     pf.add_parameters("hk_Layer_1.ref", par_type="grid", par_style="multiplier", upper_bound=10.0,
                       lower_bound=0.1)
 
@@ -6291,17 +6291,37 @@ def xsec_pars_as_obs_test(tmp_path):
     pf.mod_sys_cmds.append("mfnwt {0}".format(nam_file))
 
     pst = pf.build_pst(None)
-    pst.add_pars_as_obs(pst_path=t_d)
+    par_sigma_range = 1
+    pst.add_pars_as_obs(pst_path=t_d,par_sigma_range=par_sigma_range)
     diff = set(pst.par_names) - set(pst.obs_names)
     print(diff)
     assert len(diff) == 0
 
-    pst.write(os.path.join(t_d,"pest.pst"))
+    pst.write(os.path.join(t_d,"pest.pst"),version=2)
 
     pyemu.os_utils.run("{0} {1}".format(ies_exe_path,"pest.pst"),cwd=t_d)
     pst = pyemu.Pst(os.path.join(t_d,"pest.pst"))
     print(pst.phi)
     assert np.isclose(pst.phi, 0), pst.phi
+    par = pst.parameter_data
+    obs = pst.observation_data.loc[par.parnme,:]
+    assert len(par) == len(obs)
+    print(obs)
+    assert np.isclose(np.abs((obs["less_than"] - par["parubnd"]).sum()),0)
+    assert np.isclose(np.abs((obs["greater_than"] - par["parlbnd"]).sum()), 0)
+
+    logidx = par.partrans == "log"
+    obs.loc[logidx,"greater_than"] = np.log10(obs.loc[logidx,"greater_than"])
+    obs.loc[logidx, "less_than"] = np.log10(obs.loc[logidx, "less_than"])
+    obs["dist"] = obs.less_than - obs.greater_than
+    print(obs.dist)
+    obs.dist / par_sigma_range
+    obs.loc[logidx,"dist"] = 10**(obs.loc[logidx,"dist"])
+    weight = 1/(obs.dist/par_sigma_range)
+    print(weight)
+    print(obs.weight)
+    assert np.isclose(np.abs(weight-obs.weight).sum(),0)
+
 
 
 def draw_consistency_test(tmp_path):

autotest/pst_tests.py

@@ -112,6 +112,10 @@ def res_test(tmp_path):
 
 def pst_manip_test(tmp_path):
     import os
+
+    if os.path.exists(tmp_path):
+        shutil.rmtree(tmp_path)
+    os.makedirs(tmp_path)
     from pyemu import Pst
     pst_dir = os.path.join("pst")
     org_path = os.path.join(pst_dir, "pest.pst")
@@ -127,9 +131,46 @@ def pst_manip_test(tmp_path):
     new_pst = Pst(new_path)
     assert all(new_pst.observation_data.obsval == pst.observation_data.obsval)
     assert all(new_pst.parameter_data.parval1 == pst.parameter_data.parval1)
-
+    org_par = pst.parameter_data.copy()
+    pst.dialate_par_bounds(1.0)
+    diff = np.abs(org_par.parubnd - pst.parameter_data.parubnd).sum()
+    assert diff < 1e-7
+    diff = np.abs(org_par.parlbnd - pst.parameter_data.parlbnd).sum()
+    assert diff < 1e-7
+    pst.dialate_par_bounds(1.0)
+    diff = np.abs(org_par.parubnd - pst.parameter_data.parubnd).sum()
+    assert diff < 1e-7
+    diff = np.abs(org_par.parlbnd - pst.parameter_data.parlbnd).sum()
+    assert diff < 1e-7
+    pst.dialate_par_bounds(2.0)
+    print(pst.parameter_data.loc[:,["parubnd","parubnd_org"]])
+    pst.write(new_path)
+    pst = Pst(new_path)
+    pst.write(new_path, version=2)
+    new_pst = Pst(new_path)
+    new_par = new_pst.parameter_data
+    assert np.all(new_par.parubnd.values > org_par.parubnd.values)
+    assert np.all(new_par.parlbnd.values < org_par.parlbnd.values)
+
+    pst.dialate_par_bounds(0.5)
+    new_par = pst.parameter_data
+    #print(new_par.parubnd)
+    #print(org_par.parubnd)
+    assert np.isclose(np.abs(new_par.parubnd - org_par.parubnd).sum(),0)
+    assert np.isclose(np.abs(new_par.parlbnd - org_par.parlbnd).sum(), 0)
+
+    pst.parameter_data["parval1"] = pst.parameter_data.parubnd
+    pst.dialate_par_bounds(2.0)
+    pst.dialate_par_bounds(0.5)
+    new_par = pst.parameter_data
+    assert np.isclose(np.abs(new_par.parubnd - org_par.parubnd).sum(), 0)
+    assert np.isclose(np.abs(new_par.parlbnd - org_par.parlbnd).sum(), 0)
 
 def load_test(tmp_path):
+    import os
+    if os.path.exists(tmp_path):
+        shutil.rmtree(tmp_path)
+    os.makedirs(tmp_path)
     from pyemu import Pst
     pst_dir = setup_tmp("pst", tmp_path)
     # just testing all sorts of different pst files
@@ -177,6 +218,14 @@ def load_test(tmp_path):
                 exceptions.append(pst_file + " v2 reload fail: " + str(e))
                 continue
 
+            org_par = p.parameter_data.copy()
+            p.dialate_par_bounds(1.0)
+            diff = np.abs(org_par.parubnd - p.parameter_data.parubnd).sum()
+            assert diff < 1e-5
+            diff = np.abs(org_par.parlbnd - p.parameter_data.parlbnd).sum()
+            print(diff)
+            assert diff < 1e-5
+
 
     # with open("load_fails.txt",'w') as f:
     #    [f.write(pst_file+'\n') for pst_file in load_fails]
@@ -1533,7 +1582,9 @@ def interface_check_test():
     with this.
     """
     d = 'temp'
-    parrep_test(d)
+    #load_test(d)
+    pst_manip_test(d)
+    #parrep_test(d)
     #interface_check_test()
     # new_format_test_2()
     #write2_nan_test()

autotest/pst_tests_2.py

@@ -1098,6 +1098,11 @@ def results_mou_1_test():
         #print(df)
         assert df is not None
 
+def dialate_bound_test():
+    import pyemu
+
+
+
 if __name__ == "__main__":
     results_ies_3_test()
     results_ies_1_test()