added GPR based emulation helper

autotest/la_tests.py

@@ -741,9 +741,10 @@ def dsi_normscoretransform_test():
 
 if __name__ == "__main__":
     #dsi_normscoretransform_test()
+    ends_freyberg_test("temp")
+    #ends_freyberg_dsi_test("temp")
     #ends_freyberg_dev()
     #ends_freyberg_dsi_test("temp")
-    ends_freyberg_dsi_ztz_test('temp')
     #plot_freyberg_dsi()
     #obscomp_test()
     #alternative_dw()

autotest/pst_from_tests.py

@@ -4744,6 +4744,7 @@ def list_float_int_index_test(tmp_path):
     assert np.isclose(diff,bparval1).all(), diff.loc[~np.isclose(diff,bparval1)]
 
 
+#def mf6_freyberg_thresh_invest(setup_freyberg_mf6):
 def mf6_freyberg_thresh_test(tmp_path):
 
     import numpy as np
@@ -5168,6 +5169,7 @@ def plot_thresh(m_d):
                     prarr = np.zeros((nrow,ncol)) - 1
                     prarr[kobs.i,kobs.j] = pr_oe.loc[real,kobs.obsnme]
                     prarr[ib==0] = np.nan
+                    print(pt_oe)
                     ptarr = np.zeros((nrow, ncol)) - 1
                     ptarr[kobs.i, kobs.j] = pt_oe.loc[real, kobs.obsnme]
                     ptarr[ib == 0] = np.nan
@@ -5719,8 +5721,8 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     import flopy
     # except:
     #     return
-    # import sys
-    # sys.path.insert(0,os.path.join("..","..","pypestutils"))
+    import sys
+    sys.path.insert(0,os.path.join("..","..","pypestutils"))
 
     import pypestutils as ppu
 
@@ -5785,7 +5787,7 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     ymn = m.modelgrid.yvertices.min()
     ymx = m.modelgrid.yvertices.max()
 
-    numpp = 20
+    numpp = 30
     xvals = np.random.uniform(xmn,xmx,numpp)
     yvals = np.random.uniform(ymn, ymx, numpp)
     pp_locs = pd.DataFrame({"x":xvals,"y":yvals})
@@ -5965,9 +5967,9 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     assert cat1par.shape[0] == num_cat_arrays
     assert cat2par.shape[0] == num_cat_arrays
 
-    par.loc[cat1par, "parval1"] = 0.8
+    par.loc[cat1par, "parval1"] = 0.9
     par.loc[cat1par, "parubnd"] = 1.0
-    par.loc[cat1par, "parlbnd"] = 0.6
+    par.loc[cat1par, "parlbnd"] = 0.8
     par.loc[cat1par,"partrans"] = "none"
 
     # since the apply method only looks that first proportion, we can just fix this one
@@ -5984,7 +5986,7 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     print(pst.npar,pst.npar_adj)
 
     org_par = par.copy()
-    num_reals = 100
+    num_reals = 200
     pe = pf.draw(num_reals, use_specsim=False)
     pe.enforce()
     print(pe.shape)
@@ -6046,14 +6048,17 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     # reset away from the truth...
     pst.parameter_data.loc[:,"parval1"] = org_par.parval1.values.copy()
 
-    pst.control_data.noptmax = 1
+    pst.control_data.noptmax = 13
     pst.pestpp_options["ies_par_en"] = "prior.jcb"
-    pst.pestpp_options["ies_num_reals"] = 30
+    pst.pestpp_options["ies_num_reals"] = 200
     pst.pestpp_options["ies_subset_size"] = -10
     pst.pestpp_options["ies_no_noise"] = True
-    #pst.pestpp_options["ies_bad_phi_sigma"] = 2.0
+    pst.pestpp_options["ies_bad_phi_sigma"] = 1.5
     pst.pestpp_options["overdue_giveup_fac"] = 100.0
     pst.pestpp_options["ies_multimodal_alpha"] = 0.99
+    pst.pestpp_options["ies_n_iter_reinflate"] = 3
+
+
     #pst.pestpp_options["panther_agent_freeze_on_fail"] = True
 
     #pst.write(os.path.join(pf.new_d, "freyberg.pst"))
@@ -6092,16 +6097,17 @@ def mf6_freyberg_ppu_hyperpars_thresh_invest(tmp_path):
     # assert phidf["mean"].min() < phidf["mean"].max()
 
 
+
 if __name__ == "__main__":
     #mf6_freyberg_pp_locs_test('.')
     #mf6_subdir_test(".")
