diff --git a/autotest/la_tests.py b/autotest/la_tests.py index 0fdfa3ea0..fdc11461c 100644 --- a/autotest/la_tests.py +++ b/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() diff --git a/autotest/pst_from_tests.py b/autotest/pst_from_tests.py index 1a280a9af..b1cc8d671 100644 --- a/autotest/pst_from_tests.py +++ b/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() diff --git a/autotest/utils/constr_template/additive_par.dat b/autotest/utils/constr_template/additive_par.dat new file mode 100644 index 000000000..ca5916f8a --- /dev/null +++ b/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 diff --git a/autotest/utils/constr_template/additive_par.dat.tpl b/autotest/utils/constr_template/additive_par.dat.tpl new file mode 100644 index 000000000..e8e9cc5fa --- /dev/null +++ b/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 ~ diff --git a/autotest/utils/constr_template/constr.pst b/autotest/utils/constr_template/constr.pst new file mode 100644 index 000000000..c7f97f515 --- /dev/null +++ b/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) diff --git a/autotest/utils/constr_template/dv.dat b/autotest/utils/constr_template/dv.dat new file mode 100644 index 000000000..ac9ebec51 --- /dev/null +++ b/autotest/utils/constr_template/dv.dat @@ -0,0 +1,2 @@ +dv_0 0.5000000000 +dv_1 2.5000000000 diff --git a/autotest/utils/constr_template/dv.dat.tpl b/autotest/utils/constr_template/dv.dat.tpl new file mode 100644 index 000000000..138fe8fe4 --- /dev/null +++ b/autotest/utils/constr_template/dv.dat.tpl @@ -0,0 +1,3 @@ +ptf ~ +dv_0 ~ dv_0 ~ +dv_1 ~ dv_1 ~ diff --git a/autotest/utils/constr_template/forward_run.py b/autotest/utils/constr_template/forward_run.py new file mode 100644 index 000000000..155fbef66 --- /dev/null +++ b/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) diff --git a/autotest/utils/constr_template/obj.dat b/autotest/utils/constr_template/obj.dat new file mode 100644 index 000000000..fba583782 --- /dev/null +++ b/autotest/utils/constr_template/obj.dat @@ -0,0 +1,2 @@ +obj_1 0.5 +obj_2 7.0 diff --git a/autotest/utils/constr_template/obj.dat.ins b/autotest/utils/constr_template/obj.dat.ins new file mode 100644 index 000000000..f2e2a1b71 --- /dev/null +++ b/autotest/utils/constr_template/obj.dat.ins @@ -0,0 +1,3 @@ +pif ~ +l1 w !obj_1! +l1 w !obj_2! diff --git a/autotest/utils/hosaki_template/forward_run.py b/autotest/utils/hosaki_template/forward_run.py new file mode 100644 index 000000000..45c8d6939 --- /dev/null +++ b/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') diff --git a/autotest/utils/hosaki_template/par.csv b/autotest/utils/hosaki_template/par.csv new file mode 100644 index 000000000..f948ea6b2 --- /dev/null +++ b/autotest/utils/hosaki_template/par.csv @@ -0,0 +1,3 @@ +parnme,parval1 +par1,4.0000000000 +par2,2.0000000000 diff --git a/autotest/utils/hosaki_template/par.csv.tpl b/autotest/utils/hosaki_template/par.csv.tpl new file mode 100644 index 000000000..d6c871149 --- /dev/null +++ b/autotest/utils/hosaki_template/par.csv.tpl @@ -0,0 +1,4 @@ +ptf ~ +parnme,parval1 +par1,~ par1 ~ +par2,~ par2 ~ diff --git a/autotest/utils/hosaki_template/pest.pst b/autotest/utils/hosaki_template/pest.pst new file mode 100644 index 000000000..03982ea0c --- /dev/null +++ b/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 diff --git a/autotest/utils/hosaki_template/sim.csv b/autotest/utils/hosaki_template/sim.csv new file mode 100644 index 000000000..1ee6b2ed7 --- /dev/null +++ b/autotest/utils/hosaki_template/sim.csv @@ -0,0 +1,2 @@ +obsnme,obsval +sim,-2.345811576101292 diff --git a/autotest/utils/hosaki_template/sim.csv.ins b/autotest/utils/hosaki_template/sim.csv.ins new file mode 100644 index 000000000..534dcae96 --- /dev/null +++ b/autotest/utils/hosaki_template/sim.csv.ins @@ -0,0 +1,3 @@ +pif ~ +l1 +l1 ~,~ !sim! diff --git a/autotest/utils/zdt1_template/additive_par.dat b/autotest/utils/zdt1_template/additive_par.dat new file mode 100644 index 000000000..919332b7c --- /dev/null +++ b/autotest/utils/zdt1_template/additive_par.dat @@ -0,0 +1,2 @@ +obj1_add_par 0.00000000000000e+00 +obj2_add_par 0.00000000000000e+00 diff --git a/autotest/utils/zdt1_template/additive_par.dat.tpl b/autotest/utils/zdt1_template/additive_par.dat.tpl new file mode 100644 index 000000000..2b2a325f7 --- /dev/null +++ b/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 ~ diff --git a/autotest/utils/zdt1_template/dv.dat b/autotest/utils/zdt1_template/dv.dat new file mode 100644 index 000000000..6d7efddd6 --- /dev/null +++ b/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 diff --git a/autotest/utils/zdt1_template/dv.dat.tpl b/autotest/utils/zdt1_template/dv.dat.tpl new file mode 100644 index 000000000..e105ab03a --- /dev/null +++ b/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 ~ diff --git a/autotest/utils/zdt1_template/forward_run.py b/autotest/utils/zdt1_template/forward_run.py new file mode 100644 index 000000000..4a84cd903 --- /dev/null +++ b/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) diff --git a/autotest/utils/zdt1_template/obj.dat b/autotest/utils/zdt1_template/obj.dat new file mode 100644 index 000000000..95a3a3fd4 --- /dev/null +++ b/autotest/utils/zdt1_template/obj.dat @@ -0,0 +1,2 @@ +obj_1 5.0 +obj_2 38.416876048223 diff --git a/autotest/utils/zdt1_template/obj.dat.ins b/autotest/utils/zdt1_template/obj.dat.ins new file mode 100644 index 000000000..f2e2a1b71 --- /dev/null +++ b/autotest/utils/zdt1_template/obj.dat.ins @@ -0,0 +1,3 @@ +pif ~ +l1 w !obj_1! +l1 w !obj_2! diff --git a/autotest/utils/zdt1_template/prior.csv b/autotest/utils/zdt1_template/prior.csv new file mode 100644 index 000000000..d06a2b6d6 --- /dev/null +++ b/autotest/utils/zdt1_template/prior.csv @@ -0,0 +1,101 @@ +,dv_0,dv_1,dv_10,dv_11,dv_12,dv_13,dv_14,dv_15,dv_16,dv_17,dv_18,dv_19,dv_2,dv_20,dv_21,dv_22,dv_23,dv_24,dv_25,dv_26,dv_27,dv_28,dv_29,dv_3,dv_4,dv_5,dv_6,dv_7,dv_8,dv_9,obj2_add_par,obj1_add_par +0,0.3892082938839153,0.3375828863010191,0.5670722890739773,0.41811967939927597,0.7565874566173032,0.2453511245811651,0.6795951911292893,0.44728811160979076,0.2410500605204497,0.4672732299516216,0.19837541667560515,0.32166473648364013,0.3974635246113443,0.11288810475072303,0.2513416886268678,0.49769070682202976,0.7331072119077706,0.3614650678871929,0.8442477094182492,0.05123069484673154,0.5819633794688204,0.059304732101717716,0.524400627797456,0.4560375282936565,0.7321070585074505,0.2640676741047735,0.33164821774213066,0.5850861029940552,0.7287958975791088,0.38324982181933753,-0.3023285715257477,-0.15228415437190174 +1,0.49860477123020097,0.8597184502273669,0.4341822318447873,0.6414635592291258,0.7287795781837912,0.9198310666107179,0.41121271428093276,0.2677372144392737,0.42279032497528635,0.3876080962364676,0.40314348492772567,0.11806073626273667,0.27229469681118057,0.5421258063912178,0.5308351087448286,0.9245669502532601,0.4248429486787485,0.2331464471777333,0.17382927841531298,0.9244690726367231,0.6051735681669173,0.6103294890804305,0.3479211474224788,0.7871296709399045,0.48425157638476285,0.17452505059692486,0.8142230999713655,0.4762712565015677,0.30727471066141143,0.4074198004899751,0.04485217062225495,0.06515572152952187 +2,0.8331796906669864,0.9097659006914592,0.24241214777773729,0.47795858007778963,0.303176106448553,0.966041959329597,0.6247825885527591,0.5377395761823779,0.22747136102215126,0.29388345730011345,0.6390219158094586,1.0,0.4562491834001192,0.6752505260176712,0.3160648925928623,0.6051672752031345,0.5582536074663861,0.4247346152532194,0.5242749066097406,0.7775635755947452,0.3418311738454445,0.0,0.6024206302982346,0.6005177459869007,1.0,0.46323288098600635,0.9619152475523862,0.4124541899738293,0.5175176886645153,0.27826954166788076,0.32433118923346904,0.2085834583724079 +3,0.7540493774241512,0.6402212511157729,0.2612423276281949,0.734526252143329,0.2730599582496569,0.8471537060110217,0.644559698429711,0.4619720266347997,0.47801486339103205,0.20697673523773646,0.3984114996267316,0.5738027090766014,0.42539192653079694,0.4318783597908865,0.9864995730617958,0.0,0.0,0.25610508361965123,0.47629695746091166,0.7511145735176797,0.32527549309691856,0.5096046026871508,0.5954828908450752,0.5109657979072918,0.36899041437417884,0.7993766605015611,0.851179035593284,0.7157508929128341,0.6191142347570509,0.32909709909093265,0.27835025405168706,0.2324613314535705 +4,0.5335337291818187,0.4321697008395472,0.23065148612786962,0.441094370911613,0.4889878534722329,0.6385345543289249,0.4958514851767654,1.0,0.49535986005202765,0.4740369646268914,0.21169997281484754,0.5271579277109154,0.2449339816916956,0.5904049013436724,0.5989478966207956,0.8390022394814635,0.2986211396305482,0.4086218979906134,0.6189476936568237,0.40331589137226176,0.0,0.3227853881422331,0.44258550413831266,0.7407079913514822,0.543924778879736,0.43532078281150627,0.3321109950275374,0.6192804140923938,0.5875005838289535,0.4075891851242481,0.37236447947543017,0.017126556229116054 +5,0.5645202168079191,0.07784867926333527,0.5760659953558616,0.5358000204688023,0.7806685482896993,0.7048785059290386,0.5449642585790735,0.5762478576714267,0.8854527280349507,0.47084418384829674,0.3708225385659071,0.409757274046699,0.8758212021727222,0.418792976610261,0.3845966881581217,0.949260734224771,0.5028146675775068,0.5592855470993762,0.6716440135282976,1.0,0.45408674861357123,0.683980024594848,0.7535671131301261,0.9232162807899797,0.8947886588151286,0.6559931066422795,0.5875113960683169,0.5165763998640986,0.6679639974508111,0.2695714250534415,-0.07103281756124605,0.4235275001049192 +6,0.3341250491707543,0.4388821857511083,0.7947597130382908,0.5179701406952957,0.3605552764657207,0.7205446703057583,0.5428399404356787,0.3437243797139914,0.7560300855633973,0.652358887949782,0.2632372055408294,0.28538551967782677,0.0,0.39339102964250583,0.6845543543028759,0.7955032854467468,0.3212369918416569,0.42411765215075964,0.7483905823028056,1.0,0.16193887390725803,0.8394071943333662,0.5757218572456744,0.4955616566819399,0.8613264737014997,0.1422187121833337,0.3694052370992519,0.4923846389954708,0.5748977014489525,0.8263104653206552,-0.29753378700923455,-0.2934422021551347 +7,0.6536089055896674,0.40013145470081857,0.3946638210576323,0.7864181895086053,0.2743916974917528,0.363027379020211,0.26923047476705164,0.5995273025462396,0.0,0.2524431961407231,0.5164176631397807,0.208898214771716,0.41140876637733903,0.41008456996637765,0.26587755670511415,0.4764601487483848,0.9788839135772542,0.4847309627546741,0.7230022129428118,0.48428130991354784,0.4132155555502643,1.0,0.6935360917945338,0.5072369316934332,0.489798957126783,0.7180332878126382,0.3099646237685028,0.3702935181017991,0.27488677706010833,0.3561433289419259,0.13740784948683996,-0.17234492280400435 +8,0.7670879656944485,0.37194335900627074,0.43046617725916636,1.0,0.0710953496355366,0.36026141697821423,0.1800638308957936,0.8074115196442517,0.14818254427834382,0.5476103929154662,0.3912998226130874,0.4574895097078723,0.9206077411914755,0.7574243375102457,0.5896473668081459,0.2798564408999016,0.6009593500038238,0.6132647102577555,0.18837193090210452,0.533227664599071,0.51994698729107,0.12752581225164483,0.6734080918983045,0.47142109443784197,0.7408924100066916,0.5008344655745653,0.46737964687719324,0.5365757811887085,0.5265178061127661,0.866385689782705,0.26339449816473576,-0.5 +9,0.12451274396410578,0.28361915447657055,0.5595327019841656,0.7359335343835858,0.5240283961903607,0.3328917264348144,0.8904963667566337,0.43804036653920375,0.49347474560079746,0.12804432034940294,0.4775406852970878,1.0,0.13183882254010337,0.47263145581130506,0.0,0.25034128562368835,0.5823248305751466,0.28524163707529754,0.20737937801119483,0.5751460515131224,0.7019482672037589,0.37866697235828944,0.6608174690281619,0.02544117841825494,0.4455886115988362,0.6531079263896704,0.48054505524833924,0.9025829879459694,0.6039810837899053,0.0,0.0913804198827878,0.11870678486447396 +10,0.33830923570150107,0.3173262648378088,0.45833985704955044,0.5703004611032456,0.6714444568090626,0.4205535124345735,0.9589549146283434,0.5040506654519658,0.7327874777672351,0.4509459986412378,0.6057473357777341,0.4742985728890859,0.4633227417662734,0.4649788203998686,0.7619486026281899,0.4236944852550666,0.6521812052198597,0.597740498695332,0.3780182950154965,0.46005728508034605,0.5051043637234642,0.0,0.6149362123870789,0.3374653509406202,0.38092069537003287,0.2796945467966424,1.0,1.0,0.5476096515586885,0.3868025475090833,-0.0996402036107076,0.4180054805654091 +11,0.8058587547680536,0.6291254408101858,0.8120519320780579,0.5154230505553101,0.4784048404948046,0.2897047839899436,0.7326133730977505,0.32119371821983234,0.8034236595380826,0.5367276839860672,0.7460322366737184,0.9662301500719104,0.28519934410424325,0.52020217861093,0.2178544243011078,0.5628237511392192,0.4786405047187286,0.3349002493060763,0.41744924627635094,0.0,0.7716665777472929,0.5931180817853458,0.49527415719614526,0.5316943801598019,0.18898073470224197,0.10193895749456472,0.2833685710299606,0.2275622257146489,0.513649153391419,0.3982587404005349,0.09440198149812346,-0.24791299863093325 +12,0.0,0.7504969508679334,1.0,0.2255870016156773,0.40953629410935294,0.7627553667045804,0.5822906131831478,0.5678033278738454,0.49630284757067483,0.8345429448309301,1.0,0.6014420886369851,0.4827117055969254,0.12623367321648776,0.3158675878440177,0.5460472833693492,0.5377962199948906,0.8689087911401919,0.38544634316828463,0.7285834484138032,0.25773111647956454,0.3107255226305865,0.6523612784623033,0.6602275189187137,0.13188276385326758,0.5924623483078405,0.8785080325706076,0.5277812852580482,0.5827073075657925,0.8343971649830082,0.30208084528736334,-0.22852642593654074 +13,0.795500336423606,0.12493325776239689,0.24939498764342793,0.7416695692466972,0.5148774727429922,0.7382378962418136,0.0,0.8672708610854626,0.3795730271245865,0.5268355689394022,0.45461967296766215,0.1547570107756006,0.6239328774244969,0.25278193181219755,0.5626196122776804,0.0,0.42764891952631345,0.6987325054464677,0.7028406480056061,0.49774327608074215,0.19056911102220403,0.8536376227430016,0.5581197669881124,0.8316514917675668,0.842712505766195,0.4320832697150728,0.060357435952358374,0.31803902815984486,0.19631571965501549,0.7928048170925068,-0.33750007384831227,0.3403498718245541 +14,0.7379498301332585,0.6673380580563615,0.19773053364138998,0.7214449952399258,0.0,0.6130316372866449,0.0038527661509922817,0.04440267692772093,0.05005249471943246,0.4892493481115953,0.16404118944826823,0.48757382681962014,0.801894778654141,0.23811584602637836,0.4393526718337134,0.8545657234354698,0.44792116253508896,0.3127050930663382,0.37955684393217465,0.17636020013463982,0.32072500436200707,0.3636270468003992,0.22597901478981958,0.17041700615023475,0.4934081513520084,0.3491330425303674,0.5396887041762903,0.5503514882229776,0.4795821457169112,0.1315443087417586,-0.030573743120962248,0.14592109804995004 +15,0.748084740798811,0.6180694663765297,0.7415504203754109,0.4463027002666346,0.6368166523703016,0.7557057079459595,0.5052012761733147,0.5150055745874943,0.6992523959596326,0.3135520839589955,0.4491293429311124,0.4806465223624711,0.07812144858495401,0.36155490780009947,0.8659251586999284,0.30207985130918424,0.5931735266852775,0.20223693605751536,0.0,0.5171552835518247,0.5372403880775465,0.5701516853067787,0.1758786502244271,0.611941275174389,0.846264276458776,0.7066246467109801,0.6297159293502745,0.05182027097298614,0.4034028125908979,0.677928666797995,0.07247488416310034,-0.2563812292709337 +16,0.41236178997160844,0.38648660909516663,0.8519867158705907,0.47685522152829246,0.296511586297715,0.436729680930284,0.23951117678677764,0.8535998530817572,0.2329739410940327,0.30588594385543244,0.46816644595333434,0.7915124139088271,0.36160100465254696,0.717955624742227,0.5743803736880569,0.5559495446858493,0.19997421017753697,0.34868300036529315,0.9077771521306222,0.027104432249938193,0.22228855980997825,0.2779709188873005,1.0,0.5622174136586668,0.02934134645982983,0.834195093319928,0.5296472087734173,0.39200418661588005,0.8792580293405666,0.0,0.4859622865979145,-0.18369548915930567 +17,0.5328978607736337,0.4732779755887909,1.0,0.25137915809865585,0.8072375446161749,0.3809178103978596,0.395532284091509,0.2320100349397471,1.0,0.5968881836196044,0.36254936422246387,0.8106367585022811,0.7445466974328384,0.3855064934971697,0.11121748104119383,0.6916622130666826,0.0942884540399479,0.3254419705326168,0.08691476889728517,0.7735555026032835,0.6366101200324867,0.494192344649162,0.58956467704313,0.32380347756600814,0.35055245697952964,0.702008753291864,0.40198710791835246,0.6118577376325203,0.5663915621353168,0.2556789393833056,0.05700002409178315,0.0034468263122778586 +18,0.7199196082428148,0.4329907803443198,0.36503295348213916,0.3619353830234108,0.21786738820350382,0.3310664919136429,0.7943274109853706,0.049812105671942,0.4517376430532292,0.3261438961818016,0.1833975952307922,0.42517591977689523,0.7917388035454074,0.8204159910341093,0.8999319497432934,0.9560323782605618,0.82514606522103,0.2533378716713788,0.5602708826644748,0.45795515165558853,0.6599052824754135,0.5825568209990007,0.41168237485040926,0.7363888944911964,0.5450692166390727,0.24177121274756752,0.663220911883525,0.34248176451133466,0.6424220506167746,0.9682417835068879,-0.20264633817086003,0.2164190070306167 +19,0.0,1.0,0.0,0.37597687846561656,0.7387961630236837,0.0,0.5691226622885586,0.5582862751353072,0.6700642842531028,0.30271821124862397,0.33219705937398347,0.0033296102796739335,0.1156519002515427,0.8316038279019453,0.7808604639666346,0.43554555908892834,0.794219807362811,0.44805384745050203,0.6857637137558468,0.7617002442220826,0.40549412161419157,0.25487096988388735,0.5219939447364226,0.5414257986963572,0.6251605038738586,0.5483961739241178,0.5035768652312077,0.6604123153190306,0.294234356120059,0.894118082857759,-0.04103775198612426,-0.2332160964724896 +20,0.35705438847685234,0.8334425476599991,0.4544647592791773,1.0,0.0,0.6535906461950488,0.3966000968374976,0.6649855421789383,0.0009190046821864195,0.6464086463045409,0.5036867309034176,0.5086256757613461,0.5062099649264263,0.11839357208540147,0.5850910428142738,0.34063385873268687,0.5070706592616145,0.4025524544502572,0.5830333954543296,0.30290317498409725,0.6222196834555802,0.48252220697325027,0.1543339999941445,0.7455136146360926,0.5870409939994357,0.5184066911328433,0.4680668066186868,0.5678093173854378,0.33438425974421027,0.39868173103188614,0.5,0.08048966312177333 +21,1.0,0.13494405954391642,0.4971997582224016,0.3906536487450635,0.30441479767606644,0.5512683943488449,0.6789453740337835,0.838392585068372,0.4895285734326481,0.7720508088946789,0.7295623407041373,0.0,0.32060257913276013,0.510660463470788,0.5844151334733729,0.1700992043979347,0.415372043026878,0.5378170854252597,0.734153402033579,0.21948676094450664,0.5028225748722821,1.0,0.5799605420872622,0.4355427292922164,0.42291383506008806,0.22535822116802862,0.7177808780702303,0.9412387004840927,0.5387231015367203,0.2703199325848515,-0.02693171548902961,0.3350541933981163 +22,0.6665119312750094,0.4041207585501796,0.8171702489202188,0.6192870319477259,0.33923169897816463,0.8033122236683341,0.587739769669757,0.6147847663698122,0.40005276063390793,0.47303889438928975,0.5960674587782993,0.5744578404032484,0.46050249852716824,0.6196478030837397,0.0,0.760759954868399,0.7904866416978084,0.3645297359661244,0.223017869257309,0.2526425672207412,1.0,0.0,0.4386186827692983,1.0,0.039924920532282804,0.4634465897214007,0.7469068081662684,0.4300226590015149,0.4389377074896301,0.73800034308582,0.18861830666119922,-0.013826129412225392 +23,0.47416460426005813,0.47625775561875944,0.7091818433352897,0.6777073411105259,0.0,0.12979968220343213,0.05795752011204558,0.49262255008846284,0.9150940433684287,0.5762709091872571,0.30834631685381353,0.6794632104583898,0.4981267339284384,0.516229674102819,0.8101234012956496,0.6726346946896893,0.6874437451713379,0.1314005955918645,0.5489093461148407,0.477461441209365,0.9745563170409394,0.6238589138642666,0.006746408127842551,0.6505983554270346,0.39348271859726097,0.4107754661017027,0.36646103921609924,0.33734970664911246,0.4102171203660586,0.5089071855267274,0.39981714809653784,0.5 +24,0.9181517821538281,0.23490005342749043,0.2917935940944015,0.08483012578416332,0.7415837446835043,0.20383604715502085,0.2903404052301066,0.6602948136953651,0.4933311837062574,0.425521574829669,0.7796809535890139,0.0,0.4377694065000133,0.5878216319482059,0.738866768617195,0.47087649753668487,0.354340611688996,0.6993388911381847,0.25285348459173296,0.2782001479223396,0.6451957157402957,0.8051805894189419,0.3580353950952707,0.9188221303364488,0.20829633324058733,0.7297075823920345,0.19711704000322272,0.25424459639195907,0.39710384962019346,0.6801073697105313,-0.16109175379104312,0.09467523944389558 +25,0.5124968783135825,0.24763072903018074,0.1854693407924437,0.6616086240362039,0.2984342178073236,0.5934504364557349,0.967537085576764,0.33957649603269935,0.2738825939015097,0.44320960258456077,0.578246188366743,0.3193525448601634,0.5882465597884318,0.5219976066414322,0.4230003723313364,0.5931249634150793,0.35082021849637923,0.47755792956038723,0.5307158147278285,0.31027940276876165,0.43307178835247234,0.6968267871190456,0.23323168270946315,0.5559000060956387,0.9199421351683519,0.0,0.41847691748272864,0.7641759977174303,0.2017569782535314,0.46600029430820233,-0.14947027453895953,0.24254428655269986 +26,0.7525586532782247,0.6920740820142374,0.6036075961675319,0.9287296406910062,0.7371772625251365,0.9860661465422791,0.43634163924947134,0.42442628262970544,0.5691495591354346,0.4118432822422594,0.6545654846401769,0.6622753323199961,0.48977773754777043,0.5775590649944715,0.45480287697791955,0.8422578221614533,0.2378990228055043,0.3754731462772313,0.0,0.23970674949832377,0.1271445774380448,0.18905225992796226,0.7835135141986588,0.46748666853469917,0.4684228568589018,0.2608748065347941,0.5923640080822372,0.0037140463526625367,0.8935591658628372,0.6167268233895133,0.25656846406228845,0.15103824145432826 +27,0.42229191199606686,0.22046077435218786,0.29217188153839346,0.5831385979669094,0.6997575461693605,0.8524633166236891,0.806816532896502,0.7401695548160251,0.22779287001443327,0.5004110506048336,0.37203231977897033,0.634794002521047,0.1963348375626367,0.697823617785416,0.224314769704471,0.3065575701799929,0.4138311968206998,0.021899623361435527,0.393227571618575,0.32819070982372506,0.12720327102267076,0.23966368291478735,0.745956352479916,0.44694695676666946,0.16083948793265562,0.4828548234645529,0.6436159727574187,0.5417178892186638,0.4401265971546765,0.5627525862374141,0.2537272536957636,0.2571683044653383 +28,0.37335617403325305,0.18026481074943784,0.5395010154973237,0.24221865354776095,0.8307349997334099,0.501323923969167,0.5210445946132223,0.004229799540153623,0.441037197827855,0.5284806147834374,0.6570677646263253,0.42164495395075186,0.11858246856940413,0.4996432069495869,0.45289169784735606,0.5089272300019297,0.9116756351458498,1.0,0.1253499873315776,0.7094010686008684,0.6619613391879752,0.5358767155899419,0.052071154998980185,0.4928220169682771,0.9257096030645343,0.04018349763654322,0.49721825222970034,0.46426041098021115,0.60305840210445,0.7360881168170457,-0.13717155294816227,-0.007571715055007451 +29,0.35161030929997983,1.0,0.4914838920691488,0.08889326951711762,0.5160973163367255,0.8214054927517754,0.49448536043800473,0.5561144828020697,0.3746171063389644,0.17728029858271094,0.2984425697967503,0.6145692447908928,0.649999681192406,0.7980353577739825,0.42553287713740057,0.03995708460366848,0.7605162705324747,0.46453390282769147,0.7737564518297309,0.8078943194784431,0.7503939654596381,0.7020796875546318,0.25193316575906133,0.5464553937458336,0.16452358606690043,0.5932016981699959,0.4398912173387613,0.0,0.6354564099460203,0.6508291226426525,-0.048415847445828525,0.1865175132722409 +30,0.40644115943342113,0.708246850261206,0.0,0.3849458874133796,0.6213595822951666,0.08064467863538938,0.40693872941558434,0.6016144926706337,0.1423320512098178,0.5190541546134514,0.9884955190745447,0.362291051262872,0.35492448665270504,0.44554246205476866,0.19119585640269227,0.17560776865461197,0.558376243239899,0.5255348955564395,0.7667764793974012,0.6478547214907806,0.4904481012230273,0.2254967869205578,0.32757499007804913,0.38498113067701295,0.02920645112862691,0.7222410887886099,0.10409246940487349,0.24272477201973647,0.5273769773155075,0.7384929246187919,0.15539791397654526,-0.10344419557762848 +31,0.7363602666431073,0.20916916541789765,0.2028427068878707,0.4330807096312061,0.25333236041611584,0.5420620159626676,0.41542631875525937,0.09597834025875962,0.55571209699112,0.4991283073423938,0.012036262169093204,0.32064930199318775,0.534364160039312,0.6146568978672918,0.6536320611001705,0.6575659856246635,0.1548482309915175,0.32311433019489155,0.9596846735369342,0.3442919388480325,0.0,0.48973819431704535,0.0,0.0,0.6722023761489271,0.5837140802609024,0.9004357658310631,0.6158189958010277,0.49289342745473325,1.0,0.15753834543748071,-0.161155462390967 +32,0.4480824487266255,0.45508877869120673,0.3397499132657962,0.283149416841451,0.7435605017659075,0.1747057888363105,0.25387140155655363,0.5191882615833577,0.3707621745488931,0.2341376432495148,0.2181093870725126,0.7269417077315599,0.7760806760042114,0.6615980602222682,0.3670715059602058,0.4615039675510239,1.0,0.5694454121138188,0.5534786204935401,0.07520792498879386,0.7512310563512157,0.9220085884771159,0.05576983460694135,0.6077120372491935,0.3126534805286397,0.27928513258958665,0.0,0.5689709341765392,0.6436018903915792,0.5252399968529982,-0.24227152138263597,0.5 +33,0.03384882807285733,0.2294979338404311,0.947648907638071,0.2255524219332185,0.502099721064539,0.2445577719202094,0.2891926172601008,0.2191007898496467,0.34472730128336576,0.36503497670387236,0.3864243907548844,0.41654507696837956,0.46594884597773667,0.46462768128246323,0.42726007412192324,0.47111695627402395,0.19620195353082387,0.385807502834169,0.5006485763751674,0.007174322039380565,1.0,0.4726226452617037,0.7317229820955904,0.5398703846331434,0.595374056512329,0.3222509266033647,0.6876931991545664,0.12357728125651796,0.6642297307152725,0.21611912680589523,0.424312735679113,0.1488862134229457 +34,0.3806666348719033,0.5860521355341737,0.07252305587393215,0.5146066392089099,0.407957942056678,0.4218294978737738,0.4072144499722711,0.54694160478419,0.3325327354891533,0.43180078725416465,0.5743639113153972,0.1647388073833837,0.5492218650592217,0.49128694980982035,0.7066760426089607,0.7404438406198023,0.4333697552879297,0.4981968886927916,0.4685690739455869,0.8408933250803138,0.3636273197090415,0.5541819411998017,0.8029387849442173,0.7523714481742858,0.10188817582155713,0.5786263163373463,0.3838538121203498,0.7971283068897626,0.33350476482365654,0.44390168768836485,-0.38569021871605824,0.08175191336983373 +35,0.76814425978686,0.5029386410249297,0.43032702103395354,0.3029829633430045,1.0,0.9973245545644158,0.20072684624158504,0.5779596406241168,0.2775446316988419,0.43811104017506136,0.3624191360601696,0.30284538249434273,0.39281850744036423,0.2865337694004819,0.8925123289689704,0.0,0.2848677243040109,0.5107412684929205,0.37383046745060433,0.4368375093391221,0.8514595171039133,0.7700854199345375,0.0,0.5905541787820032,0.21826312356846922,0.34323738283458993,0.994967678238948,1.0,0.07078618952858873,0.3590865000691438,0.1952442238357629,0.1132881649402739 +36,0.0,0.43857849827197415,0.9454208388000362,0.35948650788807623,0.40812808191742916,0.9714352245638949,0.7096746828167555,0.0,0.7516503287291412,0.6661127926072481,0.2643618971511382,0.6556064918352645,0.395200540044236,0.35592315357652626,0.5382586980350895,0.4686003869342093,0.43442512581287596,1.0,0.5713609493521026,1.0,0.8969790246326204,0.3878053777661755,0.29156311870731233,0.21690025624376202,0.6487243678511606,0.0,0.20750481937409687,0.807580162066199,0.019783014723730286,0.2433422613830729,-0.5,-0.008827594799510844 +37,0.06915609412424795,0.37979356326717967,0.6112523600702258,0.5502098727126205,0.44938332077769333,0.3215402342421958,0.2663939608304208,0.4347160440209046,0.5621572117480776,0.7108968297149574,0.8002166264510318,0.6189744288459651,0.627643580842488,0.4092914131414772,0.5970745742747433,0.5166248545727649,0.2885584363530186,0.35377994552456216,0.23869811247854605,0.7970302359595054,0.46080606562105747,0.27034752102747234,0.15709097305232494,0.21878111147939255,0.543277931491675,0.1538438854752533,0.40492598561157883,0.4083099364556863,0.45488322258379343,0.15282030991339307,0.08285013333383717,-0.2397648616426061 +38,0.6173836776514297,0.3228129924318802,0.867407743062242,0.003880703774527372,0.3628181607589718,0.5740488526890826,0.6579635386762479,0.3657138709598895,0.7292576643208641,0.5298040885536022,0.20047805903250437,0.574397192182156,0.741469298051631,0.17479058210245452,0.3379452876311658,0.6539457384029509,0.8707924235793065,0.49090258521237334,0.6255569602928316,0.15204486970838477,0.5991780165497186,0.0,0.38123691191361064,0.5304965023659884,0.9121071594602492,0.5280866451148002,0.704529323366023,0.31097190279850867,0.17571621159073259,0.5692380503072701,0.213007160515199,0.11120172107873882 +39,0.6540794952996453,0.5689882336381857,0.3818778658985737,0.10561414928471924,0.6673642001771205,0.437444211194987,0.08341355641791026,0.9060514845353209,0.5019215554585107,0.5573780408325435,0.2568980441241816,0.8271322578209062,0.7402823269516202,0.7229504374365157,0.08102589589400983,0.37282927828093604,0.7128267333522847,0.5573709279779492,0.725098782892547,0.7656344996479989,0.6618096618022248,0.8588794029338616,0.270200136812361,0.5437179662084877,0.5381676810929852,0.8499530152843462,0.3696463472083177,0.47229130557759474,0.19761216445214524,0.3496996602413669,0.30406824828090334,-0.07369748304824995 +40,0.4019360520588497,0.8062170604877141,0.2598043729419975,0.5585946709630225,0.38123648763593687,0.31001559505206167,0.7633910398261003,0.3826030624819621,0.8119302483247113,0.5070982848032664,0.48890241124832207,0.4420285378324783,0.9653360545974381,0.703793567813276,0.199389226763953,0.4236644810408122,0.674835520239039,0.8006572209556388,0.20772865029450344,0.6068796996620763,1.0,0.25327481315809325,0.7030637343678442,0.36645222184345605,0.19499683357539327,0.5481521675953617,0.7816618265003932,0.7066258885968891,0.2761816312890686,0.22909531354776658,0.1321449622835776,-0.17252628809680146 +41,0.4657080374547574,1.0,0.8985826218936122,0.6329861840217392,0.19366835292679935,0.7570708875070479,0.5976678580795607,0.00432760450373626,0.39017913208943295,0.7055404549687383,0.4439792338557208,1.0,0.41759381016124497,0.7944808862524138,0.5962712942001724,0.769848000617392,0.09922003188471284,0.0,0.45082146912043025,0.21131510618640975,0.22418751071885562,0.8470234013835292,0.6357368658824247,0.2796640076727468,0.444989146061159,0.43312180980254633,1.0,0.680012632273317,0.08435054581800916,0.8484390836427358,0.3175766854115931,-0.056309084260066644 +42,0.25437203226717975,0.7405529640947847,0.7665506555430186,1.0,0.369219411369062,0.3625409642377021,1.0,0.4984421474288169,0.1300279471677378,0.3233145232611891,0.4763621114370199,0.7975211627516955,0.17081212214367747,0.3700829203780607,0.1949827836876123,0.3686888506344412,0.4124847113993505,0.7897015278343607,0.30254260508797715,0.4682852231332404,0.6633655358612972,0.852624857964591,0.2969781405496118,0.6604188100789665,0.44460750219601686,0.3286439550536137,0.0,0.565821865026941,0.48073658575589767,0.6432113671098333,0.38721835464500287,0.08281966595869816 +43,0.542421787100194,0.35040958207841655,0.6387452803268256,0.2011121281913068,0.7666313152506492,0.533482773839198,0.5340316733162248,0.6472768849927115,0.8568626604561416,0.4727246780910734,0.4951230475384573,0.2157005596627029,0.6785894106672237,0.49617251219806513,0.5072791929029781,0.0,0.27336053569128554,0.21666045395097944,0.6596015924435759,0.8488419151302343,1.0,0.3713262437258773,0.5243677168394014,0.4191788872106431,0.7763442752197072,0.3987944263032859,0.5925521371385463,0.8088131128089121,0.5593264668748867,0.5859138832065252,-0.07597882297176754,0.2546651449715935 +44,0.11934404950938254,0.15232386785979063,0.9158115748821885,0.8357498079827133,0.01667588200494885,0.7297160589326461,0.5526700385639781,0.6652334219986058,0.6483444084841872,0.0,0.8110031648078144,0.8120910793190967,0.31687796951589664,0.3591342619936134,0.14207495923604924,0.39688722585756636,0.6621909593597104,0.4475599125066919,0.6463919773455908,0.3562143317243923,0.3795205366792173,0.7564311076803361,0.6694617158935562,0.3897439718709119,0.29436028625798444,0.7125829165171607,0.3211183592286111,0.5175227972868746,0.8150661237931196,0.18250992746175781,0.3410320375694926,-0.5 +45,0.8357853976789928,0.29106050939931016,0.4163401198274355,0.18168332132191056,0.0,0.26885898100822136,0.6136585665418666,0.8482155739472277,0.4636236113719412,1.0,0.6537554255402769,0.1763623003062733,0.1291693853481522,0.3121660999351664,0.5112694338184851,0.5187208549707297,0.4374382942781911,0.3599258173178828,0.24625641072511473,0.2317985541124979,0.42829317543582746,0.7183234739453641,0.8033876894392531,0.25327543218810505,0.0,0.5576838086322431,0.35340631959545077,0.7357042731477481,0.5144644547561482,0.3771975011047435,-0.3590721487942868,-0.17596933081827934 +46,0.7717166411225942,0.6033965251992219,0.20414641140605744,0.5162689693855738,0.49737385589409533,0.4601731381587007,0.3850123553432262,0.320897265195519,0.7091920843995366,0.2900197046315636,0.4740066936072325,0.959137142752744,0.6940599382785669,0.18794151090841743,0.5122440365651137,0.0,0.6105367526783212,0.568038318651322,0.3302301376648266,0.34779285595308784,0.46435278148368675,0.5777877516317235,0.3535154924756093,1.0,0.5385372512050169,0.2902632567337933,0.8746776265723908,0.46329671923224247,0.5277329201067152,0.5257960909847993,-0.013024945805936921,0.15628892687252188 +47,0.7522479695079713,0.7596980016772817,0.12446574468913257,0.6674236806814049,0.1609973058970698,0.5830845507600383,0.5498439055922929,0.2687587859894979,0.8894636749498623,0.5431333849204025,0.41486620114849027,0.34537570197745404,0.16794130974573407,0.7896308165383514,0.522570410222444,0.9635635510622883,0.25804953414631165,0.4185361787768531,0.35306259115809646,0.6266976246359509,0.7175480133902996,0.0866785547769896,0.9846097934324168,0.7047098482988364,0.4725280963871301,0.18732920560012323,0.5546843060555487,0.4436937794237245,0.2761446811864673,0.6595859643731052,0.1394354609109112,-0.15353992008518072 +48,0.5289521756357923,0.521849613323665,0.7977140244464467,0.37096149714359816,0.6398017045973328,0.25082731944640635,0.6120789728025101,0.4994059121781213,0.5193023878686647,0.5708096365880255,0.2775213234028929,0.8262384138783723,0.45585048452926946,0.0,0.2831460972432631,0.6178639057511037,0.5073331860067065,0.5779514940276583,0.8919657366218481,0.19952442105981916,0.5890260505521827,0.3876892587805304,0.33216076289466345,0.25877033470497823,0.5278233832709357,0.6616295454356952,0.3433660388178311,1.0,0.5354378446733917,0.32611454115145466,-0.2408101444996015,-0.030232508395270623 +49,0.7722470528188823,0.46488160689150426,0.5593094522790282,0.17593356338275573,0.4175304258908066,0.48662895758171726,0.14060612858263116,0.7453979377553097,0.7434965016585147,0.3075167692419261,0.08172544973120921,0.73588939113584,0.522676888828758,0.489165053543218,0.9342520638133335,0.7172953057132184,0.0,0.19232512306478955,0.18896344543400095,0.397171785766273,0.3495558735243057,0.3133365882906135,0.21828403727120632,0.7944026101311603,0.6447874364542461,0.465867964231539,0.4406696980383443,0.36632748396491027,0.3094358001367741,0.35882168327814545,0.08379718554637505,-0.23479259873491337 +50,0.328669562139562,1.0,0.6062183624308901,0.4563357998404968,0.6673935689755826,0.03238101377071617,0.51633628045323,0.6748491460388429,1.0,0.6586830913023242,0.8044634492267271,0.2779429370411622,0.4625555020065523,0.5729415900398946,0.7066029898447274,0.32449668828835,0.503617042777779,0.7562013034524264,0.3665968576473986,0.6441549274178784,0.00820707731405923,0.2008480601324239,0.6637843410663107,0.7311985378386203,0.6278780411086975,0.6374567690530563,0.16242059389261226,0.9298869385205311,0.7294538930922985,0.2671025268900461,-0.5,0.07556750105970983 +51,0.22004250710036782,0.5500491590617412,0.4859869797087288,0.6166239712872548,0.12547649863242233,0.2796066893326739,0.5965073516628407,0.733584603049483,0.5011223948814792,0.4877501467013627,0.3426719348746994,0.47100092072690863,0.31214956112986036,0.5804485423002905,0.0,0.9712679461917902,0.26881937956071056,0.40686333651014467,0.33701349533288216,0.8167563067653226,0.8672864560104833,0.48465769316157803,0.43262247658209163,0.33788398976933287,0.6902566967842805,0.3775830738292172,0.4510983191080792,0.585187418832383,0.16407759239385156,0.6997701355807762,-0.4542614516740559,-0.000716454916878933 +52,0.670226273129089,0.45540638133451306,0.8957258202044824,0.5303294614021319,0.0,0.46410595213984235,0.6835746326617405,0.16537360417808572,0.7495295077015267,0.43855675160269875,0.5271812758471781,0.7949484224947838,0.8176558700628014,0.39629485867417147,0.6254034229232187,0.3737669375086843,0.4366349563558758,0.7024169731149235,0.4369372405789565,0.36492835463296563,0.7334639804181114,0.2344084023655187,0.7889850274651449,0.6295601438404367,0.013257457740745549,0.7140131394136943,0.237151360083884,0.49776871921604704,0.5768637268139245,0.9913860775660253,0.2540109240284344,0.5 +53,0.41252229714481986,0.3620121756091235,0.19980475471231773,0.4793246210214164,0.5607769231250923,0.47835380033264274,0.26282545454788864,0.9008647699626886,0.3405183814213174,0.4067060111839509,0.484889855164844,0.7890281947394217,0.845036806022464,0.6758408403820984,0.8919186228840696,0.4189151338140633,0.5615902441229695,0.8272768436591407,0.6463881099163977,0.6876974562663555,0.5196492735754803,0.5252031501631367,0.6247597553370913,0.582237200021855,0.6388219739532379,0.536260529720487,0.6741406745250701,0.26664250321145316,0.21274287954886584,0.5335002498317641,-0.2390563942871251,0.2804945492434957 +54,0.7539750280726325,0.18602379310972617,0.25501842106111733,0.5811497579192071,0.42286742709149094,1.0,0.7396856038059865,0.6856539637810996,1.0,0.6124538224155212,0.33928176952132827,0.0,0.37253958869678916,0.2779671332936341,0.7344398812057851,0.6453351772049489,0.31399984635929223,0.5107253502762068,0.0030555817838154864,0.45416270061554675,0.2392882588663639,0.4697396770584457,0.29313895959277014,0.8070516478588867,0.2661174029584379,0.6854317117163197,0.7522276561115689,0.031450170670522304,0.6792380882193947,0.09643020986794937,-0.2558293611775097,0.12141292900790981 +55,0.7079432773160843,0.8084460746207504,0.6371797173342149,0.5974191731108015,1.0,0.882674897351308,0.2854781322188486,0.4681278122898024,0.2797969992420546,0.4226272322954437,0.29587213533865153,0.11849394799111373,0.553367041850036,0.7531302105388005,0.6525907481379929,0.15098477234732466,0.8897592127249815,0.13648726072425071,0.7222712897325151,0.10197903746318787,0.41184631554298595,0.7650134352912523,0.4884511436554735,0.564620132665154,0.37708550457831225,0.23858230429698396,0.6234292480961756,0.5950837130366096,0.5220507368681818,0.7803365519740771,0.18988682121589953,-0.07430213903353028 +56,0.6650004452913516,0.37741030499187367,0.3201327134979229,0.42063696240251425,0.24866388985156734,0.779042084091146,0.7836571576109721,0.20759597348871994,0.37382872748919416,0.7709787162195245,0.4236315721032213,0.4547001909061683,0.2911529975934545,0.38545538910038624,0.10740329412704003,0.3979449948582884,0.25833554273647963,0.6515604905254766,0.48268866268621297,0.5147473463993274,0.4242711461614114,0.4913674744636295,0.6937326392394967,0.5278975252117245,0.37172542530386954,0.6617861252854472,0.6760668689434248,0.3711598934639797,0.6282917229966536,0.42620915670266357,0.2981616006401188,0.15249488636713532 +57,0.3430172993195735,0.5639708472006024,0.7507875556766443,0.8966538342470973,0.4909046149804963,0.40717573343896635,0.4200861653509199,0.36087449694670626,0.5002253900196805,0.7806975702298828,0.4169668333240941,0.02824076603907305,0.6681076154961838,1.0,0.5922485097897155,0.7787955048177673,0.6977313435859497,0.387706619107283,0.5696920236576792,0.057745833886607645,0.604610295962417,0.8760322336810129,0.3436758411038261,0.2311403185308546,0.5263676391801627,1.0,0.5440455188970942,0.5824816000806494,0.39323788930339254,0.5278169141480658,-0.38396771371615074,0.07402606472817112 +58,0.4512105389768403,0.0,0.714604247487129,0.3328378920232181,0.7732599795289394,0.0,0.3428180138455148,0.26666809864335195,0.7516672203940578,0.47122688608315266,0.810526782351087,0.604872703913704,0.07717889653873261,0.47144898458307216,0.5146570778734423,0.37552963202026624,0.7415098859325501,0.6073151194551418,0.5400374617861537,0.248809621572036,0.4828624896007176,0.5654128256131159,0.5539959463320202,0.3302465922901728,0.6183170084237942,0.541116835570702,0.8508747144842218,0.22584899590678753,0.1644507744409039,0.766808032975101,0.012676637662959423,-0.44591380656079566 +59,0.5274641113799028,0.7033067865633456,0.3532715680954934,0.4588827262437748,0.19931727329744525,0.3374577616550324,0.6007041112245843,0.3940671768375066,0.5946955810471668,0.46650248395825933,0.6188306915737619,0.3499773805644141,0.31932695655102883,1.0,0.48631897084804454,0.373850088461946,0.0,0.07097927988666575,0.4106653767825066,0.24651316927557887,0.30626858530190904,0.20715292226225346,0.0,0.37901120296889135,0.2531883544689001,0.43008428647400854,0.5366181327158215,0.5513648495886668,0.7794725799691187,0.552236874404002,-0.3246830441852656,0.12997708264687033 +60,0.5855214741085165,0.3017969973922196,0.8908312664822289,0.8577769998484699,0.5596081302064271,0.5478074063288815,0.4355011073019859,0.5898570513054174,0.7037714673211241,0.6918430040880598,0.8160597341033526,0.5854058474521595,0.86143098593621,0.6992695731029794,0.3081548802409446,0.12742556718490644,0.4338665568103827,0.0,0.7543977733887135,0.8332328918491663,0.5856828009639305,0.770196092711597,0.38005842208754215,0.35100702519858556,0.28016559383731154,0.761614485216112,0.4043972178949004,0.460824749070444,0.8277134130471668,0.5672079609913292,-0.17957305641426552,-0.3428926246080902 +61,0.4189699438149529,0.7511663020944568,0.41784182386973834,0.18599199988790371,1.0,0.1455484732140574,0.38220867733512964,0.8522099482660412,0.3962133734478122,1.0,0.8242784149394689,0.7874694951679417,0.3791941680241144,1.0,0.4396096937280193,0.0826073992519269,0.0,0.44934128688093744,0.48623469111754924,0.23216478009724806,0.7630732087525152,0.2360418369125558,0.2626719966803097,0.3011545459066042,0.4870178201144248,0.5077688799040879,0.8625140429491498,0.6684057782471535,0.665021889564108,0.8195968201272841,-0.36480708363342507,-0.06977695840859495 +62,0.8533094622382583,0.5835695154246467,0.5464570337598088,0.2587008408417742,0.552299329916641,0.9491159976676455,0.4327848644713381,0.21005502103766538,0.533050256258788,0.5616826423786936,0.6705069834958512,0.3934685493564,0.7166645888899552,0.20442462446276544,0.4075521605592085,0.30555133247993815,0.6414236433804532,0.5165186890253776,0.19052518776181981,0.6371837952709346,0.5550788388628707,0.41765623004429986,0.6716377735606709,0.25510190852587433,0.4638811910359897,0.5072923001318755,0.822372380277775,0.547477169502391,0.7459102520802917,0.7361582270034372,-0.0872799765450663,0.18811837629168196 +63,0.6370212824553722,0.4305910081557619,0.492051606792852,0.7986116554843987,0.395862828757934,0.28168750293494804,0.4254353636071851,0.16528031349745193,0.7435070639494037,0.5126466517017918,0.08256240219506794,0.5622370027787099,0.5652159152361261,0.6681683056850876,0.9324753005723931,0.753996739943225,0.4213632699305001,0.6700464111454734,0.6394519470773101,0.24933032177486214,0.6185965470720054,0.3935326643716802,0.9453663566386186,0.3850032783009759,0.2344780443014935,0.5290184785810802,0.4429126353038732,0.23767267086208516,0.4856630647407753,0.6911984650701086,0.04455016360011427,-0.04089112867644124 +64,0.3724625103213113,0.2758354495033989,0.3224604225680284,0.9810344885012052,0.66234525945504,0.47944015972287374,0.18975191697281213,0.3677501803850127,0.7389708005287663,0.9203917015600118,0.8450692490939112,0.5070429654279668,0.16851676881332867,0.8903927930636724,0.6838843190281472,0.8560869186757327,0.9878240351795553,0.3399408726172742,0.440007267206725,0.5760683562262088,0.0,0.9990912478818506,0.3707919390327594,0.6079150466391734,0.6056788881521775,0.510659781634102,0.3266945137507782,0.48448913782672154,0.46706206959253205,0.3886287605206773,0.43952870755079826,0.053970841263118914 +65,0.7849746613664341,0.7566205521725637,0.4420042303549989,0.5982266068711057,0.5148707593979682,0.3562345758310519,0.20112907589259404,0.2695770432541844,0.6372346811120753,0.5296753960585368,0.8759802923462614,0.7397195085552504,0.35264313823173965,0.284499174187113,0.48616816460170126,0.3035620652842637,0.4101821717588261,0.3933739886251211,0.8521392651219033,0.11830733194749177,0.04439640109485832,0.7359914881379441,0.3181296652000033,0.2936441617962065,0.48243654911964606,0.5318233045604398,0.4477272615549239,0.20240867382340022,0.548047545380242,0.37442020059998343,0.5,0.25886181123112534 +66,0.562851796259624,0.7242833195734355,0.32785183267471396,0.9668354914331967,0.5603524467294086,0.5151782236634522,0.472024694597169,0.7876471914799146,0.5555234608533115,0.5884054328811648,0.4529743234468179,0.7480133575280415,0.07394285547714391,0.09235469366531918,0.28073391732368286,0.3858214166057733,0.6548617677914047,0.6066944578874788,0.24436410955482263,0.5941867907870462,0.7176872734935571,0.0,0.5662877507910274,0.859031139064131,0.4661584021618939,0.08185593761596605,0.44883422060967904,0.2111401278616265,0.821555461081229,0.9450556137814783,0.1953321014473891,-0.0799612727482627 +67,0.5089500383081038,0.27188135344979003,0.4076617486405587,0.5176799447493615,0.641062236032227,0.25094103171169807,0.0,0.6426664813820299,0.3387748231428215,0.4771775044605186,0.3889692361576658,0.4201419162881013,0.5518066332315879,0.6055865062511536,0.41926397096106877,0.5586047610132908,0.38010563208373743,0.6142139432955496,0.3428600887932695,0.4399018887860394,0.30557401701190007,0.7107822861395177,0.7267126472375234,0.6057612134773656,0.7342039756365202,0.1704225141908045,0.0,0.917835454331554,0.46882307313081195,0.5972819825877729,0.5,-0.024870149615351705 +68,0.5866472579301244,0.5527847607821561,0.4443104462777382,0.6212350881628101,0.866964621239843,0.026317377861153346,0.8966787366923481,0.6001621538355085,0.47063758116299415,0.5597712793062329,0.11536395663531385,0.740906689681717,0.05572294777118575,0.4236976339679488,0.9417549248055417,0.2653456999776381,0.9749644763383079,0.2730791445098645,0.6787001799680408,0.680807001963658,0.3411246549286856,0.5962299266409257,0.132043013980353,0.6569901983131338,0.4372663739488553,0.5488573514712662,0.5202183354335411,0.14517031953409132,0.821491327583938,0.5094906607032107,-0.14920150760598885,-0.043253813733408045 +69,0.0,0.9333399025660762,0.2951560811912415,0.2373117169147131,0.15965795699061314,0.7068822624113702,0.19120044095763805,0.5293724465105757,0.3024939740621373,0.7028261630959195,0.3322329272366892,0.49891145829222433,0.482304300462162,0.7476756809023621,0.32542692817741903,0.452053208264027,0.3000645916728918,0.4360007397441559,0.6883372608599123,0.4029203515212782,0.29039191689446675,0.45357297411529923,0.5929369399873108,0.24768412499759523,0.0007343254937092758,0.6365797688049565,0.20401382474739188,0.2505155735075629,0.13411515618172976,0.5100674840897554,-0.011077122807735153,-0.5 +70,0.7764096477841558,0.5566232867567226,0.19662500017961249,1.0,0.5920474884143224,0.4543769135447838,0.7607962810329696,0.5669525068859843,0.3277497319288019,1.0,0.6719270576232518,0.4056854288545524,0.2553324789638559,1.0,0.508063005619778,0.6083237698975847,0.3595183800597256,0.21697421007073714,0.6598200576097463,0.9513989628847508,0.15246236831189636,0.6927088249262099,0.10372935779302755,0.26302213291370274,0.22715674015971743,0.4141653333701708,0.2311075354713366,0.6593239419912271,0.3066066820501366,0.5795744417802429,0.2768895273623977,0.36228850291816933 +71,0.4262636092559942,0.30876404893575465,0.5914675559641872,0.6362626899673093,0.17328308090114092,0.04631796081974526,0.7676363273528153,0.5508088361668038,0.8978865656057373,0.0,0.4577326767970369,0.5571456071781352,0.1932649514827824,0.6005220072270617,0.32219511494981967,0.7100555159714457,0.6605015953906863,0.6938906149410076,0.814477106179698,0.4981482820424959,0.6419092079421541,0.17614831274810405,0.4329850978512799,0.7234605080358275,0.6082670128731003,0.9201255595639155,0.6354199234170677,0.06310344235189153,0.8624518931495262,0.8048837954095084,-0.47256777474210254,0.015061400133228622 +72,0.8926611151702155,0.0,0.19127194841352663,0.625275377877066,0.7026212833010629,0.0,0.6830744909194946,0.2513675886888233,0.5628383887761244,0.6586564381257757,0.3465840042217367,0.9229687161390363,0.9101179280794125,0.4892576856997975,0.21132465776929782,0.32558671059665323,0.6719246333788154,0.5362135910062958,0.5664463479776324,0.7122559700373687,0.5025562390508914,0.8766104684190209,0.40622524998190696,0.6415922868586182,0.8000596793931578,0.5477732622491239,0.0,0.0,0.35470702804438003,0.35998144870000204,0.053071324305219686,0.3041164472570473 +73,0.45419010265023924,0.5251254200011489,0.291311939519302,0.8407634543515594,0.2088617781128218,0.4393268039282698,0.603539047574565,0.3297081184159053,0.5976179138283327,1.0,0.6312025808959824,0.8058511622276873,0.9258894937391706,0.6615272955519735,0.5546000547505676,0.3605035481417298,0.6698420946961661,0.9573390424736032,0.8544096303463141,0.4898357244339549,0.5193061019408821,0.6149612118820242,0.8465416195097013,0.07994274968889581,0.07013337002911685,0.46029549785742735,0.2262176232697914,0.37243335213667095,0.790714698105807,0.7964109209049838,0.048481369550199126,-0.41902487244393716 +74,0.5884389835480314,0.22426714118106617,0.4229603025318303,0.47686679832637213,0.3460398015854096,0.6871522657985631,0.6754665863317875,0.6128862048535578,0.23787718027298216,0.4717042807245642,0.3484649486655873,0.7889096341364592,0.17945702674803315,0.5895058146695312,0.6202485026633697,0.3908123281016713,0.0,0.0,0.5461879244656569,0.9595020217876276,0.62715748918992,0.47245660906351733,0.15659874385118,0.2907789358538445,0.5329074790882007,0.5638548874091535,0.4786888255522671,0.485075572382609,0.09657652137739847,0.4482585143318049,-0.06725354101842611,0.30410993758010163 +75,0.13449409322531358,0.4116933767075943,0.43911481491537174,0.6429746537381239,0.038842406927606044,0.20992637106276946,0.5712004217698833,0.48332982561592197,0.0,0.17638330934407864,0.5169775040324076,0.3366875208868428,0.14325980906712055,0.45066348550648466,0.4422842897593333,0.6278530018785955,0.5641673207961422,0.6025168579262402,0.3796196834935649,0.8915616095411193,0.5080895065253159,0.03387340228680685,0.256147371254307,0.1613928177284124,0.7822033701057469,0.7755054173150356,0.7835505320038326,0.32717110905168223,0.4209024651979809,0.37976182682447734,0.5,0.10007060428771795 +76,0.6161388208053726,0.6171895782979467,0.4810000628081463,0.48914394451851073,0.5232156900907154,0.840919571882605,0.4523211888165327,0.0,0.6622975772820183,0.24858294026174982,0.8073624839753115,0.4269158027348272,0.24179169844607273,0.0613123338787564,0.5960573527489551,0.507423580068971,0.5826226513295278,0.751354176661776,0.5043026760384841,0.7199183966806821,0.5769420029981731,0.48997507798917134,0.7478011363632642,0.5524425313315022,0.34892861560474925,0.3050188109325338,0.5518086012662635,0.38862790095711075,0.9091592039072317,0.8396729350374081,-0.04764717327345088,-0.11695373457547374 +77,0.2410821587552287,0.11065079848442005,0.320307018340087,0.4039736517955476,0.33821180226913317,0.7531650494008917,0.789614250236921,0.17960614131547942,0.3045955661369366,0.3205170812680558,0.6972361279630415,0.2384546913144605,0.3720922525748124,0.7292514324550028,0.6286088640181459,0.17776726774384888,0.0,0.0,0.3509536117848202,0.7758935732271766,0.3632552282174019,0.16076633892010977,1.0,0.7318849992697579,0.4916567634503535,0.3892562116784354,1.0,0.8158422558332844,0.5440418683577822,0.49915152977779653,-0.4300320456287038,-0.014146719571444035 +78,0.5800862584203833,0.7430704007494301,0.658868611059073,0.46514148803960104,0.3103072862286981,0.16754978283634353,0.06919956037864788,0.6012702791010407,0.42450297817533644,0.40285311089318643,0.29516647174736854,0.4528252708581529,0.10604427187218513,0.3733271986390453,0.4399782619228241,0.39117032084311365,0.9280093458449112,0.9560064036369134,0.6417087919554818,0.0,0.4167418380913969,0.5333126060552156,0.7626792926812778,0.05775601422019688,0.5754543713214022,0.5362049925123034,0.6508957480585762,0.6443628700051598,0.13859524842855958,0.6288630055036608,0.2671965793866974,-0.12619157204807 +79,0.3830839342514333,0.5605563425415135,0.19792343147007407,0.19198356327081245,0.5094176100488251,0.11787439164628027,0.2982397503524278,0.11581092935023224,0.3408153089738385,0.4737439768889116,0.6070106421098882,0.47572344804267164,0.33966635150933344,0.8037580342678741,0.4271383154215672,0.026864675100884028,0.9400649318109718,0.7337991108723134,0.6847426925274228,0.7704501397108925,0.4603684761392425,0.34293093831202937,1.0,0.8120955433329913,0.5121727704310906,0.10287800536853653,0.4952084924841317,0.21243746357684068,0.1443468310614821,0.6721605123638954,-0.3578545123985546,-0.4662317451430373 +80,0.16164701411523708,0.5453273368129118,1.0,0.5431336047071377,0.7951347943480784,0.18295801459296873,0.5750550365492066,0.3875035858249998,0.0,0.5946967988021741,0.37116146922011606,0.5065763993363138,0.3888407468625838,0.739101737228485,0.6334757058367494,0.8482795362427529,0.15980887882410916,0.543231292206122,0.3463687968156014,0.5887372877079532,0.23990058258269414,0.7666248473872479,0.48251927360020497,0.3862683460337661,0.8766835624221029,0.8227102783924396,0.9563862125208189,0.5843403258558691,0.5348219185417755,0.7288662961234653,0.2191975730829104,-0.12094588517339598 +81,0.5870130253739732,0.644909919817602,0.6203007296740044,0.27260440382023243,0.5854874528658305,0.41873347381716003,0.43924780605570024,0.5148862303962867,0.7728214129222679,0.470026354662507,0.6109444979549834,0.7452555848225523,0.33885246396292834,0.8279634285417539,0.08696785114808914,0.2476392239913866,0.7149529523035912,0.5519134975027595,0.2216806855462481,0.08291562123683582,0.8030388988413022,0.16885510574343016,0.5629927124747528,0.7494186205362918,0.6046975657191946,0.6625625896083926,1.0,0.8329364421154052,0.36735386792237207,0.5044739501801615,-0.4790246291764833,0.04174827024970315 +82,0.0,0.4374677773286896,0.9925062581116088,0.7451134532664903,0.6113543910339178,0.5038623086732812,0.5547925569539454,0.15394695994558572,0.4538028623846119,0.45050274986113603,0.8525265361083235,0.5277985049985882,0.6518383392915101,0.39761069045661673,0.635735539993383,0.8779498823839555,0.7975265859438292,0.42141690981370744,0.5657339999641589,0.19621999513660332,0.7476655432834375,0.2641380930250904,0.2810122295612149,0.11008452445578432,0.9299777760774904,0.4485522543784665,0.4737521632685032,0.6622383255372728,0.06949098669621889,0.5534091735865352,-0.06257986086795506,-0.12524506529765989 +83,0.20318129831440412,0.772428565110985,0.7183142463865054,0.46781605013237004,0.3220610230727281,0.8088404113488747,0.9447163166149048,0.26727792959097935,0.12770034298367922,0.13451385239805536,0.6055689663502621,0.9509615398225972,1.0,0.08884314687434625,0.7199199132521774,0.7463667210394476,0.1689510507424608,0.6159163927785317,0.2827296002771419,0.0,0.0,0.6718843869524431,0.7075389145062488,0.5412954649503233,0.42352195769271334,0.4586941150906984,0.41202707838929997,0.5182908462514131,0.5048408819520657,0.4621549061128577,-0.06721369395692647,0.2767264467883054 +84,0.5913444782866055,0.4548931395990663,0.22945158961069406,0.5217950743534246,0.0,0.31964333721829075,0.434126973039747,0.4477855964394353,0.5960051081806675,0.5468903405464732,0.5988473452413375,0.06510326052840698,1.0,0.6226285461393269,0.45339489193453625,0.6764713937463696,0.5016147854174835,1.0,0.5412111443133572,0.4535461279395134,0.3375507797390682,0.7901683314930275,0.6831727009867217,0.6928176199128854,0.4508700782092094,0.45911617993720044,0.6404781380597082,0.36069588211308545,0.8652281212558368,0.27198923094304417,0.11888108683822671,-0.031411528691290935 +85,0.5111243484603043,0.6841635107687981,0.2469571645886416,0.48919009234627053,0.7180785399370924,0.5280647853136773,0.4506286644038671,0.42181675560153997,0.4340173230381058,0.4433048087219369,0.33721904313310647,0.5105647614559925,0.9697558336304265,0.45502287528728896,0.5537562414673314,0.48155741297121823,0.7612836303934725,0.3280201538748829,0.20429920893122072,0.4911466034547927,0.46340018632551355,0.47969128161574565,0.5003506366804936,0.5369669208778461,0.4329138055966947,0.6008744461401712,0.25174590817219,0.5358676931331116,0.8851919539015736,0.5945242602043511,0.04265631479491806,-0.0994902421412142 +86,0.4949018830606719,0.3157001922701804,0.8008985496092087,0.3811789003131693,0.005412869386976216,0.5228660381276687,0.2396726209629299,0.0,0.20916353042412889,0.23772394701279115,0.4677181750907341,0.7997690491602468,0.4893340728725225,0.23553508699471964,0.4587635681242299,0.5740387059537074,0.2557903160180508,0.7708226299127922,0.40982206776064684,0.5066311971225967,0.4041747608033011,0.7427429426943354,0.3765159204751619,0.6500244685961503,0.48529541361011147,0.5722090130392191,0.45233610540109165,0.782014990927149,0.3432753604776251,0.3372237673144589,-0.4613506455230872,-0.3955887680433119 +87,0.3855300034949842,0.9753495711527034,0.21628734142570671,0.49908771552929615,0.7906179032604913,0.4084392369774662,0.6700732086882257,0.7432805644546522,0.6194897137214798,0.6651275603452652,0.7287665817526896,0.8149764361381644,0.35858198935162,0.9313370764848967,0.6274609485502317,0.9873013702779384,0.840764719078658,0.526052125484992,0.6813637512369751,0.5456764629882626,0.5699097944568309,0.9410403685321012,0.6149147960257284,0.37683968757899544,0.4231860124518507,0.28393320403882605,0.4883836189125367,0.31346309934900096,0.41039013089579557,0.7047508204540602,0.18416495530187244,0.4721842121218343 +88,0.938825620055938,0.44872411165243686,0.3038121737535815,0.42338838999104994,0.366196436218441,0.19176508188032038,0.8924952342738228,0.7048508719506402,1.0,0.22259226907662943,0.9602904703414944,0.6454776867986131,0.6238765292517277,0.2537558981245511,0.5760947975671155,0.5959483548966521,0.2807138576640953,0.0,0.8078821692427112,1.0,0.6534589648148527,0.9884045786118754,0.4346170322139931,0.5472562600821447,0.3876668312497893,0.254455963072452,0.009487569442127586,0.4072733673082581,0.6559908760886257,0.7253163547962904,-0.3314233989330351,-0.3373548971603823 +89,0.3123828253473098,0.7679189685962856,0.25732754289618004,0.3655572949312649,0.45337854545660716,0.42459669790797766,1.0,0.9136119652379532,0.7075262328319502,0.6010392163414429,0.6022062352398405,0.5582885686542022,0.897507244892782,0.5888285740878225,0.06444762122058822,0.3476170587363663,0.42392807629450663,0.4938977179617262,0.0,0.8030522793880959,0.6458464760950521,0.3993349147349139,0.6557228799902829,0.11141198237953537,0.21001660575158587,0.646800491704022,0.6891732526808051,0.2685895393631248,0.7917228644006054,0.3427430588987587,0.13573026139681937,0.2089609555875807 +90,0.5422045481994687,0.4844626329202054,0.50787389956154,0.3599820536566434,0.59951012180878,0.8179098548441082,0.6100626719268758,0.6167012854634457,0.4045806767999911,0.5093273278798466,0.05811427003224723,0.696703834272464,0.25019025228931263,0.6481286522990505,0.7302634397690546,0.7730965862800201,0.34347107548929934,0.05467314821659491,0.5879434734233636,0.7381537032235057,0.06417384086639116,0.6014655537246582,0.4345144165002432,0.5644220428979017,0.21807352705230332,0.645808327362537,0.3349052749212559,0.3510018518122391,0.605599429012144,0.45893332144382865,0.11308277899062359,0.36754994510260847 +91,0.6169557734799321,0.8383687738818464,0.4952567701803625,0.7661778215754695,0.19959853719445597,0.25417465610254886,0.03168414783737766,0.5854281179368142,0.2740005816746705,0.0,0.5650005163151752,0.189094806461151,0.4092488758798813,0.8383078506617856,0.22500419548793243,0.5372924948103439,0.6552394066906395,0.9336217878409825,0.21455528290116282,0.5029853248556954,1.0,0.720361187553128,0.6177934737662886,0.4706917558569956,0.154195695387226,0.4483126800310774,0.4008686038614082,0.23081337974419902,0.2767523992634775,0.345198634313781,0.33700474911073264,-0.5 +92,0.0,0.45531906617275986,0.1531715008874603,0.3949887561934082,0.49913413950358726,0.9343865131344913,0.14446255839422134,0.6132094437229434,0.5786058111805248,0.5570660271859157,0.7441914946521633,0.618313967455493,0.7777292837046399,0.07192788029740804,0.40951558377734365,0.4360433679025223,0.22065457628975954,0.45025105810998023,0.6402363133235385,0.3927839618919361,0.5565026412318798,0.04110094371942774,0.8679020778744633,0.12568603804537154,0.6641547551259991,0.8608994981487452,0.5011330475666346,0.807580808746582,0.4308918949012231,0.7078800642507732,-0.23145024524092608,-0.21581215285776195 +93,0.2694053190969835,0.7831797668899754,0.6321649587968877,0.4312953600592351,0.5396699439709958,0.11666135393669114,0.0,0.822074810706204,0.45128040319102475,0.7182684588740559,0.6534110762670867,0.34312476988732743,0.4726459430516131,0.43827280388100526,0.3216848315018206,0.16231723172569862,0.5144422874284493,0.12590135439484273,0.355712642710469,0.41914977892571303,0.9538861437951405,0.4612871494057248,0.6622702296899453,0.6230467333237544,0.38041988398939836,0.4734785985511577,0.5132130193173413,0.5037941166910144,0.5565477334892124,0.350312309414785,0.26918007081516154,-0.02770531517611326 +94,0.6666500084712537,0.28871948475298204,0.8730854392554007,0.3120020291697324,0.9120987304181639,0.5702745922123175,0.2638801313147978,0.7697852384582144,0.24778130386198638,0.8615189892997501,0.2665638649018289,0.4124839166350119,0.3953433643756823,0.36073237211889,0.5096763320443566,0.7239437332772848,0.6940519453321261,0.6738111449809068,0.3527437773245446,0.5062589203786457,0.9071350349461802,1.0,0.38692228016688,0.5524968127771072,0.6461147686526466,0.304347920904348,0.4632013663042951,0.801211450476998,0.8873902063644976,0.5519860893282365,-0.042665507136309525,-0.07543592262427655 +95,0.40399710411397294,0.6611952853526666,0.7806790060042248,0.7814082085038858,0.6326114166560098,0.5382000230683094,0.5182869420988676,0.0,0.38264568725317716,0.40485754665110574,0.1770080905335163,0.5635746098663962,0.39373381882135117,0.4572374493282391,0.15608388765005782,0.5738498474901829,0.5338159085559167,0.4024210230140839,0.5075833864166158,0.137044382335882,0.2859715125862241,0.3411588260417962,0.7922638425073126,0.48967716513367077,0.6255959750689413,0.6133461426649881,0.0,0.7627643858906852,1.0,0.29729124312542166,-0.03941334288891294,-0.4940103567481645 +96,0.30353151236660525,0.6626607627205008,0.5200151783308468,0.5534903469167705,0.7103860954364796,0.3525187281806021,0.4044118693250328,0.5664920474950936,0.7473317914864679,0.8051230258715045,0.6661023388808562,0.39361400811376784,0.42750076855771435,0.9100224900518472,0.39196064808855746,0.503265047130121,0.37878248714276574,0.3732750187690387,0.6585973739033126,0.6236000657678453,0.12227046018143417,0.9345725959430724,0.8745163214425666,0.36353215403926864,0.7488562878416849,0.0,0.8712443104163776,0.7585879800594992,0.8649582233936299,0.3428177309162853,-0.0948767657784465,-0.08884454501948058 +97,0.459099748023998,0.04710591451426299,0.7228347048132102,0.7698915890475564,0.705514744064665,0.31472458008069704,0.5858570233612174,0.7168109459876191,0.5420427254864731,0.29015423264670426,0.5209695778853178,0.2688803472925929,0.38774624426787296,0.5474705109644213,0.03993570815192632,0.25247591502339706,0.5178723755067391,0.7278003865278063,0.3803208760794633,0.6167930818016187,0.6679057860415283,0.3684930325843533,0.6125878023904876,0.47165014859418514,0.8196497469523513,0.43193316965754525,0.36231798183311703,0.4187031133238797,0.7233973536215068,0.48965516238804463,0.058048408050828455,0.03511985456329148 +98,0.5122677414604488,0.5013272449668944,0.42476913106912273,0.5585735892323557,0.052351602849811385,0.7142955795768947,0.015704160864150274,0.8184382742412444,0.4141167580109011,0.20054129704906054,0.9213716819075262,0.34639717895944006,0.6812587102524907,0.5900796876157719,0.3056600520934021,0.2909485929332891,0.8604453596547472,0.5972366333719288,0.8966457138535313,0.9121317963682467,0.3088850402050586,0.0,0.545783492335858,0.7055450491361126,1.0,0.5037953117288198,0.12188963247456303,0.8413696968343282,0.20414868658062485,0.3406633674832776,-0.007550907450139682,0.4292562176133695 +99,0.30027476937610575,0.1971881860058456,0.34170154891461857,0.6332741839553133,0.9129150782187159,0.36742454292946086,0.44455371167361984,0.7083453271204787,0.4195326726150156,0.49705466955871797,0.5932150144810101,0.4265571127682233,0.07509510539295028,0.44983955696942035,0.247283621689931,1.0,0.2378751117843687,0.34078771641400984,0.41385845333774757,0.48428560888076044,0.19473539749810959,0.44556911294368545,0.6253363329156638,0.4475069928061957,0.4331541889318046,0.5812785728485341,0.007666880198995973,0.7564308419295487,0.6658394680634336,0.33514120248833423,-0.5,-0.18846975808019237 diff --git a/autotest/utils/zdt1_template/zdt1.pst b/autotest/utils/zdt1_template/zdt1.pst new file mode 100644 index 000000000..e049699eb --- /dev/null +++ b/autotest/utils/zdt1_template/zdt1.pst @@ -0,0 +1,67 @@ +pcf +* control data + restart estimation + 32 2 2 0 1 + 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 + 100 1.000000E-02 3 3 1.000000E-02 3 + 0 0 0 +* singular value decomposition + 1 + 10000000 1.000000E-06 +1 +* parameter groups +decvars relative 1.0000000000E-02 0.0 switch 2.0000000000E+00 parabolic 1.0000000000E-05 5.0000000000E-01 smaller +obj_add absolute 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 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_1 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_10 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_11 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_12 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_13 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_14 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_15 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_16 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_17 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_18 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_19 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_2 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_20 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_21 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_22 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_23 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_24 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_25 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_26 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_27 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_28 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_29 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_3 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_4 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_5 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_6 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_7 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_8 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 1.0000000000E+00 0.0000000000E+00 1 +dv_9 none relative 5.0000000000E-01 0.0000000000E+00 1.0000000000E+00 decvars 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 +obj1_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 +* observation data +obj_1 5.0000000000E+00 1.0000000000E+00 less_than +obj_2 3.8416876048E+01 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 +++opt_dec_var_groups(decvars) +++mou_objectives(obj_1,obj_2) +++mou_generator(pso) +++opt_risk(0.5) +++mou_population_size(100) +++mou_save_population_every(1) diff --git a/autotest/utils_tests.py b/autotest/utils_tests.py index 39d894381..8e2c07378 100644 --- a/autotest/utils_tests.py +++ b/autotest/utils_tests.py @@ -1,12 +1,46 @@ import os import shutil import pytest +import platform # if not os.path.exists("temp"): # os.mkdir("temp") from pathlib import Path +import pandas as pd import pyemu from pst_from_tests import _get_port +ext = '' +local_bins = False # change if wanting to test with local binary exes +if local_bins: + bin_path = os.path.join("..", "..", "bin") + if "linux" in platform.system().lower(): + pass + bin_path = os.path.join(bin_path, "linux") + elif "darwin" in platform.system().lower(): + pass + bin_path = os.path.join(bin_path, "mac") + else: + bin_path = os.path.join(bin_path, "win") + ext = '.exe' +else: + bin_path = '' + if "windows" in platform.system().lower(): + ext = '.exe' + +mf_exe_path = os.path.join(bin_path, "mfnwt") +mt_exe_path = os.path.join(bin_path, "mt3dusgs") +usg_exe_path = os.path.join(bin_path, "mfusg_gsi") +mf6_exe_path = os.path.join(bin_path, "mf6") +pp_exe_path = os.path.join(bin_path, "pestpp-glm") +ies_exe_path = os.path.join(bin_path, "pestpp-ies") +mou_exe_path = os.path.join(bin_path, "pestpp-mou") +swp_exe_path = os.path.join(bin_path, "pestpp-swp") + +mf_exe_name = os.path.basename(mf_exe_path) +mf6_exe_name = os.path.basename(mf6_exe_path) + + + def add_pi_obj_func_test(tmp_path): import os @@ -2619,11 +2653,423 @@ def ppu_geostats_test(tmp_path): # plt.show() # exit() +def gpr_compare_invest(): + import numpy as np + from sklearn.gaussian_process import GaussianProcessRegressor + case = "zdt1" + use_chances = False + m_d = os.path.join(case+"_gpr_baseline") + org_d = os.path.join("utils",case+"_template") + t_d = case+"_template" + if os.path.exists(t_d): + shutil.rmtree(t_d) + shutil.copytree(org_d,t_d) + if os.path.exists(m_d): + shutil.rmtree(m_d) + + pst = pyemu.Pst(os.path.join(t_d, case+".pst")) + pst.pestpp_options["mou_generator"] = "pso" + if use_chances: + pst.pestpp_options["opt_risk"] = 0.95 + pst.pestpp_options["opt_stack_size"] = 50 + pst.pestpp_options["opt_recalc_chance_every"] = 10000 + pst.pestpp_options["opt_chance_points"] = "single" + else: + pst.pestpp_options["opt_risk"] = 0.5 + + pop_size = 60 + num_workers = 60 + noptmax_full = 30 + noptmax_inner = 10 + noptmax_outer = 5 + port = 4554 + pst.control_data.noptmax = noptmax_full + pst.pestpp_options["mou_population_size"] = pop_size + pst.pestpp_options["mou_save_population_every"] = 1 + pst.write(os.path.join(t_d, case+".pst")) + if not os.path.exists(m_d): + pyemu.os_utils.start_workers(t_d, mou_exe_path, case+".pst", num_workers, worker_root=".", + master_dir=m_d, verbose=True, port=port) + #shutil.copytree(t_d,m_d) + #pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=m_d) + # use the initial population files for training + dv_pops = [os.path.join(m_d,"{0}.0.dv_pop.csv".format(case))] + obs_pops = [f.replace("dv_","obs_") for f in dv_pops] + + pst_fname = os.path.join(m_d,case+".pst") + gpr_t_d = os.path.join(case+"_gpr_template") + pyemu.helpers.prep_for_gpr(pst_fname,dv_pops,obs_pops,gpr_t_d=gpr_t_d,nverf=int(pop_size*.1),\ + plot_fits=True,apply_standard_scalar=False,include_emulated_std_obs=True) + gpst = pyemu.Pst(os.path.join(gpr_t_d,case+".pst")) + shutil.copy2(os.path.join(m_d,case+".0.dv_pop.csv"),os.path.join(gpr_t_d,"initial_dv_pop.csv")) + gpst.pestpp_options["mou_dv_population_file"] = "initial_dv_pop.csv" + gpst.control_data.noptmax = noptmax_full + gpst.write(os.path.join(gpr_t_d,case+".pst"),version=2) + gpr_m_d = gpr_t_d.replace("template","master") + if os.path.exists(gpr_m_d): + shutil.rmtree(gpr_m_d) + pyemu.os_utils.start_workers(gpr_t_d, mou_exe_path, case+".pst", num_workers, worker_root=".", + master_dir=gpr_m_d, verbose=True, port=port) + + #o1 = pd.read_csv(os.path.join(m_d,case+".{0}.obs_pop.csv".format(max(0,pst.control_data.noptmax)))) + o1 = pd.read_csv(os.path.join(m_d,case+".pareto.archive.summary.csv")) + o1 = o1.loc[o1.generation == o1.generation.max(), :] + o1 = o1.loc[o1.is_feasible == True, :] + o1 = o1.loc[o1.nsga2_front == 1, :] + + + import matplotlib.pyplot as plt + o2 = pd.read_csv(os.path.join(gpr_m_d, case + ".{0}.obs_pop.csv".format(max(0, gpst.control_data.noptmax)))) + fig,ax = plt.subplots(1,1,figsize=(5,5)) + ax.scatter(o1.obj_1,o1.obj_2,c="r",s=10) + ax.scatter(o2.obj_1,o2.obj_2,c="0.5",s=10,alpha=0.5) + plt.tight_layout() + plt.savefig("gpr_{0}_compare_noiter.pdf".format(case)) + plt.close(fig) + + # now lets try an inner-outer scheme... + + gpst.control_data.noptmax = noptmax_inner + gpst.write(os.path.join(gpr_t_d,case+".pst"),version=2) + gpr_t_d_iter = gpr_t_d+"_outeriter{0}".format(0) + if os.path.exists(gpr_t_d_iter): + shutil.rmtree(gpr_t_d_iter) + shutil.copytree(gpr_t_d,gpr_t_d_iter) + for iouter in range(1,noptmax_outer+1): + #run the gpr emulator + gpr_m_d_iter = gpr_t_d_iter.replace("template","master") + complex_m_d_iter = t_d.replace("template", "master_complex_retrain_outeriter{0}".format(iouter)) + if os.path.exists(gpr_m_d_iter): + shutil.rmtree(gpr_m_d_iter) + pyemu.os_utils.start_workers(gpr_t_d_iter, mou_exe_path, case+".pst", num_workers, worker_root=".", + master_dir=gpr_m_d_iter, verbose=True, port=port) + o2 = pd.read_csv(os.path.join(gpr_m_d_iter,case+".{0}.obs_pop.csv".format(gpst.control_data.noptmax))) + + # now run the final dv pop thru the "complex" model + final_gpr_dvpop_fname = os.path.join(gpr_m_d_iter,case+".archive.dv_pop.csv") + assert os.path.exists(final_gpr_dvpop_fname) + complex_model_dvpop_fname = os.path.join(t_d,"gpr_outeriter{0}_dvpop.csv".format(iouter)) + if os.path.exists(complex_model_dvpop_fname): + os.remove(complex_model_dvpop_fname) + # load the gpr archive and do something clever to pick new points to eval + # with the complex model + dvpop = pd.read_csv(final_gpr_dvpop_fname,index_col=0) + if dvpop.shape[0] > pop_size: + arc_sum = pd.read_csv(os.path.join(gpr_m_d_iter,case+".pareto.archive.summary.csv")) + as_front_map = {member:front for member,front in zip(arc_sum.member,arc_sum.nsga2_front)} + as_crowd_map = {member: crowd for member, crowd in zip(arc_sum.member, arc_sum.nsga2_crowding_distance)} + as_feas_map = {member: feas for member, feas in zip(arc_sum.member, arc_sum.feasible_distance)} + as_gen_map = {member: gen for member, gen in zip(arc_sum.member, arc_sum.generation)} + + dvpop.loc[:,"front"] = dvpop.index.map(lambda x: as_front_map.get(x,np.nan)) + dvpop.loc[:, "crowd"] = dvpop.index.map(lambda x: as_crowd_map.get(x, np.nan)) + dvpop.loc[:,"feas"] = dvpop.index.map(lambda x: as_feas_map.get(x,np.nan)) + dvpop.loc[:, "gen"] = dvpop.index.map(lambda x: as_gen_map.get(x, np.nan)) + #drop members that have missing archive info + dvpop = dvpop.dropna() + if dvpop.shape[0] > pop_size: + dvpop.sort_values(by=["gen","feas","front","crowd"],ascending=[False,True,True,False],inplace=True) + dvpop = dvpop.iloc[:pop_size,:] + dvpop.drop(["gen","feas","front","crowd"],axis=1,inplace=True) + + #shutil.copy2(final_gpr_dvpop_fname,complex_model_dvpop_fname) + dvpop.to_csv(complex_model_dvpop_fname) + pst.pestpp_options["mou_dv_population_file"] = os.path.split(complex_model_dvpop_fname)[1] + pst.control_data.noptmax = -1 + pst.write(os.path.join(t_d,case+".pst"),version=2) + + pyemu.os_utils.start_workers(t_d, mou_exe_path, case+".pst", num_workers, worker_root=".", + master_dir=complex_m_d_iter, verbose=True, port=port) + + # plot the complex model results... + o2 = pd.read_csv(os.path.join(complex_m_d_iter, case + ".pareto.archive.summary.csv")) + o2 = o2.loc[o2.generation == o2.generation.max(), :] + #o2 = o2.loc[o2.is_feasible==True,:] + o2 = o2.loc[o2.nsga2_front == 1, :] + fig, ax = plt.subplots(1, 1, figsize=(5, 5)) + ax.scatter(o1.obj_1, o1.obj_2,c="r",s=10,label="full complex") + ax.scatter(o2.obj_1, o2.obj_2,c="0.5",s=10,alpha=0.5,label="mixed emulated-complex") + ax.legend(loc="upper right") + ax.set_xlim(0,10) + ax.set_ylim(0,20) + plt.tight_layout() + plt.savefig("gpr_{0}_compare_iterscheme_{1}.pdf".format(case,iouter)) + plt.close(fig) + + # now add those complex model input-output pop files to the list and retrain + # the gpr + dv_pops.append(os.path.join(complex_m_d_iter,case+".0.dv_pop.csv")) + obs_pops.append(os.path.join(complex_m_d_iter,case+".0.obs_pop.csv")) + gpr_t_d_iter = gpr_t_d+"_outeriter{0}".format(iouter) + pyemu.helpers.prep_for_gpr(pst_fname,dv_pops,obs_pops,gpr_t_d=gpr_t_d_iter,nverf=int(pop_size*.1), + plot_fits=True,apply_standard_scalar=False,include_emulated_std_obs=True) + gpst_iter = pyemu.Pst(os.path.join(gpr_t_d_iter,case+".pst")) + #aggdf = pd.read_csv(os.path.join(gpr_t_d,"gpr_aggregate_training_data.csv"),index_col=0) + #aggdf.index = ["outeriter{0}_member{1}".format(iouter,i) for i in range(aggdf.shape[0])] + restart_gpr_dvpop_fname = "gpr_restart_dvpop_outeriter{0}.csv".format(iouter) + #aggdf.to_csv(os.path.join(gpr_t_d_iter,restart_gpr_dvpop_fname)) + shutil.copy2(os.path.join(complex_m_d_iter,case+".0.dv_pop.csv"),os.path.join(gpr_t_d_iter,restart_gpr_dvpop_fname)) + gpst_iter.pestpp_options["mou_dv_population_file"] = restart_gpr_dvpop_fname + gpst_iter.control_data.noptmax = gpst.control_data.noptmax + gpst_iter.write(os.path.join(gpr_t_d_iter,case+".pst"),version=2) + + +def gpr_constr_test(): + import numpy as np + from sklearn.gaussian_process import GaussianProcessRegressor + case = "constr" + use_chances = False + m_d = os.path.join(case + "_gpr_baseline") + org_d = os.path.join("utils", case + "_template") + t_d = case + "_template" + if os.path.exists(t_d): + shutil.rmtree(t_d) + shutil.copytree(org_d, t_d) + if os.path.exists(m_d): + shutil.rmtree(m_d) + + pst = pyemu.Pst(os.path.join(t_d, case + ".pst")) + pst.pestpp_options["mou_generator"] = "pso" + if use_chances: + pst.pestpp_options["opt_risk"] = 0.95 + pst.pestpp_options["opt_stack_size"] = 50 + pst.pestpp_options["opt_recalc_chance_every"] = 10000 + pst.pestpp_options["opt_chance_points"] = "single" + else: + pst.pestpp_options["opt_risk"] = 0.5 + + pop_size = 20 + num_workers = 20 + noptmax_full = 3 + noptmax_inner = 3 + noptmax_outer = 3 + port = 4554 + pst.control_data.noptmax = -1 + pst.pestpp_options["mou_population_size"] = pop_size + pst.pestpp_options["mou_save_population_every"] = 1 + pst.write(os.path.join(t_d, case + ".pst")) + #if not os.path.exists(m_d): + # pyemu.os_utils.start_workers(t_d, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=m_d, verbose=True, port=port) + if os.path.exists(m_d): + shutil.rmtree(m_d) + shutil.copytree(t_d,m_d) + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=m_d) + # use the initial population files for training + dv_pops = [os.path.join(m_d, "{0}.0.dv_pop.csv".format(case))] + obs_pops = [f.replace("dv_", "obs_") for f in dv_pops] + + pst_fname = os.path.join(m_d, case + ".pst") + gpr_t_d = os.path.join(case + "_gpr_template") + pyemu.helpers.prep_for_gpr(pst_fname, dv_pops, obs_pops, gpr_t_d=gpr_t_d, nverf=int(pop_size * .1), \ + plot_fits=True, apply_standard_scalar=False, include_emulated_std_obs=True) + gpst = pyemu.Pst(os.path.join(gpr_t_d, case + ".pst")) + #shutil.copy2(os.path.join(m_d, case + ".0.dv_pop.csv"), os.path.join(gpr_t_d, "initial_dv_pop.csv")) + #gpst.pestpp_options["mou_dv_population_file"] = "initial_dv_pop.csv" + gpst.pestpp_options.pop("mou_dv_population_file",None) #= "initial_dv_pop.csv" + + gpst.control_data.noptmax = noptmax_full + gpst.write(os.path.join(gpr_t_d, case + ".pst"), version=2) + gpr_m_d = gpr_t_d.replace("template", "master") + if os.path.exists(gpr_m_d): + shutil.rmtree(gpr_m_d) + #pyemu.os_utils.start_workers(gpr_t_d, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=gpr_m_d, verbose=True, port=port) + shutil.copytree(gpr_t_d,gpr_m_d) + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=gpr_m_d) + + # o1 = pd.read_csv(os.path.join(m_d,case+".{0}.obs_pop.csv".format(max(0,pst.control_data.noptmax)))) + o1 = pd.read_csv(os.path.join(m_d, case + ".pareto.archive.summary.csv")) + o1 = o1.loc[o1.generation == o1.generation.max(), :] + o1 = o1.loc[o1.is_feasible == True, :] + o1 = o1.loc[o1.nsga2_front == 1, :] + + # import matplotlib.pyplot as plt + # o2 = pd.read_csv(os.path.join(gpr_m_d, case + ".{0}.obs_pop.csv".format(max(0, gpst.control_data.noptmax)))) + # fig, ax = plt.subplots(1, 1, figsize=(5, 5)) + # ax.scatter(o1.obj_1, o1.obj_2, c="r", s=10) + # ax.scatter(o2.obj_1, o2.obj_2, c="0.5", s=10, alpha=0.5) + # plt.tight_layout() + # plt.savefig("gpr_{0}_compare_noiter.pdf".format(case)) + # plt.close(fig) + + # now lets try an inner-outer scheme... + + gpst.control_data.noptmax = noptmax_inner + gpst.write(os.path.join(gpr_t_d, case + ".pst"), version=2) + gpr_t_d_iter = gpr_t_d + "_outeriter{0}".format(0) + if os.path.exists(gpr_t_d_iter): + shutil.rmtree(gpr_t_d_iter) + shutil.copytree(gpr_t_d, gpr_t_d_iter) + for iouter in range(1, noptmax_outer + 1): + # run the gpr emulator + gpr_m_d_iter = gpr_t_d_iter.replace("template", "master") + complex_m_d_iter = t_d.replace("template", "master_complex_retrain_outeriter{0}".format(iouter)) + if os.path.exists(gpr_m_d_iter): + shutil.rmtree(gpr_m_d_iter) + shutil.copytree(gpr_t_d_iter,gpr_m_d_iter) + + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=gpr_m_d_iter) + + #pyemu.os_utils.start_workers(gpr_t_d_iter, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=gpr_m_d_iter, verbose=True, port=port) + + o2 = pd.read_csv(os.path.join(gpr_m_d_iter, case + ".{0}.obs_pop.csv".format(gpst.control_data.noptmax))) + + # now run the final dv pop thru the "complex" model + final_gpr_dvpop_fname = os.path.join(gpr_m_d_iter, case + ".archive.dv_pop.csv") + assert os.path.exists(final_gpr_dvpop_fname) + complex_model_dvpop_fname = os.path.join(t_d, "gpr_outeriter{0}_dvpop.csv".format(iouter)) + if os.path.exists(complex_model_dvpop_fname): + os.remove(complex_model_dvpop_fname) + # load the gpr archive and do something clever to pick new points to eval + # with the complex model + dvpop = pd.read_csv(final_gpr_dvpop_fname, index_col=0) + if dvpop.shape[0] > pop_size: + arc_sum = pd.read_csv(os.path.join(gpr_m_d_iter, case + ".pareto.archive.summary.csv")) + as_front_map = {member: front for member, front in zip(arc_sum.member, arc_sum.nsga2_front)} + as_crowd_map = {member: crowd for member, crowd in zip(arc_sum.member, arc_sum.nsga2_crowding_distance)} + as_feas_map = {member: feas for member, feas in zip(arc_sum.member, arc_sum.feasible_distance)} + as_gen_map = {member: gen for member, gen in zip(arc_sum.member, arc_sum.generation)} + + dvpop.loc[:, "front"] = dvpop.index.map(lambda x: as_front_map.get(x, np.nan)) + dvpop.loc[:, "crowd"] = dvpop.index.map(lambda x: as_crowd_map.get(x, np.nan)) + dvpop.loc[:, "feas"] = dvpop.index.map(lambda x: as_feas_map.get(x, np.nan)) + dvpop.loc[:, "gen"] = dvpop.index.map(lambda x: as_gen_map.get(x, np.nan)) + # drop members that have missing archive info + dvpop = dvpop.dropna() + if dvpop.shape[0] > pop_size: + dvpop.sort_values(by=["gen", "feas", "front", "crowd"], ascending=[False, True, True, False], + inplace=True) + dvpop = dvpop.iloc[:pop_size, :] + dvpop.drop(["gen", "feas", "front", "crowd"], axis=1, inplace=True) + + # shutil.copy2(final_gpr_dvpop_fname,complex_model_dvpop_fname) + dvpop.to_csv(complex_model_dvpop_fname) + pst.pestpp_options["mou_dv_population_file"] = os.path.split(complex_model_dvpop_fname)[1] + pst.control_data.noptmax = -1 + pst.write(os.path.join(t_d, case + ".pst"), version=2) + if os.path.exists(complex_m_d_iter): + shutil.rmtree(complex_m_d_iter) + shutil.copytree(t_d,complex_m_d_iter) + #pyemu.os_utils.start_workers(t_d, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=complex_m_d_iter, verbose=True, port=port) + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=complex_m_d_iter) + + # plot the complex model results... + o2 = pd.read_csv(os.path.join(complex_m_d_iter, case + ".pareto.archive.summary.csv")) + o2 = o2.loc[o2.generation == o2.generation.max(), :] + # o2 = o2.loc[o2.is_feasible==True,:] + o2 = o2.loc[o2.nsga2_front == 1, :] + # fig, ax = plt.subplots(1, 1, figsize=(5, 5)) + # ax.scatter(o1.obj_1, o1.obj_2, c="r", s=10, label="full complex") + # ax.scatter(o2.obj_1, o2.obj_2, c="0.5", s=10, alpha=0.5, label="mixed emulated-complex") + # ax.legend(loc="upper right") + # ax.set_xlim(0, 10) + # ax.set_ylim(0, 20) + # plt.tight_layout() + # plt.savefig("gpr_{0}_compare_iterscheme_{1}.pdf".format(case, iouter)) + # plt.close(fig) + + # now add those complex model input-output pop files to the list and retrain + # the gpr + dv_pops.append(os.path.join(complex_m_d_iter, case + ".0.dv_pop.csv")) + obs_pops.append(os.path.join(complex_m_d_iter, case + ".0.obs_pop.csv")) + gpr_t_d_iter = gpr_t_d + "_outeriter{0}".format(iouter) + pyemu.helpers.prep_for_gpr(pst_fname, dv_pops, obs_pops, gpr_t_d=gpr_t_d_iter, nverf=int(pop_size * .1), + plot_fits=True, apply_standard_scalar=False, include_emulated_std_obs=True) + gpst_iter = pyemu.Pst(os.path.join(gpr_t_d_iter, case + ".pst")) + # aggdf = pd.read_csv(os.path.join(gpr_t_d,"gpr_aggregate_training_data.csv"),index_col=0) + # aggdf.index = ["outeriter{0}_member{1}".format(iouter,i) for i in range(aggdf.shape[0])] + #restart_gpr_dvpop_fname = "gpr_restart_dvpop_outeriter{0}.csv".format(iouter) + # aggdf.to_csv(os.path.join(gpr_t_d_iter,restart_gpr_dvpop_fname)) + #shutil.copy2(os.path.join(complex_m_d_iter, case + ".0.dv_pop.csv"), + # os.path.join(gpr_t_d_iter, restart_gpr_dvpop_fname)) + gpst_iter.pestpp_options.pop("mou_dv_population_file",None)# = restart_gpr_dvpop_fname + gpst_iter.control_data.noptmax = gpst.control_data.noptmax + gpst_iter.write(os.path.join(gpr_t_d_iter, case + ".pst"), version=2) + + psum_fname = os.path.join(complex_m_d_iter,case+".pareto.archive.summary.csv") + assert os.path.exists(psum_fname) + psum = pd.read_csv(psum_fname) + #assert 1.0 in psum.obj_1.values + #assert 1.0 in psum.obj_2.values + + +def gpr_zdt1_test(): + import numpy as np + from sklearn.gaussian_process import GaussianProcessRegressor + case = "zdt1" + use_chances = False + m_d = os.path.join(case + "_gpr_baseline") + org_d = os.path.join("utils", case + "_template") + t_d = case + "_template" + if os.path.exists(t_d): + shutil.rmtree(t_d) + shutil.copytree(org_d, t_d) + if os.path.exists(m_d): + shutil.rmtree(m_d) + + pst = pyemu.Pst(os.path.join(t_d, case + ".pst")) + pst.pestpp_options["mou_generator"] = "pso" + if use_chances: + pst.pestpp_options["opt_risk"] = 0.95 + pst.pestpp_options["opt_stack_size"] = 50 + pst.pestpp_options["opt_recalc_chance_every"] = 10000 + pst.pestpp_options["opt_chance_points"] = "single" + else: + pst.pestpp_options["opt_risk"] = 0.5 + + pop_size = 10 + num_workers = 10 + noptmax_full = 5 + + port = 4554 + pst.control_data.noptmax = -1 + pst.pestpp_options["mou_population_size"] = pop_size + pst.pestpp_options["mou_save_population_every"] = 1 + pst.write(os.path.join(t_d, case + ".pst")) + #if not os.path.exists(m_d): + # pyemu.os_utils.start_workers(t_d, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=m_d, verbose=True, port=port) + + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=t_d) + m_d = t_d + dv_pops = [os.path.join(m_d, "{0}.0.dv_pop.csv".format(case))] + obs_pops = [f.replace("dv_", "obs_") for f in dv_pops] + + pst_fname = os.path.join(m_d, case + ".pst") + gpr_t_d = os.path.join(case + "_gpr_template") + pyemu.helpers.prep_for_gpr(pst_fname, dv_pops, obs_pops, gpr_t_d=gpr_t_d, nverf=int(pop_size * .1), \ + plot_fits=True, apply_standard_scalar=False, include_emulated_std_obs=True) + gpst = pyemu.Pst(os.path.join(gpr_t_d, case + ".pst")) + shutil.copy2(os.path.join(m_d, case + ".0.dv_pop.csv"), os.path.join(gpr_t_d, "initial_dv_pop.csv")) + gpst.pestpp_options["mou_dv_population_file"] = "initial_dv_pop.csv" + gpst.control_data.noptmax = noptmax_full + gpst.write(os.path.join(gpr_t_d, case + ".pst"), version=2) + gpr_m_d = gpr_t_d.replace("template", "master") + if os.path.exists(gpr_m_d): + shutil.rmtree(gpr_m_d) + #pyemu.os_utils.start_workers(gpr_t_d, mou_exe_path, case + ".pst", num_workers, worker_root=".", + # master_dir=gpr_m_d, verbose=True, port=port) + pyemu.os_utils.run("{0} {1}.pst".format(mou_exe_path,case),cwd=gpr_t_d) + gpr_m_d = gpr_t_d + + psum_fname = os.path.join(gpr_m_d,case+".pareto.archive.summary.csv") + assert os.path.exists(psum_fname) + psum = pd.read_csv(psum_fname) + print(psum.obj_1.min()) + print(psum.obj_2.min()) + assert psum.obj_1.min() < 0.05 + if __name__ == "__main__": #ppu_geostats_test(".") - while True: - thresh_pars_test() + #gpr_compare_invest() + gpr_constr_test() + #gpr_zdt1_test() + #while True: + # thresh_pars_test() #obs_ensemble_quantile_test() #geostat_draws_test("temp") # ac_draw_test("temp") diff --git a/examples/gpr_emulation_hosaki.ipynb b/examples/gpr_emulation_hosaki.ipynb new file mode 100644 index 000000000..1037fbada --- /dev/null +++ b/examples/gpr_emulation_hosaki.ipynb @@ -0,0 +1,1033 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4a6476de-2712-42c7-98e6-c44b101a2599", + "metadata": {}, + "source": [ + "# A verbose walk thru of gaussian-process-regression emulation for global optimization\n", + "\n", + "In this notebook, we go thru much pain to demonstrate how GPR can be used for so-called \"bayesian optimization\" (BO) with a very naive in-filling strategy. Many of the steps and operations in this notebook are done only for the purposes of hoping to explain how GPR emulation works in the context of global evolutionary optimization.\n", + "\n", + "We will be using the well-known hosaki function for this analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b606d2c1-77a3-4c5a-a86f-daba87f86ffa", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import shutil\n", + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pyemu" + ] + }, + { + "cell_type": "markdown", + "id": "0be48c02-b616-4d9b-a0f9-73367807f359", + "metadata": {}, + "source": [ + "First we need to prepare for the GPR experiments. Let's get the original model directory and copy it" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a93b7992-77dd-47d7-8312-c17657798d4b", + "metadata": {}, + "outputs": [], + "source": [ + "org_d = os.path.join(\"..\",\"autotest\",\"utils\",\"hosaki_template\")\n", + "assert os.path.exists(org_d)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4c697d1-5965-4322-808c-3d7f43b6e67f", + "metadata": {}, + "outputs": [], + "source": [ + "t_d = \"hosaki_template\"\n", + "if os.path.exists(t_d):\n", + " shutil.rmtree(t_d)\n", + "shutil.copytree(org_d,t_d)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "696588ca-673c-4353-bc2b-7b58a8f28f2d", + "metadata": {}, + "outputs": [], + "source": [ + "pst = pyemu.Pst(os.path.join(t_d,\"pest.pst\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af8605ae-1a6a-4732-96c7-0e53492b9f5d", + "metadata": {}, + "outputs": [], + "source": [ + "par = pst.parameter_data\n", + "par" + ] + }, + { + "cell_type": "markdown", + "id": "6a114403-c957-42d8-bb0c-1688835b866a", + "metadata": {}, + "source": [ + "The hosaki function only has two decision variables - its very simple. But in general, GPR-based emulation does not scale to extreme dimensions (like those seem in ensemble-based DA), but well within the space of global (multi-objective) optimization, so its an obvious pairing with pestpp-mou...\n", + "\n", + "The hosaki function only has 1 objective and no constraints and our goal is to minimize this single objective:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea64aa47-fe75-494d-9c36-e28a6aebc43e", + "metadata": {}, + "outputs": [], + "source": [ + "pst.observation_data" + ] + }, + { + "cell_type": "markdown", + "id": "e34038c7-39f6-47e6-8114-aa3c75afe762", + "metadata": {}, + "source": [ + "First, lets sweep over decision variable space so that we can map the \"true\" objective function surface. You would never do this in practice but it helps us understand the GPR emulation process:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6e8e756-9513-4140-922c-76cc7a347a83", + "metadata": {}, + "outputs": [], + "source": [ + "pvals = []\n", + "#feel free to change this to a larger number (like 20-30) to get a higher resolution map\n", + "# we keep it small so that the pyemu testing on github is faster...\n", + "sweep_steps = 15\n", + "for p1 in np.linspace(par.parlbnd.iloc[0],par.parubnd.iloc[0],sweep_steps):\n", + " for p2 in np.linspace(par.parlbnd.iloc[0],par.parubnd.iloc[0],sweep_steps):\n", + " pvals.append([p1,p2])\n", + "pvals = pd.DataFrame(pvals,columns=pst.par_names)\n", + "pvals.to_csv(os.path.join(t_d,\"sweep_in.csv\"))\n", + "pst.pestpp_options[\"ies_par_en\"] = \"sweep_in.csv\"\n", + "pst.control_data.noptmax = -1\n", + "sweep_d = os.path.join(\"hosaki_sweep\")\n", + "pst.pestpp_options[\"ies_include_base\"] = False\n", + "\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n", + "port = 5544\n", + "num_workers = 50\n", + "\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=sweep_d,worker_root='.',\n", + " port=port)\n" + ] + }, + { + "cell_type": "markdown", + "id": "77d54f6a-3ee6-4604-97c5-5875ee83d56e", + "metadata": {}, + "source": [ + "And lets load up the sweep results so we can plot:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "182ae0b3-cb09-41ae-8699-c29299b6170c", + "metadata": {}, + "outputs": [], + "source": [ + "sweep_pe = pd.read_csv(os.path.join(sweep_d,\"pest.0.par.csv\"),index_col=0)\n", + "sweep_oe = pd.read_csv(os.path.join(sweep_d,\"pest.0.obs.csv\"),index_col=0)\n", + "sweep_x = sweep_pe.loc[:,pst.par_names[0]].values.reshape(sweep_steps,sweep_steps)\n", + "sweep_y = sweep_pe.loc[:,pst.par_names[1]].values.reshape(sweep_steps,sweep_steps)\n", + "sweep_z = sweep_oe.loc[:,pst.obs_names[0]].values.reshape(sweep_steps,sweep_steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30530fcc-e87c-415c-94d6-fa6461cee3a9", + "metadata": {}, + "outputs": [], + "source": [ + "def get_obj_map(ax,_sweep_x,_sweep_y,_sweep_z,label=\"objective function\",\n", + " levels=[-2,-1,0,0.5],vmin=-2,vmax=0.5,cmap=\"magma\"): \n", + " \"\"\"a simple function to plot an objective function surface\"\"\"\n", + " cb = ax.pcolormesh(_sweep_x,_sweep_y,_sweep_z,vmin=vmin,vmax=vmax,cmap=cmap)\n", + " plt.colorbar(cb,ax=ax,label=label)\n", + " ax.contour(_sweep_x,_sweep_y,_sweep_z,levels=levels,colors=\"w\")\n", + " \n", + " ax.set_aspect(\"equal\")\n", + " return ax" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1593010d-1caf-4ea6-926e-57d3b468c956", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "_ = get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.set_title(\"truth\",loc=\"left\")\n", + "plt.show()\n", + "plt.close(fig)" + ] + }, + { + "cell_type": "markdown", + "id": "f83a0143-3fed-44ef-9f94-60073fbf4a33", + "metadata": {}, + "source": [ + "So that is the true objective function surface with the global minimum at (4,2) and a local minimum at (1,2).\n", + "\n", + "The GPR BO workflow starts by first evaluating an initial population with the complex and expensive to run \"model\", which here is just the simple hosaki function. In practice, this will be an expensive and complex process based model. To make this demo more interesting, we will limit the decision variable search space to not include the global minimum (and really 10 members in the initial population is pretty crazy small - just for demo!). If we dont limit the search space and use a larger initial population, the GPR process nails it in the first go..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4db01f0d-71bd-4cdd-bd0d-d5e870861e19", + "metadata": {}, + "outputs": [], + "source": [ + "m_d = \"hosaki_model_master\"\n", + "if os.path.exists(m_d):\n", + " shutil.rmtree(m_d)\n", + "pst.pestpp_options[\"mou_population_size\"] = 10\n", + "pst.control_data.noptmax = -1\n", + "par = pst.parameter_data\n", + "par.loc[pst.par_names[0],\"parubnd\"] = 3\n", + "par.loc[pst.par_names[0],\"parval1\"] = 1.5\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n", + "num_workers = 10\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=m_d,worker_root='.',\n", + " port=port)\n" + ] + }, + { + "cell_type": "markdown", + "id": "7a4482a5-4315-425a-8d62-1558ecc38a09", + "metadata": {}, + "source": [ + "Now load the initial training input-output pairs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8875b40-bfa1-46ad-91d1-9b68338dc03e", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop_fname = os.path.join(m_d,\"pest.0.dv_pop.csv\")\n", + "training_opop_fname = os.path.join(m_d,\"pest.0.obs_pop.csv\")" + ] + }, + { + "cell_type": "markdown", + "id": "e828f158-c1a0-4ce6-b9fa-0095738b333e", + "metadata": {}, + "source": [ + "and we can plot them to visual where we have training points on the actual true surface" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d0a1063-66ee-42c9-a7d5-faf303509c69", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop = pd.read_csv(training_dvpop_fname,index_col=0)\n", + "training_opop = pd.read_csv(training_opop_fname,index_col=0)\n", + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + "ax.set_title(\"true objective function surface with training points\")\n", + "plt.show()\n", + "plt.close(fig)\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "89d437f7-1020-4644-ad30-3170a3dd37fc", + "metadata": {}, + "source": [ + "Those white triangles are the places where we have actually run the model - in practice we will never have the colormap or the contours - those are just to help us understand what is happening...and remember: that is a purposefully sparse training dataset!\n", + "\n", + "This function is where the magic happens: we use the decision variable population and resulting output/observation population to setup the GPR emulation process. Note that `pyemu.helpers.prep_for_gpr()` function setups and training a GPR emulator for each non-zero weighted observation in the control file... " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1c3ae81-26bb-41c8-bccc-6c32c55c100d", + "metadata": {}, + "outputs": [], + "source": [ + "def prep_for_gpr(dvpops,opops,gpr_t_d,noptmax=-1):\n", + " # this helper does the heavy lifting. We are including the optional GPR emulated standard deviations\n", + " # as \"observations\". These quantities are the standard deviation of the emulated objective function value\n", + " # with respect to the GPR process - this is a unique and very powerful aspect of GPR vs other\n", + " #data-driven techniques. We will use these obs to visualize the uncertainty in the \n", + " # GPR objectve function surface later...\n", + " pyemu.helpers.prep_for_gpr(os.path.join(t_d,\"pest.pst\"),dvpops,opops,gpr_t_d=gpr_t_d,plot_fits=True,\n", + " include_emulated_std_obs=True)\n", + " # now load the newly created gpr-based pst file:\n", + " gpst = pyemu.Pst(os.path.join(gpr_t_d,\"pest.pst\"))\n", + " #and copy the sweep input file that lets us run the full response surface\n", + " # in practice, you dont need this either...\n", + " shutil.copy2(os.path.join(sweep_d,\"sweep_in.csv\"),os.path.join(gpr_t_d,\"sweep_in.csv\"))\n", + " # some bits and bobs:\n", + " gpst.control_data.noptmax = noptmax\n", + " gpst.pestpp_options[\"ies_include_base\"] = False\n", + " gpst.pestpp_options[\"mou_save_population_every\"] = 1\n", + " gpst.pestpp_options.pop(\"mou_dv_population_file\",None)\n", + " par = gpst.parameter_data\n", + " par.loc[:,\"parlbnd\"] = 0.0\n", + " par.loc[:,\"parubnd\"] = 5.0\n", + " par.loc[:,\"parval1\"] = 2.5 \n", + " \n", + " gpst.write(os.path.join(gpr_t_d,\"pest.pst\"),version=2)\n", + " return gpst" + ] + }, + { + "cell_type": "markdown", + "id": "3d2f2533-8920-45cd-b33c-e082d69cb0d9", + "metadata": {}, + "source": [ + "One important consideration here: In this simple demo, we only have one objective. That means we have to be careful not to drive the emulator to a single best-fit solution (which we do by using only a few generations when we run mou), but we also need to sure that we start the emulator-driven mou runs with a fully dispersed initial population. But! In a multiobjective run, the final generation for mou is naturally disperesed because of the pareto frontier search, so in the multiobjective setting, we probably want to start the emulated mou run with the existing population to keep us from having to research all of decision variable space for the nondominated solutions. This isnt a hard-and-fast rule, but it seems to be generally applicable\n", + "\n", + "The GPR helper in pyemu accepts a list of input and output filenames to make it easier to repeatedly retrain the GPR emulators:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "251d15a4-c695-4340-81e3-8c7600f6d198", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop_fnames = [training_dvpop_fname]\n", + "training_opop_fnames = [training_opop_fname]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaa524df-9a1e-4975-8c0d-9774e077839f", + "metadata": {}, + "outputs": [], + "source": [ + "gpr_t_d = t_d + \"_gpr\"\n", + "gpst = prep_for_gpr(training_dvpop_fnames,training_opop_fnames,gpr_t_d)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "208fd361-7cec-4205-9b09-dd4af6b8a624", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.observation_data" + ] + }, + { + "cell_type": "markdown", + "id": "0befe6ee-4a00-4874-be8c-66ce6118730d", + "metadata": {}, + "source": [ + "We see the original \"sim\" observation (which is the objective function we are trying to minimize), but we also see the optional GPR-estimated standard devivation\n", + "\n", + "Let's see what was created in the new GPR-based template directory:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d733ad63-126e-47f7-9bec-b28e1ab758cd", + "metadata": {}, + "outputs": [], + "source": [ + "os.listdir(gpr_t_d)" + ] + }, + { + "cell_type": "markdown", + "id": "2c83d09e-98e0-4025-bb05-28decb6b26b2", + "metadata": {}, + "source": [ + "two very important files in that dir: the new forward run python script and a series of pickle files, one per objective and one per active constraint in the optimization problem. Essentially, we need to build a GPR-based emulator for each output that is relavent to the optimization problem. We also dont want to rebuild these emulators everytime we run the model, so we store the trained GPR emulators in pickle files and load them up as needed:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b0086a7-fefb-4172-be2b-93b8ce73cba0", + "metadata": {}, + "outputs": [], + "source": [ + "with open(os.path.join(gpr_t_d,\"forward_run.py\"),'r') as f:\n", + " for line in f:\n", + " print(line,end=\"\")" + ] + }, + { + "cell_type": "markdown", + "id": "4b60a402-d9bb-419c-92ca-df0cc8919945", + "metadata": {}, + "source": [ + "So we simply loop over all relavent model outputs that have a GPR emulator and \"emulate\" the value of the model output given the current decision variable values. easy as!" + ] + }, + { + "cell_type": "markdown", + "id": "54f493fb-a1f7-4248-9a5b-b884a7ab9a84", + "metadata": {}, + "source": [ + "Now, just for learning, lets sweep over decision variable space but evaluate the GPR emulated objective function value. You would never do this in practice, but it can be informative:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "595ea30a-df3b-48f4-9375-82dff4f77e58", + "metadata": {}, + "outputs": [], + "source": [ + "gpr_sweep_d = sweep_d+\"_gpr\"\n", + "num_workers = 50\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_sweep_d,worker_root='.',\n", + " port=port)\n" + ] + }, + { + "cell_type": "markdown", + "id": "d7ecb069-da8c-4d6a-a968-2f047ad672d9", + "metadata": {}, + "source": [ + "Just a helper function to visual the results of the GPR (re-)training:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "406fab70-18e1-43e8-8edb-ffea73348dcc", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_gpr_sweep_results(_gpr_sweep_d,_training_dvpop):\n", + " #load the gpr sweep results to viz the emulated objective function surface\n", + " sweep_gpr_pe = pd.read_csv(os.path.join(_gpr_sweep_d,\"pest.0.par.csv\"),index_col=0)\n", + " sweep_gpr_oe = pd.read_csv(os.path.join(_gpr_sweep_d,\"pest.0.obs.csv\"),index_col=0)\n", + " gpr_sweep_x = sweep_gpr_pe.loc[:,gpst.par_names[0]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_y = sweep_gpr_pe.loc[:,gpst.par_names[1]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_z = sweep_gpr_oe.loc[:,gpst.obs_names[0]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_stdev_z = sweep_gpr_oe.loc[:,gpst.obs_names[1]].values.reshape(sweep_steps,sweep_steps)\n", + " \n", + " # plot it up:\n", + " fig, axes = plt.subplots(2,2,figsize=(10,8))\n", + " axes = axes.flatten()\n", + " get_obj_map(axes[0],sweep_x,sweep_y,sweep_z)\n", + " get_obj_map(axes[1],gpr_sweep_x,gpr_sweep_y,gpr_sweep_z)\n", + " axes[1].scatter(_training_dvpop.loc[:,pst.par_names[0]],_training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + " diff = sweep_z-gpr_sweep_z\n", + " amax = np.abs(diff).max()\n", + " get_obj_map(axes[2],gpr_sweep_x,gpr_sweep_y,sweep_z-gpr_sweep_z,label=\"truth minus emulated\",levels=None,vmin=-amax,vmax=amax,cmap=\"bwr\")\n", + " get_obj_map(axes[3],gpr_sweep_x,gpr_sweep_y,gpr_sweep_stdev_z,label=\"GPR stdev\",levels=None,\n", + " vmin=gpr_sweep_stdev_z.min(),vmax=gpr_sweep_stdev_z.max(),cmap=\"jet\")\n", + " axes[2].scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + " \n", + " axes[0].set_title(\"truth\",loc=\"left\")\n", + " axes[1].set_title(\"emulated with training points\",loc=\"left\")\n", + " axes[2].set_title(\"difference with training points\",loc=\"left\")\n", + " axes[3].set_title(\"GPR standard deviation\",loc=\"left\")\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa665546-50f3-4f59-8e89-73a17760f6d9", + "metadata": {}, + "outputs": [], + "source": [ + "plot_gpr_sweep_results(gpr_sweep_d,training_dvpop)" + ] + }, + { + "cell_type": "markdown", + "id": "9668435c-ccb4-486b-8423-d864251d55d0", + "metadata": {}, + "source": [ + "OK! now we can see whats happening - the emulated objective function surface is strongly controlled by the location of the training data points, and, is this case, its not a good representation of the truth surface...yet...it should also be clear that the uncertainty in the GPR emulation is lowest near the training points but is highly uncertain as we move away from the training - just like geostatistics!\n", + "\n", + "But now lets run pestpp-mou on the GPR emulated model. This is usually quite fast, especially if the process model that is being emulated takes more than a few mins to run..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6d5d98d-3143-4833-a6de-13a3319d0ca2", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.control_data.noptmax = 8\n", + "gpst.pestpp_options[\"mou_population_size\"] = 20\n", + "gpst.write(os.path.join(gpr_t_d,\"pest.pst\"))\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "num_workers = 20\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_m_d,worker_root='.',\n", + " port=port)" + ] + }, + { + "cell_type": "markdown", + "id": "fd89d01b-d4d2-4f09-a008-c10ae9ab5bf1", + "metadata": {}, + "source": [ + "Now lets plot how the pestpp-mou population evolves across the emulated surface, but also plot it on the true surface just to help us understand what is happening" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a769de5f-58b0-432a-a9d6-418ea7ee2586", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_gpr_mou_results(_gpr_m_d,_gpr_sweep_d):\n", + " sweep_gpr_pe = pd.read_csv(os.path.join(_gpr_sweep_d,\"pest.0.par.csv\"),index_col=0)\n", + " sweep_gpr_oe = pd.read_csv(os.path.join(_gpr_sweep_d,\"pest.0.obs.csv\"),index_col=0)\n", + " gpr_sweep_x = sweep_gpr_pe.loc[:,gpst.par_names[0]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_y = sweep_gpr_pe.loc[:,gpst.par_names[1]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_z = sweep_gpr_oe.loc[:,gpst.obs_names[0]].values.reshape(sweep_steps,sweep_steps)\n", + " gpr_sweep_stdev_z = sweep_gpr_oe.loc[:,gpst.obs_names[1]].values.reshape(sweep_steps,sweep_steps)\n", + " \n", + " gpr_dvpops = [os.path.join(_gpr_m_d,f) for f in os.listdir(_gpr_m_d) if len(f.split('.')) == 4 and f.endswith(\"dv_pop.csv\") and \"archive\" not in f]\n", + " gpr_dvpops_itr = [int(f.split(\".\")[1]) for f in gpr_dvpops]\n", + " gpr_dvpops = {itr:pd.read_csv(f,index_col=0) for itr,f in zip(gpr_dvpops_itr,gpr_dvpops)}\n", + " for itr in range(max(gpr_dvpops_itr)):\n", + " fig,axes = plt.subplots(1,2,figsize=(10,4))\n", + " ax = axes[0]\n", + " get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + " ax.scatter(gpr_dvpops[itr].loc[:,gpst.par_names[0]],gpr_dvpops[itr].loc[:,gpst.par_names[1]],marker='.',c='w',s=10)\n", + " ax.set_title(\"truth generation {0}\".format(itr),loc=\"left\")\n", + " ax = axes[1]\n", + " get_obj_map(ax,gpr_sweep_x,gpr_sweep_y,gpr_sweep_z)\n", + " ax.scatter(gpr_dvpops[itr].loc[:,gpst.par_names[0]],gpr_dvpops[itr].loc[:,gpst.par_names[1]],marker='.',c='w',s=10)\n", + " ax.set_title(\"emulated generation {0}\".format(itr),loc=\"left\")\n", + " plt.show()\n", + " plt.close(fig)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a29183f3-3c7c-40de-a068-5035a5cf7bb8", + "metadata": {}, + "outputs": [], + "source": [ + "plot_gpr_mou_results(gpr_m_d,gpr_sweep_d)" + ] + }, + { + "cell_type": "markdown", + "id": "738f6e14-b896-48db-8e1e-59dcbf270245", + "metadata": {}, + "source": [ + "So just what you expected? Essentially pestpp-mou converged to the minimum of the objective function we gave it, which is the emulated objective function...at this stage the emulated objective function is a poor represetation of the truth objective function....\n", + "\n", + "Now this where some more cleverness happens: Lets take that last emulated decision variable population and actually run it thru the complex \"model\" (which in this case is just the hosaki function...). This is so that we can \"in-fill\" our GPR emulator with this new points in decision variable space. In practice, a lot more cleverness needs to happen to actually decide which points, but for this lil demo, it works..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a9a889e-fa11-40d2-a3f8-da29aa7c17d6", + "metadata": {}, + "outputs": [], + "source": [ + "gpr_dvpops = [os.path.join(gpr_m_d,f) for f in os.listdir(gpr_m_d) if len(f.split('.')) == 4 and f.endswith(\"dv_pop.csv\") and \"archive\" not in f]\n", + "gpr_dvpops_itr = [int(f.split(\".\")[1]) for f in gpr_dvpops]\n", + "gpr_dvpops = {itr:pd.read_csv(f,index_col=0) for itr,f in zip(gpr_dvpops_itr,gpr_dvpops)}\n", + "gpr_dvpops[max(gpr_dvpops_itr)].to_csv(os.path.join(t_d,\"retrain_1_dvpop.csv\"))\n", + "pst.pestpp_options[\"mou_dv_population_file\"] = \"retrain_1_dvpop.csv\"\n", + "pst.control_data.noptmax = -1\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35a5e11b-9ece-4b87-ba75-d402d113dbd7", + "metadata": {}, + "outputs": [], + "source": [ + "m_d += \"_1\"\n", + "num_workers = 10\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=m_d,worker_root='.',\n", + " port=port)" + ] + }, + { + "cell_type": "markdown", + "id": "e6fbb6ab-6457-4ed9-bba6-42bdeb9d28a8", + "metadata": {}, + "source": [ + "Now load the newly evaluated training points:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26d65204-8abb-4161-baa5-1e5c21d24f27", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop_fname1 = os.path.join(m_d,\"pest.0.dv_pop.csv\")\n", + "training_opop_fname1 = os.path.join(m_d,\"pest.0.obs_pop.csv\")" + ] + }, + { + "cell_type": "markdown", + "id": "69f6a35c-9492-4a67-9b53-0f612382a345", + "metadata": {}, + "source": [ + "And combine these new points with the points we already had (the original few places where we ran the model at the beginning):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1018f548-ec0a-446d-a3a3-43cf7723bb64", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop1 = pd.read_csv(training_dvpop_fname1,index_col=0)\n", + "training_opop1 = pd.read_csv(training_opop_fname1,index_col=0)\n", + "training_dvpop = pd.concat([training_dvpop,training_dvpop1])\n", + "training_opop = pd.concat([training_opop,training_opop1])\n", + "\n", + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + "plt.show()\n", + "plt.close(fig)" + ] + }, + { + "cell_type": "markdown", + "id": "225d90ef-2cff-4799-b745-2e9a033606f7", + "metadata": {}, + "source": [ + "See how we have combined the original points and the new points?\n", + "\n", + "Now lets re-do the GPR training process with this combined decision variable and output/observation population" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8046665-4e0b-4f0f-a197-b0be3b4f7a4c", + "metadata": {}, + "outputs": [], + "source": [ + "training_dvpop_fnames.append(training_dvpop_fname1)\n", + "training_opop_fnames.append(training_opop_fname1)\n", + "gpr_t_d = t_d + \"_gpr1\"\n", + "gpst = prep_for_gpr(training_dvpop_fnames,training_opop_fnames,gpr_t_d)" + ] + }, + { + "cell_type": "markdown", + "id": "ff2b2dac-5acb-4f27-b0a7-0b4c6e419222", + "metadata": {}, + "source": [ + "Ok, now let's also re-sweep over the emulated surface to see how much better our approximation to the true objective function surface (remember dont do this in practice! we are just doing it here to help build understanding):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1709b80c-62d5-4e6c-ad6f-9b7c565ac78f", + "metadata": {}, + "outputs": [], + "source": [ + "gpr_sweep_d = sweep_d+\"_gpr1\"\n", + "num_workers = 50\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_sweep_d,worker_root='.',\n", + " port=port)" + ] + }, + { + "cell_type": "markdown", + "id": "5b48ab2c-e132-4e2b-a387-637087a558b4", + "metadata": {}, + "source": [ + "And now visualize these new surfaces with the current training dataset:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4070d242-2405-4f4f-979a-d19777477d0b", + "metadata": {}, + "outputs": [], + "source": [ + "plot_gpr_sweep_results(gpr_sweep_d,training_dvpop)" + ] + }, + { + "cell_type": "markdown", + "id": "9d986519-b800-444f-965a-7dff9a589b9a", + "metadata": {}, + "source": [ + "Comparing this to the original emulated surface, we can see that our emulated objective function surface has improved dramatically!\n", + "\n", + "Now let's run pestpp-mou on this improved emulated model and visualize how pestpp-mou navigates this new emulated surface:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e3cef15-0ac7-4ffb-9cd4-b7106106b954", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.control_data.noptmax = 8\n", + "gpst.pestpp_options[\"mou_population_size\"] = 20\n", + "gpst.write(os.path.join(gpr_t_d,\"pest.pst\"))\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "num_workers = 20\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_m_d,worker_root='.',\n", + " port=port)\n", + "plot_gpr_mou_results(gpr_m_d,gpr_sweep_d)" + ] + }, + { + "cell_type": "markdown", + "id": "49d72679-5b94-4a04-bd1a-43f4bf812f45", + "metadata": {}, + "source": [ + "Now lets do this whole process a few times:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3e78ef6-1bf0-4c05-9688-82c578f0ab77", + "metadata": {}, + "outputs": [], + "source": [ + "i = 2\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i-1)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "gpr_dvpops = [os.path.join(gpr_m_d,f) for f in os.listdir(gpr_m_d) if len(f.split('.')) == 4 and f.endswith(\"dv_pop.csv\") and \"archive\" not in f]\n", + "gpr_dvpops_itr = [int(f.split(\".\")[1]) for f in gpr_dvpops]\n", + "gpr_dvpops = {itr:pd.read_csv(f,index_col=0) for itr,f in zip(gpr_dvpops_itr,gpr_dvpops)}\n", + "gpr_dvpops[max(gpr_dvpops_itr)].to_csv(os.path.join(t_d,\"retrain_{0}_dvpop.csv\".format(i)))\n", + "pst.pestpp_options[\"mou_dv_population_file\"] = \"retrain_{0}_dvpop.csv\".format(i)\n", + "pst.control_data.noptmax = -1\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n", + "m_d = \"hosaki_model_{0}\".format(i)\n", + "num_workers = 10\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "training_dvpop_fnamei = os.path.join(m_d,\"pest.0.dv_pop.csv\")\n", + "training_opop_fnamei = os.path.join(m_d,\"pest.0.obs_pop.csv\")\n", + "\n", + "training_dvpopi = pd.read_csv(training_dvpop_fnamei,index_col=0)\n", + "training_opopi = pd.read_csv(training_opop_fnamei,index_col=0)\n", + "training_dvpop = pd.concat([training_dvpop,training_dvpopi])\n", + "training_opop = pd.concat([training_opop,training_opopi])\n", + "\n", + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + "plt.show()\n", + "plt.close(fig)\n", + "\n", + "training_dvpop_fnames.append(training_dvpop_fnamei)\n", + "training_opop_fnames.append(training_opop_fnamei)\n", + "\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "gpst = prep_for_gpr(training_dvpop_fnames,training_opop_fnames,gpr_t_d)\n", + "gpr_sweep_d = sweep_d+\"_gpr{0}\".format(i)\n", + "num_workers = 50\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_sweep_d,worker_root='.',\n", + " port=port)\n", + "\n", + "plot_gpr_sweep_results(gpr_sweep_d,training_dvpop)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c4e058c1-f67d-4280-a9fe-e9a5a4de3e60", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.control_data.noptmax = 8\n", + "gpst.pestpp_options[\"mou_population_size\"] = 20\n", + "gpst.write(os.path.join(gpr_t_d,\"pest.pst\"))\n", + "\n", + "num_workers = 20\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "\n", + "\n", + "plot_gpr_mou_results(gpr_m_d,gpr_sweep_d)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6ec140f-024a-444d-aa7b-766d6733b26a", + "metadata": {}, + "outputs": [], + "source": [ + "i = 3\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i-1)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "gpr_dvpops = [os.path.join(gpr_m_d,f) for f in os.listdir(gpr_m_d) if len(f.split('.')) == 4 and f.endswith(\"dv_pop.csv\") and \"archive\" not in f]\n", + "gpr_dvpops_itr = [int(f.split(\".\")[1]) for f in gpr_dvpops]\n", + "gpr_dvpops = {itr:pd.read_csv(f,index_col=0) for itr,f in zip(gpr_dvpops_itr,gpr_dvpops)}\n", + "gpr_dvpops[max(gpr_dvpops_itr)].to_csv(os.path.join(t_d,\"retrain_{0}_dvpop.csv\".format(i)))\n", + "pst.pestpp_options[\"mou_dv_population_file\"] = \"retrain_{0}_dvpop.csv\".format(i)\n", + "pst.control_data.noptmax = -1\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n", + "m_d = \"hosaki_model_{0}\".format(i)\n", + "num_workers = 10\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "training_dvpop_fnamei = os.path.join(m_d,\"pest.0.dv_pop.csv\")\n", + "training_opop_fnamei = os.path.join(m_d,\"pest.0.obs_pop.csv\")\n", + "\n", + "training_dvpopi = pd.read_csv(training_dvpop_fnamei,index_col=0)\n", + "training_opopi = pd.read_csv(training_opop_fnamei,index_col=0)\n", + "training_dvpop = pd.concat([training_dvpop,training_dvpopi])\n", + "training_opop = pd.concat([training_opop,training_opopi])\n", + "\n", + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + "plt.show()\n", + "plt.close(fig)\n", + "\n", + "training_dvpop_fnames.append(training_dvpop_fnamei)\n", + "training_opop_fnames.append(training_opop_fnamei)\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "\n", + "gpst = prep_for_gpr(training_dvpop_fnames,training_opop_fnames,gpr_t_d)\n", + "gpr_sweep_d = sweep_d+\"_gpr{0}\".format(i)\n", + "num_workers = 50\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_sweep_d,worker_root='.',\n", + " port=port)\n", + "\n", + "plot_gpr_sweep_results(gpr_sweep_d,training_dvpop)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6dafb39f-5b30-47ed-9e89-9bf89fb47797", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.control_data.noptmax = 8\n", + "gpst.pestpp_options[\"mou_population_size\"] = 20\n", + "gpst.write(os.path.join(gpr_t_d,\"pest.pst\"))\n", + "\n", + "num_workers = 20\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "\n", + "\n", + "plot_gpr_mou_results(gpr_m_d,gpr_sweep_d)" + ] + }, + { + "cell_type": "markdown", + "id": "7ab00c33-17a8-48a6-ae19-7e5c63136bc1", + "metadata": {}, + "source": [ + "And one last time...and this time with more emulator-based pestpp-mou generations to polish things off since we no longer want to keep a dispersed population (since we arent doing any more in-filling and retraining)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cef02969-f1b9-431c-a1ef-cca2612d8a3b", + "metadata": {}, + "outputs": [], + "source": [ + "i = 4\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i-1)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "gpr_dvpops = [os.path.join(gpr_m_d,f) for f in os.listdir(gpr_m_d) if len(f.split('.')) == 4 and f.endswith(\"dv_pop.csv\") and \"archive\" not in f]\n", + "gpr_dvpops_itr = [int(f.split(\".\")[1]) for f in gpr_dvpops]\n", + "gpr_dvpops = {itr:pd.read_csv(f,index_col=0) for itr,f in zip(gpr_dvpops_itr,gpr_dvpops)}\n", + "gpr_dvpops[max(gpr_dvpops_itr)].to_csv(os.path.join(t_d,\"retrain_{0}_dvpop.csv\".format(i)))\n", + "pst.pestpp_options[\"mou_dv_population_file\"] = \"retrain_{0}_dvpop.csv\".format(i)\n", + "pst.control_data.noptmax = -1\n", + "pst.write(os.path.join(t_d,\"pest.pst\"))\n", + "m_d = \"hosaki_model_{0}\".format(i)\n", + "num_workers = 10\n", + "pyemu.os_utils.start_workers(t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "training_dvpop_fnamei = os.path.join(m_d,\"pest.0.dv_pop.csv\")\n", + "training_opop_fnamei = os.path.join(m_d,\"pest.0.obs_pop.csv\")\n", + "\n", + "training_dvpopi = pd.read_csv(training_dvpop_fnamei,index_col=0)\n", + "training_opopi = pd.read_csv(training_opop_fnamei,index_col=0)\n", + "training_dvpop = pd.concat([training_dvpop,training_dvpopi])\n", + "training_opop = pd.concat([training_opop,training_opopi])\n", + "\n", + "fig,ax = plt.subplots(1,1,figsize=(6,5))\n", + "get_obj_map(ax,sweep_x,sweep_y,sweep_z)\n", + "ax.scatter(training_dvpop.loc[:,pst.par_names[0]],training_dvpop.loc[:,pst.par_names[1]],marker='^',c='w',s=20)\n", + "plt.show()\n", + "plt.close(fig)\n", + "\n", + "training_dvpop_fnames.append(training_dvpop_fnamei)\n", + "training_opop_fnames.append(training_opop_fnamei)\n", + "gpr_t_d = t_d + \"_gpr{0}\".format(i)\n", + "gpr_m_d = gpr_t_d.replace(\"template\",\"master\")\n", + "gpst = prep_for_gpr(training_dvpop_fnames,training_opop_fnames,gpr_t_d)\n", + "gpr_sweep_d = sweep_d+\"_gpr{0}\".format(i)\n", + "num_workers = 50\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-ies\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_sweep_d,worker_root='.',\n", + " port=port)\n", + "\n", + "plot_gpr_sweep_results(gpr_sweep_d,training_dvpop)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2c52501-30e5-4140-879c-d7a6bbfe565e", + "metadata": {}, + "outputs": [], + "source": [ + "gpst.control_data.noptmax = 15\n", + "gpst.pestpp_options[\"mou_population_size\"] = 20\n", + "gpst.write(os.path.join(gpr_t_d,\"pest.pst\"))\n", + "\n", + "num_workers = 20\n", + "pyemu.os_utils.start_workers(gpr_t_d,\"pestpp-mou\",\"pest.pst\",\n", + " num_workers=num_workers,\n", + " master_dir=gpr_m_d,worker_root='.',\n", + " port=port)\n", + "\n", + "plot_gpr_mou_results(gpr_m_d,gpr_sweep_d)" + ] + }, + { + "cell_type": "markdown", + "id": "5aaeba17-14dc-4da9-8528-e64445ac7e41", + "metadata": {}, + "source": [ + "There you have it - after evaluating only a few populations with the complex process-based model (here that was only the hosaki function), we can effectively optimize a function with the otherwise every expensive global solvers in pestpp-mou. In general, this approach can reduce the computational demand of mou by 10X to 100X, which is pretty amazing...." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31eeeee5-f54e-4dfd-8ba7-717108066370", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d58eb64-1e67-4d01-b7c9-2d024b0806a0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyemu/prototypes/__init__.py b/pyemu/prototypes/__init__.py deleted file mode 100644 index 530cc2883..000000000 --- a/pyemu/prototypes/__init__.py +++ /dev/null @@ -1,13 +0,0 @@ -"""pyEMU: python modules for Environmental Model Uncertainty analyses. These -modules are designed to work directly and seamlessly with PEST and PEST++ model -independent interface. pyEMU can also be used to setup this interface. - -Several forms of uncertainty analyses are support including FOSM-based -analyses (pyemu.Schur and pyemu.ErrVar), Monte Carlo (including -GLUE and null-space Monte Carlo) and a prototype iterative Ensemble Smoother - -""" - -from .ensemble_method import * -from .moouu import * -from .da import * diff --git a/pyemu/pst/pst_utils.py b/pyemu/pst/pst_utils.py index eecab5364..74a02b0ac 100644 --- a/pyemu/pst/pst_utils.py +++ b/pyemu/pst/pst_utils.py @@ -916,7 +916,7 @@ def try_process_output_pst(pst): ins_file, str(e) ) ) - df = _try_run_inschek(ins_file, out_file) + #df = _try_run_inschek(ins_file, out_file) if df is not None: pst.observation_data.loc[df.index, "obsval"] = df.obsval diff --git a/pyemu/utils/helpers.py b/pyemu/utils/helpers.py index c3312a85f..c7de9f913 100644 --- a/pyemu/utils/helpers.py +++ b/pyemu/utils/helpers.py @@ -4032,6 +4032,301 @@ def get_current_prop(_cur_thresh): return thresh, prop +def prep_for_gpr(pst_fname,input_fnames,output_fnames,gpr_t_d="gpr_template",gp_kernel=None,nverf=0, + plot_fits=False,apply_standard_scalar=False, include_emulated_std_obs=False): + """helper function to setup a gaussian-process-regression (GPR) emulator for outputs of interest. This + is primarily targeted at low-dimensional settings like those encountered in PESTPP-MOU + + Parameters: + pst_fname (str): existing pest control filename + input_fnames (str | list[str]): usually a list of decision variable population files + output_fnames (str | list[str]): usually a list of observation population files that + corresponds to the simulation results associated with `input_fnames` + gpr_t_d (str): the template file dir to create that will hold the GPR emulators + gp_kernel (sklearn GaussianProcess kernel): the kernel to use. if None, a standard RBF kernel + is created and used + nverf (int): the number of input-output pairs to hold back for a simple verification test + plot_fits (bool): flag to plot the fit GPRs + apply_standard_scalar (bool): flag to apply sklearn.preprocessing.StandardScaler transform before + training/executing the emulator. Default is False + include_emulated_std_obs (bool): flag to include the estimated standard deviation in the predicted + response of each GPR emulator. If True, additional obserations are added to the GPR pest interface + , one for each nominated observation quantity. Can be very useful for designing in-filling strategies + + Returns: + None + + Note: + requires scikit-learn + + """ + + from sklearn.gaussian_process import GaussianProcessRegressor + from sklearn.gaussian_process.kernels import RBF, Matern, ConstantKernel + from sklearn.preprocessing import StandardScaler + from sklearn.pipeline import Pipeline + + import pickle + import inspect + pst = pyemu.Pst(pst_fname) + + # work out input variable names + input_groups = pst.pestpp_options.get("opt_dec_var_groups",None) + par = pst.parameter_data + if input_groups is None: + print("using all adjustable parameters as inputs") + input_names = pst.adj_par_names + else: + input_groups = set([i.strip() for i in input_groups.lower().strip().split(",")]) + print("input groups:",input_groups) + adj_par = par.loc[pst.adj_par_names,:].copy() + adj_par = adj_par.loc[adj_par.pargp.apply(lambda x: x in input_groups),:] + input_names = adj_par.parnme.tolist() + print("input names:",input_names) + + #work out constraints and objectives + ineq_names = pst.less_than_obs_constraints.tolist() + ineq_names.extend(pst.greater_than_obs_constraints.tolist()) + obs = pst.observation_data + objs = pst.pestpp_options.get("mou_objectives",None) + constraints = [] + + if objs is None: + print("'mou_objectives' not found in ++ options, using all ineq tagged non-zero weighted obs as objectives") + objs = ineq_names + else: + objs = objs.lower().strip().split(',') + constraints = [n for n in ineq_names if n not in objs] + + print("objectives:",objs) + print("constraints:",constraints) + output_names = objs + output_names.extend(constraints) + + print("loading input and output files") + if isinstance(input_fnames,str): + input_fnames = [input_fnames] + if isinstance(output_fnames,str): + output_fnames = [output_fnames] + if len(output_fnames) != len(input_fnames): + raise Exception("len(input_fnames) != len(output_fnames)") + + if os.path.exists(gpr_t_d): + shutil.rmtree(gpr_t_d) + os.makedirs(gpr_t_d) + + dfs = [] + for input_fname,output_fname in zip(input_fnames,output_fnames): + if input_fname.lower().endswith(".csv"): + input_df = pd.read_csv(os.path.join(input_fname),index_col=0) + elif input_fname.lower().endswith(".jcb"): + input_df = pyemu.ParameterEnsemble.from_binary(pst=pst,filename=input_fname)._df + else: + raise Exception("unrecognized input_fname extension:'{0}', looking for csv or jcb".\ + format(input_fname.lower())) + + if output_fname.lower().endswith(".csv"): + output_df = pd.read_csv(os.path.join(output_fname),index_col=0) + elif output_fname.lower().endswith(".jcb"): + output_df = pyemu.ObservationEnsemble.from_binary(pst=pst,filename=output_fname)._df + else: + raise Exception("unrecognized output_fname extension:'{0}', looking for csv or jcb".\ + format(output_fname.lower())) + + if input_df.shape[0] != output_df.shape[0]: + raise Exception("input rows != output rows for {0} and {1}".\ + format(input_fname,output_fname)) + input_df = input_df.loc[:,input_names] + assert input_df.shape == input_df.dropna().shape + + output_df = output_df.loc[:, output_names] + assert output_df.shape == output_df.dropna().shape + + input_df.loc[:,output_names] = output_df.values + dfs.append(input_df) + print("...loaded",input_fname,output_fname) + + df = pd.concat(dfs) + assert df.shape == df.dropna().shape + df.to_csv(os.path.join(gpr_t_d,"gpr_aggregate_training_data.csv")) + print("aggregated training dataset shape",df.shape,"saved to",pst_fname + ".aggresults.csv") + + + if gp_kernel is None: + #gp_kernel = ConstantKernel(constant_value=1.0,constant_value_bounds=(1e-8,1e8)) *\ + # RBF(length_scale=1000.0, length_scale_bounds=(1e-8, 1e8)) + gp_kernel = Matern(length_scale=100.0, length_scale_bounds=(1e-4, 1e4), nu=0.5) + + for hp in gp_kernel.hyperparameters: + print(hp) + + cut = df.shape[0] - nverf + X_train = df.loc[:, input_names].values.copy()[:cut, :] + X_verf = df.loc[:, input_names].values.copy()[cut:, :] + + model_fnames = [] + if plot_fits: + import matplotlib.pyplot as plt + from matplotlib.backends.backend_pdf import PdfPages + pdf = PdfPages(os.path.join(gpr_t_d,"gpr_fits.pdf")) + for output_name in output_names: + + y_verf = df.loc[:,output_name].values.copy()[cut:] + y_train = df.loc[:, output_name].values.copy()[:cut] + print("training GPR for {0} with {1} data points".format(output_name,y_train.shape[0])) + gaussian_process = GaussianProcessRegressor(kernel=gp_kernel, n_restarts_optimizer=20) + if not apply_standard_scalar: + print("WARNING: not applying StandardScalar transformation - user beware!") + pipeline = Pipeline([("gpr",gaussian_process)]) + else: + pipeline = Pipeline([("std_scalar", StandardScaler()), ("gpr", gaussian_process)]) + pipeline.fit(X_train, y_train) + print(output_name,"optimized kernel:",pipeline["gpr"].kernel_) + if plot_fits: + print("...plotting fits for",output_name) + predmean,predstd = pipeline.predict(df.loc[:, input_names].values.copy(), return_std=True) + df.loc[:,"predmean"] = predmean + df.loc[:,"predstd"] = predstd + isverf = np.zeros_like(predmean) + isverf[cut:] = 1 + df.loc[:,"isverf"] = isverf + for input_name in input_names: + fig,ax = plt.subplots(1,1,figsize=(6,6)) + #print(X_train[:,ii]) + #print(y_train) + ddf = df[[input_name,output_name,"predmean","predstd","isverf"]].copy() + ddf.sort_values(by=input_name,inplace=True) + ax.scatter(ddf.loc[ddf.isverf==0,input_name],ddf.loc[ddf.isverf==0,output_name],marker=".",c="r",label="training") + if nverf > 0: + ax.scatter(ddf.loc[ddf.isverf==1,input_name],ddf.loc[ddf.isverf==1,output_name],marker="^",c="c",label="verf") + ax.plot(ddf[input_name],ddf["predmean"],"k--",label="GPR mean") + ax.fill_between(ddf[input_name],ddf["predmean"] - (2*ddf["predstd"]),ddf["predmean"]+(2*ddf["predstd"]),alpha=0.5,fc='0.5',label="+/- 95%") + ax.set_title("input:{0}, output:{1}".format(input_name,output_name),loc="left") + ax.legend() + plt.tight_layout() + pdf.savefig() + + plt.close(fig) + + + + model_fname = os.path.split(pst_fname)[1]+"."+output_name+".pkl" + if os.path.exists(os.path.join(gpr_t_d,model_fname)): + print("WARNING: model_fname '{0}' exists, overwriting...".format(model_fname)) + with open(os.path.join(gpr_t_d,model_fname),'wb') as f: + pickle.dump(pipeline,f) + + model_fnames.append(model_fname) + if nverf > 0: + pred_mean,pred_std = pipeline.predict(X_verf,return_std=True) + vdf = pd.DataFrame({"y_verf":y_verf,"y_pred":pred_mean,"y_pred_std":pred_std}) + verf_fname = os.path.join(gpr_t_d,"{0}_gpr_verf.csv".format(output_name)) + vdf.to_csv(verf_fname) + print("saved ",output_fname,"verfication csv to",verf_fname) + mabs = np.abs(vdf.y_verf - vdf.y_pred).mean() + print("...mean abs error",mabs) + if plot_fits: + pdf.close() + mdf = pd.DataFrame({"output_name":output_names,"model_fname":model_fnames},index=output_names) + minfo_fname = os.path.join(gpr_t_d,"gprmodel_info.csv") + mdf.to_csv(minfo_fname) + print("GPR model info saved to",minfo_fname) + + #write a template file + tpl_fname = os.path.join(gpr_t_d,"gpr_input.csv.tpl") + with open(tpl_fname,'w') as f: + f.write("ptf ~\nparnme,parval1\n") + for input_name in input_names: + f.write("{0},~ {0} ~\n".format(input_name)) + other_pars = list(set(pst.par_names)-set(input_names)) + aux_tpl_fname = None + + if len(other_pars) > 0: + + aux_tpl_fname = os.path.join(gpr_t_d,"aux_par.csv.tpl") + print("writing aux par tpl file: ",aux_tpl_fname) + with open(aux_tpl_fname,'w') as f: + f.write("ptf ~\n") + for input_name in other_pars: + f.write("{0},~ {0} ~\n".format(input_name)) + #write an ins file + ins_fname = os.path.join(gpr_t_d,"gpr_output.csv.ins") + with open(ins_fname,'w') as f: + f.write("pif ~\nl1\n") + for output_name in output_names: + if include_emulated_std_obs: + f.write("l1 ~,~ !{0}! ~,~ !{0}_gprstd!\n".format(output_name)) + else: + f.write("l1 ~,~ !{0}!\n".format(output_name)) + tpl_list = [tpl_fname] + if aux_tpl_fname is not None: + tpl_list.append(aux_tpl_fname) + input_list = [f.replace(".tpl","") for f in tpl_list] + gpst = pyemu.Pst.from_io_files(tpl_list,input_list, + [ins_fname],[ins_fname.replace(".ins","")],pst_path=".") + par_names = pst.par_names + assert len(set(par_names).symmetric_difference(set(gpst.par_names))) == 0 + for col in pst.parameter_data.columns: + # this gross thing is to avoid a future error warning in pandas - + # why is it getting so strict?! isnt python duck-typed? + if col in gpst.parameter_data.columns and\ + gpst.parameter_data.dtypes[col] != pst.parameter_data.dtypes[col]: + gpst.parameter_data[col] = gpst.parameter_data[col].astype(pst.parameter_data.dtypes[col]) + gpst.parameter_data.loc[par_names,col] = pst.parameter_data.loc[par_names,col].values + + for col in pst.observation_data.columns: + # this gross thing is to avoid a future error warning in pandas - + # why is it getting so strict?! isnt python duck-typed? + if col in gpst.observation_data.columns and \ + gpst.observation_data.dtypes[col] != pst.observation_data.dtypes[col]: + gpst.observation_data[col] = gpst.obsveration_data[col].astype(pst.observation_data.dtypes[col]) + gpst.observation_data.loc[output_names,col] = pst.observation_data.loc[output_names,col].values + if include_emulated_std_obs: + stdobs = [o for o in gpst.obs_names if o.endswith("_gprstd")] + assert len(stdobs) > 0 + gpst.observation_data.loc[stdobs,"weight"] = 0.0 + gpst.pestpp_options = pst.pestpp_options + gpst.prior_information = pst.prior_information.copy() + lines = [line[4:] for line in inspect.getsource(gpr_forward_run).split("\n")][1:] + + with open(os.path.join(gpr_t_d, "forward_run.py"), 'w') as f: + + for line in lines: + f.write(line + "\n") + + gpst.control_data.noptmax = 0 + gpst.model_command = "python forward_run.py" + gpst_fname = os.path.split(pst_fname)[1] + gpst.write(os.path.join(gpr_t_d,gpst_fname),version=2) + print("saved gpr pst:",gpst_fname,"in gpr_t_d",gpr_t_d) + try: + pyemu.os_utils.run("pestpp-mou {0}".format(gpst_fname),cwd=gpr_t_d) + except Exception as e: + print("WARNING: pestpp-mou test run failed: {0}".format(str(e))) + gpst.control_data.noptmax = pst.control_data.noptmax + gpst.write(os.path.join(gpr_t_d, gpst_fname), version=2) + + +def gpr_forward_run(): + """the function to evaluate a set of inputs thru the GPR emulators.\ + This function gets added programmatically to the forward run process""" + import os + import pandas as pd + import numpy as np + import pickle + from sklearn.gaussian_process import GaussianProcessRegressor + input_df = pd.read_csv("gpr_input.csv",index_col=0) + mdf = pd.read_csv(os.path.join("gprmodel_info.csv"),index_col=0) + mdf.loc[:,"sim"] = np.nan + mdf.loc[:,"sim_std"] = np.nan + for output_name,model_fname in zip(mdf.output_name,mdf.model_fname): + guassian_process = pickle.load(open(model_fname,'rb')) + sim = guassian_process.predict(np.atleast_2d(input_df.parval1.values),return_std=True) + mdf.loc[output_name,"sim"] = sim[0] + mdf.loc[output_name,"sim_std"] = sim[1] + + mdf.loc[:,["output_name","sim","sim_std"]].to_csv("gpr_output.csv",index=False) + def randrealgen_optimized(nreal, tol=1e-7, max_samples=1000000): """