diff --git a/Fig 8.png b/Fig 8.png index 5c47846..d2e13ef 100755 Binary files a/Fig 8.png and b/Fig 8.png differ diff --git a/Figure_8.m b/Figure_8.m deleted file mode 100755 index 937380c..0000000 --- a/Figure_8.m +++ /dev/null @@ -1,408 +0,0 @@ - -clear all; clc - -% view the CONUS map to make sure the gscd conus regression works - -lat = load('lat_2d'); -lon = load('lon_2d'); -data = readtable('./CONUS_clustering_id.csv'); -data = sortrows(data,[4,5]); -cluster_id = data.clustering_7; -class_maxtri = nan(224, 464); -lat_index = data.Lat_index; -lon_index = data.Lon_index; - -sim = readtable('../gscd_conus_regression/qmean.csv'); -sim = sortrows(sim,[2,3]); -sim = sim.V8; -class_maxtri = nan(224, 464); -for i = 1:length(cluster_id) - class_maxtri(lat_index(i),lon_index(i)) = sim(i); -end - -subplot(2,1,1) -pcolor(lon, lat, class_maxtri); hold on; shading flat; -caxis([0,5]) - - - - -obs = readtable('../gscd_CONUS_land_cell.csv'); -obs = sortrows(obs, [2,3]); -obs = obs.qmean_mm_yr/365; -class_maxtri = nan(224, 464); -for i = 1:length(cluster_id) - class_maxtri(lat_index(i),lon_index(i)) = obs(i); -end - -subplot(2,1,2) -pcolor(lon, lat, class_maxtri); hold on; shading flat; -caxis([0,5]) - - - - - - -%% find the CAMELS basin parameter constrain - -clear all; clc; - -temp = readtable('../../ensemble_sim/ensemble_q10'); -out_regional = nan(464, 7); %1-gauge id, 2-lat, 3-lon, 4-cluster id, 5-q10_mm_d, 6-q90_mm_d, 7-qmean_mm_d -out_regional(:,1) = temp.Var1; -out_regional(:,2) = temp.Var2; -out_regional(:,3) = temp.Var3; -out_regional(:,4) = temp.Var4; -out_regional = sortrows(out_regional, 1); -temp = sortrows(temp, 1); - -para_qmean = load('camel_regional_id_1000_qmean.txt'); -par1000 = load('./sucessful_par_id_1000'); % id 1-1500 -par1307 = load('./sucessful_par_id_1307'); % id 1-1500 -index1000 = ismember(par1307, par1000); - - -out_para_basin_r_p = nan(464, 1000); %id 1-1500 -temp = readtable('../../ensemble_sim/ensemble_qmean'); -temp = sortrows(temp, 1); -for i = 1:464 - temp1 = table2array(temp(i, 5:end)); - temp1 = temp1(index1000); - temp2 = nan(10,1); - for j = 1:7 - if out_regional(i,4)==j - par_id = para_qmean(j,2:11); - for k = 1:10 - out_para_basin_r_p(i,k) = par1000(par_id(k)); - end - end - end -end - - - - - - -temp = readtable('../gscd_camel_regression/qmean.csv'); -var(:,9) = temp.V7; -var(:,10) = temp.V8; -var = sortrows(var,1); - -temp = readtable('../../ensemble_sim/ensemble_q10'); -temp = sortrows(temp,1); -cluster_id = table2array(temp(:,4)); - -out_para_basin_hybrid = nan(464, 1000); %id 1-1500 -temp = readtable('../../ensemble_sim/ensemble_qmean'); -temp = sortrows(temp,1); -for i = 1:464 - temp1 = table2array(temp(i, 5:end)); - temp1 = temp1(index1000); - low = var(i,9); - up = var(i,10); - indices = find(temp1 > low & temp1 < up); % 1- 350 - - % 1) use full constrain members - if length(indices)>=10 - temp2 = nan(length(indices), 1); - for j = 1:length(indices) - out_para_basin_hybrid(i,j) = par1000(indices(j)); - end - end - - % 2) if no constrain, use the top 10 members - if length(indices)==0 - temp2 = nan(10,1); - - for k = 1:7 - if cluster_id(i)==k - par_id = para_qmean(k,2:21); - for j = 1:10 - out_para_basin_hybrid(i,j) = par1000(par_id(j)); - end - end - end - end - - % 3) mixed with constain but no 10 parameters - if length(indices)>0 && length(indices)<10 - temp2 = nan(10, 1); - for j = 1:length(indices) - temp2(j,1) = temp1(indices(j)); - out_para_basin_hybrid(i,j) = par1000(indices(j)); - end - - need_num = 10 - length(indices); - count1 = length(indices)+1; - - - for k = 1:7 - if cluster_id(i)==k - par_id = para_qmean(k,2:21); - - for j = 1:20 - if ismember(par_id(j), indices) - continue - else - out_para_basin_hybrid(i,count1) = par1000(par_id(j)); -% temp2(count1,1) = temp1(par_id(j)); - count1 = count1 + 1; - if count1 == 11 - break - end - end - end - end - end - - end -end - - - - - - - - - - - - - - - -%% find GSCD regression constrain - -para = load('../../conus_ensemble_sim/parameter_ensemble_LHS700'); -para_id_350 = load('../../conus_ensemble_sim/succ_id_350'); % 1-700 -para_id_400 = load('../../conus_ensemble_sim/succ_id_400'); % 1-700 - -ensemble_qmean_400 = load('../../conus_ensemble_sim/Vol_Annual'); -ensemble_qmean_350 = nan(350, 50629); - -count = 1; -for i = 1:400 - isValueInArray = ismember(para_id_400(i), para_id_350); - if isValueInArray - ensemble_qmean_350(count,:) = ensemble_qmean_400(i,:); - count = count + 1; - end -end -clear ensemble_qmean_400; - -para_qmean = load('conus_regional_id_350_qmean.txt'); - -US = load('us_coor.txt'); -SL = load('sl_coor.txt'); -scatter_size = 20; label_size = 10; legend_size = 10; tick_size = 10; colorbar_size= 10;title_size = 10; line_width = 1; text_size = 10; - - -lat = load('lat_2d'); -lon = load('lon_2d'); -data = readtable('./CONUS_clustering_id.csv'); -data = sortrows(data,[4,5]); - -cluster_id = data.clustering_7; - -class_maxtri = nan(224, 464); -lat_index = data.Lat_index; -lon_index = data.Lon_index; - - - - -% find GSCD regression -temp = readtable('../gscd_conus_regression/qmean.csv'); -var = nan(50629, 4); %1-lat, 2-lon, 3-low, 4-up V4 V5 (95 CI), V6 V7 (95 PI) -var(:,1) = temp.V2; -var(:,2) = temp.V3; -var(:,3) = temp.V6; -var(:,4) = temp.V7; -var = sortrows(var,[1,2]); - -out_para_regression = nan(50629, 350); % para id from 1-700 -out_sim_regression = nan(50629, 350); - -for i = 1:50629 - low = var(i,3); - up = var(i,4); - indices = find(ensemble_qmean_350(:,i) > low & ensemble_qmean_350(:,i) < up); % 1- 350 - if length(indices)>=10 - for j = 1:length(indices) - out_sim_regression(i,j) = ensemble_qmean_350(indices(j),i); - out_para_regression(i,j) = para_id_350(indices(j)); - end - end - - if length(indices)==0 - temp2 = nan(10,1); - - for k = 1:7 - if cluster_id(i)==k - par_id = para_qmean(k,2:21); - for j = 1:10 - out_para_regression(i,j) = para_id_350(par_id(j)); - end - end - end - end - - - if length(indices)>0 && length(indices)<10 - temp2 = nan(10, 1); - for j = 1:length(indices) - out_para_regression(i,j) = para_id_350(indices(j)); - end - - need_num = 10 - length(indices); - count1 = length(indices)+1; - - - for k = 1:7 - if cluster_id(i)==k - par_id = para_qmean(k,2:21); - - for j = 1:20 - if ismember(par_id(j), indices) - continue - else - out_para_regression(i,count1) = para_id_350(par_id(j)); -% temp2(count1,1) = temp1(par_id(j)); - count1 = count1 + 1; - if count1 == 11 - break - end - end - end - end - end - - end - - -end - - - - - - - - - - - - -%% choose parameter - -US = load('us_coor.txt'); -SL = load('sl_coor.txt'); -scatter_size = 20; label_size = 11; legend_size = 11; tick_size = 11; colorbar_size= 11;title_size = 10; line_width = 1; text_size = 9; - -para_index = 1; - -cmin = 0.02; -cmax = 5; - -default_value = 0.5; - -figure; - -para = load('../../parameter_ensemble_LHS1500'); - -para_uniform = para(par1000,para_index); - -clc - - - - - -% c1 (1144000), c6 (05507600) . c5(08195000), c7 (02363000) c4 (09312600) -gauge_id = 1144000; -metric_daily = nan(464,1); -for i = 1:464 - - if out_regional(i,1)==gauge_id - - - - - temp = out_para_basin_hybrid(i,:); temp = temp(~isnan(temp)); - temp1 = nan(length(temp),1); - - data_pdf = nan(size(temp1,1), 15); - - for k = 1:15 - for j = 1:length(temp) - temp1(j,1) = para(temp(j),k); - end - - - - data_pdf(:,k) = temp1; - - end - - end -end - -% normalize the value -min_p = [0.02, 1, 0.0005, 0.01, 0.8, 0.1, 0.1, 0.8, 180, 0.1, 1.4, 0.05, 0.01, 10, 0.5]; -max_p = [5, 2, 0.07, 0.02, 1.2, 5, 5, 1.2, 220, 0.4, 9.5, 2, 0.5, 60, 1]; - -n_data = data_pdf; -for i = 1:15 - n_data(:,i) = (data_pdf(:,i) - min_p(i))/(max_p(i)-min_p(i)); -end - - - -for i = 1:size(n_data,1) - - plot(1:15, n_data(i,:), 'LineWidth', 0.5); hold on; - xlim([1,15]) - ylim([0,1]) - ylabel('Normalized Range', 'FontSize', label_size); - xlabel('Minimum Values', 'FontSize', label_size); - - if i == 40 - for j = 2:14 - xline(j, 'k', 'LineWidth', 2); - end - end -end - -xticks(1:15); -set(gca(),'XTickLabel',{sprintf('fff\\newline0.02'), sprintf('N_{bf}\\newline1'), sprintf('K_{bf}\\newline0.0005'), sprintf('S_{y}\\newline0.01'),... - sprintf('B\\newline0.8'), sprintf('ψ_{sat}\\newline0.1'), sprintf('k_{sat}\\newline0.1'), sprintf('Ɵ_{sat}\\newline0.8'), ... - sprintf('N_{melt}\\newline180'), sprintf('k_{acc}\\newline0.1'), sprintf('p_{sno}\\newline1.4'), sprintf('p_{lip}\\newline0.05'),... - sprintf('f_{wet}\\newline0.01'), sprintf('d_{max}\\newline10'), sprintf('Ɵ_{ini}\\newline0.5')}); - -ax1 = gca; -ax2=axes('Position',ax1.Position,'XAxisLocation','top','xlim',[1 15], 'ylim',[0 1], 'YTick', [], 'Color', 'none'); -xticks(ax2, 1:1:15); -% xlabel('Minimum Values', 'FontSize', label_size); -set(ax2,'XTickLabel',{'5', '2', '0.07', '0.02', '1.2', '5', '5', '1.2', '220', '0.4', '9.5', '2', '0.5', '60', '1'}); -xlabel(ax2, 'Maximum Values', 'FontSize', label_size); - - - - - - - -%% output the plot - -fig = gcf; -fig.PaperUnits = 'inches'; -% fig.PaperPosition = [0 0 10 5]; - -fig.PaperPositionMode = 'auto'; -fig.PaperPosition = [0 0 10 6]; -print('./figure', '-dpng', '-r300') - - diff --git a/README.md b/README.md index d0b2564..acfe18c 100644 --- a/README.md +++ b/README.md @@ -46,15 +46,17 @@ Clone the CLM5 repository to set up the CLM5 model, you will need to download th ## Reproduce My Figures | Figure Numbers | Script Name | Description | Figure | |:--------------:|:-----------:|:-----------:|:------:| -| 3 | Figure_3.m | Describe the regional mean daily FDC | | -| 4 | Figure_4.m | Describe the relative bias using 3 regionalization methods | | -| 5 | Figure_5.m | Describe the Qmean PDF of the behavioral parameter for 7 sites | | -| 6 | Figure_6.m | Describe the default and behavior parameters over the CONUS | | -| 7 | Figure_7.m | Describe the behavioral parameters for 7 sites | | -| 8 | Figure_8.m | Describe the 15 behavioral parameters for one site | | -| 9 | Figure_9.m | Describe the ensemble daily FDC for one site using one, two, three constraints | | -| 10 | Figure_10.m | Describe the Qmean prediction over the CONUS | | -| 11 | Figure_11.m | Describe the Q10 prediction over the CONUS | | -| 12 | Figure_12.m | Describe the Q90 prediction over the CONUS | | +| 1 | | Regionalization strategy | | +| 2 | | CAMELS basin and grid cell clustering | | +| 3 | Figure_3.m | Describe the regional mean daily FDC | | +| 4 | Figure_4.m | Describe the relative bias using 3 regionalization methods | | +| 5 | Figure_5.m | Describe the Qmean PDF of the behavioral parameter for 7 sites | | +| 6 | Figure_6.m | Describe the default and behavior parameters over the CONUS | | +| 7 | Figure_7.m | Describe the behavioral parameters for 7 sites | | +| 8 | figure\_8\_plot\_james\_params.py, figure\_8\_plot\_parallel.py | Describe the 15 behavioral parameters for one site | | +| 9 | Figure_9.m | Describe the ensemble daily FDC for one site using one, two, three constraints | | +| 10 | Figure_10.m | Describe the Qmean prediction over the CONUS | | +| 11 | Figure_11.m | Describe the Q10 prediction over the CONUS | | +| 12 | Figure_12.m | Describe the Q90 prediction over the CONUS | | | S1 | Figure_S1.m | Describe the Q10 PDF of the behavioral parameter for 7 sites | | | S2 | Figure_S2.m | Describe the Q90 PDF of the behavioral parameter for 7 sites | | diff --git a/figure_8_plot_james_params.py b/figure_8_plot_james_params.py new file mode 100644 index 0000000..a039cdc --- /dev/null +++ b/figure_8_plot_james_params.py @@ -0,0 +1,28 @@ +import numpy as np +import pandas as pd +from plot_parallel import custom_parallel_coordinates + +axis_labels = [r"$fff$", r"$N_{bf}$", r"$K_{bf}$", r"$S_{y}$", r"$B$", r"$k_{sat}$", r"$\theta_{sat}$", + r"$N_{melt}$", r"$k_{acc}$", r"$p_{sno}$", r"$p_{lip}$", r"$N_{bf}$", r"$f_{wet}$", r"$d_{max}$", r"$\theta_{ini}$"] + +util_names = ['OWASA', 'Durham', 'Cary', 'Raleigh', 'Chatham', 'Pittsboro', 'Regional'] +objs = ['REL', 'RF', 'NPC', 'PFC', 'WCC', 'UC'] +dvs = ['RT', 'TT', 'InfT'] + +# set parallel plot function parameters +fontsize = 10 +figsize = (10, 4) + +data_file = 'data.csv' + +data_df = pd.read_csv(data_file, index_col=None) +num_axis = data_df.shape[1] +axis_abbrevs = data_df.columns +colors = ['#DC851F', '#48A9A6','#355544'] +fig_filename = f'parallel_params.pdf' +fig_title = f'Behavioral hydrological parameters (C1-Northeast Basin)' +custom_parallel_coordinates(data_df, columns_axes=axis_abbrevs, axis_labels=axis_labels, + zorder_by=0, ideal_direction='upwards', zorder_direction='ascending', + alpha_base=0.6, lw_base=1.5, fontsize=fontsize, figsize=figsize, + minmaxs=['max']*num_axis, figtitle=fig_title, + save_fig_filename = fig_filename, kde_plot=True) \ No newline at end of file diff --git a/figure_8_plot_parallel.py b/figure_8_plot_parallel.py new file mode 100644 index 0000000..6ffea45 --- /dev/null +++ b/figure_8_plot_parallel.py @@ -0,0 +1,359 @@ +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.collections import PatchCollection +from matplotlib.patches import Rectangle +from matplotlib.lines import Line2D +from pandas.plotting import parallel_coordinates +import statsmodels.api as sm + +### function to normalize data based on direction of preference and whether each objective is minimized or maximized +### -> output dataframe will have values ranging from 0 (which maps to bottom of figure) to 1 (which maps to top) +def reorganize_objs(objs, columns_axes, ideal_direction, minmaxs): + ''' + Function to normalize data based on direction of preference and whether each objective is minimized or maximized. + Output dataframe will have values ranging from 0 (which maps to bottom of figure) to 1 (which maps to top) + Args: + objs (pd.DataFrame): dataframe of objective values + columns_axes (list): list of the subset of column names to be used as parallel axes + ideal_direction (str): 'upwards' or 'downwards' - direction of preference for each axis + minmaxs (list): list of 'max' or 'min' for each axis, indicating whether each axis should be maximized or minimized + Returns: + objs_reorg (pd.DataFrame): reorganized dataframe of objective values + top_values (pd.Series): series of top values for each axis + bottom_values (pd.Series): series of bottom values for each axis + ''' + ### if min/max directions not given for each axis, assume all should be maximized + if minmaxs is None: + minmaxs = ['max']*len(columns_axes) + + ### get subset of dataframe columns that will be shown as parallel axes + objs_reorg = objs.copy() + + ### reorganize & normalize data to go from 0 (bottom of figure) to 1 (top of figure), + ### based on direction of preference for figure and individual axes + if ideal_direction == 'downwards': + top_values = objs_reorg.min(axis=0) + bottom_values = objs_reorg.max(axis=0) + for i, minmax in enumerate(minmaxs): + if minmax == 'min': + # changed here + objs_reorg.iloc[:, i] = (objs.iloc[:, i].max(axis=0) - objs.iloc[:, i]) / \ + (objs.iloc[:, i].max(axis=0) - objs.iloc[:, i].min(axis=0)) + + else: + bottom_values[i], top_values[i] = top_values[i], bottom_values[i] + objs_reorg.iloc[:, -1] = (objs.iloc[:, -1] - objs.iloc[:, -1].min(axis=0)) / \ + (objs.iloc[:, -1].max(axis=0) - objs.iloc[:, -1].min(axis=0)) + + elif ideal_direction == 'upwards': + #top_values = objs_reorg.max(axis=0) + top_values = [1.0]*len(columns_axes) + #bottom_values = objs_reorg.min(axis=0) + bottom_values = [0.0]*len(columns_axes) + + for i, minmax in enumerate(minmaxs): + if minmax == 'max': + continue + ''' + objs_reorg.iloc[:, i] = (objs.iloc[:, i] - objs.iloc[:, i].min(axis=0)) / \ + (objs.iloc[:, i].max(axis=0) - objs.iloc[:, i].min(axis=0)) + ''' + else: + bottom_values[i], top_values[i] = top_values[i], bottom_values[i] + objs_reorg.iloc[:, i] = (objs.iloc[:, i].max(axis=0) - objs.iloc[:, i]) / \ + (objs.iloc[:, i].max(axis=0) - objs.iloc[:, i].min(axis=0)) + + objs_reorg = objs_reorg[columns_axes] + + return objs_reorg, top_values, bottom_values + +def get_color(value, color_by_continuous, color_palette_continuous, + color_by_categorical, color_dict_categorical): + ''' + Function to get color based on continuous color map or categorical map + Args: + value (float or str): value to be colored + color_by_continuous (int): index of column to be colored by continuous color map + color_palette_continuous (str): name of continuous color map + color_by_categorical (int): index of column to be colored by categorical color map + color_dict_categorical (dict): dictionary of categorical color map + + Returns: + color (str): color to be used for given value + ''' + if color_by_continuous is not None: + #color = colormaps.get_cmap(color_palette_continuous)(value) + color = cm.get_cmap(color_palette_continuous)(value) + elif color_by_categorical is not None: + color = color_dict_categorical[value] + return color + +def get_zorder(norm_value, zorder_num_classes, zorder_direction): + ''' + Function to get zorder value for ordering lines on plot. + Works by binning a given axis' values and mapping to discrete classes. + + Args: + norm_value (float): normalized value of objective + zorder_num_classes (int): number of classes to bin values into + zorder_direction (str): 'ascending' or 'descending' - direction of preference for zorder + Returns: + zorder (int): zorder value for ordering lines on plot + + ''' + xgrid = np.arange(0, 1.001, 1/zorder_num_classes) + if zorder_direction == 'ascending': + #return 4 + np.sum(norm_value > xgrid) + return np.sum(norm_value > xgrid) + elif zorder_direction == 'descending': + #return 4 + np.sum(norm_value < xgrid) + return 4 + np.sum(norm_value < xgrid) + +### customizable parallel coordinates plot +def custom_parallel_coordinates(objs, columns_axes=None, axis_labels=None, + ideal_direction='upwards', minmaxs=None, + color_by_continuous=None, color_palette_continuous=None, + color_by_categorical=None, color_palette_categorical=None, + colorbar_ticks_continuous=None, color_dict_categorical=None, + zorder_by=None, zorder_num_classes=10, zorder_direction='ascending', + alpha_base=0.8, brushing_dict=None, alpha_brush=0.1, + lw_base=1.5, fontsize=14, figtitle = None, + figsize=(11,6), save_fig_filename=None, + single_solution = False, single_solution_idx = 0, single_solution_color = 'red', + many_solutions = False, many_solutions_idx = None, many_solutions_color = None, + kde_plot=False): + ''' + Function to create a customizable parallel coordinates plot. + + Args: + objs (pd.DataFrame): dataframe of objective values + columns_axes (list): list of the subset of column names to be used as parallel axes + axis_labels (list): list of axis labels + ideal_direction (str): 'upwards' or 'downwards' - direction of preference for each axis + minmaxs (list): list of 'max' or 'min' for each axis, indicating whether each axis should be maximized or minimized + color_by_continuous (array): the array of values to be colored by continuous color map + color_palette_continuous (str): name of continuous color map + color_by_categorical (array): index of column to be colored by categorical color map + color_palette_categorical (str): name of categorical color map + colorbar_ticks_continuous (list): list of ticks for continuous color map + color_dict_categorical (dict): dictionary of categorical color map + zorder_by (int): index of column to be used for ordering lines on plot + zorder_num_classes (int): number of classes to bin values into + zorder_direction (str): 'ascending' or 'descending' - direction of preference for zorder + alpha_base (float): transparency of lines + brushing_dict (dict): dictionary of brushing criteria + alpha_brush (float): transparency of brushed lines + lw_base (float): line width + fontsize (int): font size + figtitle (str): title of figure + figsize (tuple): figure size + save_fig_filename (str): filename to save figure to + + ''' + + # get all column names of the original dataframe + #all_columns = objs.columns + + ### verify that all inputs take supported values + assert ideal_direction in ['upwards','downwards'] + assert zorder_direction in ['ascending', 'descending'] + if minmaxs is not None: + for minmax in minmaxs: + assert minmax in ['max','min'] + + assert color_by_continuous is None or color_by_categorical is None + + if columns_axes is None: + columns_axes = objs.columns + if axis_labels is None: + axis_labels = columns_axes + + ### create figure + fig,ax = plt.subplots(1,1,figsize=figsize, gridspec_kw={'hspace':0.1, 'wspace':0.1}) + + ### reorganize & normalize objective data that you want to plot + objs_reorg, tops, bottoms = reorganize_objs(objs, columns_axes, ideal_direction, minmaxs) + objs_reorg = objs_reorg.fillna(1) + + ### apply any brushing criteria + if brushing_dict is not None: + satisfice = np.zeros(objs.shape[0]) == 0. + + ### iteratively apply all brushing criteria to get satisficing set of solutions + for col_idx, (threshold, operator) in brushing_dict.items(): + if operator == '<': + satisfice = np.logical_and(satisfice, objs.iloc[:,col_idx] < threshold) + elif operator == '<=': + satisfice = np.logical_and(satisfice, objs.iloc[:,col_idx] <= threshold) + elif operator == '>': + satisfice = np.logical_and(satisfice, objs.iloc[:,col_idx] > threshold) + elif operator == '>=': + satisfice = np.logical_and(satisfice, objs.iloc[:,col_idx] >= threshold) + elif operator == '==': + satisfice = np.logical_and(satisfice, objs.iloc[:,col_idx] == threshold) + + + ### loop over all solutions/rows & plot on parallel axis plot + satisficing_solutions= [] + satisficing_np = None + if brushing_dict is not None: + satisficing_solutions = [i for i,e in enumerate(satisfice) if e == True] + satisficing_np = np.array(satisficing_solutions) + np.savetxt('satisficing.csv', satisficing_np, delimiter=',') + + for i in range(objs_reorg.shape[0]): + if color_by_continuous is not None: + color = get_color(objs_reorg[columns_axes[color_by_continuous]].iloc[i], + color_by_continuous, color_palette_continuous, + color_by_categorical, color_dict_categorical) + elif color_by_categorical is not None: + color = get_color(objs[color_by_categorical].iloc[i], + color_by_continuous, color_palette_categorical, + color_by_categorical, color_dict_categorical) + + else: + color = 'silver' + + ### order lines according to ascending or descending values of one of the objectives + if zorder_by is None: + zorder = 4 + else: + zorder = get_zorder(objs_reorg[columns_axes[zorder_by]].iloc[i], + zorder_num_classes, zorder_direction) + + ### apply any brushing? + if brushing_dict is not None: + if satisfice.iloc[i]: + alpha = alpha_base + lw = lw_base + #satisficing_solutions.append(i) + else: + alpha = alpha_brush + lw = 1.5 + zorder = 2 + elif single_solution: + if i == single_solution_idx: + alpha = alpha_base + lw = lw_base + color = single_solution_color + else: + alpha = alpha_brush + lw = 1.5 + zorder = 2 + else: + alpha = alpha_base + lw = lw_base + ### loop over objective/column pairs & plot lines between parallel axes + #for j in range(objs_reorg.shape[1]-1): + for j in range(0,len(columns_axes)-1): + y = [objs_reorg.iloc[i, j], objs_reorg.iloc[i, j+1]] + x = [j, j+1] + ax.plot(x, y, c=color, alpha=alpha, zorder=zorder, lw=lw) + + ### add top/bottom ranges + for j in range(len(columns_axes)-1): + if j < len(columns_axes): + ax.annotate(str(tops[j]), [j, 1.01], ha='center', va='bottom', + zorder=zorder, fontsize=fontsize) + if j == len(columns_axes): + ax.annotate(str(bottoms[j]) + '+', [j, -0.01], ha='center', va='top', + zorder=zorder, fontsize=fontsize) + else: + ax.annotate(str(bottoms[j]), [j, -0.01], ha='center', va='top', + zorder=zorder, fontsize=fontsize) + ax.plot([j,j], [0,1], c='lightgrey', zorder=1, linewidth=1.5) + + if single_solution: + for j in range(len(columns_axes)-1): + ax.plot([j, j+1], [objs_reorg.iloc[single_solution_idx, j], objs_reorg.iloc[single_solution_idx, j+1]], + c=single_solution_color, zorder=8, lw=3.5) + + if many_solutions: + for i in range(len(many_solutions_idx)): + idx = many_solutions_idx[i] + for j in range(len(columns_axes)-1): + ax.plot([j, j+1], [objs_reorg.iloc[idx, j], objs_reorg.iloc[idx, j+1]], + c=many_solutions_color[i], zorder=8, lw=4) + kde= sm.nonparametric.KDEUnivariate(objs_reorg.iloc[j, :]) + if j == len(columns_axes)-2: + ax.plot([j, j+1], [objs_reorg.iloc[idx, j], objs_reorg.iloc[idx, j+1]], + c=many_solutions_color[i], zorder=8, lw=4, label=f'Solution {idx}') + # plot a kde plot on the y axis + kde = sm.nonparametric.KDEUnivariate(objs_reorg.iloc[idx, :]) + ax.legend(loc='lower center', bbox_to_anchor=(0.5, 0), ncol=len(many_solutions_idx), frameon=False) + + if kde_plot: + bw = 0.025 + for j in range(len(columns_axes)): + y = np.arange(0, 1, 0.01) + x = [] + + kde = sm.nonparametric.KDEUnivariate(objs_reorg.iloc[:, j]) + kde.fit(bw=bw) + + for yy in y: + xx = kde.evaluate(yy) * 0.08 + if np.isnan(xx): + x.append(0) + else: + x.append(xx[0]) + x = np.array(x) + ax.fill_betweenx(y, x+j, j, x, color='lightseagreen', alpha=0.7, zorder=41, fc='turquoise', ec='k') + #ax.fill_betweenx(y, len(columns_axes)+1, len(columns_axes)+2, x, color='lightseagreen', alpha=0.7, zorder=0, fc='turquoise', ec='k') + + ### other aesthetics + ax.set_xticks([]) + ax.set_yticks([]) + + for spine in ['top','bottom','left','right']: + ax.spines[spine].set_visible(False) + + if ideal_direction == 'upwards': + ax.arrow(-0.15,0.1,0,0.7, head_width=0.08, head_length=0.05, color='k', lw=1.5) + elif ideal_direction == 'downwards': + ax.arrow(-0.15,0.9,0,-0.7, head_width=0.08, head_length=0.05, color='k', lw=1.5) + ax.annotate('Direction of preference', xy=(-0.3,0.5), ha='center', va='center', + rotation=90, fontsize=fontsize) + + #ax.set_xlim(-0.4, len(columns_axes)-2) + ax.set_ylim(-0.4, 1.1) + + for i,l in enumerate(axis_labels): + ax.annotate(l, xy=(i,-0.1), ha='center', va='top', fontsize=fontsize) + ax.patch.set_alpha(0) + + + ### colorbar for continuous legend + if color_by_continuous is not None: + mappable = cm.ScalarMappable(cmap=color_palette_continuous) + mappable.set_clim(vmin=0, vmax=1) + cb = plt.colorbar(mappable, ax=ax, orientation='horizontal', shrink=0.4, + label='Robustness', pad=0.03, + alpha=alpha_base) + if colorbar_ticks_continuous is not None: + _ = cb.ax.set_xticks(colorbar_ticks_continuous, colorbar_ticks_continuous, + fontsize=fontsize) + _ = cb.ax.set_xlabel(cb.ax.get_xlabel(), fontsize=fontsize) + + ### categorical legend + ''' + elif color_by_categorical is not None: + leg = [] + for label,color in color_dict_categorical.items(): + leg.append(Line2D([0], [0], color=color, lw=3, + alpha=alpha_base, label=label)) + _ = ax.legend(handles=leg, loc='lower center', + ncol=max(3, len(color_dict_categorical)), + bbox_to_anchor=[0.5,-0.07], frameon=False, fontsize=fontsize) + ''' + + if figtitle is not None: + plt.title(figtitle) + + ### save figure + if save_fig_filename is not None: + plt.savefig(save_fig_filename, bbox_inches='tight', dpi=300) + + \ No newline at end of file