-    mf6_freyberg_ppu_hyperpars_invest(".")
-    # mf6_freyberg_ppu_hyperpars_thresh_invest(".")
+    #mf6_freyberg_ppu_hyperpars_invest(".")
+    #mf6_freyberg_ppu_hyperpars_thresh_invest(".")
     #while True:
     #    mf6_freyberg_thresh_test(".")
     plot_thresh("master_thresh_nonstat")
     #mf6_freyberg_thresh_test(".")
-    plot_thresh("master_thresh_nonstat")
+    #plot_thresh("master_thresh_nonstat")
     #plot_thresh("master_thresh_nonstat_nim")
 
     # invest()

autotest/utils/constr_template/additive_par.dat

@@ -0,0 +1,4 @@
+obj1_add_par 0.00000000000000e+00
+obj2_add_par 0.00000000000000e+00
+constr1_add_par 0.00000000000000000e+00
+constr2_add_par 0.00000000000000000e+00

autotest/utils/constr_template/additive_par.dat.tpl

@@ -0,0 +1,5 @@
+ptf ~
+obj1_add_par ~   obj1_add_par   ~
+obj2_add_par ~   obj2_add_par   ~
+constr1_add_par ~   constr1_add_par   ~
+constr2_add_par ~   constr2_add_par   ~

autotest/utils/constr_template/constr.pst

@@ -0,0 +1,41 @@
+pcf
+* control data
+             restart          estimation
+         6         2         2         2         2
+         2         1              single               point         1
+   2.000000E+01  -3.000000E+00   3.000000E-01   1.000000E-02        -7
+   1.000000E+01   1.000000E+01   1.000000E-03
+   1.000000E-01
+         0   1.000000E-02         3         3   1.000000E-02         3
+         0         0         0
+* singular value decomposition
+         1
+  10000000    1.000000E-06
+1
+* parameter groups
+obj_add                absolute              1.0000000000E-02      0.0 switch                2.0000000000E+00      parabolic             1.0000000000E-05      5.0000000000E-01      smaller              
+decvars                relative              1.0000000000E-02      0.0 switch                2.0000000000E+00      parabolic             1.0000000000E-05      5.0000000000E-01      smaller              
+* parameter data
+dv_0                  none                  relative              5.0000000000E-01      1.0000000000E-01      1.0000000000E+00      decvars               1.0000000000E+00      0.0000000000E+00      1          
+dv_1                  none                  relative              2.5000000000E+00      0.0000000000E+00      5.0000000000E+00      decvars               1.0000000000E+00      0.0000000000E+00      1          
+constr1_add_par       none                  relative              0.0000000000E+00      -5.0000000000E-01     5.0000000000E-01      obj_add               1.0000000000E+00      0.0000000000E+00      1          
+constr2_add_par       none                  relative              0.0000000000E+00      -5.0000000000E-01     5.0000000000E-01      obj_add               1.0000000000E+00      0.0000000000E+00      1          
+obj1_add_par          none                  relative              0.0000000000E+00      -5.0000000000E-01     5.0000000000E-01      obj_add               1.0000000000E+00      0.0000000000E+00      1          
+obj2_add_par          none                  relative              0.0000000000E+00      -5.0000000000E-01     5.0000000000E-01      obj_add               1.0000000000E+00      0.0000000000E+00      1          
+* observation groups
+less_than            
+greater_than         
+* observation data
+obj_1                 5.0000000000E-01      1.0000000000E+00      less_than            
+obj_2                 3.0000000000E+00      1.0000000000E+00      less_than            
+* model command line
+python forward_run.py
+* model input/output
+./dv.dat.tpl ./dv.dat
+./additive_par.dat.tpl ./additive_par.dat
+./obj.dat.ins ./obj.dat
+* prior information
+const_1               9.0 * dv_0 + 1.0 * dv_1 = 6.0 1.0000000000E+00     greater_than         
+const_2               9.0 * dv_0 - 1.0 * dv_1 = 1.0 1.0000000000E+00     greater_than         
+++opt_dec_var_groups(decvars)
+++mou_objectives(obj_1,obj_2)

autotest/utils/constr_template/dv.dat

@@ -0,0 +1,2 @@
+dv_0 0.5000000000
+dv_1 2.5000000000

autotest/utils/constr_template/dv.dat.tpl

@@ -0,0 +1,3 @@
+ptf ~
+dv_0 ~ dv_0     ~
+dv_1 ~ dv_1     ~

autotest/utils/constr_template/forward_run.py

