From 080453deff43f3f170f10cd04134d54995e572db Mon Sep 17 00:00:00 2001 From: Brioch Date: Wed, 11 Nov 2020 22:37:27 +1300 Subject: [PATCH 1/3] Adding support for headers in direct pars and strings in index_cols --- autotest/pst_from_tests.py | 19 +++++++++++++++++++ pyemu/utils/helpers.py | 37 +++++++++++++++++++++++++++++++++---- pyemu/utils/pst_from.py | 28 ++++++++++++++++++++++------ 3 files changed, 74 insertions(+), 10 deletions(-) diff --git a/autotest/pst_from_tests.py b/autotest/pst_from_tests.py index 712251cd6..984cc6b16 100644 --- a/autotest/pst_from_tests.py +++ b/autotest/pst_from_tests.py @@ -1240,6 +1240,24 @@ def mf6_freyberg_direct_test(): list_files = ["freyberg6.wel_stress_period_data_{0}.txt".format(t) for t in range(1, m.nper + 1)] + # make dummy versions with headers + for fl in list_files[0:2]: + with open(os.path.join(template_ws, fl), 'r') as fr: + lines = [line for line in fr] + with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: + fw.write("k i j flux \n") + for line in lines: + fw.write(line) + + fl = "freyberg6.wel_stress_period_data_3.txt" + with open(os.path.join(template_ws, fl), 'r') as fr: + lines = [line for line in fr] + with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: + fw.write("well k i j flux \n") + for i, line in enumerate(lines): + fw.write(f"well{i}" + line) + + list_files.sort() for list_file in list_files: kper = int(list_file.split(".")[1].split('_')[-1]) - 1 @@ -1278,6 +1296,7 @@ def mf6_freyberg_direct_test(): # test par mults are working b_d = os.getcwd() os.chdir(pf.new_d) + pst.write_input_files() try: pyemu.helpers.apply_list_and_array_pars( arr_par_file="mult2model_info.csv", chunk_len=1) diff --git a/pyemu/utils/helpers.py b/pyemu/utils/helpers.py index a92b89f05..bb8ec8f64 100644 --- a/pyemu/utils/helpers.py +++ b/pyemu/utils/helpers.py @@ -3978,6 +3978,7 @@ def apply_genericlist_pars(df): else: storehead = [] # work out if headers are used for index_cols + # big assumption here that int type index cols will not be written as headers index_col_eg = df_mf.index_cols.iloc[-1][0] if isinstance(index_col_eg, str): # TODO: add test for model file with headers @@ -4024,12 +4025,16 @@ def apply_genericlist_pars(df): else: mlts = pd.read_csv(mlt.mlt_file) # get mult index to align with org_data, - # mult idxs will always be written zero based + # mult idxs will always be written zero based if int # if original model files is not zero based need to add 1 add1 = int(mlt.zero_based == False) mlts.index = pd.MultiIndex.from_tuples( - mlts.sidx.apply(lambda x: tuple(add1 + np.array(literal_eval(x)))), - names=mlt.index_cols, + mlts.sidx.apply(lambda x: + [add1 + int(xx.strip()) + if xx.strip().isdigit() + else xx.strip('\'\"') + for xx in x.strip('()').split(',')]), + names=mlt.index_cols ) if mlts.index.nlevels < 2: # just in case only one index col is used mlts.index = mlts.index.get_level_values(0) @@ -4350,7 +4355,8 @@ def build_jac_test_csv(pst, num_steps, par_names=None, forward=True): return df -def _write_df_tpl(filename, df, sep=",", tpl_marker="~", **kwargs): +def _write_df_tpl(filename, df, sep=",", tpl_marker="~", + headerlines=None, **kwargs): """function write a pandas dataframe to a template file. """ @@ -4360,9 +4366,32 @@ def _write_df_tpl(filename, df, sep=",", tpl_marker="~", **kwargs): with open(filename, "w") as f: f.write("ptf {0}\n".format(tpl_marker)) f.flush() + if headerlines is not None: + _add_headerlines(f, headerlines) df.to_csv(f, sep=sep, mode="a", **kwargs) +def _add_headerlines(f, headerlines): + lc = 0 + for key in sorted(headerlines.keys()): + if key > lc: + lc += 1 + continue + # TODO if we want to preserve mid-table comments, + # these lines might help - will also need to + # pass comment_char through so it can be + # used by the apply methods + # to = key - lc + # df.iloc[fr:to].to_csv( + # fp, sep=',', mode='a', header=hheader, # todo - presence of header may cause an issue with this + # **kwargs) + # lc += to - fr + # fr = to + f.write(headerlines[key]) + f.flush() + lc += 1 + + def setup_fake_forward_run( pst, new_pst_name, org_cwd=".", bak_suffix="._bak", new_cwd="." ): diff --git a/pyemu/utils/pst_from.py b/pyemu/utils/pst_from.py index 2c70b5475..1ebb7700e 100644 --- a/pyemu/utils/pst_from.py +++ b/pyemu/utils/pst_from.py @@ -661,6 +661,7 @@ def _par_prep( ) = self._prep_arg_list_lengths( filenames, fmts, seps, skip_rows, index_cols, use_cols ) + storehead = None if index_cols is not None: for filename, sep, fmt, skip in zip(filenames, seps, fmts, skip_rows): file_path = os.path.join(self.new_d, filename) @@ -811,7 +812,8 @@ def _par_prep( file_dict[fnames[j]].shape, ) ) - return index_cols, use_cols, file_dict, fmt_dict, sep_dict, skip_dict + return (index_cols, use_cols, file_dict, fmt_dict, sep_dict, skip_dict, + storehead) def _next_count(self, prefix): if prefix not in self._prefix_count: @@ -1286,6 +1288,8 @@ def add_parameters( ) + " already used for 'direct' parameterization" ) + else: + self.direct_org_files.append(filenames[0]) # Default par data columns used for pst par_data_cols = pyemu.pst_utils.pst_config["par_fieldnames"] self.logger.log( @@ -1348,6 +1352,7 @@ def add_parameters( fmt_dict, sep_dict, skip_dict, + headerlines, # needed for direct pars (need to be in tpl) ) = self._par_prep( filenames, idx_cols, @@ -1474,6 +1479,7 @@ def add_parameters( zero_based=self.zero_based, input_filename=in_fileabs, par_style=par_style, + headerlines=headerlines, # only needed for direct pars ) assert np.mod(len(df), len(use_cols)) == 0.0, ( "Parameter dataframe wrong shape for number of cols {0}" @@ -2061,6 +2067,7 @@ def write_list_tpl( zero_based=True, input_filename=None, par_style="multiplier", + headerlines=None, ): """ Write template files for a list style input. @@ -2124,6 +2131,7 @@ def write_list_tpl( ij_in_idx=ij_in_idx, xy_in_idx=xy_in_idx, zero_based=zero_based, + headerlines=headerlines, ) else: df_tpl = _get_tpl_or_ins_df( @@ -2253,6 +2261,7 @@ def _write_direct_df_tpl( xy_in_idx=None, zero_based=True, gpname=None, + headerlines=None ): """ @@ -2314,7 +2323,8 @@ def _write_direct_df_tpl( inames = ["{0}".format(i) for i in range(len(index_cols))] if not zero_based: - df_ti.loc[:, "sidx"] = df_ti.sidx.apply(lambda x: tuple(xx - 1 for xx in x)) + df_ti.loc[:, "sidx"] = df_ti.sidx.apply( + lambda x: tuple(xx - 1 if isinstance(xx,int) else xx for xx in x)) df_ti.loc[:, "idx_strs"] = df_ti.sidx.apply( lambda x: j.join([fmt.format(iname, xx) for xx, iname in zip(x, inames)]) ).str.replace(" ", "") @@ -2430,8 +2440,13 @@ def _write_direct_df_tpl( direct_tpl_df.loc[:, use_col] = ( df_ti.loc[:, use_col].apply(lambda x: "~ {0} ~".format(x)).values ) - - pyemu.helpers._write_df_tpl(tpl_filename, direct_tpl_df, index=False, header=False) + # nasty assumption to distiguish if we need to write headers + if isinstance(direct_tpl_df.columns[0], str): + header = True + else: + header = False + pyemu.helpers._write_df_tpl(tpl_filename, direct_tpl_df, index=False, + header=header, headerlines=headerlines) return df_ti @@ -2529,8 +2544,9 @@ def _get_tpl_or_ins_df( else: inames = ["{0}".format(i) for i in range(len(index_cols))] - if not zero_based: - df_ti.loc[:, "sidx"] = df_ti.sidx.apply(lambda x: tuple(xx - 1 for xx in x)) + if not zero_based: # only if indices are ints (trying to support strings as par ids) + df_ti.loc[:, "sidx"] = df_ti.sidx.apply( + lambda x: tuple(xx - 1 if isinstance(xx, int) else xx for xx in x)) df_ti.loc[:, "idx_strs"] = df_ti.sidx.apply( lambda x: j.join([fmt.format(iname, xx) for xx, iname in zip(x, inames)]) From 91e66d38d0b64a08728bf03185f9d6552272fd4a Mon Sep 17 00:00:00 2001 From: Brioch Date: Wed, 11 Nov 2020 23:43:45 +1300 Subject: [PATCH 2/3] zero-based adjust tweeks for str in index_cols --- pyemu/utils/helpers.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyemu/utils/helpers.py b/pyemu/utils/helpers.py index 7d52959e5..cea48ac66 100644 --- a/pyemu/utils/helpers.py +++ b/pyemu/utils/helpers.py @@ -4006,10 +4006,11 @@ def apply_genericlist_pars(df): add1 = int(mlt.zero_based == False) mlts.index = pd.MultiIndex.from_tuples( mlts.sidx.apply(lambda x: - [add1 + int(xx.strip()) + [add1 + int(xx) if xx.strip().isdigit() - else xx.strip('\'\"') - for xx in x.strip('()').split(',')]), + else xx.strip('\'\" ') + for xx in x.strip('()').split(',') + if xx]), names=mlt.index_cols ) if mlts.index.nlevels < 2: # just in case only one index col is used From 9358a87fd4d97284347521d98f57cd420cd38587 Mon Sep 17 00:00:00 2001 From: Brioch Date: Thu, 12 Nov 2020 16:42:15 +1300 Subject: [PATCH 3/3] more index_col and direct par tidy Opened up method for constant direct pars Added more test and checks on direct pars and varied tabular file types (headers etc) --- autotest/pst_from_tests.py | 351 ++++++++++++++++++++++++++++++++++--- pyemu/utils/pst_from.py | 65 ++++--- 2 files changed, 369 insertions(+), 47 deletions(-) diff --git a/autotest/pst_from_tests.py b/autotest/pst_from_tests.py index c8163bfad..7391b0bc4 100644 --- a/autotest/pst_from_tests.py +++ b/autotest/pst_from_tests.py @@ -1201,11 +1201,13 @@ def mf6_freyberg_direct_test(): # pf.add_observations('freyberg.hds.dat', insfile='freyberg.hds.dat.ins2', # index_cols='obsnme', use_cols='obsval', prefix='hds') - df = pd.read_csv(os.path.join(tmp_model_ws, "sfr.csv"), index_col=0) + # Add stream flow observation + # df = pd.read_csv(os.path.join(tmp_model_ws, "sfr.csv"), index_col=0) pf.add_observations("sfr.csv", insfile="sfr.csv.ins", index_cols="time", use_cols=["GAGE_1","HEADWATER","TAILWATER"],ofile_sep=",") + # Setup geostruct for spatial pars v = pyemu.geostats.ExpVario(contribution=1.0, a=1000) - gr_gs = pyemu.geostats.GeoStruct(variograms=v,transform="log") + gr_gs = pyemu.geostats.GeoStruct(variograms=v, transform="log") rch_temporal_gs = pyemu.geostats.GeoStruct(variograms=pyemu.geostats.ExpVario(contribution=1.0, a=60)) pf.extra_py_imports.append('flopy') ib = m.dis.idomain[0].array @@ -1213,17 +1215,18 @@ def mf6_freyberg_direct_test(): "rch_recharge": [.5, 1.5]} dts = pd.to_datetime("1-1-2018") + pd.to_timedelta(np.cumsum(sim.tdis.perioddata.array["perlen"]), unit="d") print(dts) + # setup from array style pars for tag, bnd in tags.items(): lb, ub = bnd[0], bnd[1] arr_files = [f for f in os.listdir(tmp_model_ws) if tag in f and f.endswith(".txt")] if "rch" in tag: for arr_file in arr_files: + # indy direct grid pars for each array type file recharge_files = ["recharge_1.txt","recharge_2.txt","recharge_3.txt"] pf.add_parameters(filenames=arr_file, par_type="grid", par_name_base="rch_gr", pargp="rch_gr", zone_array=ib, upper_bound=1.0e-3, lower_bound=1.0e-7, par_style="direct") - - for arr_file in arr_files: + # additional constant mults kper = int(arr_file.split('.')[1].split('_')[-1]) - 1 pf.add_parameters(filenames=arr_file, par_type="constant", par_name_base=arr_file.split('.')[1] + "_cn", @@ -1232,16 +1235,17 @@ def mf6_freyberg_direct_test(): datetime=dts[kper]) else: for arr_file in arr_files: + # grid mults pure and simple pf.add_parameters(filenames=arr_file, par_type="grid", par_name_base=arr_file.split('.')[1] + "_gr", pargp=arr_file.split('.')[1] + "_gr", zone_array=ib, upper_bound=ub, lower_bound=lb, geostruct=gr_gs) - + # Add a variety of list style pars list_files = ["freyberg6.wel_stress_period_data_{0}.txt".format(t) for t in range(1, m.nper + 1)] # make dummy versions with headers - for fl in list_files[0:2]: + for fl in list_files[0:2]: # this adds a header to well file with open(os.path.join(template_ws, fl), 'r') as fr: lines = [line for line in fr] with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: @@ -1249,13 +1253,14 @@ def mf6_freyberg_direct_test(): for line in lines: fw.write(line) - fl = "freyberg6.wel_stress_period_data_3.txt" - with open(os.path.join(template_ws, fl), 'r') as fr: - lines = [line for line in fr] - with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: - fw.write("well k i j flux \n") - for i, line in enumerate(lines): - fw.write(f"well{i}" + line) + # fl = "freyberg6.wel_stress_period_data_3.txt" # Add extra string col_id + for fl in list_files[2:4]: + with open(os.path.join(template_ws, fl), 'r') as fr: + lines = [line for line in fr] + with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: + fw.write("well k i j flux \n") + for i, line in enumerate(lines): + fw.write(f"well{i}" + line) list_files.sort() @@ -1269,16 +1274,49 @@ def mf6_freyberg_direct_test(): # add temporally indep, but spatially correlated wel flux pars pf.add_parameters(filenames=list_file, par_type="grid", par_name_base="wel_grid_{0}".format(kper), pargp="wel_{0}".format(kper), index_cols=[0, 1, 2], use_cols=[3], - upper_bound=0.0, lower_bound=-1000, geostruct=gr_gs,par_style="direct", + upper_bound=0.0, lower_bound=-1000, geostruct=gr_gs, par_style="direct", transform="none") + # Adding dummy list pars with different file structures + list_file = "new_freyberg6.wel_stress_period_data_1.txt" # with a header + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nwell_mlt', + pargp='nwell_mult', index_cols=['k', 'i', 'j'], use_cols='flux', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, + transform="none") + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nwell_grid', + pargp='nwell', index_cols=['k', 'i', 'j'], use_cols='flux', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none") + # with skip instead + list_file = "new_freyberg6.wel_stress_period_data_2.txt" + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nwell_grid', + pargp='nwell', index_cols=[0, 1, 2], use_cols=3, + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none", mfile_skip=1) + + list_file = "new_freyberg6.wel_stress_period_data_3.txt" + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nwell_mlt', + pargp='nwell_mult', index_cols=['well', 'k', 'i', 'j'], use_cols='flux', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, + transform="none") + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nwell_grid', + pargp='nwell', index_cols=['well', 'k', 'i', 'j'], use_cols='flux', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none") + # with skip instead + list_file = "new_freyberg6.wel_stress_period_data_4.txt" + pf.add_parameters(filenames=list_file, par_type="grid", + par_name_base='nwell_grid', pargp='nwell', + index_cols=[0, 1, 2, 3], # or... {'well': 0, 'k': 1, 'i': 2, 'j': 3}, + use_cols=4, upper_bound=10, lower_bound=-10, + geostruct=gr_gs, par_style="direct", transform="none", + mfile_skip=1) list_file = "freyberg6.ghb_stress_period_data_1.txt" - pf.add_parameters(filenames=list_file, par_type="grid", par_name_base=["ghb_stage","ghb_cond"], + pf.add_parameters(filenames=list_file, par_type="constant", par_name_base=["ghb_stage","ghb_cond"], pargp=["ghb_stage","ghb_cond"], index_cols=[0, 1, 2], use_cols=[3,4], - upper_bound=[35,150], lower_bound=[32,50], geostruct=gr_gs, par_style="direct", + upper_bound=[35,150], lower_bound=[32,50], par_style="direct", transform="none") - # add model run command pf.mod_sys_cmds.append("mf6") print(pf.mult_files) @@ -1303,6 +1341,17 @@ def mf6_freyberg_direct_test(): except Exception as e: os.chdir(b_d) raise Exception(str(e)) + # TODO Some checks on resultant par files... + list_files = [f for f in os.listdir('.') + if f.startswith('new_') and f.endswith('txt')] + # check on that those dummy pars compare to the model versions. + for f in list_files: + n_df = pd.read_csv(f, sep="\s+") + o_df = pd.read_csv(f.strip('new_'), sep="\s+", header=None) + o_df.columns = ['k', 'i', 'j', 'flux'] + assert n_df.loc[:, o_df.columns].eq(o_df).all().all(), ( + "Something broke with alternative style model files" + ) os.chdir(b_d) num_reals = 100 @@ -1352,6 +1401,7 @@ def mf6_freyberg_direct_test(): if np.abs(arr.max() - rch_val) > 1.0e-6 or np.abs(arr.min() - rch_val) > 1.0e-6: raise Exception("recharge too diff") + def mf6_freyberg_varying_idomain(): import numpy as np import pandas as pd @@ -1531,16 +1581,269 @@ def xsec_test(): assert pst.phi < 1.0e-7 +def mf6_freyberg_short_direct_test(): + + import numpy as np + import pandas as pd + pd.set_option('display.max_rows', 500) + pd.set_option('display.max_columns', 500) + pd.set_option('display.width', 1000) + try: + import flopy + except: + return + + org_model_ws = os.path.join('..', 'examples', 'freyberg_mf6') + tmp_model_ws = "temp_pst_from_direct" + if os.path.exists(tmp_model_ws): + shutil.rmtree(tmp_model_ws) + os.mkdir(tmp_model_ws) + sim = flopy.mf6.MFSimulation.load(sim_ws=org_model_ws) + # sim.set_all_data_external() + sim.simulation_data.mfpath.set_sim_path(tmp_model_ws) + # sim.set_all_data_external() + m = sim.get_model("freyberg6") + sim.set_all_data_external() + sim.write_simulation() + + # to by pass the issues with flopy + # shutil.copytree(org_model_ws,tmp_model_ws) + # sim = flopy.mf6.MFSimulation.load(sim_ws=org_model_ws) + # m = sim.get_model("freyberg6") + + # SETUP pest stuff... + os_utils.run("{0} ".format("mf6"), cwd=tmp_model_ws) + + template_ws = "new_temp_direct" + sr = m.modelgrid + # set up PstFrom object + pf = PstFrom(original_d=tmp_model_ws, new_d=template_ws, + remove_existing=True, + longnames=False, spatial_reference=sr, + zero_based=False, start_datetime="1-1-2018") + # obs + # using tabular style model output + # (generated by pyemu.gw_utils.setup_hds_obs()) + # pf.add_observations('freyberg.hds.dat', insfile='freyberg.hds.dat.ins2', + # index_cols='obsnme', use_cols='obsval', prefix='hds') + + # Add stream flow observation + # df = pd.read_csv(os.path.join(tmp_model_ws, "sfr.csv"), index_col=0) + pf.add_observations("sfr.csv", insfile="sfr.csv.ins", index_cols="time", + use_cols=["GAGE_1","HEADWATER","TAILWATER"],ofile_sep=",") + # Setup geostruct for spatial pars + v = pyemu.geostats.ExpVario(contribution=1.0, a=1000) + gr_gs = pyemu.geostats.GeoStruct(variograms=v, transform="log") + rch_temporal_gs = pyemu.geostats.GeoStruct(variograms=pyemu.geostats.ExpVario(contribution=1.0, a=60)) + pf.extra_py_imports.append('flopy') + ib = m.dis.idomain[0].array + tags = {"npf_k_": [0.1, 10.], "npf_k33_": [.1, 10], "sto_ss": [.1, 10], "sto_sy": [.9, 1.1], + "rch_recharge": [.5, 1.5]} + dts = pd.to_datetime("1-1-2018") + pd.to_timedelta(np.cumsum(sim.tdis.perioddata.array["perlen"]), unit="d") + print(dts) + # setup from array style pars + for tag, bnd in tags.items(): + lb, ub = bnd[0], bnd[1] + arr_files = [f for f in os.listdir(tmp_model_ws) if tag in f and f.endswith(".txt")] + if "rch" in tag: + for arr_file in arr_files: + nmebase = arr_file.split('.')[1].replace('layer','').replace('_','').replace("npf",'').replace("sto",'').replace("recharge",'') + # indy direct grid pars for each array type file + recharge_files = ["recharge_1.txt","recharge_2.txt","recharge_3.txt"] + pf.add_parameters(filenames=arr_file, par_type="grid", par_name_base="rchg", + pargp="rch_gr", zone_array=ib, upper_bound=1.0e-3, lower_bound=1.0e-7, + par_style="direct") + # additional constant mults + kper = int(arr_file.split('.')[1].split('_')[-1]) - 1 + pf.add_parameters(filenames=arr_file, par_type="constant", + par_name_base=nmebase + "cn", + pargp="rch_const", zone_array=ib, upper_bound=ub, lower_bound=lb, + geostruct=rch_temporal_gs, + datetime=dts[kper]) + else: + for arr_file in arr_files: + nmebase = arr_file.split('.')[1].replace( + 'layer', '').replace('_','').replace("npf",'').replace( + "sto",'').replace("recharge",'') + # grid mults pure and simple + pf.add_parameters(filenames=arr_file, par_type="grid", par_name_base=nmebase, + pargp=arr_file.split('.')[1] + "_gr", zone_array=ib, upper_bound=ub, + lower_bound=lb, + geostruct=gr_gs) + + # Add a variety of list style pars + list_files = ["freyberg6.wel_stress_period_data_{0}.txt".format(t) + for t in range(1, m.nper + 1)] + # make dummy versions with headers + for fl in list_files[0:2]: # this adds a header to well file + with open(os.path.join(template_ws, fl), 'r') as fr: + lines = [line for line in fr] + with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: + fw.write("k i j flx \n") + for line in lines: + fw.write(line) + + # fl = "freyberg6.wel_stress_period_data_3.txt" # Add extra string col_id + for fl in list_files[2:4]: + with open(os.path.join(template_ws, fl), 'r') as fr: + lines = [line for line in fr] + with open(os.path.join(template_ws, f"new_{fl}"), 'w') as fw: + fw.write("well k i j flx \n") + for i, line in enumerate(lines): + fw.write(f"w{i}" + line) + + + list_files.sort() + for list_file in list_files: + kper = int(list_file.split(".")[1].split('_')[-1]) - 1 + #add spatially constant, but temporally correlated wel flux pars + pf.add_parameters(filenames=list_file, par_type="constant", par_name_base="wel{0}".format(kper), + pargp="twel_mlt_{0}".format(kper), index_cols=[0, 1, 2], use_cols=[3], + upper_bound=1.5, lower_bound=0.5, datetime=dts[kper], geostruct=rch_temporal_gs) + + # add temporally indep, but spatially correlated wel flux pars + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base="wel{0}".format(kper), + pargp="wel_{0}_direct".format(kper), index_cols=[0, 1, 2], use_cols=[3], + upper_bound=0.0, lower_bound=-1000, geostruct=gr_gs, par_style="direct", + transform="none") + # Adding dummy list pars with different file structures + list_file = "new_freyberg6.wel_stress_period_data_1.txt" # with a header + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nw', + pargp='nwell_mult', index_cols=['k', 'i', 'j'], use_cols='flx', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, + transform="none") + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nw', + pargp='nwell', index_cols=['k', 'i', 'j'], use_cols='flx', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none") + # with skip instead + list_file = "new_freyberg6.wel_stress_period_data_2.txt" + pf.add_parameters(filenames=list_file, par_type="grid", par_name_base='nw', + pargp='nwell', index_cols=[0, 1, 2], use_cols=3, + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none", mfile_skip=1) + + list_file = "new_freyberg6.wel_stress_period_data_3.txt" + pf.add_parameters(filenames=list_file, par_type="grid", + pargp='nwell_mult', index_cols=['well', 'k', 'i', 'j'], use_cols='flx', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, + transform="none") + pf.add_parameters(filenames=list_file, par_type="grid", + pargp='nwell', index_cols=['well', 'k', 'i', 'j'], use_cols='flx', + upper_bound=10, lower_bound=-10, geostruct=gr_gs, par_style="direct", + transform="none") + # with skip instead + list_file = "new_freyberg6.wel_stress_period_data_4.txt" + pf.add_parameters(filenames=list_file, par_type="grid", + par_name_base='nw', pargp='nwell', + index_cols=[0, 1, 2, 3], # or... {'well': 0, 'k': 1, 'i': 2, 'j': 3}, + use_cols=4, upper_bound=10, lower_bound=-10, + geostruct=gr_gs, par_style="direct", transform="none", + mfile_skip=1) + + list_file = "freyberg6.ghb_stress_period_data_1.txt" + pf.add_parameters(filenames=list_file, par_type="constant", par_name_base=["ghbst","ghbc"], + pargp=["ghb_stage","ghb_cond"], index_cols=[0, 1, 2], use_cols=[3,4], + upper_bound=[35,150], lower_bound=[32,50], par_style="direct", + transform="none") + + # add model run command + pf.mod_sys_cmds.append("mf6") + print(pf.mult_files) + print(pf.org_files) + + # build pest + pst = pf.build_pst('freyberg.pst') + cov = pf.build_prior(fmt="non") + cov.to_coo("prior.jcb") + pst.try_parse_name_metadata() + df = pd.read_csv(os.path.join(tmp_model_ws, "heads.csv"), index_col=0) + pf.add_observations("heads.csv", insfile="heads.csv.ins", index_cols="time", use_cols=list(df.columns.values), + prefix="hds", rebuild_pst=True) + + # test par mults are working + b_d = os.getcwd() + os.chdir(pf.new_d) + pst.write_input_files() + try: + pyemu.helpers.apply_list_and_array_pars( + arr_par_file="mult2model_info.csv", chunk_len=1) + except Exception as e: + os.chdir(b_d) + raise Exception(str(e)) + # TODO Some checks on resultant par files... + list_files = [f for f in os.listdir('.') + if f.startswith('new_') and f.endswith('txt')] + # check on that those dummy pars compare to the model versions. + for f in list_files: + n_df = pd.read_csv(f, sep="\s+") + o_df = pd.read_csv(f.strip('new_'), sep="\s+", header=None) + o_df.columns = ['k', 'i', 'j', 'flx'] + assert np.allclose(o_df.values, + n_df.loc[:, o_df.columns].values, + rtol=1e-4), ( + f"Something broke with alternative style model file: {f}" + ) + os.chdir(b_d) + + num_reals = 100 + pe = pf.draw(num_reals, use_specsim=True) + pe.to_binary(os.path.join(template_ws, "prior.jcb")) + assert pe.shape[1] == pst.npar_adj, "{0} vs {1}".format(pe.shape[0], pst.npar_adj) + assert pe.shape[0] == num_reals + + pst.control_data.noptmax = 0 + pst.pestpp_options["additional_ins_delimiters"] = "," + + pst.write(os.path.join(pf.new_d, "freyberg.pst")) + pyemu.os_utils.run("{0} freyberg.pst".format(ies_exe_path), cwd=pf.new_d) + + res_file = os.path.join(pf.new_d, "freyberg.base.rei") + assert os.path.exists(res_file), res_file + pst.set_res(res_file) + print(pst.phi) + assert pst.phi < 0.1, pst.phi + + + # turn direct recharge to min and direct wel to min and + # check that the model results are consistent + par = pst.parameter_data + rch_par = par.loc[par.pargp.apply( + lambda x: "rch_gr" in x and "direct" in x), "parnme"] + wel_par = par.loc[par.pargp.apply( + lambda x: "wel_" in x and "direct" in x), "parnme"] + par.loc[rch_par,"parval1"] = par.loc[rch_par, "parlbnd"] + # this should set wells to zero since they are negative values in the control file + par.loc[wel_par,"parval1"] = par.loc[wel_par, "parubnd"] + pst.write(os.path.join(pf.new_d, "freyberg.pst")) + pyemu.os_utils.run("{0} freyberg.pst".format(ies_exe_path), cwd=pf.new_d) + lst = flopy.utils.Mf6ListBudget(os.path.join(pf.new_d, "freyberg6.lst")) + flx, cum = lst.get_dataframes(diff=True) + wel_tot = flx.wel.apply(np.abs).sum() + print(flx.wel) + assert wel_tot < 1.0e-6, wel_tot + + # shortpars so not going to be able to get ij easily + # rch_files = [f for f in os.listdir(pf.new_d) + # if ".rch_recharge" in f and f.endswith(".txt")] + # rch_val = par.loc[rch_par,"parval1"][0] + # i, j = par.loc[rch_par, ["i", 'j']].astype(int).values.T + # for rch_file in rch_files: + # arr = np.loadtxt(os.path.join(pf.new_d, rch_file))[i, j] + # print(rch_file, rch_val, arr.mean(), arr.max(), arr.min()) + # if np.abs(arr.max() - rch_val) > 1.0e-6 or np.abs(arr.min() - rch_val) > 1.0e-6: + # raise Exception("recharge too diff") if __name__ == "__main__": - # freyberg_test() - # freyberg_prior_build_test() - #mf6_freyberg_test() - #mf6_freyberg_shortnames_test() - # mf6_freyberg_da_test() - #mf6_freyberg_direct_test() - #mf6_freyberg_varying_idomain() + freyberg_test() + freyberg_prior_build_test() + mf6_freyberg_test() + mf6_freyberg_shortnames_test() + mf6_freyberg_da_test() + mf6_freyberg_direct_test() + mf6_freyberg_short_direct_test() + mf6_freyberg_varying_idomain() xsec_test() diff --git a/pyemu/utils/pst_from.py b/pyemu/utils/pst_from.py index efe645560..dd1437dcb 100644 --- a/pyemu/utils/pst_from.py +++ b/pyemu/utils/pst_from.py @@ -17,7 +17,7 @@ # the tolerable percent difference (100 * (max - min)/mean) # used when checking that constant and zone type parameters are in fact constant (within # a given zone) -# DIRECT_PAR_PERCENT_DIFF_TOL = 1.0 +DIRECT_PAR_PERCENT_DIFF_TOL = 1.0 def _get_datetime_from_str(sdt): @@ -120,7 +120,8 @@ def __init__( self._prefix_count = {} self.get_xy = None - + self.add_pars_callcount = 0 + self.ijwarned = {} self.initialize_spatial_reference() self._setup_dirs() @@ -228,17 +229,29 @@ def parse_kij_args(self, args, kwargs): if ij_id is not None: i, j = [args[ij] for ij in ij_id] else: + if not self.ijwarned[self.add_pars_callcount]: + self.logger.warn( + ( + "get_xy() warning: position of i and j in index_cols " + "not specified, assume (i,j) are final two entries in " + "index_cols." + ) + + ) + self.ijwarned[self.add_pars_callcount] = True # assume i and j are the final two entries in index_cols i, j = args[-2], args[-1] else: - self.logger.warn( - ( - "get_xy() warning: need locational information " - "(e.g. i,j) to generate xy, " - "insufficient index cols passed to interpret: {}" - "" - ).format(str(args)) - ) + if not self.ijwarned[self.add_pars_callcount]: + self.logger.warn( + ( + "get_xy() warning: need locational information " + "(e.g. i,j) to generate xy, " + "insufficient index cols passed to interpret: {}" + "" + ).format(str(args)) + ) + self.ijwarned[self.add_pars_callcount] = True i, j = None, None return i, j @@ -1205,10 +1218,15 @@ def add_parameters( par_name_base (`str`): basename for parameters that are set up index_cols (`list`-like): if not None, will attempt to parameterize expecting a tabular-style model input file. `index_cols` - defines the unique columns used to set up pars. - May be passed as a dictionary using the keys `i` and `j` - to allow columns that relate to model rows and columns to be - identified and processed to x,y. + defines the unique columns used to set up pars. If passed as a + list of `str`, stings are expected to denote the columns + headers in tabular-style parameter files; if `i` and `j` in + list, these columns will be used to define spatial position for + spatial correlations (if required). WARNING: If passed as list + of `int`, `i` and `j` will be assumed to be in last two entries + in the list. Can be passed as a dictionary using the keys + `i` and `j` to explicitly speficy the columns that relate to + model rows and columns to be identified and processed to x,y. use_cols (`list`-like or `int`): for tabular-style model input file, defines the columns to be parameterised pargp (`str`): Parameter group to assign pars to. This is PESTs @@ -1262,7 +1280,8 @@ def add_parameters( # TODO support passing par_file (i,j)/(x,y) directly where information # is not contained in model parameter file - e.g. no i,j columns - + self.add_pars_callcount += 1 + self.ijwarned[self.add_pars_callcount] = False # this keeps denormal values for creeping into the model input arrays if ult_ubound is None: ult_ubound = self.ult_ubound_fill @@ -1302,9 +1321,9 @@ def add_parameters( ) ) if geostruct is not None: - if geostruct.sill != 1.0 and par_style != "multiplier": + if geostruct.sill != 1.0: # and par_style != "multiplier": #TODO !=? self.logger.warn( - "geostruct sill != 1.0 for 'multiplier' style parameters" + "geostruct sill != 1.0" # for 'multiplier' style parameters" ) if geostruct.transform != transform: self.logger.warn( @@ -2303,6 +2322,7 @@ def _write_direct_df_tpl( pandas.DataFrame with paranme and pargroup define for each `use_col` """ + # TODO much of this duplicates what is in _get_tpl_or_ins_df() -- could posssibly be consolidated # work out the union of indices across all dfs sidx = set() @@ -2320,7 +2340,7 @@ def _write_direct_df_tpl( else: inames = ["idx{0}".format(i) for i in range(len(index_cols))] else: - fmt = "{1|3}" + fmt = "{1:3}" j = "" if isinstance(index_cols[0], str): inames = index_cols @@ -2408,7 +2428,7 @@ def _write_direct_df_tpl( df_ti.loc[:, use_col] += suffix # todo check that values are constant within zones and - # assign parval1 + # assign parval1 raise NotImplementedError( "list-based direct zone-type parameters not implemented" ) @@ -2837,13 +2857,12 @@ def grid_namer(i, j): def _check_diff(org_arr, input_filename, zval=None): - percent_diff = 100.0 * np.abs( - np.nanmax(org_arr) - np.nanmin(org_arr) / np.nanmean(org_arr) - ) + percent_diff = 100.0 * np.abs((np.nanmax(org_arr) - np.nanmin(org_arr)) + / np.nanmean(org_arr)) if percent_diff > DIRECT_PAR_PERCENT_DIFF_TOL: message = "_check_diff() error: direct par for file '{0}'".format( input_filename - ) + "exceeds tolerance for percent difference: {2} > {3}".format( + ) + "exceeds tolerance for percent difference: {0} > {1}".format( percent_diff, DIRECT_PAR_PERCENT_DIFF_TOL ) if zval is not None: