diff --git a/UQPyL/__init__.py b/UQPyL/__init__.py index 1c43282e..7e5b6b42 100644 --- a/UQPyL/__init__.py +++ b/UQPyL/__init__.py @@ -1,6 +1,6 @@ from . import problems, surrogates, optimization, sensibility, DoE, utility -__version__ = "2.0.4" +__version__ = "2.0.5" __author__ = "wmtSky" __all__=[ diff --git a/UQPyL/optimization/asmo.py b/UQPyL/optimization/asmo.py index 78906d38..fdf96f0c 100644 --- a/UQPyL/optimization/asmo.py +++ b/UQPyL/optimization/asmo.py @@ -113,7 +113,7 @@ def run(self, oneStep=False): history_best_decs[FE]=BestX history_best_objs[FE]=BestY - Result={'best_dec':BestX, 'best_obj':BestY, 'history_best_decs': history_best_decs, 'history_best_objs':history_best_objs ,'FE':FE} + Result={'best_dec':BestX, 'best_obj':BestY, 'history_best_decs': history_best_decs, 'history_best_objs':history_best_objs ,'FEs':FE} return Result diff --git a/UQPyL/sensibility/morris.py b/UQPyL/sensibility/morris.py index efa1fdd0..462dd729 100644 --- a/UQPyL/sensibility/morris.py +++ b/UQPyL/sensibility/morris.py @@ -94,7 +94,7 @@ def analyze(self, X: Optional[np.ndarray]=None, Y: Optional[np.ndarray]=None, ve N=int(X.shape[0]/self.num_levels) - for i in range(N): + for i in range(num_trajectory): X_sub=X[i*(n_input+1):(i+1)*(n_input+1), :] Y_sub=Y[i*(n_input+1):(i+1)*(n_input+1), :] diff --git a/examples/draw_p.py b/examples/draw_p.py new file mode 100644 index 00000000..9be672d8 --- /dev/null +++ b/examples/draw_p.py @@ -0,0 +1,33 @@ +import os +os.chdir('./examples') + +import numpy as np +import matplotlib.pyplot as plt + +values=np.loadtxt('./sensibility.txt', dtype=np.float32) +with open('./si_names.txt', 'r') as f: + name=f.readlines() + name=[n.strip() for n in name] +colors=['#03AED2', '#FFC55A', '#75A47F', '#D2649A'] +labels=['Sobol\'', 'FAST', 'Morris', 'MARS'] +x=np.arange(len(name)) +bar_width=0.2 +for i in range(4): + plt.bar(x + i * bar_width, values[:, i], bar_width, label=labels[i], color=colors[i]) +for i in range(len(name) - 1): + plt.axvline(x=i+1-bar_width, color='black', linestyle='--', linewidth=0.75) + + +plt.xlim(-0.5, 25.8) +plt.xticks(x+0.3, [f'P{i+1}' for i in range(len(name))], fontsize=14) +plt.yticks(fontsize=14) +plt.tick_params(axis='both', which='both', length=0) +plt.xlabel('Parameters', fontsize=14) +plt.ylabel('Sensibility', fontsize=14) +plt.legend(loc='upper right', fancybox=True, shadow=True, ncol=4, fontsize=14) +ax = plt.gca() +ax.spines['top'].set_linewidth(2) +ax.spines['bottom'].set_linewidth(2) +ax.spines['left'].set_linewidth(2) +ax.spines['right'].set_linewidth(2) +plt.show() \ No newline at end of file diff --git a/examples/example_optimization.py b/examples/example_optimization.py index eea6fc8e..b7e6b5df 100644 --- a/examples/example_optimization.py +++ b/examples/example_optimization.py @@ -8,35 +8,32 @@ - 6. MOEAD/D - Multi-objective - 7. MO_ASMO - Multi-objective - Surrogate ''' -#tmp +#--------------tmp-------------------# import sys sys.path.append(".") -from scipy.io import loadmat -print(sys.path) import os os.chdir('./examples') -# +#-----------tmp--------------------# import numpy as np ################1. Genetic Algorithm (GA) - Single objective################ -# print('################1. Genetic Algorithm (GA) - Single objective################') -# from UQPyL.optimization import GA -# from UQPyL.problems import Sphere +print('################1. Genetic Algorithm (GA) - Single objective################') +from UQPyL.optimization import GA +from UQPyL.problems import Sphere + +problem=Sphere(n_input=30, ub=100, lb=-100) +ga=GA(problem, n_samples=50) +res=ga.run() +print('Best objective:', res['best_obj']) +print('Best decisions:', res['best_dec']) +print('FE:', res['FEs']) -# problem=Sphere(n_input=30, ub=100, lb=-100) -# ga=GA(problem, n_samples=50) -# res=ga.run() -# print('Best objective:', res['best_obj']) -# print('Best decisions:', res['best_dec']) -# print('FE:', res['FEs']) - # import matplotlib.pyplot as plt # FEs_objs=res['history_best_objs'] # FEs=list(FEs_objs.keys()) # objs=list(FEs_objs.values()) # plt.plot(FEs, np.log10(objs)) - ################2. SCE-UA - Single objective################ # print('################2. SCE-UA - Single objective################') # from UQPyL.optimization import SCE_UA @@ -83,7 +80,7 @@ problem=Sphere(n_input=30, ub=100, lb=-100) rbf=RBF(kernel=Cubic()) asmo=ASMO(problem, rbf, n_init=50) -res=asmo.run(maxFE=1000) +res=asmo.run() print('Best objective:', res['best_obj']) print('Best decisions:', res['best_dec']) print('FE:', res['FEs']) diff --git a/examples/parallel_swat1(1).py b/examples/parallel_swat1(1).py new file mode 100644 index 00000000..74b6ce01 --- /dev/null +++ b/examples/parallel_swat1(1).py @@ -0,0 +1,380 @@ +import re +import os +import itertools +import pandas as pd +import numpy as np +from datetime import datetime, timedelta +import tempfile +import shutil +from UQPyL.DoE import LHS +from UQPyL.problems import Problem +from UQPyL.utility.metrics import r2_score +from datetime import datetime, timedelta +import subprocess +import multiprocessing + +def _get_default_paras(file_path, name): + pattern=r"(\s*)(\-?\d+\.?\d*)(\s*.*{}.*)".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for line in lines: + match = re.search(pattern, line) + if match: + return float(match.group(2)) + +def _get_default_paras_for_sol(file_path, name): + pattern=r".*{}.*".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for line in lines: + match = re.search(pattern, line) + if match: + numbers=re.findall(r"\d+\.\d+", line) + return " ".join(numbers) + +def _set_paras_for_sol(file_path, name, value, mode, origin_value=None): + pattern=r"\s*{}.*".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(pattern, line) + if match: + numbers=re.findall(r"\d+\.\d+", line) + for j, number in enumerate(numbers): + if mode=='r': + value=float(origin_value[j])*(1+value) + elif mode=='a': + value=float(origin_value[j])+value + line=line.replace(str(number), "{:.2f}".format(value)) + lines[i]=line + break + f.seek(0) + f.writelines(lines) + f.truncate() + +def _set_paras(file_path, name, value, mode, origin_value=None): + + pattern=r"(\s*)(\-?\d+\.?\d*)(\s*.*{}.*)".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(pattern, line) + + if match: + if mode=='r': + value=origin_value*(1+value) + elif mode=='a': + value=origin_value+value + + new_text = re.sub(pattern, match.group(1)+"{}".format(value)+match.group(3), line) + lines[i]=new_text + break + f.seek(0) + f.writelines(lines) + f.truncate() +############################### + + +def sub_process(x, swat_cups_indices, lock, swat_cups, log_queue): + result=-1 + try: + with lock: + if not swat_cups_indices: + log_queue.put("No more swat_cup") + return None + + index=swat_cups_indices.pop() + + swat_cup=swat_cups[index] + result=swat_cup.evaluate(x.reshape(1,-1)) + finally: + with lock: + swat_cups_indices.append(index) + log_queue.put("swat_cup {} finished".format(index)) + + return result + +def main(swat_cups, X): + + + with multiprocessing.Manager() as manager: + lock=manager.Lock() + swat_cups_indices=manager.list(range(len(swat_cups))) + log_queue=manager.Queue() + + try: + with multiprocessing.Pool(processes=12) as pool: + results=pool.starmap(sub_process, [(row, swat_cups_indices, lock, swat_cups, log_queue) for row in X]) + + finally: + pass + + while not log_queue.empty(): + print(log_queue.get()) + + return results + +class SWAT_CUP(Problem): + + HRU_suffix=["chm", "gw", "hru", "mgt", "sdr", "sep", "sol"] + Watershed_suffiex=["pnd", "rte", "sub", "swq", "wgn", "wus"] + n_output=1 + + def __init__(self, work_path, para_file_name, observe_file_name, swat_exe_path, rch_id): + self.work_path=work_path + self.para_file_name=para_file_name + self.observe_file_name=observe_file_name + self.exe_path=os.path.join(self.work_path, swat_exe_path) + self.rch_id=rch_id + + self._initial() + self._recond_default_values() #assign self.lb self.ub self.n_input + self._get_observed_data() + + def evaluate(self, X): + n=X.shape[0] + Y=np.zeros((n,1)) + for i in range(n): + self._set_values(X[i]) + try: + subprocess.run(self.exe_path, cwd=self.work_path,stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) + simulation_data=self._get_simulation_data() + rch_flows=simulation_data.query('RCH=={}'.format(self.rch_id))['FLOW_OUTcms'].to_numpy() + Y[i,0]=-r2_score(self.observed_data, rch_flows[self.begin_calibration-1:self.end_calibration]) + except: + Y[i,0]=-np.inf + return Y + + def _get_simulation_data(self): + file_name="output.rch" + file_path=os.path.join(self.work_path, file_name) + begin_ind=0 + with open(file_path, "r") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(r".*FLOW_OUTcms.*", line) + if match: + begin_ind=i + break + colspecs = [ + (7, 11), # RCH + (22, 26), # MON + (52, 62), # FLOW_OUTcms + ] + + simulation_data = pd.read_fwf(file_path, colspecs=colspecs, header=None, names=['RCH', 'MON', 'FLOW_OUTcms'], skiprows=begin_ind+1) + + return simulation_data + + def _get_observed_data(self): + file_path=os.path.join(self.work_path, self.observe_file_name) + data=[] + with open(file_path, "r") as f: + lines=f.readlines() + pattern = re.compile(r'(\d+)\s+FLOW_OUT_(\d+)_(\d+)\s+(\d+\.?\d*)') + for line in lines: + match = pattern.match(line) + if match: + index, day, year = map(int, match.groups()[:-1]) + value = float(match.groups()[-1]) + data.append([index, day, year, value]) + observed_data = pd.DataFrame(data, columns=['index', 'day', 'year', 'value']) + observed_data = observed_data.astype({'index': int, 'day': int, 'year': int, 'value': float}).set_index('index') + values_array = observed_data['value'].to_numpy(dtype=float) + + first_record = observed_data.iloc[0] + self.begin_calibration=(datetime(int(first_record['year']), 1, 1) + timedelta(int(first_record['day']) - 1)-self.begin_date).days + + last_record = observed_data.iloc[-1] + self.end_calibration=(datetime(int(last_record['year']), 1, 1) + timedelta(int(last_record['day']) - 1)-self.begin_date).days + + self.observed_data=values_array + + def _set_values(self, paras_values): + + paras_values=paras_values.ravel() + for i in range(self.n_input): + para_name=self.paras_list[i] + mode=self.paras_infos.loc[para_name, 'mode'] + value=paras_values[i] + file_suffix=self.paras_files.loc[para_name, 'file'] + + if file_suffix in self.HRU_suffix: + for sub in self.total_sub_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + if file_suffix=='sol': + _set_paras_for_sol(file_path, para_name, value, mode, str.split(self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0])) + else: + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + elif file_suffix in self.Watershed_suffiex: + for sub in self.watershed_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, sub+"."+file_suffix) + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + elif file_suffix=="bsn": + file_name="basins."+file_suffix + file_path=os.path.join(self.work_path, file_name) + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + + def _recond_default_values(self): + + para_file_name=os.path.join(self.work_path, self.para_file_name) + paras_infos=pd.read_csv(para_file_name, sep=' ', names=['Parameter', 'mode', 'low_bound', 'up_bound'], index_col='Parameter') + self.lb= paras_infos['low_bound'].values + self.ub= paras_infos['up_bound'].values + paras_list=paras_infos.index.tolist() + self.x_labels=paras_list + num_paras=paras_infos.shape[0] + self.n_input=num_paras + ##generate default_values + default_values_path=os.path.join(self.work_path, 'default_values.xlsx') + if not os.path.exists(default_values_path): + default_values=pd.DataFrame(columns=['para_name', 'file_name', 'value']) + + for i in range(num_paras): + para_name=paras_list[i] + file_suffix=self.paras_files.loc[para_name, 'file'] + + if file_suffix in self.HRU_suffix: + for sub in self.total_sub_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + if file_suffix=='sol': + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras_for_sol(file_path, para_name)] + else: + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] + elif file_suffix in self.Watershed_suffiex: + for sub in self.watershed_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] + elif file_suffix=="bsn": + file_path=os.path.join(self.work_path, "basins."+file_suffix) + default_values.loc[len(default_values)]=[para_name, "basins."+file_suffix, _get_default_paras(file_path, para_name)] + + + default_values.to_excel(default_values_path, index=True) + else: + default_values=pd.read_excel(default_values_path, index_col=0) + self.default_values=default_values + self.paras_infos=paras_infos + self.paras_list=paras_list + def _initial(self): + Watershed={} + + NYSKIP=-1 + num_years=-1 + begin_year=-1 + begin_day=-1 + end_day=-1 + #read control file fig.cio + with open(os.path.join(self.work_path, "file.cio"), "r") as f: + lines=f.readlines() + for line in lines: + match1 = re.search(r"(\s*)(\d+)(\s*.*NBYR.*)", line) + if num_years<0 and match1: + num_years=int(match1.group(2)) + + match2 = re.search(r"(\s*)(\d+)(\s*.*IYR.*)", line) + if begin_year<0 and match2: + begin_year=int(match2.group(2)) + + match3 = re.search(r"(\s*)(\d+)(\s*.*IDAF.*)", line) + if begin_day<0 and match3: + begin_day=int(match3.group(2)) + + match4 = re.search(r"(\s*)(\d+)(\s*.*IDAL.*)", line) + if end_day<0 and match4: + end_day=int(match4.group(2)) + + match5 = re.search(r"(\s*)(\d+)(\s*.*NYSKIP.*)", line) + if NYSKIP<0 and match5: + NYSKIP=int(match5.group(2)) + + self.begin_date=datetime(begin_year+NYSKIP, 1, 1) + timedelta(begin_day - 1) + self.end_date=datetime(begin_year+num_years-1, 1, 1) + timedelta(end_day - 1) + self.simulation_days=(self.end_date-self.begin_date).days+1 + + #read control file fig.fig + with open(os.path.join(self.work_path, "fig.fig"), "r") as f: + lines=f.readlines() + for line in lines: + match = re.search(r'(\d+)\.sub', line) + if match: + Watershed[match.group(1)]=[] + + # #read sub files + # def read_sub_files(file_path, id, queue): + # file_name=id+".sub" + # res=[] + # with open(os.path.join(file_path, file_name), "r") as f: + # lines=f.readlines() + # for line in lines: + # match = re.search(r'(\d+)\.mgt', line) + # if match: + # res.append(match.group(1)) + # queue.put((id, res)) + # result_queue=queue.Queue() + # threads=[] + # for sub in Watershed: + # thread=threading.Thread(target=read_sub_files, args=(self.work_path, sub, result_queue)) + # thread.start() + + for sub in Watershed: + file_name=sub+".sub" + with open(os.path.join(self.work_path, file_name), "r") as f: + lines=f.readlines() + for line in lines: + match = re.search(r'(\d+)\.mgt', line) + if match: + Watershed[sub].append(match.group(1)) + + self.Watershed=Watershed + self.total_sub_list = list(itertools.chain.from_iterable(Watershed.values())) + self.watershed_list=list(Watershed.keys()) + #read para_file + file_path=os.path.join(self.work_path, 'SWAT_paras_files.txt') + self.paras_files = pd.read_csv(file_path, sep=' ', names=['parameter', 'file'], index_col='parameter') + +if __name__=="__main__": + + from UQPyL.sensibility import Sobol, Morris + from UQPyL.DoE import LHS + work_temp_dir=tempfile.mkdtemp() + + n_parallels=12 + + work_temp_dirs=[os.path.join(work_temp_dir, "instance{}".format(i)) for i in range(n_parallels)] + + for work_temp in work_temp_dirs: + shutil.copytree("D:\SiHuRiver\model\FuTIanSi001\Scenarios\Test\TxtInOut", work_temp) + + swat_cups=[SWAT_CUP(work_path=work_temp_dir, + para_file_name="paras_infos.txt", + observe_file_name="observed.txt", + swat_exe_path="swat_64rel.exe", + rch_id=40) for work_temp_dir in work_temp_dirs] + for i in range(10): + lhs=LHS(criterion='classic') + num=1000 + X=lhs.sample(nt=num, nx=swat_cups[1].n_input) + # sobol=Sobol(problem=swat_cups[0], cal_second_order=False) + # X=sobol.sample(128) + # mor=Morris(problem=swat_cups[0], num_levels=4) + # X=mor.sample(num_trajectory=128) + X=X*(swat_cups[1].ub-swat_cups[1].lb)+swat_cups[1].lb + res=np.array(main(swat_cups, X)).ravel().reshape(-1,1) + np.savetxt("res_Y_{}_{}.txt".format(num, i), res) + np.savetxt("res_X_{}_{}.txt".format(num, i), X) + for temp_dir in work_temp_dirs: + shutil.rmtree(temp_dir) + + + + + + + + + diff --git a/examples/sensibility.txt b/examples/sensibility.txt new file mode 100644 index 00000000..336385d1 --- /dev/null +++ b/examples/sensibility.txt @@ -0,0 +1,26 @@ +0.200532725 0.576757356 0.259118928 0.7882 +0 0.003816324 0 0 +0.031582953 0.004924289 0.011140644 0 +0.00304414 0.004062538 0 0 +0.03348554 0.005293611 0.063790081 0 +0 0.004924289 0 0 +0.03348554 0.004062538 0.032686761 0 +0.004 0.006 0.03 0.005 +0.006 0.004 0.01 0 +0.027 0.004 0.05 0 +0.027777778 0.036685953 0.24175762 0.0449 +0.157914764 0.023390373 0.000169654 0.0553 +0.255707763 0.081743198 0.097890629 0.0311 +0.027777778 0.004308753 0.000622066 0 +0.002663623 0.004801182 0.01176271 0 +0.0304414 0.004554967 0.011310298 0 +0.00456621 0.004554967 0.011423401 0 +0.022450533 0.026714268 0.001244133 0.0165 +0.01369863 0.165209898 0.119436747 0.0568 +0.023592085 0.004062538 0.027484024 0.0013 +0.003805175 0.004 0.001244133 0 +0.027016743 0.004801182 0.001074478 0 +0.003805175 0.004801182 0.000508963 0.0014 +0.027777778 0.003939431 0.000622066 0 +0.003805175 0.004185646 0.010970989 0 +0.027777778 0.004185646 0.000452412 0 \ No newline at end of file diff --git a/examples/si_names.txt b/examples/si_names.txt new file mode 100644 index 00000000..05c5a9ea --- /dev/null +++ b/examples/si_names.txt @@ -0,0 +1,26 @@ +CN2 +GW_DELAY +ALPHA_BF +GWQMN +GW_REVAP +REVAPMN +RCHRG_DP +SOL_AWC +SOL_K +SOL_ALB +CH_N2 +CH_K2 +ALPHA_BNK +TLAPS +SLSUBBSN +HRU_SLP +OV_N +CANMX +ESCO +EPCO +SFTMP +SMTMP +SMFMX +SMFMN +TIMP +SURLAG \ No newline at end of file diff --git a/examples/temp.py b/examples/temp.py new file mode 100644 index 00000000..4ceb91e4 --- /dev/null +++ b/examples/temp.py @@ -0,0 +1,528 @@ +import re +import os +import itertools +import pandas as pd +import numpy as np +import sys +from datetime import datetime, timedelta +sys.path.append(".") +from scipy.io import loadmat +os.chdir('./examples') +from UQPyL.utility.metrics import r_square + +def _get_default_paras(file_path, name): + pattern=r"(\s*)(\-?\d+\.?\d*)(\s*.*{}.*)".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for line in lines: + match = re.search(pattern, line) + if match: + return float(match.group(2)) + +def _get_default_paras_for_sol(file_path, name): + pattern=r".*{}.*".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for line in lines: + match = re.search(pattern, line) + if match: + numbers=re.findall(r"\d+\.\d+", line) + return " ".join(numbers) + +def _set_paras_for_sol(file_path, name, value, mode, origin_value=None): + pattern=r"\s*{}.*".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(pattern, line) + if match: + numbers=re.findall(r"\d+\.\d+", line) + for j, number in enumerate(numbers): + if mode=='r': + value=float(origin_value[j])*(1+value) + elif mode=='a': + value=float(origin_value[j])+value + line=line.replace(str(number), "{:.2f}".format(value)) + lines[i]=line + break + f.seek(0) + f.writelines(lines) + f.truncate() + +def _set_paras(file_path, name, value, mode, origin_value=None): + + pattern=r"(\s*)(\-?\d+\.?\d*)(\s*.*{}.*)".format(name) + with open(file_path, "r+") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(pattern, line) + + if match: + if mode=='r': + value=origin_value*(1+value) + elif mode=='a': + value=origin_value+value + + new_text = re.sub(pattern, match.group(1)+"{}".format(value)+match.group(3), line) + lines[i]=new_text + # origin_value=match.group(2) + # if re.fullmatch(r'\-?\d+', origin_value) is not None and re.fullmatch(r'\-?\d+', str(value)) is not None: + # if mode=='r': + # value=origin_value*(1+value) + # new_text = re.sub(pattern, match.group(1)+"{}".format(value)+match.group(3), line) + # lines[i]=new_text + # else: + # raise ValueError("Please check the type of input value for {} in {}.".format(name, file_path)) + break + f.seek(0) + f.writelines(lines) + f.truncate() +# #user setting +# # work_path="D:\SiHuRiver\model\FuTIanSi001\Scenarios\Test\TxtInOut" +# # para_file_name="paras_infos.txt" + +# # begin_year=2009 +# # begin_day=1 +# # end_year=2016 +# # end_day=366 + +# # begin_data=datetime(begin_year, 1, 1) + timedelta(begin_day - 1) +# # end_data=datetime(end_year, 1, 1) + timedelta(end_day - 1) +# # days_between=(end_data-begin_data).days+1 +# # begin_days_calibration=(datetime(2011,1,1)-begin_data).days-1 +# # end_days_calibration=(datetime(2014,12,31)-begin_data).days + +# #initial +# # paras_files = pd.read_csv('SWAT_paras_files.txt', sep=' ', names=['parameter', 'file'], index_col='parameter') +# # HRU_suffix=["chm", "gw", "hru", "mgt", "sdr", "sep", "sol"] +# # Watershed_suffiex=["pnd", "rte", "sub", "swq", "wgn", "wus"] + +# # Watershed={} +# # with open(os.path.join(work_path, "fig.fig"), "r") as f: +# # lines=f.readlines() +# # for line in lines: +# # match = re.search(r'(\d+)\.sub', line) +# # if match: +# # Watershed[match.group(1)]=[] + +# # for sub in Watershed: +# # file_name=sub+".sub" +# # with open(os.path.join(work_path, file_name), "r") as f: +# # lines=f.readlines() +# # for line in lines: +# # match = re.search(r'(\d+)\.mgt', line) +# # if match: +# # Watershed[sub].append(match.group(1)) + +# # total_sub_list = list(itertools.chain.from_iterable(Watershed.values())) +# # watershed_list=list(Watershed.keys()) + +# # paras_infos=pd.read_csv(para_file_name, sep=' ', names=['Parameter', 'mode', 'low_bound', 'up_bound'], index_col='Parameter') +# # low_bound= paras_infos['low_bound'].values +# # up_bound= paras_infos['up_bound'].values +# # paras_list=paras_infos.index.tolist() +# # num_paras=paras_infos.shape[0] +# # ##generate default_values +# # default_values=pd.DataFrame(columns=['para_name', 'file_name', 'value']) + +# # for i in range(num_paras): +# # para_name=paras_list[i] +# # file_suffix=paras_files.loc[para_name, 'file'] + +# # if file_suffix in HRU_suffix: +# # for sub in total_sub_list: +# # file_name=sub+"."+file_suffix +# # file_path=os.path.join(work_path, file_name) +# # if file_suffix=='sol': +# # default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras_for_sol(file_path, para_name)] +# # else: +# # default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] +# # elif file_suffix in Watershed_suffiex: +# # for sub in watershed_list: +# # file_name=sub+"."+file_suffix +# # file_path=os.path.join(work_path, file_name) +# # default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] +# # elif file_suffix=="bsn": +# # file_path=os.path.join(work_path, "basins."+file_suffix) +# # default_values.loc[len(default_values)]=[para_name, "basins."+file_suffix, _get_default_paras(file_path, para_name)] + +# # default_values.to_excel('default_values.xlsx', index=True) + +# #setting values +# # paras_values=np.random.random((1, num_paras))*(up_bound-low_bound)+low_bound +# # paras_values=paras_values.ravel() +# # for i in range(num_paras): +# # para_name=paras_list[i] +# # mode=paras_infos.loc[para_name, 'mode'] +# # value=paras_values[i] +# # file_suffix=paras_files.loc[para_name, 'file'] + +# # if file_suffix in HRU_suffix: +# # for sub in total_sub_list: +# # file_name=sub+"."+file_suffix +# # file_path=os.path.join(work_path, file_name) +# # if file_suffix=='sol': +# # _set_paras_for_sol(file_path, para_name, value, mode, str.split(default_values.loc[(default_values['para_name'] == para_name) & (default_values['file_name'] == file_name), 'value'].values[0])) +# # else: +# # _set_paras(file_path, para_name, value, mode, default_values.loc[(default_values['para_name'] == para_name) & (default_values['file_name'] == file_name), 'value'].values[0]) +# # elif file_suffix in Watershed_suffiex: +# # for sub in watershed_list: +# # file_name=sub+"."+file_suffix +# # file_path=os.path.join(work_path, sub+"."+file_suffix) +# # _set_paras(file_path, para_name, value, mode, default_values.loc[(default_values['para_name'] == para_name) & (default_values['file_name'] == file_name), 'value'].values[0]) +# # elif file_suffix=="bsn": +# # file_name="basins."+file_suffix +# # file_path=os.path.join(work_path, file_name) +# # _set_paras(file_path, para_name, value, mode, default_values.loc[(default_values['para_name'] == para_name) & (default_values['file_name'] == file_name), 'value'].values[0]) +# #observe +# file_name="observed.txt" +# file_path=os.path.join(work_path, file_name) + +# data=[] +# with open(file_path, "r") as f: +# lines=f.readlines() +# pattern = re.compile(r'(\d+)\s+FLOW_OUT_(\d+)_(\d+)\s+(\d+)') +# for line in lines: +# match = pattern.match(line) +# if match: +# index, day, year, value = map(int, match.groups()) +# data.append([index, day, year, value]) +# observed_data = pd.DataFrame(data, columns=['index', 'day', 'year', 'value']).set_index('index') +# values_array = observed_data['value'].to_numpy() + +# ############################evaluation################## + +# def evaluate(): + +# file_name="output.rch" +# file_path=os.path.join(work_path, file_name) +# begin_ind=0 +# with open(file_path, "r") as f: +# lines=f.readlines() +# for i, line in enumerate(lines): +# match = re.search(r".*FLOW_OUTcms.*", line) +# if match: +# begin_ind=i +# break +# colspecs = [ +# (7, 11), # RCH +# (22, 26), # MON +# (52, 62), # FLOW_OUTcms +# ] +# # name_list=['RCH', 'MON', 'FLOW_OUTcms'] + +# df = pd.read_fwf(file_path, colspecs=colspecs, header=None, names=['RCH', 'MON', 'FLOW_OUTcms'], skiprows=begin_ind+1) +# df_rch_40 = df.query('RCH==40') +# flow_outcms_np = df_rch_40['FLOW_OUTcms'].to_numpy() + +# r2=r_square(values_array[begin_days_calibration:end_days_calibration], flow_outcms_np[begin_days_calibration:end_days_calibration]) + +# return r2 + + +############################### +from UQPyL.problems import Problem +from UQPyL.utility.metrics import r_square +from datetime import datetime, timedelta +import subprocess +class SWAT_CUP(Problem): + HRU_suffix=["chm", "gw", "hru", "mgt", "sdr", "sep", "sol"] + Watershed_suffiex=["pnd", "rte", "sub", "swq", "wgn", "wus"] + n_output=1 + def __init__(self, work_path, para_file_name, observe_file_name, swat_exe_path, rch_id): + self.work_path=work_path + self.para_file_name=para_file_name + self.observe_file_name=observe_file_name + self.exe_path=os.path.join(self.work_path, swat_exe_path) + self.rch_id=rch_id + + self._initial() + self._recond_default_values() #assign self.lb self.ub self.n_input + self._get_observed_data() + + def evaluate(self, X): + n=X.shape[0] + Y=np.zeros((n,1)) + for i in range(n): + self._set_values(X[i]) + try: + subprocess.run(self.exe_path, cwd=self.work_path,stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) + simulation_data=self._get_simulation_data() + rch_flows=simulation_data.query('RCH=={}'.format(self.rch_id))['FLOW_OUTcms'].to_numpy() + Y[i,0]=-r_square(self.observed_data, rch_flows[self.begin_calibration-1:self.end_calibration]) + except: + Y[i,0]=-np.inf + return Y + + def _get_simulation_data(self): + file_name="output.rch" + file_path=os.path.join(self.work_path, file_name) + begin_ind=0 + with open(file_path, "r") as f: + lines=f.readlines() + for i, line in enumerate(lines): + match = re.search(r".*FLOW_OUTcms.*", line) + if match: + begin_ind=i + break + colspecs = [ + (7, 11), # RCH + (22, 26), # MON + (52, 62), # FLOW_OUTcms + ] + # name_list=['RCH', 'MON', 'FLOW_OUTcms'] + + simulation_data = pd.read_fwf(file_path, colspecs=colspecs, header=None, names=['RCH', 'MON', 'FLOW_OUTcms'], skiprows=begin_ind+1) + + return simulation_data + + def _get_observed_data(self): + file_path=os.path.join(self.work_path, self.observe_file_name) + data=[] + with open(file_path, "r") as f: + lines=f.readlines() + pattern = re.compile(r'(\d+)\s+FLOW_OUT_(\d+)_(\d+)\s+(\d+\.?\d*)') + for line in lines: + match = pattern.match(line) + if match: + index, day, year = map(int, match.groups()[:-1]) + value = float(match.groups()[-1]) + data.append([index, day, year, value]) + observed_data = pd.DataFrame(data, columns=['index', 'day', 'year', 'value']) + observed_data = observed_data.astype({'index': int, 'day': int, 'year': int, 'value': float}).set_index('index') + values_array = observed_data['value'].to_numpy(dtype=float) + + first_record = observed_data.iloc[0] + self.begin_calibration=(datetime(int(first_record['year']), 1, 1) + timedelta(int(first_record['day']) - 1)-self.begin_date).days + + last_record = observed_data.iloc[-1] + self.end_calibration=(datetime(int(last_record['year']), 1, 1) + timedelta(int(last_record['day']) - 1)-self.begin_date).days + + self.observed_data=values_array + + def _set_values(self, paras_values): + + paras_values=paras_values.ravel() + for i in range(self.n_input): + para_name=self.paras_list[i] + mode=self.paras_infos.loc[para_name, 'mode'] + value=paras_values[i] + file_suffix=self.paras_files.loc[para_name, 'file'] + + if file_suffix in self.HRU_suffix: + for sub in self.total_sub_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + if file_suffix=='sol': + _set_paras_for_sol(file_path, para_name, value, mode, str.split(self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0])) + else: + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + elif file_suffix in self.Watershed_suffiex: + for sub in self.watershed_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, sub+"."+file_suffix) + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + elif file_suffix=="bsn": + file_name="basins."+file_suffix + file_path=os.path.join(self.work_path, file_name) + _set_paras(file_path, para_name, value, mode, self.default_values.loc[(self.default_values['para_name'] == para_name) & (self.default_values['file_name'] == file_name), 'value'].values[0]) + + def _recond_default_values(self): + + para_file_name=os.path.join(self.work_path, self.para_file_name) + paras_infos=pd.read_csv(para_file_name, sep=' ', names=['Parameter', 'mode', 'low_bound', 'up_bound'], index_col='Parameter') + self.lb= paras_infos['low_bound'].values + self.ub= paras_infos['up_bound'].values + paras_list=paras_infos.index.tolist() + self.x_labels=paras_list + num_paras=paras_infos.shape[0] + self.n_input=num_paras + ##generate default_values + default_values_path=os.path.join(self.work_path, 'default_values.xlsx') + if not os.path.exists(default_values_path): + default_values=pd.DataFrame(columns=['para_name', 'file_name', 'value']) + + for i in range(num_paras): + para_name=paras_list[i] + file_suffix=self.paras_files.loc[para_name, 'file'] + + if file_suffix in self.HRU_suffix: + for sub in self.total_sub_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + if file_suffix=='sol': + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras_for_sol(file_path, para_name)] + else: + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] + elif file_suffix in self.Watershed_suffiex: + for sub in self.watershed_list: + file_name=sub+"."+file_suffix + file_path=os.path.join(self.work_path, file_name) + default_values.loc[len(default_values)]=[para_name, file_name, _get_default_paras(file_path, para_name)] + elif file_suffix=="bsn": + file_path=os.path.join(self.work_path, "basins."+file_suffix) + default_values.loc[len(default_values)]=[para_name, "basins."+file_suffix, _get_default_paras(file_path, para_name)] + + + default_values.to_excel(default_values_path, index=True) + else: + default_values=pd.read_excel(default_values_path, index_col=0) + self.default_values=default_values + self.paras_infos=paras_infos + self.paras_list=paras_list + def _initial(self): + Watershed={} + + NYSKIP=0 + #read control file fig.cio + with open(os.path.join(self.work_path, "file.cio"), "r") as f: + lines=f.readlines() + for line in lines: + match1 = re.search(r"(\s*)(\d+)(\s*.*NBYR.*)", line) + if match1: + num_years=int(match1.group(2)) + + match2 = re.search(r"(\s*)(\d+)(\s*.*IYR.*)", line) + if match2: + begin_year=int(match2.group(2)) + + match3 = re.search(r"(\s*)(\d+)(\s*.*IDAF.*)", line) + if match3: + begin_day=int(match3.group(2)) + + match4 = re.search(r"(\s*)(\d+)(\s*.*IDAL.*)", line) + if match4: + end_day=int(match4.group(2)) + + match5 = re.search(r"(\s*)(\d+)(\s*.*NYSKIP.*)", line) + if match5: + NYSKIP=int(match5.group(2)) + + self.begin_date=datetime(begin_year+NYSKIP, 1, 1) + timedelta(begin_day - 1) + self.end_date=datetime(begin_year+num_years-1, 1, 1) + timedelta(end_day - 1) + self.simulation_days=(self.end_date-self.begin_date).days+1 + + #read control file fig.fig + with open(os.path.join(self.work_path, "fig.fig"), "r") as f: + lines=f.readlines() + for line in lines: + match = re.search(r'(\d+)\.sub', line) + if match: + Watershed[match.group(1)]=[] + + #read sub files + for sub in Watershed: + file_name=sub+".sub" + with open(os.path.join(self.work_path, file_name), "r") as f: + lines=f.readlines() + for line in lines: + match = re.search(r'(\d+)\.mgt', line) + if match: + Watershed[sub].append(match.group(1)) + + self.Watershed=Watershed + self.total_sub_list = list(itertools.chain.from_iterable(Watershed.values())) + self.watershed_list=list(Watershed.keys()) + #read para_file + file_path=os.path.join(self.work_path, 'SWAT_paras_files.txt') + self.paras_files = pd.read_csv(file_path, sep=' ', names=['parameter', 'file'], index_col='parameter') + +# import shutil +# import tempfile +# import multiprocessing +# import copy +# class SWAT_CUP_Parallel(): +# def __init__(self, work_path, para_file_name, observe_file_name, swat_exe_path, rch_id, n_parallel=4): +# self.n_parallel=n_parallel +# self.work_temp_dir = tempfile.mkdtemp() +# self.original_work_path=work_path +# self.swat_exe_path=swat_exe_path + +# swat_cup=SWAT_CUP(work_path=work_path, para_file_name=para_file_name, observe_file_name=observe_file_name, swat_exe_path=swat_exe_path, rch_id=rch_id) + +# swat_cup_parallels= [copy.deepcopy(swat_cup) for _ in range(self.n_parallel)] + +# for i in range(self.n_parallel): +# swat_cup_parallels[i].work_path=os.path.join(self.work_temp_dir, f'instance_{i}') + +# def copy_and_run(self, instance_id, inputs_queue): +# temp_dir=self.work_temp_dir +# temp_folder_path = os.path.join(temp_dir, f'instance_{instance_id}') +# shutil.copytree(self.original_work_path, temp_folder_path) + +# outputs=[] +# while not inputs_queue.empty(): +# try: +# intput_index, input_data = inputs_queue.get_nowait() +# Y=SWAT_CUP_Parallel[instance_id].evaluate(input_data.reshape(1,-1)) +# outputs.append(1) +# except Exception as e: +# break + +# shutil.rmtree(temp_dir) +# return outputs + +# def run(self, X): +# manager = multiprocessing.Manager() +# inputs_queue = manager.Queue() + + +# results=[] +# for input_index, input_data in enumerate(X): +# inputs_queue.put((input_index, input_data)) + +# with multiprocessing.Pool(processes=self.num_instances) as pool: +# async_results = [pool.apply_async(self.copy_and_run_instance, args=(i, inputs_queue)) for i in range(self.num_instances)] +# pool.close() +# pool.join() + +# for async_result in async_results: +# try: +# result = async_result.get() +# results.extend(result) +# except Exception as e: +# print(f"Error in processing: {e}") + +# a=1 +# print(results) +# return results +# if __name__ == '__main__': +# swat_pro=SWAT_CUP(work_path="H:\SiHuRiver\model\FuTIanSi001\Scenarios\Test\TxtInOut", +# para_file_name="paras_infos.txt", +# observe_file_name="observed.txt", +# swat_exe_path="swat_64rel.exe", +# rch_id=40) +# problem=swat_pro +# from UQPyL.DoE import LHS +# lhs=LHS(problem=problem) +# X=lhs.sample(12, problem.n_input) + +# parallel=SWAT_CUP_Parallel(work_path="H:\SiHuRiver\model\FuTIanSi001\Scenarios\Test\TxtInOut", +# para_file_name="paras_infos.txt", +# observe_file_name="observed.txt", +# swat_exe_path="swat_64rel.exe", +# rch_id=40, n_parallel=3) +# parallel.run(X) + + +# ##########################main pocess########################## +swat_pro=SWAT_CUP(work_path="H:\SiHuRiver\model\FuTIanSi001\Scenarios\Test\TxtInOut", + para_file_name="paras_infos.txt", + observe_file_name="observed.txt", + swat_exe_path="swat_64rel.exe", + rch_id=40) +problem=swat_pro +# print("################1.Sobol################") +# from UQPyL.sensibility import Sobol +# sobol_method=Sobol(problem=problem, cal_second_order=False) #Using Sobol Sequence and saltelli_sequence +# X=sobol_method.sample(128) +# Y=problem.evaluate(X) +# Si=sobol_method.analyze(X, Y, verbose=True) +from UQPyL.sensibility import Morris +# from UQPyL.sensibility import Morris +morris_method=Morris(problem=problem, num_levels=4) #Using Morris Sampler +X=np.loadtxt('./res_X.txt') +Y=np.loadtxt('./res_Y.txt') +# # X=morris_method.sample(128) +# # Y=problem.evaluate(X) +Si=morris_method.analyze(X, Y, verbose=True) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index d0eec1b4..e0d70fc4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,7 @@ build-backend = "setuptools.build_meta" [project] name = "UQPyL" -version = "2.0.4" +version = "2.0.5" authors = [ {name = "wmtSky", email = "wmtsky@hhu.edu.cn"} ] diff --git a/setup.py b/setup.py index 2488578a..46a6a54e 100644 --- a/setup.py +++ b/setup.py @@ -41,7 +41,7 @@ setup( name="UQPyL", author="wmtSky", - version="2.0.4", + version="2.0.5", author_email="wmtsky@hhu.edu.cn", # ... 其他常规setup参数 ... ext_modules=extensions, # 如果有自定义的编译行为