@@ -0,0 +1,29 @@
+import os
+import numpy as np
+import pandas as pd
+def constr(x):
+    return (x[0],(1 + x[1]) / x[0]),[]
+
+def helper(func):
+    pdf = pd.read_csv("dv.dat",delim_whitespace=True,index_col=0, header=None, names=["parnme","parval1"])
+    #obj1,obj2 = func(pdf.values)
+    objs,constrs = func(pdf.values)
+
+    if os.path.exists("additive_par.dat"):
+        obj1,obj2 = objs[0],objs[1]
+        cdf = pd.read_csv("additive_par.dat", delim_whitespace=True, index_col=0,header=None, names=["parnme","parval1"])
+        obj1[0] += cdf.parval1.values[0]
+        obj2[0] += cdf.parval1.values[1]
+        for i in range(len(constrs)):
+            constrs[i] += cdf.parval1.values[i+2]
+
+    with open("obj.dat",'w') as f:
+        for i,obj in enumerate(objs):
+            f.write("obj_{0} {1}\n".format(i+1,float(obj)))
+        #f.write("obj_2 {0}\n".format(float(obj2)))
+        for i,constr in enumerate(constrs):
+            f.write("constr_{0} {1}\n".format(i+1,float(constr)))
+    return objs,constrs
+
+if __name__ == '__main__':
+    helper(constr)

autotest/utils/constr_template/obj.dat

@@ -0,0 +1,2 @@
+obj_1 0.5
+obj_2 7.0

autotest/utils/constr_template/obj.dat.ins

@@ -0,0 +1,3 @@
+pif ~
+l1 w !obj_1!
+l1 w !obj_2!

autotest/utils/hosaki_template/forward_run.py

@@ -0,0 +1,12 @@
+import pandas as pd
+import numpy as np
+pdf = pd.read_csv('par.csv')
+pval1 = float(pdf.iloc[0,1])
+pval2 = float(pdf.iloc[1,1])
+term1 = (pval2**2)*np.e**(-pval2)
+term2 = 1. - (8.*(pval1)) + (7.*(pval1**2))
+term3 = (-(7./3.)*(pval1**3)) + (0.25*(pval1**4))
+sim = (term2 + term3) * term1
+with open('sim.csv','w') as f:
+    f.write('obsnme,obsval\n')
+    f.write('sim,'+str(sim)+'\n')

autotest/utils/hosaki_template/par.csv

@@ -0,0 +1,3 @@
+parnme,parval1
+par1,4.0000000000
+par2,2.0000000000

autotest/utils/hosaki_template/par.csv.tpl

@@ -0,0 +1,4 @@
+ptf ~
+parnme,parval1
+par1,~   par1   ~
+par2,~   par2   ~

autotest/utils/hosaki_template/pest.pst

@@ -0,0 +1,28 @@
+pcf
+* control data
+             restart          estimation
+         2         1         1         0         1
+         1         1              single               point         1
+   2.000000E+01  -3.000000E+00   3.000000E-01   1.000000E-02        -7
+   1.000000E+01   1.000000E+01   1.000000E-03
+   1.000000E-01
+        25   1.000000E-02         3      1000   1.000000E-02         3
+         0         0         0
+* singular value decomposition
+         1
+  10000000    1.000000E-06
+1
+* parameter groups
+pargp                  relative              1.0000000000E-02      0.0 switch                2.0000000000E+00      parabolic             1.0000000000E-05      5.0000000000E-01      smaller              
+* parameter data
+par1                  none                  relative              2.5000000000E+00      0.0000000000E+00      5.0000000000E+00      pargp                 1.0000000000E+00      0.0000000000E+00      1          
+par2                  none                  relative              2.5000000000E+00      0.0000000000E+00      5.0000000000E+00      pargp                 1.0000000000E+00      0.0000000000E+00      1          
+* observation groups
+obgnme               
+* observation data
+sim                   -2.3458115761E+00     1.0000000000E+02      less_than              
+* model command line
+python forward_run.py
+* model input/output
+./par.csv.tpl ./par.csv
+./sim.csv.ins ./sim.csv

autotest/utils/hosaki_template/sim.csv

@@ -0,0 +1,2 @@
+obsnme,obsval
+sim,-2.345811576101292

autotest/utils/hosaki_template/sim.csv.ins

@@ -0,0 +1,3 @@
+pif ~
+l1
+l1 ~,~ !sim!

autotest/utils/zdt1_template/additive_par.dat

@@ -0,0 +1,2 @@
+obj1_add_par 0.00000000000000e+00
+obj2_add_par 0.00000000000000e+00

