diff --git a/pyAPEX/Utility/sensitive_0.PAR b/pyAPEX/Utility/sensitive_0.PAR deleted file mode 100644 index b0f1cca..0000000 --- a/pyAPEX/Utility/sensitive_0.PAR +++ /dev/null @@ -1 +0,0 @@ -2, 4, 7, 8, 14, 15, 17, 20, 23, 34, 42, 46, 50, 69, 70, 72 \ No newline at end of file diff --git a/pyAPEX/calibration.py b/pyAPEX/calibration.py deleted file mode 100644 index e1ddfe1..0000000 --- a/pyAPEX/calibration.py +++ /dev/null @@ -1,59 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Dec 21 07:42:35 2022 - -@author: Mahesh.Maskey -""" - - -import os -import os.path -from pathlib import Path -from configobj import ConfigObj -from Utility.apex_utility import print_progress_bar -from Utility.apex_utility import copy_best_output -from Utility.apex_utility import read_summary -from pyCALAPEX import calAPEX -print('/014') - -src_dir = Path(os.path.dirname(os.path.realpath(__file__))) -config = ConfigObj(str(src_dir / 'runtime.ini')) -site = 'Farm_1' -crops = None -scenario = 'grazing' -obs_attribute='runoff' -config['Site'] = site -config['Scenario'] = scenario -out_dir = 'Output/'+ obs_attribute -cluster_dir= f'/home/mahesh.maskey/SWMRU_APEX/OklahomaWRE/APEX_V2_output/' -# cluster_dir = '../Output' -sen_obj = calAPEX(src_dir, config=config, cluster_dir=cluster_dir, site=site, scenario=scenario, - obs_attribute=obs_attribute, mod_attribute='WYLD', out_dir=out_dir) -print('Congratulation! Calibration done') - -print('Importing the outputs within the criteria') - - -if scenario == 'non_grazing': - scenario = 'pyAPEX_n' -else: - scenario = 'pyAPEX_g' - -cluster_dir= f'/home/mahesh.maskey/SWMRU_APEX/OklahomaWRE/APEX_V2_output/{site}/{scenario}/Output' -id_best_runs, local_dir = read_summary(file_name='best_stats.csv', obs_attribute='runoff', in_dir=None) - - -k = 0 -print_progress_bar(k, len(id_best_runs), prefix='Complete', suffix='', decimals=1, length=50, fill='█') -for id_run in id_best_runs: - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='001RUN', iteration = id_run, extension= 'OUT', crop = None) - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='daily_outlet', iteration = id_run, extension= 'csv', crop = None) - if crops is None: - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='daily_basin', iteration = id_run, extension= 'csv', crop = None) - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='annual', iteration = id_run, extension= 'csv', crop = None) - else: - for crop in crops: - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='daily_basin', iteration = id_run, extension= 'csv', crop = crop) - copy_best_output(input_dir=cluster_dir, save_dir=local_dir, attribute='annual', iteration = id_run, extension= 'csv', crop = crop) - print_progress_bar(k, len(id_best_runs), prefix=f'{k}: {id_run}', suffix='Complete', decimals=1, length=50, fill='█') - k = k + 1 diff --git a/pyAPEX/load_module.sh b/pyAPEX/load_module.sh deleted file mode 100644 index 69af272..0000000 --- a/pyAPEX/load_module.sh +++ /dev/null @@ -1,16 +0,0 @@ -#!/bin/bash - -#SBATCH --account=swmru_apex -#SBATCH --job-name="python_module" -#SBATCH --partition atlas -#SBATCH --ntasks 1 -#SBATCH --cpus-per-task=1 -#SBATCH --time 20:00:00 -##SBATCH --mail-user=emailAddress -##SBATCH --mail-type=BEGIN,END,FAIL - -module load miniconda singularity -source activate -conda activate /project/swmru_apex/py_env - -python calibration.py diff --git a/pyAPEX/pyAPEXSCU.py b/pyAPEX/pyAPEXSCU.py deleted file mode 100644 index cfde46d..0000000 --- a/pyAPEX/pyAPEXSCU.py +++ /dev/null @@ -1,654 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Sep 16 13:12:23 2022 - -@author: Mahesh.Maskey, Brian Stucky (USDA) -""" -import os -import pandas as pd -import numpy as np -import subprocess -import random -import fortranformat as ff -from Utility.apex_utility import read_sensitive_params, copy_rename_file, get_scanario_name -from Utility.apex_utility import modify_list -from Utility.apex_utility import get_acy -from Utility.apex_utility import get_daily_sad -from Utility.apex_utility import get_daily_dps -from Utility.apex_utility import get_daily_dws -from Utility.apex_utility import read_param_file -from Utility.overwrite_param import overwrite_param -from Utility.apex_utility import savedata_rel1 -from Utility.apex_utility import organize2save -from Utility.apex_utility import do_validate_fill -from Utility.apex_utility import calculate_nutrients -from Utility.easypy import easypy as ep -from Utility.pyAPEXpost import pyAPEXpost as ap -from datetime import datetime - - -class simAPEX: - def __init__( - self, config, src_dir, winepath, in_obj, model_mode, scale=None, - isall=True - ): - """ - self.config: A configuration object. - src_dir: The location of the source repository for managing APEX runs. - winepath: Path to an executable Wine container image. - in_obj: An inAPEX object. - """ - self.config = config - self.src_dir = src_dir - self.winepath = winepath - self.run_name = self.config['run_name'] - self.file_limits = self.src_dir / 'Utility' / self.config['file_limits'] - - self.dir_calibration = in_obj.dir_calibration - self.dir_sensitivity = os.path.join(os.path.dirname(in_obj.dir_sensitivity), 'OutputSensitivity') - self.dir_uncertainty = os.path.join(os.path.dirname(in_obj.dir_uncertainty), 'OutputUncertainty') - in_obj.dir_sensitivity, in_obj.dir_uncertainty = self.dir_sensitivity, self.dir_uncertainty - self.senstitivty_out_file = os.path.join(self.dir_sensitivity, 'Output_summary.txt') - self.uncertainty_outflie = os.path.join(self.dir_uncertainty, 'Output_summary.txt') - - self.get_range() - if (model_mode == "calibration"): - n_discrete = int(self.config['n_discrete']) - n_simul = int(self.config['n_simulation']) - self.generate_param_set(n_discrete, isall) - else: - dir_res = self.config['dir_calibrate_res'] - file_pem = self.config['file_pem'] - self.file_pem = dir_res + '/' + file_pem - # reads the performance statistics file to get the best params - self.get_stats() - file_parm = self.config['file_param'] - # reads the parameter file from the calibration results - self.file_parm = dir_res + '/' + file_parm - self.read_params() - if (model_mode == "sensitivity"): - # file for basic info - if os.path.exists(self.dir_sensitivity) is False: - os.makedirs(self.dir_sensitivity) - f = open(self.senstitivty_out_file, 'w') - f.close() - self.get_best4sa(scale='daily') - self.generate_sensitive_params(isall) - elif (model_mode == "uncertainty"): - if os.path.exists(self.dir_uncertainty) is False: - os.makedirs(self.dir_uncertainty) - f = open(self.uncertainty_outflie, 'w') - f.close() - self.get_params4ua(scale) - self.generate_uncertaintity_params(isall) - - self.recommended_parameters = self.recc_params - self.n_params = len(self.recc_params) - self.pem = ['CODAD', 'RMSEAD', 'NRMSEAD', 'NSEAD', 'PBIASAD', 'IOAAD', 'OF1AD', 'OF2AD', - 'CODDC', 'RMSEDC', 'NRMSEDC', 'NSEDC', 'PBIASDC', 'IOADC', 'OF1DC', 'OF2DC', - 'CODDV', 'RMSEDV', 'NRMSEDV', 'NSEDV', 'PBIASDV', 'IOADV', 'OF1DV', 'OF2DV', - 'CODAM', 'RMSEAM', 'NRMSEAM', 'NSEAM', 'PBIASAM', 'IOAAM', 'OF1AM', 'OF2AM', - 'CODMC', 'RMSEMC', 'NRMSEMC', 'NSEMC', 'PBIASMC', 'IOAMC', 'OF1MC', 'OF2MC', - 'CODMV', 'RMSEMV', 'NRMSEMV', 'NSEMV', 'PBIASMV', 'IOAMV', 'OF1MV', 'OF2MV', - 'CODAY', 'RMSEAY', 'NRMSEAY', 'NSEAY', 'PBIASAY', 'IOAAY', 'OF1AY', 'OF2AY', - 'CODYC', 'RMSEYC', 'NRMSEYC', 'NSEYC', 'PBIASYC', 'IOAYC', 'OF1YC', 'OF2YC', - 'CODYV', 'RMSEYV', 'NRMSEYV', 'NSEYV', 'PBIASYV', 'IOAYV', 'OF1YV', 'OF2YV'] - if (model_mode == "calibration"): - - id_start = int(self.config['n_start']) - n_simul = int(self.config['n_simulation']) - else: - id_start = int(self.config['n_start']) - id_start = 0 - n_simul = self.parameters_matrix.shape[0] - - self.assign_output_df() - if (model_mode == "calibration"): - df_PEM_runoff, df_PEM_sediment_YSD, df_PEM_sediment_USLE = pd.DataFrame(columns=self.pem), pd.DataFrame( - columns=self.pem), pd.DataFrame(columns=self.pem) - df_PEM_sediment_MUSL, df_PEM_sediment_REMX, df_PEM_sediment_MUSS = pd.DataFrame( - columns=self.pem), pd.DataFrame(columns=self.pem), pd.DataFrame(columns=self.pem) - df_PEM_sediment_MUST, df_PEM_sediment_RUS2, df_PEM_sediment_RUSL = pd.DataFrame( - columns=self.pem), pd.DataFrame(columns=self.pem), pd.DataFrame(columns=self.pem) - self.df_p_set = pd.DataFrame() - else: - if 'RunId' in self.df_calpem_sets: - df_PEM_runoff = pd.DataFrame(self.df_calpem_sets).T - df_PEM_runoff.index = df_PEM_runoff.RunId.values - df_PEM_runoff = df_PEM_runoff.drop(['RunId'], axis=1) - else: - df_PEM_runoff = pd.DataFrame(self.df_calpem_sets).T - df_PEM_sediment_YSD, df_PEM_sediment_USLE = pd.DataFrame(columns=self.pem), pd.DataFrame(columns=self.pem) - df_PEM_sediment_MUSL, df_PEM_sediment_REMX, df_PEM_sediment_MUSS = pd.DataFrame( - columns=self.pem), pd.DataFrame(columns=self.pem), pd.DataFrame(columns=self.pem) - df_PEM_sediment_MUST, df_PEM_sediment_RUS2, df_PEM_sediment_RUSL = pd.DataFrame( - columns=self.pem), pd.DataFrame(columns=self.pem), pd.DataFrame(columns=self.pem) - - self.df_p_set = self.pbest - # start simulation - t0 = datetime.now() - for i in range(id_start, n_simul): - try: - t1 = datetime.now() - if (model_mode == "calibration"): - if isall: - self.pick_param(allparam=True, i=i) - else: - self.pick_param(allparam=False, i=i) - # params = np.array([v.replace(',', '') for v in params], dtype=np.float64) - else: - self.p = self.parameters_matrix[i, :] - self.p = overwrite_param(in_obj.param_file, in_obj.simparam_file, self.p) - modify_list(in_obj.list_file, in_obj.simparam_file) - - t2 = datetime.now() - print('Calling APEXgraze') - if os.name == 'nt': - p = subprocess.run(['APEXgraze.exe'], capture_output=True, text=True) - print(p.stdout) - else: - p = subprocess.run([self.winepath, 'APEXgraze.exe'], capture_output=True, text=True) - print(p.stderr) - print(p.stdout) - # Read run name from APEXRUN.DAT file - curr_directory = os.getcwd() - self.scenario_name = get_scanario_name(curr_directory) - #Saving standard output file together with runs for final summary - copy_rename_file(curr_directory, self.scenario_name, itr_id=i, extension='OUT', in_obj=in_obj, - model_mode=model_mode) - - # rename run_name_[iteration].out e.g., 001RUN_0000106.out - print( - f'Completed simulation {(i + 1)}/{n_simul} in {round((datetime.now() - t2).total_seconds(), 3)} seconds') - modify_list(in_obj.list_file, in_obj.param_file) - df_SAD = get_daily_sad(self.run_name) - df_SAD = calculate_nutrients(df_SAD) - df_DPS = get_daily_dps(self.run_name) - df_DWS = get_daily_dws(self.run_name) - df_annual = get_acy(self.run_name) - is_pesticide = self.config['is_pesticide'] - warm_years = int(self.config['warm_years']) - calib_years = int(self.config['calib_years']) - print('Saving model performance statistics') - df_PEM_runoff = do_validate_fill(model_mode, df_SAD, df_PEM_runoff, in_obj, 'WYLD', 'runoff', i, - warm_years, calib_years) - if is_pesticide: - df_PEM_sediment_YSD = do_validate_fill(model_mode, df_DPS, df_PEM_sediment_YSD, in_obj, 'YSD', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_USLE = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_USLE, in_obj, 'USLE', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_MUSL = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_MUSL, in_obj, 'MUSL', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_REMX = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_REMX, in_obj, 'REMX', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_MUSS = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_MUSS, in_obj, 'MUSS', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_MUST = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_MUST, in_obj, 'MUST', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_RUS2 = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_RUS2, in_obj, 'RUS2', - 'sediment', i, warm_years, calib_years) - df_PEM_sediment_RUSL = do_validate_fill(model_mode, df_SAD, df_PEM_sediment_RUSL, in_obj, 'RUSL', - 'sediment', i, warm_years, calib_years) - file_outlet = 'daily_outlet_' + str(i + 1).zfill(7) + '.csv' - df_outlet = df_DWS[['Y', 'M', 'D', 'RFV', 'WYLD', 'TMX', 'TMN', 'PET', 'Q', 'CN', - 'SSF', 'PRK', 'IRGA', 'USLE', 'MUSL', 'REMX', - 'MUSS', 'MUST', 'RUS2', 'RUSL']] - if is_pesticide: - df_outlet['YSD'] = df_DPS['YSD'] - savedata_rel1(df_outlet, file_outlet, model_mode, in_obj) - - file_basin = 'daily_basin_' + str(i + 1).zfill(7) + '.csv' - df_basin = df_SAD[['Y', 'M', 'D', 'CPNM', 'LAI', 'BIOM', 'STL', - 'STD', 'STDL', 'PRCP', 'WYLD', 'TMX', 'TMN', - 'PET', 'ET', 'Q', 'CN', 'SSF', 'PRK', 'QDR', - 'IRGA', 'USLE', 'MUSL', 'REMX', 'MUSS', 'MUST', - 'RUS2', 'RUSL', 'YN', 'YP', 'QN', 'QP', 'QDRN', - 'QPRP', 'SSFN', 'RSFN', 'QRFN', 'QRFP', 'QDRP', - 'DPRK', 'TN', 'TP']] - savedata_rel1(df_basin, file_basin, model_mode, in_obj) - - file_annual = 'annual_' + str(i + 1).zfill(7) + '.csv' - df_year = df_annual[['YR', 'CPNM', 'YLDG', 'YLDF', 'BIOM', 'WS', - 'NS', 'PS', 'TS', 'AS', 'SS']] - savedata_rel1(df_year, file_annual, model_mode, in_obj) - - crops = df_basin.CPNM.unique() - if len(crops) > 1: - for cp in crops: - file_basin = 'daily_basin_' + str(i + 1).zfill(7) + '_' + cp + '.csv' - df_basin_c = df_basin[df_basin['CPNM'] == cp] - savedata_rel1(df_basin_c, file_basin, model_mode, in_obj) - file_annual = 'annual_' + str(i + 1).zfill(7) + '_' + cp + '.csv' - df_year_c = df_year[df_year['CPNM'] == cp] - savedata_rel1(df_year_c, file_annual, model_mode, in_obj) - - df_p = pd.DataFrame(self.p) - df_p = df_p.T - df_p.columns = self.param_list - if 'RunId' in self.df_p_set: - self.df_p_set.index = self.df_p_set.RunId.values - self.df_p_set = self.df_p_set.drop(['RunId'], axis=1) - self.df_p_set = organize2save(self.df_p_set, df_p, i, axis=0) - print('Saving parmeters in APEXPARM.DAT') - savedata_rel1(self.df_p_set, 'APEXPARM', model_mode, in_obj) - print('Parameters') - print(self.p) - print('---------------------------------------------------------------') - print(f'Completed run no. {i + 1} in {round((datetime.now() - t2).total_seconds(), 3)} seconds') - print('---------------------------------------------------------------\n') - print( - f'Completed simulation {(i + 1)}/{n_simul} in {round((datetime.now() - t1).total_seconds(), 3)} seconds') - except Exception as e: - print(e) - print(f'error occurs in simulation {i}') - continue - print(f'\nCompleted {n_simul - id_start} runs in {round((datetime.now() - t0).total_seconds(), 3)} seconds') - - def get_best4sa(self, scale): - criteria = [float(self.config['COD_criteria']), float(self.config['NSE_criteria']), - float(self.config['PBAIS_criteria'])] - f = open(self.senstitivty_out_file, 'a') - f.writelines('---------Calibration and validation resukt summary') - f.writelines('\n---------------------------------------------------------------') - f.writelines('Calibratiion criteria according to Moriasi et a. (20015)\n') - f.writelines(f'Nash Sutcliffe efficiency, NSE {criteria[1]}\n') - f.writelines(f'Percent Bias, PBIAS: {criteria[2]}\n') - f.writelines(f'Coefficient of determination, COD: {criteria[0]}\n') - f.writelines('\n---------------------------------------------------------------') - stats_out = ap.compile_stats(self.stats, criteria, scale) - if scale == 'daily': - best_value, id_run, best_stats = ap.find_best(stats_out[1], metric='OF2DC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2DC', stats='min') - elif scale == 'monthly': - best_value, id_run, best_stats = ap.find_best(stats_out[2], metric='OF2MC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2MC', stats='min') - else: - best_value, id_run, best_stats = ap.find_best(stats_out[3], metric='OF2YC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2YC', stats='min') - df_pem_criteria = stats_out[0] - df_pem_criteria.to_csv(f'{self.dir_sensitivity}/selected_Statistcs_runoff.csv') - run_vec = df_pem_criteria.RunId - df_stats_daily = stats_out[1] - df_stats_monthly = stats_out[2] - df_stats_yearly = stats_out[3] - df_stats = pd.concat([df_stats_daily, df_stats_monthly, df_stats_yearly], axis=0) - self.stats_citerion = df_stats - self.best_stat = best_stats - self.best_obj_value = best_value - self.id_best = id_run - - f.writelines('Number of solutions: {:7d}'.format(df_pem_criteria.shape[0])) - f.writelines('Best Objective function value: {:.4f}'.format(best_value)) - f.writelines('\nBest iterstion {:7d}'.format(id_run)) - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nModel performnce for the best set at calibration\n') - f.writelines('\nCOD: {:0.3f}'.format(best_stats[1])) - f.writelines('\nRMSE: {:0.3f}'.format(best_stats[2])) - f.writelines('\nnRMSE: {:0.3f}'.format(best_stats[3])) - f.writelines('\nNSE: {:0.3f}'.format(best_stats[4])) - f.writelines('\nPBIAS: {:0.3f}'.format(best_stats[5])) - f.writelines('\nIOA: {:0.3f}'.format(best_stats[6])) - f.writelines('\nObjValue1: {:0.3f}'.format(best_stats[7])) - f.writelines('\nObjValue2: {:0.3f}'.format(best_stats[8])) - f.writelines('\nModel performnce for the best set at validation\n') - f.writelines('\nCOD: {:0.3f}'.format(best_stats[9])) - f.writelines('\nRMSE: {:0.3f}'.format(best_stats[10])) - f.writelines('\nnRMSE: {:0.3f}'.format(best_stats[11])) - f.writelines('\nNSE: {:0.3f}'.format(best_stats[12])) - f.writelines('\nPBIAS: {:0.3f}'.format(best_stats[13])) - f.writelines('\nIOA: {:0.3f}'.format(best_stats[14])) - f.writelines('\nObjValue1: {:0.3f}'.format(best_stats[15])) - f.writelines('\nObjValue2: {:0.3f}'.format(best_stats[16])) - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nRange of model performnce for the runs within criteria') - f.writelines('\nNotes: Last two three charecters refer to sets and time scales:') - f.writelines('\n*AD: all daily set; *DC: daily set for calibration, *DV: daily sets for validation') - f.writelines('\n*D: daily, *M: monthly, *Y: yearly, C: calibration, V: validation') - f.writelines('\n---------------------------------------------------------------') - for col in stats_out[4].columns[1:-2]: - f.writelines(f'\n{col}: {(round(stats_out[4][col].min(), 3))}-{(round(stats_out[4][col].max(), 3))}') - pbest = self.param[self.param['RunId'] == id_run] - self.pbest = pbest - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nBest model parameter set from calibration') - for col in pbest.columns[1:]: - f.writelines(f'\n{col}: {(round(pbest[col].min(), 8))}-{(round(pbest[col].max(), 8))}') - df_param_criteria = ap.get_param_bests(self.file_parm, run_vec, scale) - df_param_criteria.to_csv(f'{self.dir_sensitivity}/selected_APEXPARM.csv') - f.writelines('\n---------------------------------------------------------------') - f.writelines('\n Range of model parameters for the runs within criteria') - f.writelines('\n---------------------------------------------------------------') - for col in df_param_criteria.columns[1:]: - f.writelines( - f'\n{col}: {(round(df_param_criteria[col].min(), 8))}-{(round(df_param_criteria[col].max(), 8))}') - f.close() - return self - - def get_stats(self): - df_pem = pd.read_csv(self.file_pem) - df_pem.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - self.stats = df_pem - return self - - def read_params(self): - df_param = pd.read_csv(self.file_parm) - df_param.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - self.param = df_param - return self - - def get_params4ua(self, scale): - criteria = [float(self.config['COD_criteria']), float(self.config['NSE_criteria']), - float(self.config['PBAIS_criteria'])] - f = open(self.uncertainty_outflie, 'a') - f.writelines('---------Calibration and validation resukt summary') - f.writelines('\n---------------------------------------------------------------') - f.writelines('Calibratiion criteria according to Moriasi et a. (20015)\n') - f.writelines(f'Nash Sutcliffe efficiency, NSE {criteria[1]}\n') - f.writelines(f'Percent Bias, PBIAS: {criteria[2]}\n') - f.writelines(f'Coefficient of determination, COD: {criteria[0]}\n') - f.writelines('\n---------------------------------------------------------------') - stats_out = ap.compile_stats(self.stats, criteria, scale) - if scale == 'daily': - best_value, id_run, best_stats = ap.find_best(stats_out[1], metric='OF2DC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2DC', stats='min') - df_stats4ua = stats_out[1] - elif scale == 'monthly': - best_value, id_run, best_stats = ap.find_best(stats_out[2], metric='OF2MC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2MC', stats='min') - df_stats4ua = stats_out[2] - else: - best_value, id_run, best_stats = ap.find_best(stats_out[3], metric='OF2YC', stats='min') - _, _, self.df_calpem_sets = ap.find_best(stats_out[4], metric='OF2YC', stats='min') - df_stats4ua = stats_out[3] - df_stats_daily = stats_out[1] - df_stats_daily.columns = ['RunId', 'CODC', 'RMSEC', 'NRMSEC', 'NSEC', 'PBIASC', - 'IOAC', 'OF1C', 'OF2C', 'CODV', 'RMSEV', 'NRMSEV', 'NSEV', - 'PBIASV', 'IOAV', 'OF1V', 'OF2V'] - df_stats_daily['SCALE'] = 'daily' - df_stats_monthly = stats_out[2] - df_stats_monthly.columns = ['RunId', 'CODC', 'RMSEC', 'NRMSEC', 'NSEC', 'PBIASC', - 'IOAC', 'OF1C', 'OF2C', 'CODV', 'RMSEV', 'NRMSEV', 'NSEV', - 'PBIASV', 'IOAV', 'OF1V', 'OF2V'] - df_stats_monthly['SCALE'] = 'monthly' - df_stats_yearly = stats_out[3] - df_stats_yearly.columns = ['RunId', 'CODC', 'RMSEC', 'NRMSEC', 'NSEC', 'PBIASC', - 'IOAC', 'OF1C', 'OF2C', 'CODV', 'RMSEV', 'NRMSEV', 'NSEV', - 'PBIASV', 'IOAV', 'OF1V', 'OF2V'] - df_stats_yearly['SCALE'] = 'yearly' - df_stats = pd.concat([df_stats_daily, df_stats_monthly, df_stats_yearly], axis=0) - self.stats_citerion = df_stats - self.best_stat = best_stats - self.best_obj_value = best_value - self.id_best = id_run - f.writelines('Number of solutions: {:7d}'.format(df_stats4ua.shape[0])) - f.writelines('Best Objective function value: {:.4f}'.format(best_value)) - f.writelines('\nBest iterstion {:7d}'.format(id_run)) - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nModel performnce for the best set at calibration\n') - f.writelines('\nCOD: {:0.3f}'.format(best_stats[1])) - f.writelines('\nRMSE: {:0.3f}'.format(best_stats[2])) - f.writelines('\nnRMSE: {:0.3f}'.format(best_stats[3])) - f.writelines('\nNSE: {:0.3f}'.format(best_stats[4])) - f.writelines('\nPBIAS: {:0.3f}'.format(best_stats[5])) - f.writelines('\nIOA: {:0.3f}'.format(best_stats[6])) - f.writelines('\nObjValue1: {:0.3f}'.format(best_stats[7])) - f.writelines('\nObjValue2: {:0.3f}'.format(best_stats[8])) - f.writelines('\nModel performnce for the best set at validation\n') - f.writelines('\nCOD: {:0.3f}'.format(best_stats[9])) - f.writelines('\nRMSE: {:0.3f}'.format(best_stats[10])) - f.writelines('\nnRMSE: {:0.3f}'.format(best_stats[11])) - f.writelines('\nNSE: {:0.3f}'.format(best_stats[12])) - f.writelines('\nPBIAS: {:0.3f}'.format(best_stats[13])) - f.writelines('\nIOA: {:0.3f}'.format(best_stats[14])) - f.writelines('\nObjValue1: {:0.3f}'.format(best_stats[15])) - f.writelines('\nObjValue2: {:0.3f}'.format(best_stats[16])) - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nRange of model performnce for the runs within criteria') - f.writelines('\nNotes: Last two three charecters refer to sets and time scales:') - f.writelines('\n*AD: all daily set; *DC: daily set for calibration, *DV: daily sets for validation') - f.writelines('\n*D: daily, *M: monthly, *Y: yearly, C: calibration, V: validation') - f.writelines('\n---------------------------------------------------------------') - self.stats = df_stats4ua - id_sets = df_stats4ua.RunId.values - for col in stats_out[4].columns[1:-2]: - f.writelines(f'\n{col}: {(round(stats_out[4][col].min(), 3))}-{(round(stats_out[4][col].max(), 3))}') - pbest = self.param[self.param['RunId'] == id_run] - self.pbest = pbest - f.writelines('\n---------------------------------------------------------------') - f.writelines('\nBest model parameter set from calibration') - for col in pbest.columns[1:]: - f.writelines(f'\n{col}: {(round(pbest[col].min(), 8))}-{(round(pbest[col].max(), 8))}') - df_ua_params = ap.get_param_bests(self.file_parm, id_sets, scale) - df_ua_params.to_csv(f'{self.dir_uncertainty}/selected_APEXPARM.csv') - f.writelines('\n---------------------------------------------------------------') - f.writelines('\n Range of model parameters for the runs within criteria') - f.writelines('\n---------------------------------------------------------------') - for col in df_ua_params.columns[1:]: - f.writelines(f'\n{col}: {(round(df_ua_params[col].min(), 8))}-{(round(df_ua_params[col].max(), 8))}') - f.close() - self.ua_params = df_ua_params - return self - - def get_range(self): - # import csv file containing range of parameters with recommended values for specific project - df_param_limit = pd.read_csv(self.file_limits, index_col=0, encoding="ISO-8859-1") - self.param_discription = df_param_limit.columns.to_list() - self.param_list = df_param_limit.iloc[0, :].to_list() - array_param_list = df_param_limit.iloc[1:, :].to_numpy() - self.array_param_list = array_param_list.astype(np.float64) - mat_param_list = np.asmatrix(array_param_list) - self.mat_param_list = np.squeeze(np.asarray(mat_param_list)) - return self - - def generate_param_set(self, n_sim, isall): - self.recc_params = self.mat_param_list[2, :] - if isall: - minp = self.mat_param_list[0, :] - self.min_params = np.array([v.replace(',', '') for v in minp], dtype=np.float64) - maxp = self.mat_param_list[1, :] - self.max_params = np.array([v.replace(',', '') for v in maxp], dtype=np.float64) - n_params = len(self.max_params) - else: - id_sensitive = read_sensitive_params(self.src_dir) - mat_sensitive_limt = self.mat_param_list[:, id_sensitive] - minp = mat_sensitive_limt[0, :] - self.min_params = np.array([v.replace(',', '') for v in minp], dtype=np.float64) - maxp = mat_sensitive_limt[1, :] - self.max_params = np.array([v.replace(',', '') for v in maxp], dtype=np.float64) - recc_sensitive = mat_sensitive_limt[2, :] - n_params = len(recc_sensitive) - - diff = self.max_params - self.min_params - inc = diff / n_sim - mat_params = np.zeros((n_sim + 1, n_params)) - if isall: - mat_params[0, :] = self.recc_params - else: - mat_params[0, :] = recc_sensitive - - for i in range(n_sim): - for j in range(n_params): - mat_params[i + 1, j] = self.min_params[j] + i * inc[j] - self.parameters_matrix = mat_params - return self - - def pick_param(self, i, allparam=True): - nset, nparam = self.parameters_matrix.shape - p = np.zeros((nparam)) - random.seed(i) - id_rands = ep.ranom_intmatrix(nparam, 1, limit=(0, nset - 1)) - for ip in range(nparam): - p[ip] = self.parameters_matrix[int(id_rands[ip]), ip] - if allparam: - p = p - else: - p_all = self.recc_params - id_sensitive = read_sensitive_params(self.src_dir) - p_all[id_sensitive] = p - p = p_all - self.p = p - return self.p - - def generate_sensitive_params(self, isall_try, p=None): - maxp = int(self.config['max_range']) - inc = float(self.config['increment']) - minp = -maxp - deltas = np.arange(minp, maxp + 1, inc) - nset = len(deltas) - p = self.pbest.T.values[1:].ravel() - if p is None: - recc_params = np.array(self.mat_param_list[2, :]) - else: - recc_params = p - self.recc_params = recc_params - lb_params = self.mat_param_list[0, :].astype(float) - ub_params = self.mat_param_list[1, :].astype(float) - diff = ub_params - lb_params - nparams = len(recc_params) - if isall_try: - # Changing all the parameters with in the range at once - mat_params = ep.nanmatrix(nset, nparams) - for i in range(nparams): - if i <= 69: - mat_params[:, i] = recc_params[i] - else: - for j in range(nset): - p_sen = recc_params[i] + diff[i] * deltas[j] * 0.01 - if p_sen <= lb_params[i]: - p_sen = lb_params[i] - elif p_sen >= ub_params[i]: - p_sen = ub_params[i] - else: - p_sen = p_sen - mat_params[j, i] = p_sen - # Increasing or decreasing all the parameters individually - id_params = np.arange(70, nparams)[:-4] - for id_ in id_params: - mat_id = ep.nanmatrix(nset, nparams) - mat_id[:, :] = recc_params - p_vec = [] - for j in range(nset): - p_sen = recc_params[id_] + diff[id_] * deltas * 0.01 - if p_sen <= lb_params[id_]: - p_sen = lb_params[id_] - elif p_sen >= ub_params[id_]: - p_sen = ub_params[id_] - else: - p_sen = p_sen - p_vec.append(p_sen) - mat_id[:, id_] = p_vec - mat_params = np.concatenate((mat_params, mat_id), axis=0) - del mat_id - else: - mat_params = ep.nanmatrix(nset, nparams) - idsensitive = read_sensitive_params(self.src_dir) - recc_params = p - mat_params[:, :] = recc_params - min_params = lb_params[idsensitive] - max_params = ub_params[idsensitive] - diff = max_params - min_params - # Increasing or decreasing all the selected parameters with in the range at once - for i in range(len(idsensitive)): - for j in range(nset): - p_sen = recc_params[idsensitive[i]] + diff[i] * deltas[j] * 0.01 - if p_sen <= min_params[i]: - p_sen = min_params[i] - elif p_sen >= max_params[i]: - p_sen = max_params[i] - else: - p_sen = p_sen - mat_params[j, idsensitive[i]] = p_sen - del i - # Increasing or decreasing all the selected parameters individually - for i in range(len(idsensitive)): - mat_id = ep.nanmatrix(nset, nparams) - mat_id[:, :] = recc_params - p_vec = [] - for j in range(nset): - p_sen = recc_params[idsensitive[i]] + diff[i] * deltas[j] * 0.01 - if p_sen <= min_params[i]: - p_sen = min_params[i] - elif p_sen >= max_params[i]: - p_sen = max_params[i] - else: - p_sen = p_sen - p_vec.append(p_sen) - mat_id[:, idsensitive[i]] = p_vec - mat_params = np.concatenate((mat_params, mat_id), axis=0) - del mat_id - self.parameters_matrix = mat_params - # from numpy import savetxt - # file_text = os.path.join(os.path.dirname(self.senstitivty_out_file), 'APEXPARM_SA.txt') - # savetxt(file_text, mat_params, delimiter=',') - return self - - def generate_uncertaintity_params(self, isall_try): - maxp = int(self.config['max_range_uncertaintity']) - inc = float(self.config['increment_uncertainty']) - minp = -maxp - deltas = np.arange(minp, maxp + inc, inc) - nset = len(deltas) - p = self.pbest.T.values[1:].ravel() - self.recc_params = p - df_ua_params = self.ua_params - ua_params = df_ua_params.iloc[:, 1:] - mu_vec = ua_params.mean(axis=0).values - std_vec = ua_params.std(axis=0).values - lb_params = self.mat_param_list[0, :].astype(float) - ub_params = self.mat_param_list[1, :].astype(float) - nparams = len(mu_vec) - if isall_try: - # Changing all the parameters with in the range at once - mat_params = ep.nanmatrix(nset, nparams) - for i in range(nparams): - if i <= 69: - mat_params[:, i] = mu_vec[i] - else: - for j in range(nset): - p_un = mu_vec[i] + std_vec[i] * deltas[j] - if p_un <= lb_params[i]: - p_un = lb_params[i] - elif p_un >= ub_params[i]: - p_un = ub_params[i] - else: - p_un = p_un - mat_params[j, i] = p_un - else: - mat_params = ep.nanmatrix(nset, nparams) - idsensitive = read_sensitive_params(self.src_dir) - mat_params[:, :] = p - min_params = lb_params[idsensitive] - max_params = ub_params[idsensitive] - # Increasing or decreasing all the selected parameters with in the range at once - for i in range(len(idsensitive)): - for j in range(nset): - p_un = mu_vec[idsensitive[i]] + std_vec[idsensitive[i]] * deltas[j] - if p_un <= min_params[i]: - p_un = min_params[i] - elif p_un >= max_params[i]: - p_un = max_params[i] - else: - p_un = p_un - mat_params[j, idsensitive[i]] = p_un - self.parameters_matrix = mat_params - return self - - def assign_output_df(self): - self.df_WYLD, self.df_LAI, self.df_BIOM = pd.DataFrame(), pd.DataFrame(), pd.DataFrame() - self.df_PRCP, self.df_ET, self.df_PET = pd.DataFrame(), pd.DataFrame(), pd.DataFrame() - self.df_YSD, self.df_USLE = pd.DataFrame(), pd.DataFrame() - self.df_MUSL, self.df_REMX = pd.DataFrame(), pd.DataFrame() - self.df_MUSS, self.df_MUST = pd.DataFrame(), pd.DataFrame() - self.df_RUS2, self.df_RUSL = pd.DataFrame(), pd.DataFrame() - self.df_STL, self.df_STD, self.df_STDL = pd.DataFrame(), pd.DataFrame(), pd.DataFrame() - self.df_YN, self.df_YP = pd.DataFrame(), pd.DataFrame() - self.df_QN, self.df_QP = pd.DataFrame(), pd.DataFrame() - self.df_QDRN, self.df_QDRP = pd.DataFrame(), pd.DataFrame() - self.df_QRFN, self.df_QRFP = pd.DataFrame(), pd.DataFrame() - self.df_SSFN, self.df_RSFN = pd.DataFrame(), pd.DataFrame() - self.df_TN, self.df_TP = pd.DataFrame(), pd.DataFrame() - self.df_DPRK = pd.DataFrame() - self.df_YLDG, self.df_YLDF = pd.DataFrame(), pd.DataFrame() - self.df_WS, self.df_TS = pd.DataFrame(), pd.DataFrame() - self.df_NS, self.df_PS = pd.DataFrame(), pd.DataFrame() - return self diff --git a/pyAPEX/pyAPEXin.py b/pyAPEX/pyAPEXin.py deleted file mode 100644 index 37af5df..0000000 --- a/pyAPEX/pyAPEXin.py +++ /dev/null @@ -1,67 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 15 15:16:44 2022 - -@author: Mahesh.Maskey, Brian Stucky (USDA) -""" - -import pandas as pd -import numpy as np -import fortranformat as ff -from Utility.apex_utility import read_param_file -from Utility.apex_utility import txt2list - -class inAPEX: - def __init__(self, config, prog_dir): - """ - config: A run configuration object. - prog_dir: The location of the APEX simulation softare. - """ - # Specify the paths for all required input files. - self.control_file = prog_dir / 'APEXCONT.DAT' - self.list_file = prog_dir / 'APEXFILE.DAT' - self.param_file = prog_dir / 'APEXPARM.DAT' - self.simparam_file = prog_dir / 'simAPEXPARM.DAT' - self.file_observe = prog_dir / config['file_observe'] - self.sim_id_range = f'_{(config["n_start"] + 1):07}-{(config["n_simulation"]):07}' - self.dir_calibration = prog_dir.parent / (config['dir_calibration'] + self.sim_id_range) - self.dir_sensitivity = prog_dir.parent / (config['dir_sensitivity'] + self.sim_id_range) - self.dir_uncertainty = prog_dir.parent / (config['dir_uncertainty'] + self.sim_id_range) - - self.get_control_period() - # Search subarea file in APEXFILE list - file_list = txt2list(self.list_file) - df_file_list = pd.DataFrame(file_list) - list_area = df_file_list.iloc[np.where(df_file_list.iloc[:, 0]=='FSUBA')[0][0], 1] - area_list = txt2list(prog_dir / list_area) - self.file_subarea = prog_dir / area_list[0][1] - self.watershed_area = self.get_watershed_area() - self.ndates = len(self.vec_date) - self.nmonth = len(self.vec_month) - self.nyears = len(self.vec_year) - - def get_control_period(self): - # Extracts the simulation period, start date from APEXCONT.DAT file - # and computes the date vector in three time scales: days, month, and year - lines = read_param_file(self.control_file ) - read_format = ff.FortranRecordReader('(20I6)') - line_read = read_format.read(lines[0]) - self.n_years = line_read[0] - self.startyear = line_read[1] - self.startmonth = line_read[2] - self.startday = line_read[3] - self.lastyear = self.startyear+self.n_years-1 - self.startdate = pd.to_datetime(str(self.startyear) + '/' +str(self.startmonth) + '/' + str(self.startday)) - self.stopdate = pd.to_datetime(str(self.lastyear) + '/12/31') - self.vec_date = pd.date_range(self.startdate, self.stopdate, freq='d') - self.vec_month = pd.date_range(self.startdate, self.stopdate, freq='M') - self.vec_year = pd.date_range(self.startdate, self.stopdate, freq='Y') - return self - - def get_watershed_area(self): - # Reads the subarea file and extracts the watershed area in hactare. - lines = read_param_file(self.file_subarea) - read_format = ff.FortranRecordReader('(10F8.3)') - line_read = read_format.read(lines[3]) - return line_read[0] - diff --git a/pyAPEX/pyCALAPEX.py b/pyAPEX/pyCALAPEX.py deleted file mode 100644 index 43d7940..0000000 --- a/pyAPEX/pyCALAPEX.py +++ /dev/null @@ -1,416 +0,0 @@ -# ## Calibration of APEX model for surface runoff data - -# ### Import libraries - -import os -import pandas as pd -import numpy as np -import warnings -from configobj import ConfigObj -from Utility.pyAPEXpost import pyAPEXpost as ap - -warnings.filterwarnings('ignore') -config = ConfigObj('runtime.ini') - -# ### File locations and input parameters -# Main input -# cluster_dir = 'G:/PostDocResearch/USDA-ARS/Project/OklahomaWRE/Cluster_data/' -# mod_attribute = ['WYLD', 'USLE', 'MUSL', 'REMX', 'MUSS', 'RUS2', 'RUSL', 'YSD'] -class calAPEX: - def __init__(self, src_dir, config, cluster_dir, site, scenario, obs_attribute, mod_attribute, out_dir): - - """ - - - - """ - self.stats_daily = None - self.config = config - self.obs_attribute = obs_attribute - self.site = site - self.scenario = scenario - self.mod_attribute = mod_attribute - self.src_dir = src_dir - self.cluster_dir = cluster_dir - self.metrics = ['OF', 'NSE', 'PBIAS', 'COD'] - self.cal_dir = config['dir_calibration'] - self.out_dir = out_dir - if not os.path.isdir(self.out_dir): - os.makedirs(self.out_dir) - self.get_pe_files() # Set the locations of modeled and process data - # #### Set the criteria - if obs_attribute == 'runoff': - self.criteria = [float(config['COD_criteria']), float(config['NSE_criteria']), - float(config['PBAIS_criteria'])] - elif obs_attribute == 'sediment': - self.criteria = [float(config['COD_criteria_sediment']), float(config['NSE_criteria_sediment']), - float(config['PBAIS_criteria_sediment'])] - elif obs_attribute == 'Soil_erosion': - self.criteria = [float(config['COD_criteria_sediment']), float(config['NSE_criteria_sediment']), - float(config['PBAIS_criteria_sediment'])] - self.get_stats() # ### Read statistics - self.daily_stats() - self.monthly_stats() - self.yearly_stats() - # combine all the stats in long form - self.df_stats = pd.concat([self.stats_daily, self.stats_monthly, self.stats_yearly], axis=0) - self.count_stats() - self.df_stats.to_csv(f'{self.out_dir}/best_stats.csv') - # Merge stata summer and export - # try: - self.extract_objective_function(scale='daily') - self.extract_objective_function(scale='monthly') - self.extract_objective_function(scale='yearly') - # except: - # self.get_stats_by_metric(scale='daily') - # self.get_stats_by_metric(scale='monthly') - # self.get_stats_by_metric(scale='yearly') - # self.extract_objective_function(scale='daily') - # self.extract_objective_function(scale='monthly') - # self.extract_objective_function(scale='yearly') - self.df_stats = pd.concat([self.best_stats_daily, self.best_stats_monthly, self.best_stats_yearly], axis=0) - self.df_stats.to_csv(f'{self.out_dir}/summary_stats.csv') - # Import parameters - self.import_parameter() - # concate parameter summary and export - self.df_param_best = pd.concat( - [self.df_param_best_daily, self.df_param_best_monthly, self.df_param_best_yearly], axis=0) - self.df_param_best.to_csv(f'{self.out_dir}/summary_APEXPARM.csv') - # ### Import modeled data - self.assign_read_dir() - self.import_modeled_data(scale='daily', croplist=None) - self.import_modeled_data(scale='monthly', croplist=None) - self.import_modeled_data(scale='yearly', croplist=None) - # Read observed data - self.df_observed = ap.get_measure(data_dir='Program', file_name='calibration_data.csv') - - # Seperate data for model and simulation - ## Finalize model data for outlet - # for objective function - df_model_outlet_daily, df_calibration_daily = self.final_data(scale='daily', metric='OF', location='outlet') - df_model_outlet_monthly, df_calibration_monthly = self.final_data(scale='monthly', metric='OF', - location='outlet') - df_model_outlet_yearly, df_calibration_yearly = self.final_data(scale='yearly', metric='OF', location='outlet') - self.df_model_out_OF = pd.concat([df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly], - axis=0) - del df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly - self.df_model_out_OF.to_csv(f'{self.out_dir}/model_result_outlet_OF.csv') - self.df_model_calibration = pd.concat([df_calibration_daily, df_calibration_monthly, df_calibration_yearly], - axis=0) - self.df_model_calibration.to_csv(f'{self.out_dir}/model_calibration_OF.csv') - df_model_basin_daily = self.final_data(scale='daily', metric='OF', location='basin') - df_model_basin_monthly = self.final_data(scale='monthly', metric='OF', location='basin') - df_model_basin_yearly = self.final_data(scale='yearly', metric='OF', location='basin') - self.df_model_basin_OF = pd.concat([df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly], - axis=0) - self.df_model_basin_OF.to_csv(f'{self.out_dir}/model_result_basin_OF.csv') - del df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly - df_model_annual_daily = self.final_data(scale='daily', metric='OF', location='annual') - df_model_annual_monthly = self.final_data(scale='monthly', metric='OF', location='annual') - df_model_annual_yearly = self.final_data(scale='yearly', metric='OF', location='annual') - self.df_model_annual_OF = pd.concat([df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly], - axis=0) - self.df_model_annual_OF.to_csv(f'{self.out_dir}/model_result_annual_OF.csv') - del df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly - - # for Nash-Sutcliffe efficiency - df_model_outlet_daily, df_calibration_daily = self.final_data(scale='daily', metric='NSE', location='outlet') - df_model_outlet_monthly, df_calibration_monthly = self.final_data(scale='monthly', metric='NSE', - location='outlet') - df_model_outlet_yearly, df_calibration_yearly = self.final_data(scale='yearly', metric='NSE', location='outlet') - self.df_model_out_NSE = pd.concat([df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly], - axis=0) - del df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly - self.df_model_out_NSE.to_csv(f'{self.out_dir}/model_result_outlet_NSE.csv') - self.df_model_calibration = pd.concat([df_calibration_daily, df_calibration_monthly, df_calibration_yearly], - axis=0) - self.df_model_calibration.to_csv(f'{self.out_dir}/model_calibration_NSE.csv') - df_model_basin_daily = self.final_data(scale='daily', metric='NSE', location='basin') - df_model_basin_monthly = self.final_data(scale='monthly', metric='NSE', location='basin') - df_model_basin_yearly = self.final_data(scale='yearly', metric='NSE', location='basin') - self.df_model_basin_NSE = pd.concat([df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly], - axis=0) - self.df_model_basin_NSE.to_csv(f'{self.out_dir}/model_result_basin_NSE.csv') - del df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly - df_model_annual_daily = self.final_data(scale='daily', metric='NSE', location='annual') - df_model_annual_monthly = self.final_data(scale='monthly', metric='NSE', location='annual') - df_model_annual_yearly = self.final_data(scale='yearly', metric='NSE', location='annual') - self.df_model_annual_NSE = pd.concat([df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly], - axis=0) - self.df_model_annual_NSE.to_csv(f'{self.out_dir}/model_result_annual_NSE.csv') - del df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly - - # for Coefficient of determination - df_model_outlet_daily, df_calibration_daily = self.final_data(scale='daily', metric='COD', location='outlet') - df_model_outlet_monthly, df_calibration_monthly = self.final_data(scale='monthly', metric='COD', - location='outlet') - df_model_outlet_yearly, df_calibration_yearly = self.final_data(scale='yearly', metric='COD', location='outlet') - self.df_model_out_COD = pd.concat([df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly], - axis=0) - del df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly - self.df_model_out_COD.to_csv(f'{self.out_dir}/model_result_outlet_COD.csv') - self.df_model_calibration = pd.concat([df_calibration_daily, df_calibration_monthly, df_calibration_yearly], - axis=0) - self.df_model_calibration.to_csv(f'{self.out_dir}/model_calibration_COD.csv') - df_model_basin_daily = self.final_data(scale='daily', metric='COD', location='basin') - df_model_basin_monthly = self.final_data(scale='monthly', metric='COD', location='basin') - df_model_basin_yearly = self.final_data(scale='yearly', metric='COD', location='basin') - self.df_model_basin_COD = pd.concat([df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly], - axis=0) - self.df_model_basin_COD.to_csv(f'{self.out_dir}/model_result_basin_COD.csv') - del df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly - df_model_annual_daily = self.final_data(scale='daily', metric='COD', location='annual') - df_model_annual_monthly = self.final_data(scale='monthly', metric='COD', location='annual') - df_model_annual_yearly = self.final_data(scale='yearly', metric='COD', location='annual') - self.df_model_annual_COD = pd.concat([df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly], - axis=0) - self.df_model_annual_COD.to_csv(f'{self.out_dir}/model_result_annual_COD.csv') - del df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly - - # for Perecnt bias - df_model_outlet_daily, df_calibration_daily = self.final_data(scale='daily', metric='PBIAS', location='outlet') - df_model_outlet_monthly, df_calibration_monthly = self.final_data(scale='monthly', metric='PBIAS', - location='outlet') - df_model_outlet_yearly, df_calibration_yearly = self.final_data(scale='yearly', metric='PBIAS', - location='outlet') - self.df_model_out_PBIAS = pd.concat([df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly], - axis=0) - del df_model_outlet_daily, df_model_outlet_monthly, df_model_outlet_yearly - self.df_model_out_PBIAS.to_csv(f'{self.out_dir}/model_result_outlet_PBIAS.csv') - self.df_model_calibration = pd.concat([df_calibration_daily, df_calibration_monthly, df_calibration_yearly], - axis=0) - self.df_model_calibration.to_csv(f'{self.out_dir}/model_calibration_PBIAS.csv') - df_model_basin_daily = self.final_data(scale='daily', metric='PBIAS', location='basin') - df_model_basin_monthly = self.final_data(scale='monthly', metric='PBIAS', location='basin') - df_model_basin_yearly = self.final_data(scale='yearly', metric='PBIAS', location='basin') - self.df_model_basin_PBIAS = pd.concat([df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly], - axis=0) - self.df_model_basin_PBIAS.to_csv(f'{self.out_dir}/model_result_basin_PBIAS.csv') - del df_model_basin_daily, df_model_basin_monthly, df_model_basin_yearly - df_model_annual_daily = self.final_data(scale='daily', metric='PBIAS', location='annual') - df_model_annual_monthly = self.final_data(scale='monthly', metric='PBIAS', location='annual') - df_model_annual_yearly = self.final_data(scale='yearly', metric='PBIAS', location='annual') - self.df_model_annual_PBIAS = pd.concat([df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly], - axis=0) - self.df_model_annual_PBIAS.to_csv(f'{self.out_dir}/model_result_annual_PBIAS.csv') - del df_model_annual_daily, df_model_annual_monthly, df_model_annual_yearly - - def get_pe_files(self): - if (self.mod_attribute == 'WYLD'): - file_pe = os.path.join(self.cal_dir, 'Statistics_runoff.csv') - elif (self.mod_attribute == 'YSD'): - file_pe = os.path.join(self.cal_dir, 'Statistics_sediment.csv') - else: - file_pe = os.path.join(self.cal_dir, f'Statistics_Soil_erosion_{self.mod_attribute}.csv') - file_parameter = os.path.join(self.cal_dir, 'APEXPARM.csv') - self.file_pem = file_pe - self.file_param = file_parameter - return self - - def assign_read_dir(self): - if self.scenario == 'non_grazing': - scenario1 = 'pyAPEX_n' - else: - scenario1 = 'pyAPEX_g' - self.read_dir = os.path.join(self.cluster_dir, - self.site, scenario1, 'Output') - return self - - def create_output_dir(self): - self.out_dir = os.path.join(self.out_dir, self.obs_attribute) - if os.path.exists(self.out_dir) is False: - os.makedirs(self.out_dir) - return self - - def count_stats(self): - criteria = self.criteria - nCOD_daily = np.sum(self.df_pem.CODDC >= criteria[0]) - nNSE_daily = np.sum(self.df_pem.NSEDC >= criteria[1]) - nPBIAS_daily = np.sum(self.df_pem.PBIASDC.abs() <= criteria[2]) - nCOD_monthly = np.sum(self.df_pem.CODMC >= criteria[0]) - nNSE_monthly = np.sum(self.df_pem.NSEMC >= criteria[1]) - nPBIAS_monthly = np.sum(self.df_pem.PBIASMC.abs() <= criteria[2]) - nCOD_yearly = np.sum(self.df_pem.CODYC >= criteria[0]) - nNSE_yearly = np.sum(self.df_pem.NSEYC >= criteria[1]) - nPBIAS_yearly = np.sum(self.df_pem.PBIASYC.abs() <= criteria[2]) - self.dict_count_stats = {'CODD': nCOD_daily, 'NSED': nNSE_daily, 'PBIASD': nPBIAS_daily, - 'CODM': nCOD_monthly, 'NSEM': nNSE_monthly, 'PBIASM': nPBIAS_monthly, - 'CODY': nCOD_yearly, 'NSEY': nNSE_yearly, 'PBIASY': nPBIAS_yearly} - return self - - def get_stats(self): - df_pem = pd.read_csv(self.file_pem) - df_pem.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - self.df_pem = df_pem - return self - - def daily_stats(self): - # get best stats at daily scale - stats_daily = ap.compile_stats(df=self.df_pem, criteria=self.criteria, scale='daily') - self.stats_daily = stats_daily[0] - self.daily2daily = stats_daily[1] - self.daily2monthly = stats_daily[2] - self.daily2yearly = stats_daily[3] - self.daily_all = stats_daily[4] - return self - - def get_stats_by_metric(self, scale='daily'): - # get best stats at daily scale - stats_COD = ap.compile_stats_by_metrics(df=self.df_pem, criteria=self.criteria, scale=scale, metric='COD') - stats_NSE = ap.compile_stats_by_metrics(df=self.df_pem, criteria=self.criteria, scale=scale, metric='NSE') - stats_PBIAS = ap.compile_stats_by_metrics(df=self.df_pem, criteria=self.criteria, scale=scale, metric='PBIAS') - if scale == 'daily': - self.stats_daily_COD = stats_COD[0] - self.daily2daily_COD = stats_COD[1] - self.daily2monthly_COD = stats_COD[2] - self.daily2yearly_COD = stats_COD[3] - self.daily_all_COD = stats_COD[4] - self.stats_daily_NSE = stats_NSE[0] - self.daily2daily_NSE = stats_NSE[1] - self.daily2monthly_NSE = stats_NSE[2] - self.daily2yearly_NSE = stats_NSE[3] - self.daily_all_NSE = stats_NSE[4] - self.stats_daily_PBIAS = stats_PBIAS[0] - self.daily2daily_PBIAS = stats_PBIAS[1] - self.daily2monthly_PBIAS = stats_PBIAS[2] - self.daily2yearly_PBIAS = stats_PBIAS[3] - self.daily_all_PBIAS = stats_PBIAS[4] - if scale == 'monthly': - self.stats_monthly_COD = stats_COD[0] - self.monthly2daily_COD = stats_COD[1] - self.monthly2monthly_COD = stats_COD[2] - self.monthly2yearly_COD = stats_COD[3] - self.monthly_all_COD = stats_COD[4] - self.stats_monthly_NSE = stats_NSE[0] - self.monthly2daily_NSE = stats_NSE[1] - self.monthly2monthly_NSE = stats_NSE[2] - self.monthly2yearly_NSE = stats_NSE[3] - self.monthly_all_NSE = stats_NSE[4] - self.stats_monthly_PBIAS = stats_PBIAS[0] - self.monthly2daily_PBIAS = stats_PBIAS[1] - self.monthly2monthly_PBIAS = stats_PBIAS[2] - self.monthly2yearly_PBIAS = stats_PBIAS[3] - self.monthly_all_PBIAS = stats_PBIAS[4] - if scale == 'yearly': - self.stats_yearly_COD = stats_COD[0] - self.yearly2daily_COD = stats_COD[1] - self.yearly2monthly_COD = stats_COD[2] - self.yearly2yearly_COD = stats_COD[3] - self.yearly_all_COD = stats_COD[4] - self.stats_yearly_NSE = stats_NSE[0] - self.yearly2daily_NSE = stats_NSE[1] - self.yearly2monthly_NSE = stats_NSE[2] - self.yearly2yearly_NSE = stats_NSE[3] - self.yearly_all_NSE = stats_NSE[4] - self.stats_yearly_PBIAS = stats_PBIAS[0] - self.yearly2daily_PBIAS = stats_PBIAS[1] - self.yearly2monthly_PBIAS = stats_PBIAS[2] - self.yearly2yearly_PBIAS = stats_PBIAS[3] - self.yearly_all_PBIAS = stats_PBIAS[4] - return self - - def monthly_stats(self): - # get best stats at monthly scale - stats_monthly = ap.compile_stats(df=self.df_pem, criteria=self.criteria, scale='monthly') - self.stats_monthly = stats_monthly[0] - self.monthly2daily = stats_monthly[1] - self.monthly2monthly = stats_monthly[2] - self.monthly2yearly = stats_monthly[3] - self.monthly_all = stats_monthly[4] - return self - - def yearly_stats(self): - # get best stats at yearly scale - stats_yearly = ap.compile_stats(df=self.df_pem, criteria=self.criteria, scale='daily') - self.stats_yearly = stats_yearly[0] - self.yearly2daily = stats_yearly[1] - self.yearly2monthly = stats_yearly[2] - self.yearly2yearly = stats_yearly[3] - self.yearly_all = stats_yearly[4] - return self - - def extract_objective_function(self, scale): - if scale == 'daily': - df = self.daily2daily - metrics = ['OF2DC', 'NSEDC', 'PBIASDC', 'CODDC'] - elif scale == 'monthly': - df = self.monthly2monthly - metrics = ['OF2MC', 'NSEMC', 'PBIASMC', 'CODMC'] - else: - df = self.yearly2yearly - metrics = ['OF2YC', 'NSEYC', 'PBIASYC', 'CODYC'] - stats = ['min', 'max', 'min', 'max'] - best_stats, ids_best = ap.summarize_stats(df, scale, metrics, stats) - if scale == 'daily': - self.best_stats_daily, self.ids_best_daily = best_stats, ids_best - elif scale == 'monthly': - self.best_stats_monthly, self.ids_best_monthly = best_stats, ids_best - else: - self.best_stats_yearly, self.ids_best_yearly = best_stats, ids_best - return self - - def import_parameter(self): - # daily - self.df_param_best_daily = ap.get_best_params(file_name=self.file_param, - ids_bests=self.ids_best_daily, - scale='daily') - # monthly - self.df_param_best_monthly = ap.get_best_params(file_name=self.file_param, - ids_bests=self.ids_best_monthly, - scale='monthly') - # yearly - self.df_param_best_yearly = ap.get_best_params(file_name=self.file_param, - ids_bests=self.ids_best_yearly, - scale='yearly') - return self - - def import_modeled_data(self, scale, croplist=None): - if scale == 'daily': - ids_best = self.ids_best_daily - elif scale == 'monthly': - ids_best = self.ids_best_monthly - else: - ids_best = self.ids_best_yearly - - outlet, basin, annual = ap.import_save(dir_data=self.read_dir, - attribute=self.obs_attribute, - croplist=croplist, - ids=ids_best, - scale=scale, - dir_save=self.out_dir) - if scale == 'daily': - self.outlet_daily, self.basin_daily, self.annual_daily = outlet, basin, annual - elif scale == 'monthly': - self.outlet_monthly, self.basin_monthly, self.annual_monthly = outlet, basin, annual - else: - self.outlet_yearly, self.basin_yearly, self.annual_yearly = outlet, basin, annual - return self - - def final_data(self, scale, metric, location): - if location == 'outlet': - outlet_result = ap.finalize_outlet_result(self.out_dir, - self.df_observed, - self.obs_attribute, - self.mod_attribute, - scale, metric) - df_model_daily = outlet_result[0] - df_calibration_daily = outlet_result[1] - return df_model_daily, df_calibration_daily - elif location == 'basin': - df_outlet_daily, _ = self.final_data(scale, metric, location='outlet') - df_model_daily = ap.finalize_basin_result(self.out_dir, - self.df_observed, - df_outlet_daily, - self.obs_attribute, - self.mod_attribute, - scale, metric) - return df_model_daily - elif location == 'annual': - df_outlet_daily, _ = self.final_data(scale, metric, location='outlet') - df_model_yearly = ap.finalize_annual_result(self.out_dir, - self.df_observed, - df_outlet_daily, - self.obs_attribute, - self.mod_attribute, - scale, metric) - return df_model_yearly diff --git a/pyAPEX/pySAAPEX.py b/pyAPEX/pySAAPEX.py deleted file mode 100644 index 8e0627d..0000000 --- a/pyAPEX/pySAAPEX.py +++ /dev/null @@ -1,410 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Dec 13 09:44:44 2022 - -@author: Mahesh.Maskey -""" - -import os -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt -import warnings -from configobj import ConfigObj -import seaborn as sns -from Utility.pyAPEXpost import print_progress_bar -from Utility.pyAPEXpost import pyAPEXpost as ap -from Utility.apex_utility import read_sensitive_params -from Utility.apex_utility import split_data -from Utility.easypy import easypy as ep -from Utility.apex_utility import plot_2_bar - -warnings.filterwarnings('ignore') -config = ConfigObj('runtime.ini') -print('\014') - - -# noinspection PyAttributeOutsideInit,PyTypeChecker -class senanaAPEX: - # noinspection PyTypeChecker - def __init__(self, src_dir, config_obj, out_dir, attribute, metric='OF'): - - """ - - - - """ - - self.config = config_obj - self.src_dir = src_dir - self.file_limits = self.src_dir / 'Utility' / self.config['file_limits'] - self.attribute = attribute - self.folder = config_obj['dir_sensitivity'] - self.get_range() - self.df_obs = ap.get_measure(data_dir='Program', file_name='calibration_data.csv') - id_sensitive = read_sensitive_params(self.src_dir) - for i in range(len(id_sensitive)): - id_sensitive[i] = id_sensitive[i] + 1 - param_range = self.param_range[:, id_sensitive] - self.list_bound = list(param_range[:2, :].T) - self.get_pe_files() - self.get_params(id_sensitive, is_all=False) - metric_list = ['RunId', 'CODDC', 'RMSEDC', 'NRMSEDC', 'NSEDC', 'PBIASDC', 'OF2DC'] - self.read_pem(metric_list) - self.params.insert(self.params.shape[1], 'PARAM', 'all') - self.pem4criteria.insert(self.pem4criteria.shape[1], 'PARAM', 'all') - df_pem = self.pem4criteria - df_params = self.params - self.site = config_obj['Site'] - self.scenario = config_obj['Scenario'] - if not os.path.isdir(out_dir): - os.makedirs(out_dir) - # reading criteria set in config - cod = float(config_obj['COD_criteria']) - nse = float(config_obj['NSE_criteria']) - p_bias = float(config_obj['PBAIS_criteria']) - - df_sen_summary = pd.DataFrame() - df_sens = pd.DataFrame() - self.metric = metric - maxp = int(config_obj['max_range']) - # inc = float(config_obj['increment']) - min_p = -maxp - deltas = np.linspace(min_p, maxp, 200) - - n_param = len(self.param_list) - n_step = len(deltas) - n_simulate = len(deltas) * (n_param + 1) - id_vec = np.arange(n_step, n_simulate, n_step) - for i in range(len(id_vec)): - df_pem.PARAM[id_vec[i]:id_vec[i] + n_step] = self.param_list[i] - df_params.PARAM[id_vec[i]:id_vec[i] + n_step] = self.param_list[i] - df_count = pd.DataFrame() - for metric in ['OF', 'NSE', 'PBIAS', 'COD']: - # calculate standard deviation - df_x = df_params.copy() - df_x = df_x[df_x.PARAM == 'all'] - df_x = df_x.drop(['RunId', 'PARAM'], axis=1) - df_y = df_pem[df_pem.PARAM == 'all'] - df_y = df_y[metric] - print(f'------Calculating total Standarized Regression Coefficient for {metric}----\n') - src_total = ap.standarizedRegressionCoefficientTotal(df_x, df_y, intercept=False) - df_param_mat = pd.DataFrame() - df_metric = pd.DataFrame() - df_src = pd.DataFrame() - df_pem_combined = df_pem[df_pem.PARAM == 'all'] - df_pem_criteria = df_pem_combined[ - (df_pem_combined.COD > cod) & (df_pem_combined.NSE > nse) & (df_pem_combined.PBIAS < p_bias)] - df_count1 = pd.DataFrame({'COD': len(df_pem_combined.COD >= cod), 'NSE': len(df_pem_combined.NSE >= nse), - 'PBIAS': len(df_pem_combined.PBIAS <= p_bias), 'TOTAL': df_pem_criteria.shape[0], - 'METRIC': metric}, index=['all']) - for j in range(len(self.param_list)): - param = self.param_list[j] - param_name = f'PARAM [{id_sensitive[j] - 70}]' - df_pm_j = df_pem[df_pem.PARAM == param] - df_pem_criteria = df_pm_j[(df_pm_j.COD > cod) & (df_pm_j.NSE > nse) & (df_pm_j.PBIAS < p_bias)] - df_count = pd.DataFrame({'COD': len(df_pm_j.COD >= cod), 'NSE': len(df_pm_j.NSE >= nse), - 'PBIAS': len(df_pm_j.PBIAS <= p_bias), 'TOTAL': df_pem_criteria.shape[0], - 'METRIC': metric}, index=[param_name]) - df_count = pd.concat([df_count, df_count], axis=0) - df_pa_j = df_params[df_params.PARAM == param] - df_param_individual = pd.DataFrame(df_pa_j[param].values) - df_param_mat = pd.concat([df_param_mat, df_param_individual], axis=1) - df_metric_individual = pd.DataFrame(df_pm_j[metric].values) - df_metric = pd.concat([df_metric, df_metric_individual], axis=1) - df_sen = pd.DataFrame({'PARAM': df_pa_j[self.param_list[j]], - 'METRIC': df_pm_j[metric]}) - df_sen.METRIC = ep.get_outlier_na(df_sen.METRIC) - df_sen.insert(1, 'PARAM_CHANGE', ep.find_change_percent(df_sen.PARAM.values, self.p_best[param])) - df_sen.insert(3, 'CHANGE_METRIC', ep.find_change_percent(df_sen.METRIC.values, self.pe_best[metric])) - df_sen.insert(4, 'ABSOLUTE_CHANGE', deltas) - df_sen.insert(5, 'PARAM_NAME', param) - df_sen.insert(6, 'NAME', param_name) - df_sen_non_nan = df_sen.dropna() - r2, p_value = ep.corr_test(df_sen_non_nan.PARAM_CHANGE, df_sen_non_nan.CHANGE_METRIC) - x = df_params[df_params.PARAM == param][param].values - y = df_pem[df_pem.PARAM == param][metric].values - si = ap.sensitivity_index(x, y) - df_j = pd.DataFrame({'Nmet': df_pem_criteria.shape[0], - 'minCOD': df_pem_criteria.COD.min(), - 'maxCOD': df_pem_criteria.COD.max(), - 'nCOD': len(df_pem_criteria.COD > cod), - 'minNSE': df_pem_criteria.NSE.min(), - 'maxNSE': df_pem_criteria.NSE.max(), - 'nNSE': len(df_pem_criteria.NSE > nse), - 'minPBIAS': df_pem_criteria.PBIAS.min(), - 'maxPBIAS': df_pem_criteria.PBIAS.max(), - 'nPBIAS': len(df_pem_criteria.PBIAS > p_bias), - 'MINp': df_pa_j[self.param_list[j]].min(), - 'MAXp': df_pa_j[self.param_list[j]].max(), - 'Rsquared': r2, - 'p-value': p_value, - 'SensitivityIndex': si}, - index=[param_name]) - # calculate individual standardized Regression coefficient - if j == 0: - print(f'------Calculating individual Standarized Regression Coefficient for {metric}----\n') - src = ap.standarizedRegressionCoefficient(df_pa_j[param].values, df_pm_j[metric].values, intercept=True) - df_src_ = pd.DataFrame({'SRC': src, 'Metric': metric}, index=[param_name]) - df_src = pd.concat([df_src, df_src_], axis=0) - df_j['PARAM'] = param - df_sen_summary = pd.concat([df_sen_summary, df_j], axis=0) - df_sens = pd.concat([df_sens, df_sen_non_nan], axis=0) - df_count = pd.concat([df_count, df_count1], axis=0) - df_param_mat.columns = self.param_list - df_param_all = pd.concat([df_params[df_params.PARAM == 'all'][self.param_list], df_param_mat], axis=0) - df_metric.columns = self.param_list - df_count.to_csv(f'{out_dir}/Param_within_{metric}.csv') - # prepare data for heat map - param_list = [] - for i, param in enumerate(self.param_list): - param_name = f'PARAM [{id_sensitive[i] - 70}]' - param_list.append(param_name) - - if metric == 'OF': - c_bar_title = f'Percent change in objective function value' - else: - c_bar_title = f'Percent change in {metric} value' - - df_plot_data = (df_metric - self.pe_best[metric]) * 100 / self.pe_best[metric] - df_plot_data.columns = param_list - df_plot_data.index = [round(item, 2) for item in deltas] - df_plot_data.to_csv(f'{out_dir}/Heatmap-data_{metric}.csv') - self.Heatmap_data = df_plot_data - fig, ax = plt.subplots(figsize=(8, 8), tight_layout=True) - sns.heatmap(df_plot_data, cmap='bwr', - vmin=df_plot_data.min().min(), - vmax=df_plot_data.max().max(), - cbar_kws={'label': c_bar_title}, ax=ax) - ax.set_ylabel('Percent change in parameters') - plt.savefig(f'{out_dir}/Heatmap_{metric}.png', dpi=600, bbox_inches="tight") - - # Sensitivity indices - df_src.insert(1, 'SRC_Total', src_total.values, True) - - y = df_pem[metric].values - list_bound = df_sen_summary[['MINp', 'MAXp']] - list_bound = list_bound.values.tolist() - print(f'------Calculating SOBOL and FAST indices {metric}--------- \n') - df_sobol_first, df_sobol_total = ap.getSOBOL_index(param_list, list_bound, y) # SOBOL indices - df_fast_first, df_fast_total = ap.get_FAST_index(param_list, list_bound, y) # FAST indices - - df_src_first = pd.DataFrame( - {'Sensitivity Index': df_src.SRC.values, 'Order': 'First', 'Method': 'SRC', 'PARAM': param_list}, - index=param_list) - df_src_total = pd.DataFrame( - {'Sensitivity Index': df_src.SRC_Total.values, 'Order': 'Total', 'Method': 'SRC', 'PARAM': param_list}, - index=param_list) - df_si_first = pd.concat([df_sobol_first, df_src_first, df_fast_first], axis=0) - df_si_total = pd.concat([df_sobol_total, df_src_total, df_fast_total], axis=0) - df_sobol = pd.concat([df_sobol_first, df_sobol_total], axis=0) - df_fast = pd.concat([df_fast_first, df_fast_total], axis=0) - df_SOBOL = pd.DataFrame({'First': df_sobol_first['Sensitivity Index'].values, - 'Total': df_sobol_total['Sensitivity Index'].values}, index=df_sobol_total.index) - df_FAST = pd.DataFrame({'First': df_fast_first['Sensitivity Index'].values, - 'Total': df_fast_total['Sensitivity Index'].values}, index=df_fast_total.index) - df_SRC = pd.DataFrame({'First': df_src_first['Sensitivity Index'].values, - 'Total': df_src_total['Sensitivity Index'].values}, index=df_src_total.index) - df_src = pd.concat([df_src_first, df_src_total], axis=0) - # plots - fig, axes = plt.subplots(1, 3, figsize=(12, 10), sharex=False, sharey=True, tight_layout=True) - sns.barplot(data=df_sobol, y='PARAM', x='Sensitivity Index', hue='Order', ci=None, ax=axes[0]) - # handles, labels = ax1.get_legend_handles_labels() - axes[0].grid(True) - axes[0].set_ylabel('Parameters') - axes[0].set_xscale('log') - axes[0].set_title('SOBOL index') - axes[0].get_legend().remove() - sns.barplot(data=df_fast, y='PARAM', x='Sensitivity Index', hue='Order', ci=None, ax=axes[1]) - # handles, labels = ax1.get_legend_handles_labels() - axes[1].grid(True) - axes[1].set_ylabel('Parameters') - axes[1].set_xscale('log') - axes[1].set_title('FAST Index') - axes[1].get_legend().remove() - sns.barplot(data=df_src, y='PARAM', x='Sensitivity Index', hue='Order', ci=None, ax=axes[2]) - # handles, labels = ax1.get_legend_handles_labels() - axes[2].grid(True) - axes[2].set_ylabel('Parameters') - axes[2].set_xscale('log') - axes[2].set_title('Standardized regression Coefficient') - axes[2].get_legend().remove() - axes[2].legend(loc='upper center', bbox_to_anchor=(0.5, -0.05), ncol=2) - plt.tight_layout() - plt.savefig(f'{out_dir}/Index_{metric}.png', dpi=600, bbox_inches="tight") - - fig, axes = plt.subplots(3, 1, figsize=(20, 12), sharex=False, sharey=False, tight_layout=True) - axes[0] = plot_2_bar(df_SOBOL, labels=['First', 'Total'], ax=axes[0]) - axes[1] = plot_2_bar(df_FAST, labels=['First', 'Total'], ax=axes[1]) - axes[2] = plot_2_bar(df_SRC, labels=['First', 'Total'], ax=axes[2]) - plt.savefig(f'{out_dir}/Index_{metric}_V1.png', dpi=600, bbox_inches="tight") - - df_param_all.to_csv(f'{out_dir}/Param_{metric}.csv') - df_si_first.to_csv(f'{out_dir}/Index_first_{metric}.csv') - df_si_total.to_csv(f'{out_dir}/Index_Total_{metric}.csv') - df_sen_summary.to_csv(f'{out_dir}/Stat_Summary_{metric}.csv') - - # import major output - print('------- importing results ------\n') - df_wyld = self.import_output(location='outlet', variable='WYLD') - df_wyld.to_csv(f'{out_dir}/Daily_WYLD.csv') - df_rus2 = self.import_output(location='outlet', variable='RUS2') - df_rus2.to_csv(f'{out_dir}/Daily_RUS2.csv') - df_ysd = self.import_output(location='outlet', variable='YSD') - df_ysd.to_csv(f'{out_dir}/Daily_YSD.csv') - df_lai = self.import_output(location='basin', variable='LAI') - df_lai.to_csv(f'{out_dir}/Daily_LAI.csv') - df_biom = self.import_output(location='basin', variable='BIOM') - df_biom.to_csv(f'{out_dir}/Daily_BIOM.csv') - df_prcp = self.import_output(location='basin', variable='PRCP') - df_prcp.to_csv(f'{out_dir}/Daily_PRCP.csv') - df_stl = self.import_output(location='basin', variable='STL') - df_stl.to_csv(f'{out_dir}/Daily_STL.csv') - df_std = self.import_output(location='basin', variable='STD') - df_std.to_csv(f'{out_dir}/Daily_STD.csv') - df_stdl = self.import_output(location='basin', variable='STDL') - df_stdl.to_csv(f'{out_dir}/Daily_STDL.csv') - df_et = self.import_output(location='basin', variable='ET') - df_et.to_csv(f'{out_dir}/Daily_ET.csv') - df_pet = self.import_output(location='basin', variable='PET') - df_pet.to_csv(f'{out_dir}/Daily_PET.csv') - df_dprk = self.import_output(location='basin', variable='DPRK') - df_dprk.to_csv(f'{out_dir}/Daily_DPRK.csv') - df_tp = self.import_output(location='basin', variable='TP') - df_tp.to_csv(f'{out_dir}/Daily_TP.csv') - df_tn = self.import_output(location='basin', variable='TN') - df_tn.to_csv(f'{out_dir}/Daily_TN.csv') - df_yldg = self.import_output(location='annual', variable='YLDG') - df_yldg.to_csv(f'{out_dir}/Annual_YLDG.csv') - df_yldf = self.import_output(location='annual', variable='YLDF') - df_yldf.to_csv(f'{out_dir}/Annual_YLDF.csv') - df_biom = self.import_output(location='annual', variable='BIOM') - df_biom.to_csv(f'{out_dir}/Annual_BIOM.csv') - df_ws = self.import_output(location='annual', variable='WS') - df_ws.to_csv(f'{out_dir}/Annual_WS.csv') - df_ns = self.import_output(location='annual', variable='NS') - df_ns.to_csv(f'{out_dir}/Annual_NS.csv') - df_ps = self.import_output(location='annual', variable='PS') - df_ps.to_csv(f'{out_dir}/Annual_PS.csv') - df_ts = self.import_output(location='annual', variable='TS') - df_ts.to_csv(f'{out_dir}/Annual_TS.csv') - - def get_range(self): - # import csv file containing range of parameters with recommended values for specific project - df_param_limit = pd.read_csv(self.file_limits, index_col=0, encoding="ISO-8859-1") - array_param_list = df_param_limit.iloc[1:, :].to_numpy() - array_param_list = array_param_list.astype(np.float) - mat_param_list = np.asmatrix(array_param_list) - mat_param_list = np.squeeze(np.asarray(mat_param_list)) - self.param_range = mat_param_list - return self - - def get_pe_files(self): - if self.attribute == 'WYLD': - file_pe = os.path.join(self.folder, 'Statistics_runoff.csv') - elif self.attribute == 'YSD': - file_pe = os.path.join(self.folder, 'Statistics_sediment.csv') - else: - file_pe = os.path.join(self.folder, f'Statistics_{self.attribute}.csv') - file_parameter = os.path.join(self.folder, 'APEXPARM.csv') - self.file_pem = file_pe - self.file_param = file_parameter - return self - - def get_params(self, id_sensitive, is_all): - df_param = pd.DataFrame() - df_params = pd.read_csv(self.file_param) - df_params.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - self.param_list = df_params.columns - self.params_all = df_params - if is_all: - self.params = df_params.copy() - else: - self.param_list = self.param_list[id_sensitive] - df_param = df_params[self.param_list] - df_param.insert(0, 'RunId', df_params.RunId) - self.p_best = df_param.iloc[0, :] - self.params = df_param.iloc[1:, ] - return self - - def read_pem(self, metric_list): - # read statistics - df_pem = ap.get_stats(self.file_pem) - df_pem_obs = df_pem[metric_list] - df_pem_obs.columns = ['RunId', 'COD', 'RMSE', 'NRMSE', 'NSE', 'PBIAS', 'OF'] - pe_best = df_pem_obs.iloc[0, :] - df_pem_obs = df_pem_obs.iloc[1:, ] - self.pem = df_pem - self.pem4criteria = df_pem_obs - self.pe_best = pe_best - return self - - def assign_dates(self, df_mod, n_warm, n_calib_year): - start_year = df_mod.Year.values[0] - cal_start = start_year + n_warm - cal_end = cal_start + n_calib_year - val_end = self.df_obs.Year[-1] - val_end_date = self.df_obs.Date[-1] - return start_year, cal_start, cal_end, val_end, val_end_date - - # noinspection PyBroadException - def import_output(self, location, variable): - print(f'\nExporting {variable} data') - stage_vec = [] - n_runs = self.pem4criteria.shape[0] - n_warm = int(self.config['warm_years']) - n_calib_year = int(self.config['calib_years']) - df_out = pd.DataFrame() - print_progress_bar(0, n_runs + 1, prefix='Progress', suffix='Complete', length=50, fill='█') - for i in range(1, n_runs + 1): - if location == 'outlet': - file_outlet = os.path.join(self.folder, f'daily_outlet_{i:07}.csv.csv') - data = pd.read_csv(file_outlet) - data.index = data.Date - data.index = pd.to_datetime(data.index) - data = data.drop(['Date', 'Y', 'M', 'D'], axis=1) - data.insert(0, 'Year', data.index.year, True) - try: - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - except: - continue - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end_date] - elif location == 'basin': - file_basin = os.path.join(self.folder, f'daily_basin_{i:07}.csv.csv') - data = pd.read_csv(file_basin) - data.index = data.Date - data.index = pd.to_datetime(data.index) - data = data.drop(['Date', 'Y', 'M', 'D'], axis=1) - data.insert(0, 'Year', data.index.year, True) - try: - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - except: - continue - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end_date] - else: - file_annual = os.path.join(self.folder, f'annual_{i:07}.csv.csv') - data = pd.read_csv(file_annual) - data.index = data.YR - data = data.drop(['Unnamed: 0', 'YR'], axis=1) - data.insert(0, 'Year', data.index, True) - try: - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - except: - continue - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end] - df_data_cal.insert(df_data_cal.shape[1], 'Stage', 'Calibration', True) - df_data_val.insert(df_data_val.shape[1], 'Stage', 'Validation', True) - df_data = pd.concat([df_data_cal, df_data_val], axis=0) - stage_vec = df_data.Stage.values - df_data = pd.DataFrame(df_data[variable]) - df_data.columns = [f'{i}'] - df_out = pd.concat([df_out, df_data], axis=1) - print_progress_bar(i, n_runs + 1, prefix=f'Progress: {i}', suffix='Complete', length=50, fill='█') - df_out.insert(0, 'Stage', stage_vec, True) - return df_out diff --git a/pyAPEX/pyUAAPEX.py b/pyAPEX/pyUAAPEX.py deleted file mode 100644 index f967f99..0000000 --- a/pyAPEX/pyUAAPEX.py +++ /dev/null @@ -1,263 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Sat Dec 31 09:00:28 2022 - -@author: Mahesh.Maskey -""" - -import os -import pandas as pd -from Utility.apex_utility import print_progress_bar -from Utility.apex_utility import split_data -from Utility.pyAPEXpost import pyAPEXpost as ap - -print('\014') - - -# noinspection PyTypeChecker,PyAttributeOutsideInit,PyUnboundLocalVariable,PyBroadException -class unaAPEX: - def __init__(self, src_dir, config, scale, obs_attribute, mod_attribute, out_dir, is_full=True): - - """ - - - - """ - self.config = config - self.obs_attribute = obs_attribute - self.mod_attribute = mod_attribute - self.site = config['Site'] - self.scenario = config['Scenario'] - self.src_dir = src_dir - self.criteria = [float(config['COD_criteria']), float(config['NSE_criteria']), float(config['PBAIS_criteria'])] - self.metrics = ['OF', 'NSE', 'PBIAS', 'COD'] - self.cal_dir = config['dir_calibration'] - self.una_dir = config['dir_uncertainty'] - self.out_dir = out_dir - self.get_pe_files() - self.get_params() - self.read_pem(scale) - self.df_obs = ap.get_measure(data_dir='Program', file_name='calibration_data.csv') - if not os.path.isdir(self.out_dir): - os.makedirs(self.out_dir) - self.outlet_daily = self.import_best(location='daily_outlet', crops=None) - self.basin_daily = self.import_best(location='daily_basin', crops=None) - self.annual = self.import_best(location='annual', crops=None) - self.compile_stats_param(scale) - - # import major output - print('-------Importing model results-------------------\n') - df_WYLD = self.import_output(location='outlet', variable='WYLD', is_full=is_full) - df_WYLD.to_csv(f'{out_dir}/Daily_WYLD.csv') - df_RUS2 = self.import_output(location='outlet', variable='RUS2', is_full=is_full) - df_RUS2.to_csv(f'{out_dir}/Daily_RUS2.csv') - df_YSD = self.import_output(location='outlet', variable='YSD', is_full=is_full) - df_YSD.to_csv(f'{out_dir}/Daily_YSD.csv') - df_LAI = self.import_output(location='basin', variable='LAI', is_full=is_full) - df_LAI.to_csv(f'{out_dir}/Daily_LAI.csv') - df_BIOM = self.import_output(location='basin', variable='BIOM', is_full=is_full) - df_BIOM.to_csv(f'{out_dir}/Daily_BIOM.csv') - df_PRCP = self.import_output(location='basin', variable='PRCP', is_full=is_full) - df_PRCP.to_csv(f'{out_dir}/Daily_PRCP.csv') - df_STL = self.import_output(location='basin', variable='STL', is_full=is_full) - df_STL.to_csv(f'{out_dir}/Daily_STL.csv') - df_STD = self.import_output(location='basin', variable='STD', is_full=is_full) - df_STD.to_csv(f'{out_dir}/Daily_STD.csv') - df_STDL = self.import_output(location='basin', variable='STDL', is_full=is_full) - df_STDL.to_csv(f'{out_dir}/Daily_STDL.csv') - df_ET = self.import_output(location='basin', variable='ET', is_full=is_full) - df_ET.to_csv(f'{out_dir}/Daily_ET.csv') - df_PET = self.import_output(location='basin', variable='PET', is_full=is_full) - df_PET.to_csv(f'{out_dir}/Daily_PET.csv') - df_DPRK = self.import_output(location='basin', variable='DPRK', is_full=is_full) - df_DPRK.to_csv(f'{out_dir}/Daily_DPRK.csv') - dfTP = self.import_output(location='basin', variable='TP', is_full=is_full) - dfTP.to_csv(f'{out_dir}/Daily_TP.csv') - df_TN = self.import_output(location='basin', variable='TN', is_full=is_full) - df_TN.to_csv(f'{out_dir}/Daily_TN.csv') - df_YLDG = self.import_output(location='annual', variable='YLDG', is_full=is_full) - df_YLDG.to_csv(f'{out_dir}/Annual_YLDG.csv') - df_YLDF = self.import_output(location='annual', variable='YLDF', is_full=is_full) - df_YLDF.to_csv(f'{out_dir}/Annual_YLDF.csv') - df_BIOM = self.import_output(location='annual', variable='BIOM', is_full=is_full) - df_BIOM.to_csv(f'{out_dir}/Annual_BIOM.csv') - df_WS = self.import_output(location='annual', variable='WS', is_full=is_full) - df_WS.to_csv(f'{out_dir}/Annual_WS.csv') - df_NS = self.import_output(location='annual', variable='NS', is_full=is_full) - df_NS.to_csv(f'{out_dir}/Annual_NS.csv') - df_PS = self.import_output(location='annual', variable='PS', is_full=is_full) - df_PS.to_csv(f'{out_dir}/Annual_PS.csv') - df_TS = self.import_output(location='annual', variable='TS', is_full=is_full) - df_TS.to_csv(f'{out_dir}/Annual_TS.csv') - - def get_pe_files(self): - if self.mod_attribute == 'WYLD': - file_pe = os.path.join(self.una_dir, 'Statistics_runoff.csv') - elif self.mod_attribute == 'YSD': - file_pe = os.path.join(self.una_dir, 'Statistics_sediment.csv') - else: - file_pe = os.path.join(self.una_dir, f'Statistics_{self.attribute}.csv') - file_parameter = os.path.join(self.una_dir, 'APEXPARM.csv') - self.file_pem = file_pe - self.file_param = file_parameter - return self - - def get_params(self): - df_params = pd.read_csv(self.file_param) - df_params.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - self.param_list = df_params.columns - self.params_all = df_params - self.p_best = df_params.iloc[0, :] - self.params = df_params.iloc[1:, ] - return self - - # noinspection PyUnboundLocalVariable - def read_pem(self, scale): - # read statistics - df_pem = pd.read_csv(self.file_pem) - df_pem.rename(columns={'Unnamed: 0': 'RunId'}, inplace=True) - id_best = int(df_pem.iloc[0, 0]) - df_pem = df_pem.iloc[1:, :] - pe_best = df_pem.iloc[0, :] - if scale == 'daily': - COD, NSE, PBIAS = 'CODDC', 'NSEDC', 'PBIASDC' - elif scale == 'monthly': - COD, NSE, PBIAS = 'CODMC', 'NSEMC', 'PBIASMC' - elif scale == 'yearly': - COD, NSE, PBIAS = 'CODYC', 'LYNSEY', 'PBIASYC' - df_pem_within = df_pem[ - (df_pem[COD] > self.criteria[0]) & (df_pem[NSE] > self.criteria[1]) & (df_pem[PBIAS] < self.criteria[2])] - self.pem = df_pem - self.id_best = id_best - self.pe_best = pe_best - self.pem4criteria = df_pem_within - return self - - def import_best(self, location, crops): - run_id = self.id_best - if crops is None: - file_read = f'{location}_{run_id:07}.csv' - file_path = os.path.join(self.cal_dir, self.obs_attribute, file_read) - df_data = pd.read_csv(file_path) - if location == 'annual': - df_data.index = df_data.YR - df_data = df_data.drop(['Unnamed: 0', 'YR'], axis=1) - else: - df_data.Date = pd.to_datetime(df_data.Date) - df_data.index = df_data.Date - df_data = df_data.drop('Date', axis=1) - df_data.to_csv(os.path.join(self.out_dir, f'{location}_best.csv')) - df_data.to_csv(os.path.join(self.out_dir, file_read)) - return df_data - else: - file_read = f'{location}_{run_id:07}.csv' - file_path = os.path.join(self.cal_dir, self.obs_attribute, file_read) - df_data = pd.read_csv(file_path) - if location == 'annual': - df_data.index = df_data.YR - df_data = df_data.drop(['Unnamed: 0', 'YR'], axis=1) - else: - df_data.Date = pd.to_datetime(df_data.Date) - df_data.index = df_data.Date - df_data = df_data.drop('Date', axis=1) - df_data.to_csv(os.path.join(self.out_dir, f'{location}_best.csv')) - df_data.to_csv(os.path.join(self.out_dir, file_read)) - df_list = [df_data] - for crop in crops: - file_read = f'{location}_{run_id:07}_{crop}.csv' - file_path = os.path.join(self.cal_dir, self.obs_attribute, file_read) - df_data = pd.read_csv(file_path) - df_data.index = df_data.Date - df_data = df_data.drop('Date', axis=1) - df_data.to_csv(os.path.join(self.out_dir, f'{location}_best.csv')) - df_data.to_csv(os.path.join(self.out_dir, file_read)) - df_list.append(df_data) - return df_list - - def compile_stats_param(self, scale): - out_compile = ap.compile_stats(self.pem, self.criteria, scale) - self.df_stats_daily = out_compile[0] - self.df_daily_daily = out_compile[1] - self.df_daily_monthly = out_compile[2] - self.df_daily_yearly = out_compile[3] - self.ids_best_daily = self.df_stats_daily.RunId.values.astype(int) - self.df_best_params = pd.DataFrame() - for i in range(len(self.ids_best_daily)): - df = self.params[self.params['RunId'] == self.ids_best_daily[i]] - self.df_best_params = pd.concat([self.df_best_params, df], axis=0) - return self - - def assign_dates(self, df_mod, n_warm, n_calib_year): - start_year = df_mod.Year.values[0] - cal_start = start_year + n_warm - cal_end = cal_start + n_calib_year - val_end = self.df_obs.Year[-1] - val_end_date = self.df_obs.Date[-1] - return start_year, cal_start, cal_end, val_end, val_end_date - - def import_output(self, location, variable, is_full): - print(f'\nExporting {variable} data') - if is_full: - n_runs = self.params.shape[0] - else: - n_runs = self.df_best_params.shape[0] - n_warm = int(self.config['warm_years']) - n_calib_year = int(self.config['calib_years']) - df_out = pd.DataFrame() - print_progress_bar(0, n_runs, prefix='Progress', suffix='Complete', length=50, fill='█') - for k in range(n_runs): - if is_full: - i = k + 1 - else: - i = self.ids_best_daily[k] - try: - if location == 'outlet': - file_outlet = os.path.join(self.una_dir, f'daily_outlet_{i:07}.csv.csv') - data = pd.read_csv(file_outlet) - data.index = data.Date - data.index = pd.to_datetime(data.index) - data = data.drop(['Date', 'Y', 'M', 'D'], axis=1) - data.insert(0, 'Year', data.index.year, True) - try: - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - except: - continue - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end_date] - elif location == 'basin': - file_basin = os.path.join(self.una_dir, f'daily_basin_{i:07}.csv.csv') - data = pd.read_csv(file_basin) - data.index = data.Date - data.index = pd.to_datetime(data.index) - data = data.drop(['Date', 'Y', 'M', 'D'], axis=1) - data.insert(0, 'Year', data.index.year, True) - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end_date] - else: - file_annual = os.path.join(self.una_dir, f'annual_{i:07}.csv.csv') - data = pd.read_csv(file_annual) - data.index = data.YR - data = data.drop(['Unnamed: 0', 'YR'], axis=1) - data.insert(0, 'Year', data.index, True) - start_year, cal_start, cal_end, val_end, val_end_date = self.assign_dates(data, n_warm, - n_calib_year) - df_data_cal, df_data_val = split_data(start_year, data, n_warm, n_calib_year) - df_data_val = df_data_val[df_data_val.index <= val_end] - df_data_cal.insert(df_data_cal.shape[1], 'Stage', 'Calibration', True) - df_data_val.insert(df_data_val.shape[1], 'Stage', 'Validation', True) - except Exception as e: - print(e) - print(f'error occurs in simulation {i}') - continue - df_data = pd.concat([df_data_cal, df_data_val], axis=0) - stage_vec = df_data.Stage.values - df_data = pd.DataFrame(df_data[variable]) - df_data.columns = [f'{i}'] - df_out = pd.concat([df_out, df_data], axis=1) - - print_progress_bar(k, n_runs, prefix=f'Progress: {i}', suffix='Complete', length=50, fill='█') - df_out.insert(0, 'Stage', stage_vec, True) - return df_out diff --git a/pyAPEX/sensitivity_analysis.py b/pyAPEX/sensitivity_analysis.py deleted file mode 100644 index a4434d2..0000000 --- a/pyAPEX/sensitivity_analysis.py +++ /dev/null @@ -1,24 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Dec 21 07:42:35 2022 - -@author: Mahesh.Maskey -""" - - -import os -import os.path -from pathlib import Path -from configobj import ConfigObj -from pySAAPEX import senanaAPEX -print('\014') - - -src_dir = Path(os.path.dirname(os.path.realpath(__file__))) -config = ConfigObj(str(src_dir / 'runtime.ini')) -site = 'Farm_1' -scenario = 'grazing' -config['Site'] = site -config['Scenario'] = scenario -out_dir = f'../../../post_analysis/sensitivity_analysis/{site}/{scenario}' -sen_obj = senanaAPEX(src_dir, config, out_dir, attribute='runoff', metric='OF') \ No newline at end of file diff --git a/pyAPEX/task_worker.py b/pyAPEX/task_worker.py deleted file mode 100644 index 4feef96..0000000 --- a/pyAPEX/task_worker.py +++ /dev/null @@ -1,115 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Manages APEX model runs in the context of a Slurm task that is running -concurrently with (potentially many) other model run tasks. Expects to receive -an ID range that uniquely identifies these model runs from the model runs in -all other Slurm tasks. - -@author: Brian Stucky (USDA), Mahesh.Maskey -python task_worker.py --ntasks=1 --nsims=20000 --taskidmin=1 --taskid=1 --simidstart=1 --outputdir=path/to/output/ --winepath=None --id_mode=0 -""" -import os -import os.path -from pathlib import Path -import shutil -import argparse -from pyAPEXin import inAPEX -# from pyAPEXSCU import simAPEX -from pyAPEXSCU import simAPEX -from configobj import ConfigObj -import numpy as np - -print('\014') -parser = argparse.ArgumentParser() -parser.add_argument( - '--ntasks', type=int, required=True, - help='The total number of worker tasks in the modeling run.' -) -parser.add_argument( - '--nsims', type=int, required=True, - help='The total number of simulations to run across all worker tasks.' -) -parser.add_argument( - '--taskidmin', type=int, required=True, help='The Start ID of all worker tasks.' -) -parser.add_argument( - '--taskid', type=int, required=True, help='The ID of this worker task.' -) -parser.add_argument( - '--simidstart', type=int, required=True, - help='The start simulation id to run.' -) -parser.add_argument( - '--outputdir', type=str, required=True, - help='The model run and output location.' -) -parser.add_argument( - '--winepath', type=str, required=True, - help='The location of an executable Wine container image.' -) -parser.add_argument( - '--id_mode', type=int, required=True, - help='Type of simulation. : 0: Calibration, 1:Sensitivity, 2: Uncertainty' -) - -args = parser.parse_args() - -# Get the location of the simulation software and output dir. -src_dir = Path(os.path.dirname(os.path.realpath(__file__))) -outputdir = Path(args.outputdir) - -# Calculate the simulation ID range for this worker task. -if args.nsims < args.ntasks: - raise Exception( - f'Incompatible simulation and task counts: {args.nsims}, {args.ntasks}.' - ) - -if args.ntasks > 999: - raise Exception( - f'Task count too large.' - ) -#Select the model mode 0: calibration, 1: sensitivity, 2: uncertainty -#id_mode = 2 -# evenly distribute the sim ids to each task -sim_per_task = np.ceil(args.nsims / args.ntasks).astype(int) -id_mode = args.id_mode -z = np.zeros((sim_per_task, args.ntasks)) # zero pads -sim_ids = np.arange(args.simidstart, args.simidstart + args.nsims) -sim_ids.resize(z.shape, refcheck=False) -# sim_ids = np.zeros(z.shape) - -id = args.simidstart -for c in range(args.ntasks): - for r in range(sim_per_task): - if sim_ids[r, c] > 0: - sim_ids[r, c] = id - id += 1 - -sim_id_min = sim_ids[0, args.taskid - args.taskidmin] - 1 -sim_id_max = sim_ids[:, args.taskid - args.taskidmin].max().astype(int) - -# Load and adjust the configuration settings. -config = ConfigObj(str(src_dir / 'runtime.ini')) -model_mode = ['calibration', 'sensitivity', 'uncertainty'] -config['n_start'] = sim_id_min -config['n_simulation'] = sim_id_max - -# Make a unique copy of the simulation code for this task, if needed. -if id_mode == 0: - task_progdir = outputdir / f'Program_{args.taskid:03}' - if not(task_progdir.is_dir()): - shutil.copytree(src_dir / 'Program', task_progdir) -else: - task_progdir = src_dir/'Program' - -# Run the simulations. -# isall_try = True -isall_try = False -curr_directory = os.getcwd() -os.chdir(task_progdir) -input_APEX = inAPEX(config, task_progdir) -sim_APEX = simAPEX( - config, src_dir, args.winepath, input_APEX, model_mode[id_mode], - scale='daily', isall=isall_try -) -os.chdir(curr_directory) diff --git a/pyAPEX/uncertainty_analysis.py b/pyAPEX/uncertainty_analysis.py deleted file mode 100644 index 8c9226c..0000000 --- a/pyAPEX/uncertainty_analysis.py +++ /dev/null @@ -1,25 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Dec 21 07:42:35 2022 - -@author: Mahesh.Maskey -""" - -import os.path -from pathlib import Path -from configobj import ConfigObj -from pyUAAPEX import unaAPEX - -print('\014') - -src_dir = Path(os.path.dirname(os.path.realpath(__file__))) -config = ConfigObj(str(src_dir / 'runtime.ini')) -site = 'Farm_1' -scenario = 'grazing' -config['Site'] = site -config['Scenario'] = scenario -obs_attribute = 'runoff' -out_dir = f'../../../post_analysis/Uncertainty_analysis/{site}/{scenario}/{obs_attribute}' -scale = 'daily' -mod_attribute = 'WYLD' -una_obj = unaAPEX(src_dir, config, scale, obs_attribute, mod_attribute, out_dir)