autotest/utils/zdt1_template/additive_par.dat.tpl

@@ -0,0 +1,3 @@
+ptf ~
+obj1_add_par ~   obj1_add_par   ~
+obj2_add_par ~   obj2_add_par   ~

autotest/utils/zdt1_template/dv.dat

@@ -0,0 +1,30 @@
+dv_0 0.5000000000
+dv_1 0.5000000000
+dv_2 0.5000000000
+dv_3 0.5000000000
+dv_4 0.5000000000
+dv_5 0.5000000000
+dv_6 0.5000000000
+dv_7 0.5000000000
+dv_8 0.5000000000
+dv_9 0.5000000000
+dv_10 0.50000000000
+dv_11 0.50000000000
+dv_12 0.50000000000
+dv_13 0.50000000000
+dv_14 0.50000000000
+dv_15 0.50000000000
+dv_16 0.50000000000
+dv_17 0.50000000000
+dv_18 0.50000000000
+dv_19 0.50000000000
+dv_20 0.50000000000
+dv_21 0.50000000000
+dv_22 0.50000000000
+dv_23 0.50000000000
+dv_24 0.50000000000
+dv_25 0.50000000000
+dv_26 0.50000000000
+dv_27 0.50000000000
+dv_28 0.50000000000
+dv_29 0.50000000000

autotest/utils/zdt1_template/dv.dat.tpl

@@ -0,0 +1,31 @@
+ptf ~
+dv_0 ~ dv_0     ~
+dv_1 ~ dv_1     ~
+dv_2 ~ dv_2     ~
+dv_3 ~ dv_3     ~
+dv_4 ~ dv_4     ~
+dv_5 ~ dv_5     ~
+dv_6 ~ dv_6     ~
+dv_7 ~ dv_7     ~
+dv_8 ~ dv_8     ~
+dv_9 ~ dv_9     ~
+dv_10 ~ dv_10     ~
+dv_11 ~ dv_11     ~
+dv_12 ~ dv_12     ~
+dv_13 ~ dv_13     ~
+dv_14 ~ dv_14     ~
+dv_15 ~ dv_15     ~
+dv_16 ~ dv_16     ~
+dv_17 ~ dv_17     ~
+dv_18 ~ dv_18     ~
+dv_19 ~ dv_19     ~
+dv_20 ~ dv_20     ~
+dv_21 ~ dv_21     ~
+dv_22 ~ dv_22     ~
+dv_23 ~ dv_23     ~
+dv_24 ~ dv_24     ~
+dv_25 ~ dv_25     ~
+dv_26 ~ dv_26     ~
+dv_27 ~ dv_27     ~
+dv_28 ~ dv_28     ~
+dv_29 ~ dv_29     ~

autotest/utils/zdt1_template/forward_run.py

@@ -0,0 +1,30 @@
+import os
+import numpy as np
+import pandas as pd
+def zdt1(x):
+    g = 1 + 9 * np.sum(x[1:]) / (len(x) - 1)
+    return (x[0]*10, g * (1 - np.sqrt(x[0] / g))*10),[]
+
+def helper(func):
+    pdf = pd.read_csv("dv.dat",delim_whitespace=True,index_col=0, header=None, names=["parnme","parval1"])
+    #obj1,obj2 = func(pdf.values)
+    objs,constrs = func(pdf.values)
+
+    if os.path.exists("additive_par.dat"):
+        obj1,obj2 = objs[0],objs[1]
+        cdf = pd.read_csv("additive_par.dat", delim_whitespace=True, index_col=0,header=None, names=["parnme","parval1"])
+        obj1[0] += cdf.parval1.values[0]
+        obj2[0] += cdf.parval1.values[1]
+        for i in range(len(constrs)):
+            constrs[i] += cdf.parval1.values[i+2]
+
+    with open("obj.dat",'w') as f:
+        for i,obj in enumerate(objs):
+            f.write("obj_{0} {1}\n".format(i+1,float(obj)))
+        #f.write("obj_2 {0}\n".format(float(obj2)))
+        for i,constr in enumerate(constrs):
+            f.write("constr_{0} {1}\n".format(i+1,float(constr)))
+    return objs,constrs
+
+if __name__ == '__main__':
+    helper(zdt1)

autotest/utils/zdt1_template/obj.dat

@@ -0,0 +1,2 @@
+obj_1 5.0
+obj_2 38.416876048223

autotest/utils/zdt1_template/obj.dat.ins

@@ -0,0 +1,3 @@
+pif ~
+l1 w !obj_1!
+l1 w !obj_2!