This notebook reproduces the plots found in the paper "Beschreibung von Trennoperationen mit mehrdimensionalen Partikeleigenschaftsverteilungen" published in Chemie Ingenieur Technik. The resulting plots are not necessarily the final results as used in the article, but sometimes they are only parts of the final figures that were later stitched together in InkScape. All of the produced plots are vector graphics (.svg), but change the graph_type variable in the fourth code cell below to change the type of file to your liking.

The article will most likely have been published open access under a Create Commons CC BY license, meaning that you can copy, distribute and transmit the article, adapt the article and its contents accordingly, as long as the original authors, or rather the original journal article, are attributed. But check the final article to be safe.

Have fun!

Thomas Buchwald, June 2022

Preamble

Packages

import numpy as np
import pandas as pd
from scipy import stats
## Gaussian filters for drawing equality lines in 2D histogram plots
from scipy.ndimage.filters import gaussian_filter
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm

sns.set(font_scale=1.2, style="ticks")

## package warnings will not be shown, as not to interrupt the notebook's flow ;-)
import warnings
warnings.filterwarnings('ignore')

# Import CMasher to register colormaps
import cmasher as cmr

The colormap of a 2D histogram should be accessible to the greatest audience possible. The colormaps displayed below will result in a smooth greyscale gradient even when no color is perceived or the document is printed in greyscale. The matplotlib colormaps that fulfill this criterion are shown below.

from IPython.core.display import Image, display

colormap_string = "https://matplotlib.org/3.5.0/_images/sphx_glr_colormaps_001.png"
display(Image(url=colormap_string))

graph_type="SVG"

Unstack Function

# basically np.vstack() in reverse
def unstack(a, axis=0):
    return np.moveaxis(a, axis, 0)

Specification of Color Map

colormap = "cividis_r"

from matplotlib import cm
from matplotlib.colors import ListedColormap,LinearSegmentedColormap

cmap_modified = cm.get_cmap(colormap, 256)
ncmap = ListedColormap(cmap_modified(np.linspace(0., 1.0, 256)))

cmr.view_cmap(ncmap, show_grayscale=True)

# Take 5 colors from the colormap in [0.0, 1.0] range in RGB
colors = cmr.take_cmap_colors(colormap, 5, cmap_range=(0.0, 1.0), return_fmt='rgb')

colors

['#FEE838', '#BBAD6D', '#7C7B78', '#434E6C', '#00224E']

colormap_five = cmr.get_sub_cmap(colormap, 0.0, 1.0, N=5)
cmr.view_cmap(colormap_five, show_grayscale=True)

# Take 10 colors from the colormap in [0.0, 1.0] range in RGB
colors_10 = cmr.take_cmap_colors(colormap, 10, cmap_range=(0.0, 1.0), return_fmt='rgb')

colormap_five = cmr.get_sub_cmap(colormap, 0.0, 1.0, N=10)
cmr.view_cmap(colormap_five, show_grayscale=True)

Plot Kwargs

kwargs_sep = {"gridspec_kw" : {'width_ratios': [3, 0.8], 'height_ratios': [0.8, 3]},
              "figsize" : (5.5, 5), "sharey" : 'row', "sharex" : 'col'}

kwargs_CDF = {"gridspec_kw" : {'width_ratios': [3, 0.8], 'height_ratios': [0.8, 3]},
              "figsize" : (5.5, 5), "sharey" : 'row', "sharex" : 'col'}

Data Import

Excel to Pickled DataFrame

The dataset is imported from a huge Excel table. After an initial import, the data is "pickled", so it doesn't need to be converted again. Therefore you don't need to run the next few cells – the import has already been done, the dataset resides in Data/data.pkl.

xls = pd.ExcelFile("Data/data.xlsx")xls.sheet_namesdf = pd.DataFrame() for sheet_name in tqdm(xls.sheet_names): df_temp = pd.read_excel("Data/data.xlsx", sheet_name=sheet_name) df = pd.concat([df, df_temp])xls.sheet_names[9]df.head()df.drop([0], axis=0, inplace=True)df.drop(['Filter0', 'Filter1', 'Filter2', 'Filter3', 'Filter4', 'Filter5', 'Filter6'], axis=1, inplace=True)df.columnsdf = df.astype(float)df.to_pickle("Data/data.pkl")

Reading Pickled DataFrame

df = pd.read_pickle("Data/data.pkl")

df.head()

	Id	Img Id	Da	Dp	FWidth	FLength	FThickness	ELength	EThickness	EWidth	...	Surface Area	L/W Ratio	W/L Ratio	W/T Ratio	T/W Ratio	CHull Surface Area	Sieve	Ellipticity	Fiber Length	Fiber Width
1	101093.0	21988.0	29.826	72.672	34.009	81.818	34.009	94.058	23.118	23.118	...	2794.715	2.406	0.416	1.0	1.0	8419.495	34.009	4.069	79.538	8.784
2	102570.0	22154.0	29.709	82.837	18.083	123.691	18.083	124.014	10.082	10.082	...	2772.921	6.840	0.146	1.0	1.0	6510.668	18.083	12.300	103.526	6.696
3	37450.0	15290.0	29.668	55.372	30.395	72.132	30.395	21.292	9.739	9.739	...	2765.162	2.373	0.421	1.0	1.0	4384.862	30.395	2.186	0.000	0.000
4	248683.0	41088.0	29.326	67.100	32.673	77.247	32.673	88.978	19.259	19.259	...	2701.841	2.364	0.423	1.0	1.0	7266.333	32.673	4.620	57.612	11.724
5	74808.0	19227.0	29.078	68.382	28.131	89.354	28.131	108.358	17.780	17.780	...	2656.284	3.176	0.315	1.0	1.0	6623.341	28.131	6.094	0.000	0.000

5 rows × 37 columns

df.columns

Index(['Id', 'Img Id', 'Da', 'Dp', 'FWidth', 'FLength', 'FThickness',
       'ELength', 'EThickness', 'EWidth', 'Volume', 'Area', 'Perimeter',
       'CHull  Area', 'CHull Perimeter', 'Sphericity', 'L/T Ratio',
       'T/L Aspect Ratio', 'Compactness', 'Roundness', 'Ellipse Ratio',
       'Circularity', 'Solidity', 'Concavity', 'Convexity', 'Extent', 'hash',
       'Surface Area', 'L/W Ratio', 'W/L Ratio', 'W/T Ratio', 'T/W Ratio',
       'CHull Surface Area', 'Sieve', 'Ellipticity', 'Fiber Length',
       'Fiber Width'],
      dtype='object')

plt.scatter(df["Da"], df["Sphericity"]);

len(df)

280632

Due to the limitations of the dynamic imaging system, only the discrete particle information > 150 µm is actually usable.

df = df[df.Da > 151.]
len(df)

56566

And then also some of the columns:

df.drop(['Id', 'Img Id', 'ELength', 'EThickness', 'EWidth', 
         'CHull Perimeter', 'CHull  Area', 'Compactness', 'Ellipse Ratio', 
         'Solidity', 'Concavity', 'Extent', 'hash', 'L/T Ratio',
         'L/W Ratio', 'W/T Ratio', 'T/W Ratio', 'CHull Surface Area',
         'Ellipticity', 'Fiber Length', 'Fiber Width'], axis=1, inplace=True)

df.head()

	Da	Dp	FWidth	FLength	FThickness	Volume	Area	Perimeter	Sphericity	T/L Aspect Ratio	Roundness	Circularity	Convexity	Surface Area	W/L Ratio	Sieve
1	170.000	183.797	143.392	236.617	126.947	2255233.588	22697.987	577.416	0.925	0.537	0.516	0.855	1.000	90791.948	0.606	135.169
2	169.996	186.346	164.407	248.286	97.430	2082394.586	22696.903	585.422	0.912	0.392	0.469	0.832	0.999	90787.612	0.662	130.918
3	169.996	178.306	179.598	228.420	143.185	3075613.350	22696.814	560.164	0.953	0.627	0.554	0.909	0.999	90787.254	0.786	161.391
4	169.987	190.220	135.451	267.719	107.553	2042122.471	22694.413	597.593	0.894	0.402	0.403	0.799	0.999	90777.650	0.506	121.502
5	169.979	178.485	166.651	231.977	126.724	2565142.118	22692.295	560.726	0.952	0.546	0.537	0.907	0.999	90769.179	0.718	146.688

Correlation Matrix

sns.set(font_scale=.8, style="ticks")

corrMatrix = df.corr()
f, ax = plt.subplots(figsize=(10, 9))
ax = sns.heatmap(corrMatrix, annot=True)
plt.savefig("Plots/correlation_matrix" + "." + graph_type, format=graph_type)

sns.set(font_scale=1.2, style="ticks")

Parameter Choice

# this first set of variables is used to access the dataset in data.pkl
x_name = "Dp"
y_name = "W/L Ratio"

# this second set of variables is used in the graphs for labeling
xname = "$x_\mathrm{V}$ in µm" 
yname = "$b/l$ [$\mathrm{-}$]"

edgecolor = "k" # black
linewidth = 0

x = df[x_name]
y = df[y_name]

x_min = x.min()
y_min = y.min()
x_max = x.max() * 0.65
y_max = y.max()

The classes variable actually defines the number of edges, resulting in $n - 1$ classes.

classes = 7 # - 1 = 6

Choice of Class Edges

Linear-scaled x-axis:

xedges = np.linspace(x_min, x_max, classes)
yedges = np.linspace(y_min, y_max, classes)

xedges, yedges

(array([162.696   , 213.977475, 265.25895 , 316.540425, 367.8219  ,
        419.103375, 470.38485 ]),
 array([0.231 , 0.3535, 0.476 , 0.5985, 0.721 , 0.8435, 0.966 ]))

If both $x$ and $y$ would scale linearly, the $\Delta Q (x, y)$ and $q(x, y)$ plots would look identical, as the ground area of each bar would be the same. By scaling the particle size classes logarithmically, as one should usually do, this potential source of confusion is circumvented.

Log-scaled x-axis:

xedges = np.logspace(np.log10(x_min), np.log10(x_max), classes)
yedges = np.linspace(y_min, y_max, classes)

xedges, yedges

(array([162.696     , 194.1882512 , 231.77629999, 276.64007945,
        330.18791635, 394.1007403 , 470.38485   ]),
 array([0.231 , 0.3535, 0.476 , 0.5985, 0.721 , 0.8435, 0.966 ]))

Separation

Roger's widely used classification function is defined in the following with a separation criterion that includes both particle size and width/length ratio. Note that there is not necessarily a physical meaning to df["Dp"] * df["W/L Ratio"] ** 1.8, it just reuslts in an interesting separation function in 2D space.

Reference:

Rogers (1982) A classification function for vibrating screens
https://doi.org/10.1016/0032-5910(82)80015-6

def rogers(x, xT=100, a=0, alpha=1, boolean=True, loc=.0, scale=0.15):
    # separation percentage according to Roger's model function
    ratio = x / xT
    sep_chance = (1 - a) / ( 1 + ratio * np.exp( alpha * ( 1 - ratio ** 3 ) ) ) + a
    # an additional element of chance, can be adapted to make the separaton function sharper
    if boolean == True:
        boolean = bool(np.round_(sep_chance + np.random.normal(loc=loc, scale=scale)))
        return boolean
    else:
        return sep_chance

xT = 130
a = 0.05
alpha = 30
scale = 0.15

kwargs = {"xT": xT, "a": a, "alpha": alpha, "scale": scale}
df["shape"] = df["Dp"] * df["W/L Ratio"] ** 1.8
df["sorted"] = df["shape"].apply(rogers, **kwargs)

df["sorted"].sum() / len(df)

0.4631404023618428

CDF

Bivariate Marginal CDFs, Number- and Volume-Based

In order for the CDF graphs (i.e. the resulting SVGs) to be manageable, only 20% of the data points are used:

df[...][::5]

The above states that every fifth point of a DataFrame is used.

len(df), len(df[::5])

(56566, 11314)

x = df[x_name][::5]
y = df[y_name][::5]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_CDF)
ax1.axis("off")

ax0.plot(x.sort_values(), np.arange(1, len(x) + 1, 1) / len(x), color="black")
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$Q_0(x_\mathrm{V})$")

ax2.scatter(x, y, s=2, alpha=1.0, linewidth=0, color=colors_10[-2])
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(np.arange(1, len(y) + 1, 1) / len(y), y.sort_values(), color="black")
ax3.set_xlim(0, 1)
ax3.set_ylim(y_min, y_max)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$Q_0(b/l)$")

plt.tight_layout()
plt.savefig("Plots/cdf_Q0_discrete" + "." + graph_type, format=graph_type)

The volume-based graph is scaled by the information given in the "Volume" column of the DataFrame.

df_sample = df[[x_name, y_name, "Volume"]][::5]
df_x = df_sample[[x_name, "Volume"]].sort_values(by=x_name)
df_y = df_sample[[y_name, "Volume"]].sort_values(by=y_name)

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_CDF)
ax1.axis("off")

ax0.plot(df_x[x_name], df_x["Volume"].cumsum() / df_x["Volume"].sum(), color="black")
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$Q_3(x_\mathrm{V})$")
 
scaler = 0.0000018
ax2.scatter(df_sample[x_name], df_sample[y_name], s=df_sample["Volume"]*scaler, alpha=0.2, linewidth=0, color=colors_10[-2])
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(df_y["Volume"].cumsum() / df_y["Volume"].sum(), df_y[y_name], color="black")
ax3.set_xlim(0, 1)
ax3.set_ylim(y_min, y_max)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$Q_3(b/l)$")

plt.tight_layout()
plt.savefig("Plots/cdf_Q3_discrete" + "." + graph_type, format=graph_type)

The mean values of the classes are calculated solely for printing the respective class values.

xmeans = [0.5 * (xedges[i] + xedges[i+1]) for i in range(len(xedges) - 1)]
ymeans = [0.5 * (yedges[i] + yedges[i+1]) for i in range(len(yedges) - 1)]

Q_0, xedges, yedges = np.histogram2d(df[x_name], df[y_name],
                                          weights=None, bins=[xedges, yedges], density=False)
Q_0 = Q_0 / Q_0.sum()
Q_0 = Q_0.T
Q_0 = Q_0.cumsum(axis=1)
Q_0 = Q_0.cumsum(axis=0)

x = df[x_name][::5]
y = df[y_name][::5]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_CDF)
ax1.axis("off")

ax0.plot(x.sort_values(), np.arange(1, len(x) + 1, 1) / len(x), color="black")
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$Q_0(x_\mathrm{V})$")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, Q_0, cmap=ncmap, alpha=0.8, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(Q_0):
    if z > 0.6:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.2f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(np.arange(1, len(y) + 1, 1) / len(y), y.sort_values(), color="black")
ax3.set_xlim(0, 1)
ax3.set_ylim(y_min, y_max)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$Q_0(b/l)$")

plt.tight_layout()
plt.savefig("Plots/cdf_Q0_classes" + "." + graph_type, format=graph_type)

Q_3, xedges, yedges = np.histogram2d(df[x_name], df[y_name],
                                          weights=df["Volume"], bins=[xedges, yedges], density=False)
Q_3 = Q_3 / Q_3.sum()
Q_3 = Q_3.T
Q_3 = Q_3.cumsum(axis=1)
Q_3 = Q_3.cumsum(axis=0)

df_sample = df[[x_name, y_name, "Volume"]][::5]
df_x = df_sample[[x_name, "Volume"]].sort_values(by=x_name)
df_y = df_sample[[y_name, "Volume"]].sort_values(by=y_name)

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_CDF)
ax1.axis("off")

ax0.plot(df_x[x_name], df_x["Volume"].cumsum() / df_x["Volume"].sum(), color="black")
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$Q_3(x_\mathrm{V})$")
 
X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, Q_3, cmap=ncmap, alpha=0.8, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(Q_3):
    if z > 0.6:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.2f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(df_y["Volume"].cumsum() / df_y["Volume"].sum(), df_y[y_name], color="black")
ax3.set_xlim(0, 1)
ax3.set_ylim(y_min, y_max)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$Q_3(b/l)$")

plt.tight_layout()
plt.savefig("Plots/cdf_Q3_classes" + "." + graph_type, format=graph_type)

Number-Based Bar Plot

Although the bar plots don't have numbered $x$ and $y$ axes, the data contained in them is the same as in the surface plots, as they all stem from the same dataset.

H_supply, xedges, yedges = np.histogram2d(df[x_name], df[y_name],
                                          weights=None, bins=[xedges, yedges], density=False)
# Histogram does not follow Cartesian convention,
# therefore transpose H for visualization purposes.
H_supply = H_supply.T

xmeans = [0.5 * (xedges[i] + xedges[i+1]) for i in range(len(xedges) - 1)]
ymeans = [0.5 * (yedges[i] + yedges[i+1]) for i in range(len(yedges) - 1)]

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX, YY = np.meshgrid(xmeans, ymeans)
X, Y = XX.ravel(), YY.ravel()
Z = H_supply.ravel()

bar_colors = ["gold"] * 36

fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_subplot(111, projection='3d')

top = Z
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=0.6, shade=False, color=bar_colors, edgecolor="black", lw=1.5)

ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$N$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

xticks = xedges
xlabels = ["$x_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1)
yticks = [*YY.T[0], YY.T[0][-1] + depth]
ylabels = ["$y_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major')
ax1.tick_params(axis='z', which='major', pad=10)

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(False)
ax1.yaxis.set_rotate_label(False)

plt.tight_layout()
plt.savefig("Plots/bar_plot_number" + "." + graph_type, format=graph_type)

Volume-Based Bar Plot, Weighted

H_supply, xedges, yedges = np.histogram2d(df[x_name], df[y_name],
                                          weights=df["Volume"], bins=[xedges, yedges], density=False)
H_supply = H_supply.T

xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX, YY = np.meshgrid(xmeans, ymeans)
X, Y = XX.ravel(), YY.ravel()
Z = H_supply.ravel()

bar_colors = ["gold"] * 36

fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_subplot(111, projection='3d')

top = Z / Z.cumsum().max()
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=0.6, shade=False, color=bar_colors, edgecolor="black", lw=1.5)

ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$\Delta Q_3(x, y)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

xticks = xedges
xlabels = ["$x_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1)
yticks = [*YY.T[0], YY.T[0][-1] + depth]
ylabels = ["$y_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major')
ax1.tick_params(axis='z', which='major', pad=10)

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(False)
ax1.yaxis.set_rotate_label(False)

plt.tight_layout()
plt.savefig("Plots/bar_plot_DeltaQ" + "." + graph_type, format=graph_type)

CDF Plots

2D with Marginals

H_sum = H_supply.cumsum(axis=1)
H_sum = H_sum.cumsum(axis=0)

xmeans = [0.5 * (xedges[i] + xedges[i+1]) for i in range(len(xedges) - 1)]
ymeans = [0.5 * (yedges[i] + yedges[i+1]) for i in range(len(yedges) - 1)]

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

df_x = df[[x_name, "Volume"]][::5].sort_values(by=x_name)
df_y = df[[y_name, "Volume"]][::5].sort_values(by=y_name)

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

## because the Q3 should be derived from the classes, the Q3 plots only have points at the edges of the classes
ax0.plot(xedges, np.interp(xedges, df_x[x_name], df_x["Volume"].cumsum() / df_x["Volume"].sum()), color="black")
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$Q_3(x_\mathrm{V})$")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_sum, alpha=.8, cmap=ncmap, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(H_sum / H_sum.max()):
    if z > 0.6:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.2f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(np.interp(yedges, df_y[y_name], df_y["Volume"].cumsum() / df_y["Volume"].sum()), yedges, color="black")
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, 1)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$Q_3(b/l)$")
ax3.set_ylabel("")

ax0.yaxis.labelpad=20
ax2.yaxis.labelpad=20

plt.tight_layout()
plt.savefig("Plots/CDF_2D" + "." + graph_type, format=graph_type)

3D Bar Plot

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX, YY = np.meshgrid(xmeans, ymeans)
X, Y = XX.ravel(), YY.ravel()

Z = H_sum.ravel() / H_sum.max()

bar_colors = ["gold"] * 36

fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_subplot(111, projection='3d')

top = Z
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=0.6, shade=False, color=bar_colors, edgecolor="black", lw=1.5)
ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$Q_0(x, y)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

xticks = xedges
xlabels = ["$x_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1)
yticks = [*YY.T[0], YY.T[0][-1] + depth]
ylabels = ["$y_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major')
ax1.tick_params(axis='z', which='major', pad=10)

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(False)
ax1.yaxis.set_rotate_label(False)

plt.tight_layout()
plt.savefig("Plots/CDF_bar_plot" + "." + graph_type, format=graph_type)

3D Surface Plot

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:]
ymeans = yedges[:]
XX, YY = np.meshgrid(xmeans, ymeans)
X, Y = XX.ravel(), YY.ravel()

Z = H_sum / H_sum.max()

H_array = np.zeros([7, 7])
H_array[1:, 1:] = Z

fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_subplot(111, projection='3d')

ax1.plot_surface(XX, YY, H_array, alpha=0.8, cmap=ncmap, edgecolors="k", lw=0.8, antialiased=True)
ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$Q_3(x, y)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

xticks = [*XX[0]]
xlabels = ["$x_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1)
yticks = [*YY.T[0]]
ylabels = ["$y_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major')
ax1.tick_params(axis='z', which='major', pad=10)

#ax1.contour(XX, YY, H_array, zdir='x', offset=xedges[0], cmap=ncmap, levels=[0.1, 0.5, 0.9])
#ax1.contour(XX, YY, H_array, zdir='y', offset=yedges[-1], cmap=ncmap, levels=[0.1, 0.5, 0.9])

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(False)
ax1.yaxis.set_rotate_label(False)

plt.tight_layout()
plt.savefig("Plots/CDF_surface_plot" + "." + graph_type, format=graph_type)

PDF

pdf_x_max = 0.0115
pdf_y_max = 5.0
pdf_color_cut = 15 

## edges of cmap definition
vmin = 0.0
vmax = 30

Supply Fraction

weights = df["Volume"]

H_supply, xedges, yedges = np.histogram2d(df[x_name], df[y_name], 
                                          weights=weights, bins=[xedges, yedges], density=True)
H_supply = H_supply.T * 1000

xmeans = [0.5 * (xedges[i] + xedges[i+1]) for i in range(len(xedges) - 1)]
ymeans = [0.5 * (yedges[i] + yedges[i+1]) for i in range(len(yedges) - 1)]

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df[x_name]
y = df[y_name]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

ax0.hist(x, density=True, bins=xedges, color=colors[3])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, pdf_x_max)
ax0.set_ylabel("$q_3(x_\mathrm{V})$\nin 1/µm")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_supply, alpha=.8, cmap=ncmap, vmin=vmin, vmax=vmax, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(H_supply):
    if z > pdf_color_cut:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.1f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.hist(y, density=True, bins=yedges, orientation="horizontal", color=colors[3])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, pdf_y_max)
ax3.set_xticks([0.0, 2, 4])
ax3.set_xlabel("$q_3(b/l)$ [–]")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/PDF_supply" + "." + graph_type, format=graph_type)

The graph below is only created here to dmeonstrate how a KDE can be used to draw the density distribution function.

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df[x_name]
y = df[y_name]

bins_x = np.linspace(x_min, x_max, 100)
bins_y = np.linspace(y_min, y_max, 100)

kernel_x = stats.gaussian_kde(x, bw_method="scott")
kernel_y = stats.gaussian_kde(y, bw_method="scott")

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

ax0.plot(bins_x, kernel_x(bins_x), lw=2, color=colors_10[-2])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, pdf_x_max)
ax0.set_ylabel("$q_3(x_\mathrm{V})$\nin 1/µm")

ax2.scatter(df_sample[x_name], df_sample[y_name], s=df_sample["Volume"]*scaler, alpha=0.2, linewidth=0, color=colors_10[-2])
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.plot(kernel_y(bins_y), bins_y, lw=2, color=colors_10[-2])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, pdf_y_max)
ax3.set_xticks([0.0, 2, 4])
ax3.set_xlabel("$q_3(b/l)$ [–]")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/PDF_supply_discrete" + "." + graph_type, format=graph_type)

Coarse Fraction

df_coarse = df[df["sorted"] == True].copy()

weights = df_coarse["Volume"]

H_coarse, xedges, yedges = np.histogram2d(df_coarse[x_name], df_coarse[y_name], 
                                          weights=weights, bins=[xedges, yedges], density=True)
H_coarse = H_coarse.T * 1000

H_coarse.max()

29.846588977371656

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df_coarse[::10][x_name]
y = df_coarse[::10][y_name]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

ax0.hist(x, density=True, bins=xedges, color=colors[3])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, pdf_x_max)
ax0.set_ylabel("$q_3(x_\mathrm{V})$\nin 1/µm")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_coarse, alpha=0.8, cmap=ncmap, vmin=vmin, vmax=vmax, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(H_coarse):
    if z > pdf_color_cut:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.1f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.hist(y, density=True, bins=yedges, orientation="horizontal", color=colors[3])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, pdf_y_max)
ax3.set_xticks([0.0, 2, 4])
ax3.set_xlabel("$q_3(b/l)$ [–]")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/PDF_coarse" + "." + graph_type, format=graph_type)

sep_fraction = df_coarse["Volume"].sum() / df["Volume"].sum()
sep_fraction

0.5356982683469017

Fine Fraction

df_fines = df[df["sorted"] == False].copy()

weights = df_fines["Volume"]

H_fines, xedges, yedges = np.histogram2d(df_fines[x_name], df_fines[y_name],
                                         weights=weights, bins=[xedges, yedges], density=True)
H_fines = H_fines.T * 1000

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df_fines[::10][x_name]
y = df_fines[::10][y_name]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

ax0.hist(x, density=True, bins=xedges, color=colors[3])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, pdf_x_max)
ax0.set_ylabel("$q_3(x_\mathrm{V})$\nin 1/µm")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_fines, alpha=0.8, cmap=ncmap, vmin=vmin, vmax=vmax, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(H_fines):
    if z > pdf_color_cut:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.1f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")
    
ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.hist(y, density=True, bins=yedges, orientation="horizontal", color=colors[3])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, pdf_y_max)
ax3.set_xticks([0.0, 2, 4])
ax3.set_xlabel("$q_3(b/l)$ [–]")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/PDF_fines" + "." + graph_type, format=graph_type)

Separation Function

H_sep = H_coarse / H_supply * sep_fraction
H_sep = np.nan_to_num(H_sep, nan=0.0)
H_sep

array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 3.50564918e-03, 2.56717364e-03, 1.27427594e-03,
        0.00000000e+00, 4.66526414e-03],
       [0.00000000e+00, 1.55656731e-03, 9.61594207e-04, 6.70193878e-04,
        2.69042442e-01, 6.52467504e-01],
       [1.75413785e-03, 1.53252562e-03, 2.75720295e-01, 7.97328877e-01,
        1.00909457e+00, 1.01119286e+00],
       [1.88291297e-01, 6.96917053e-01, 1.00843606e+00, 1.01016112e+00,
        1.00876190e+00, 1.01119286e+00],
       [9.79875757e-01, 1.01082113e+00, 1.00984707e+00, 1.01017864e+00,
        1.01119286e+00, 1.01119286e+00]])

T_x = []
for digit in np.arange(0, len(xmeans), 1):
    coarse = df_coarse[(df_coarse[x_name] > xedges[digit]) & (df_coarse[x_name] < xedges[digit + 1])][x_name].count()
    feed = df[(df[x_name] > xedges[digit]) & (df[x_name] < xedges[digit + 1])][x_name].count()
    T_x.append(coarse / feed * sep_fraction)

T_y = []
for digit in np.arange(0, len(ymeans), 1):
    coarse = df_coarse[(df_coarse[y_name] > yedges[digit]) & (df_coarse[y_name] < yedges[digit + 1])][y_name].count()
    feed = df[(df[y_name] > yedges[digit]) & (df[y_name] < yedges[digit + 1])][y_name].count()
    T_y.append(coarse / feed * sep_fraction)

np.max(T_x), np.max(T_y)

(0.3912419461885943, 0.5321012535954117)

X, Y = np.mgrid[x_min:x_max:6j, y_min:y_max:6j]
positions = np.vstack([X.ravel(), Y.ravel()])

H_sep

array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 3.50564918e-03, 2.56717364e-03, 1.27427594e-03,
        0.00000000e+00, 4.66526414e-03],
       [0.00000000e+00, 1.55656731e-03, 9.61594207e-04, 6.70193878e-04,
        2.69042442e-01, 6.52467504e-01],
       [1.75413785e-03, 1.53252562e-03, 2.75720295e-01, 7.97328877e-01,
        1.00909457e+00, 1.01119286e+00],
       [1.88291297e-01, 6.96917053e-01, 1.00843606e+00, 1.01016112e+00,
        1.00876190e+00, 1.01119286e+00],
       [9.79875757e-01, 1.01082113e+00, 1.00984707e+00, 1.01017864e+00,
        1.01119286e+00, 1.01119286e+00]])

H_sep_new = []

for list_ in H_sep:
    list_ = [value if value < 1.0 else 1.0 for value in list_]
    H_sep_new.append(list_)
    
H_sep_new = np.array(H_sep_new)
H_sep_new

array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 3.50564918e-03, 2.56717364e-03, 1.27427594e-03,
        0.00000000e+00, 4.66526414e-03],
       [0.00000000e+00, 1.55656731e-03, 9.61594207e-04, 6.70193878e-04,
        2.69042442e-01, 6.52467504e-01],
       [1.75413785e-03, 1.53252562e-03, 2.75720295e-01, 7.97328877e-01,
        1.00000000e+00, 1.00000000e+00],
       [1.88291297e-01, 6.96917053e-01, 1.00000000e+00, 1.00000000e+00,
        1.00000000e+00, 1.00000000e+00],
       [9.79875757e-01, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
        1.00000000e+00, 1.00000000e+00]])

# for 1d density plots
kwargs_1d_kde = {"bw_method": "scott"}
kwargs_1d_plot = {"c": "black"}

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

x_range = np.linspace(x_min, x_max, 7)
ax0.plot(xmeans, T_x, **kwargs_1d_plot)
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 1)
ax0.set_ylabel("$T(x_\mathrm{V})$")

XX, YY = np.meshgrid(xedges, yedges)
im = ax2.pcolormesh(XX, YY, H_sep, alpha=0.8, cmap=ncmap, edgecolors=edgecolor, lw=linewidth)
ax2.contour(X, Y, gaussian_filter(H_sep_new.T, 0.4), [0.1, 0.5, 0.9], colors=["black", "lightgrey", "white"], linewidths=[1, 3, 1], linestyles="dashed")
for (i, j), z in np.ndenumerate(H_sep_new):
    if z > 0.6:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.2f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

y_range = np.linspace(y_min, y_max, 7)
ax3.plot(T_y, ymeans, **kwargs_1d_plot)
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, 1)
ax3.set_xticks([0.0, 0.5, 1.0])
ax3.set_xlabel("$T(b/l)$")
ax3.set_ylabel("")

ax0.yaxis.labelpad=20
ax2.yaxis.labelpad=20

plt.tight_layout()
plt.savefig("Plots/sep_func" + "." + graph_type, format=graph_type)

Separate Colorbars

fig, (ax1) = plt.subplots(1, 1, figsize=(2,3.5))
ax1.axis("off")
im = ax1.pcolormesh(XX, YY, Z, alpha=.8, cmap=ncmap)
fig.gca().set_visible(False)
cbar = fig.colorbar(im, ticks=[0.001, 0.2, 0.4, 0.6, 0.8, 1.0])
cbar.ax.set_yticklabels(["0.0", "0.2", "0.4", "0.6", "0.8", "1.0"])
cbar.ax.set_title('$T(x_\mathrm{V}, b/l)$\n[–]\n')
plt.tight_layout()
plt.savefig("Plots/colorbar_T" + "." + graph_type, format=graph_type)

df_coarse = df[df["sorted"] == True].copy()

H_bar = np.zeros(H_coarse.shape)
H_bar[0] = [3.1, 3.1, 3.1, 2.5, 2.6, 0.0]

fig, (ax1) = plt.subplots(1, 1, figsize=(2,3.8))
ax1.axis("off")
X, Y = np.meshgrid(xedges, yedges)
im = ax1.pcolormesh(X, Y, H_supply, alpha=0.8, cmap=ncmap, vmin=vmin, vmax=vmax, edgecolors=edgecolor, lw=linewidth)
fig.gca().set_visible(False)
cbar = fig.colorbar(im)
cbar.ax.set_title('$q_3(x_\mathrm{V}, b/l)$\nin 1/mm\n')
plt.tight_layout()
plt.savefig("Plots/colorbar_q_small" + "." + graph_type, format=graph_type)

df_coarse = df[df["sorted"] == True].copy()

H_bar = np.zeros(H_coarse.shape)
H_bar[0] = [1.0, 0.0, 1.0, 0.0, 1.0, 0.0]

fig, (ax1) = plt.subplots(1, 1, figsize=(2,3.8))
ax1.axis("off")
X, Y = np.meshgrid(xedges, yedges)
im = ax1.pcolormesh(X, Y, H_bar, alpha=0.8, cmap=ncmap)
fig.gca().set_visible(False)
cbar = fig.colorbar(im)
cbar.ax.set_title('$Q_3(x_\mathrm{V}, b/l)$\n[–]\n')
plt.tight_layout()
plt.savefig("Plots/colorbar_Q_large" + "." + graph_type, format=graph_type)

Probability Density Example

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX_, YY_ = np.meshgrid(xmeans, ymeans)
X, Y = XX_.ravel(), YY_.ravel()

Z_ = H_supply.ravel()

bar_colors = ["gold"] * 36 #[colors[0]] * 3
bar_colors[21] = "royalblue"

fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_subplot(111, projection='3d')

top = Z_
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=0.6, shade=False, color=bar_colors, edgecolor="black", lw=1.5)

ax1.set_xticks([200, 300, 400])
ax1.set_yticks([0.2, 0.4, 0.6, 0.8])
ax1.set_zticks([0.00, 5., 10, 15, 20])
ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$q_3(x, y)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

xticks = xedges
xlabels = ["$x_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1)
yticks = [*YY_.T[0], YY_.T[0][-1] + depth]
ylabels = ["$y_{}$".format("{" + str(i) + "}") for i in np.arange(1, 8, 1)]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major', pad=4)
ax1.tick_params(axis='z', which='major', pad=10)

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(False)
ax1.yaxis.set_rotate_label(False)

plt.tight_layout()
plt.savefig("Plots/PDF_bar_example" + "." + graph_type, format=graph_type)

width = width[20]
width

44.86377945081384

bar_colors = ["gold"] * 3 #[colors[0]] * 3
bar_colors[1] = "royalblue" #colors[3]

X_ = [X[0], X[21], X[35]]
Y_ = [Y[0], Y[21], Y[35]]
top = [Z_[0], Z_[21] / 2, Z_[35]]

fig = plt.figure(figsize=(12, 11))
ax1 = fig.add_subplot(111, projection='3d')

bottom = np.zeros_like(top)
width = width
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X_, Y_, bottom, width, depth, top, alpha=0.5, shade=False, color=bar_colors, edgecolor="black", lw=2.0)

ax1.set_zlim(0, 10)
ax1.set_xticks([200, 300, 400])
ax1.set_yticks([0.2, 0.4, 0.6, 0.8])
ax1.set_zticks([0.00, 0.005, 0.01, 0.015, 0.02])
ax1.set_xlabel(x_name)
ax1.set_ylabel(y_name)
ax1.set_zlabel("$q_0(x)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18
ax1.tick_params(axis='z', which='major', pad=10)

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.axis("off")

plt.tight_layout()
plt.savefig("Plots/PDF_bar_example_single_bar" + "." + graph_type, format=graph_type)

Definitions

fig = plt.figure(figsize=(5.2, 4.5))
ax1 = fig.add_subplot(111)

xticks = [5, 11.5, 20, 31, 45]
yticks = [0, 10, 20, 30, 40]

ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.set_xlim(xticks[0], xticks[-1])
ax1.set_ylim(yticks[0], yticks[-1])

# Delta x_i
ax1.arrow(xticks[1], yticks[0] + 3, xticks[2] - xticks[1], 0, 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.arrow(xticks[2], yticks[0] + 3, xticks[1] - xticks[2], 0, 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.annotate("$\Delta x_i$", (xticks[1] + (xticks[2] - xticks[1]) * 0.35, yticks[0] + 4), fontsize="small")

# Delta x_i+1
ax1.arrow(xticks[2], yticks[0] + 3, xticks[3] - xticks[2], 0, 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.arrow(xticks[3], yticks[0] + 3, xticks[2] - xticks[3], 0, 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.annotate("$\Delta x_{i+1}$", (xticks[2] + (xticks[3] - xticks[2]) * 0.25, yticks[0] + 4), fontsize="small")

# Delta y_i
ax1.arrow(xticks[3] + 3.5, yticks[1], 0, yticks[2] - yticks[1],
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.arrow(xticks[3] + 3.5, yticks[2], 0, yticks[1] - yticks[2], 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.annotate("$\Delta y_j$", (xticks[3] + 4.5, yticks[1] + (yticks[2] - yticks[1]) * 0.45), fontsize="small")

# Delta y_i+1
ax1.arrow(xticks[3] + 3.5, yticks[2], 0, yticks[3] - yticks[2],
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.arrow(xticks[3] + 3.5, yticks[3], 0, yticks[2] - yticks[3], 
          head_width=1.2, head_length=1.5, fc='k', ec='k', shape="full", lw=1., length_includes_head=True)
ax1.annotate("$\Delta y_{j+1}$", (xticks[3] + 4.5, yticks[2] + (yticks[3] - yticks[2]) * 0.45), fontsize="small")

# Delta Q_r,i's
ax1.annotate("$\Delta Q_{r,i,j}$", (xticks[1] + 0.9, yticks[1] + 2.5), fontsize="small", rotation=45)
ax1.annotate("$\Delta Q_{r,i+1,j}$", (xticks[2] + 1.7, yticks[1] + 2.0), fontsize="small", rotation=45)
ax1.annotate("$\Delta Q_{r,i,j+1}$", (xticks[1] + 0.5, yticks[2] + 2.25), fontsize="small", rotation=45)
ax1.annotate("$\Delta Q_{r,i+1,j+1}$", (xticks[2] + 0.9, yticks[2] + 1.2), fontsize="small", rotation=45)

for value in xticks:
    ax1.axvline(value, c="k", lw=1.)
for value in yticks:
    ax1.axhline(value, c="k", lw=1.)
    
xlabels = ["...", "$x_{i-1}$", "$x_{i}$", "$x_{i+1}$", "..."]
ax1.set_xticks(xticks, labels=xlabels)
ax1.tick_params(axis='x', which='major', pad=1, colors="black")
ylabels = ["...", "$y_{j-1}$", "$y_{j}$", "$y_{j+1}$", "..."]
ax1.set_yticks(yticks, labels=ylabels)
ax1.tick_params(axis='y', which='major')

# x_mean
ax2 = ax1.twiny()
ax2.set_xlim(xticks[0], xticks[-1])
ax2.set_ylim(yticks[0], yticks[-1])

xticks_ = [8.25, 15.75, 25.5, 38]
xlabels_ = ["$\overline{x}_{i-1}$", "$\overline{x}_{i}$", "$\overline{x}_{i+1}$", "$\overline{x}_{i+2}$"]
ax2.set_xticks(xticks_, labels=xlabels_, fontsize="small")
ax2.tick_params(axis='x', which='major', pad=1)

# y_mean
ax3 = ax1.twinx()
ax3.set_xlim(xticks[0], xticks[-1])
ax3.set_ylim(yticks[0], yticks[-1])

yticks_ = [5, 15, 25, 35]
ylabels_ = ["$\overline{y}_{j-1}$", "$\overline{y}_{j}$", "$\overline{y}_{j+1}$", "$\overline{y}_{j+2}$"]
ax3.set_yticks(yticks_, labels=ylabels_, fontsize="small")
ax3.tick_params(axis='y', which='major')

plt.tight_layout()
plt.savefig("Plots/definitions" + "." + graph_type, format=graph_type)

Graphical Abstract

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX_, YY_ = np.meshgrid(xmeans, ymeans)
X, Y = XX_.ravel(), YY_.ravel()

Z_ = H_supply.ravel()

import matplotlib

norm = matplotlib.colors.Normalize(vmin=vmin, vmax=20, clip=True)
mapper = cm.ScalarMappable(norm=norm, cmap=ncmap)

bar_colors = []
for value in Z_:
    bar_colors.append(mapper.to_rgba(value))

fig = plt.figure(figsize=(4.3, 4))
ax1 = fig.add_subplot(111, projection='3d')

top = Z_
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=1.0, shade=False, color=bar_colors, edgecolor="black", lw=1.0)

ax1.set_xticks([200, 300, 400])
ax1.set_yticks([0.2, 0.4, 0.6, 0.8])
ax1.set_zticks([0.00, 5., 10, 15, 20])
ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$q_3(x, y)$")
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(True)
ax1.yaxis.set_rotate_label(True)

plt.tight_layout()
plt.savefig("Plots/graphical_abstract_1" + "." + graph_type, format=graph_type)

kwargs_sep_ = {"gridspec_kw" : {'width_ratios': [3, 0.8], 'height_ratios': [0.8, 3]},
              "figsize" : (4.3, 4), "sharey" : 'row', "sharex" : 'col'}

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df[x_name]
y = df[y_name]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep_)
ax1.axis("off")

ax0.hist(x, density=True, bins=xedges, color=colors[3])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, 0.01)
ax0.set_ylabel("$q_3(x)$")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_supply, alpha=1, cmap=ncmap, vmin=vmin, vmax=20, edgecolors="k", lw=1.0)

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel("$x$")
ax2.set_ylabel("$y$")

ax3.hist(y, density=True, bins=yedges, orientation="horizontal", color=colors[3])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, 3.5)
ax3.set_xticks([0.0, 1.5, 3])
ax3.set_xlabel("$q_3(y)$")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/graphical_abstract_2" + "." + graph_type, format=graph_type)

fig, (ax1) = plt.subplots(1, 1, figsize=(2,3.8))
ax1.axis("off")
X, Y = np.meshgrid(xedges, yedges)
im = ax1.pcolormesh(X, Y, H_supply, alpha=1.0, cmap=ncmap, vmin=vmin, vmax=20, edgecolors=edgecolor, lw=linewidth)
fig.gca().set_visible(False)
cbar = fig.colorbar(im)
cbar.ax.set_title('$q_3(x, y)$\n')
plt.tight_layout()
plt.savefig("Plots/graphical_abstract_3" + "." + graph_type, format=graph_type)

Graphs again for Definition of Density Graph

xmeans = [0.5 * (xedges[i] + xedges[i+1]) for i in range(len(xedges) - 1)]
ymeans = [0.5 * (yedges[i] + yedges[i+1]) for i in range(len(yedges) - 1)]

# for 1d density plots
kwargs_1d_kde = {"bw_method": 0.2}
kwargs_1d_plot = {"c": "black"}

x = df[x_name]
y = df[y_name]

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, **kwargs_sep)
ax1.axis("off")

ax0.hist(x, density=True, bins=xedges, color=colors[3])
ax0.set_xlim(x_min, x_max)
ax0.set_ylim(0, pdf_x_max)
ax0.set_ylabel("$q_3(x_\mathrm{V})$\nin 1/µm")

X, Y = np.meshgrid(xedges, yedges)
ax2.pcolormesh(X, Y, H_supply, alpha=.8, cmap=ncmap, vmin=vmin, vmax=20, edgecolors=edgecolor, lw=linewidth)
for (i, j), z in np.ndenumerate(H_supply):
    if z > pdf_color_cut:
        color="white"
    else:
        color="black"
    ax2.text(xmeans[j], ymeans[i], '{:0.1f}'.format(z), ha='center', va='center', color=color, fontsize="x-small")

ax2.set_xlim(x_min, x_max)
ax2.set_ylim(y_min, y_max)
ax2.set_xlabel(xname)
ax2.set_ylabel(yname)

ax3.hist(y, density=True, bins=yedges, orientation="horizontal", color=colors[3])
ax3.set_ylim(y_min, y_max)
ax3.set_xlim(0, pdf_y_max)
ax3.set_xticks([0.0, 2, 4])
ax3.set_xlabel("$q_3(b/l)$ [–]")
ax3.set_ylabel("")

plt.tight_layout()
plt.savefig("Plots/PDF_definition_2" + "." + graph_type, format=graph_type)

## the bar plot needs the lower edge as reference for each bar,
## so the list of edges for the 2D contour plot is simply stripped by the last value
xmeans = xedges[:-1]
ymeans = yedges[:-1]
XX_, YY_ = np.meshgrid(xmeans, ymeans)
X, Y = XX_.ravel(), YY_.ravel()

fig = plt.figure(figsize=(5.5, 5))
ax1 = fig.add_subplot(111, projection='3d')

top = Z_
bottom = np.zeros_like(top)
width = np.array(list(np.diff(xedges)) * 6)
depth = ymeans[1] - ymeans[0]

ax1.bar3d(X, Y, bottom, width, depth, top, alpha=0.8, shade=False, color=bar_colors, edgecolor="black", lw=1.0)

ax1.set_xticks([200, 300, 400])
ax1.set_yticks([0.2, 0.4, 0.6, 0.8])
ax1.set_zticks([0.00, 5., 10, 15, 20])
ax1.set_xlabel(xname)
ax1.set_ylabel(yname)
ax1.set_zlabel('$q_3(x_\mathrm{V}, b/l)$\nin 1/mm')
ax1.xaxis.labelpad=10
ax1.yaxis.labelpad=10
ax1.zaxis.labelpad=18

ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))

ax1.xaxis.set_rotate_label(True)
ax1.yaxis.set_rotate_label(True)

plt.tight_layout()
plt.savefig("Plots/PDF_definition_1" + "." + graph_type, format=graph_type)

fig, (ax1) = plt.subplots(1, 1, figsize=(2,3.8))
ax1.axis("off")
X, Y = np.meshgrid(xedges, yedges)
im = ax1.pcolormesh(X, Y, H_supply, alpha=1.0, cmap=ncmap, vmin=vmin, vmax=20, edgecolors=edgecolor, lw=linewidth)
fig.gca().set_visible(False)
cbar = fig.colorbar(im)
cbar.ax.set_title('$q_3(x_\mathrm{V}, b/l)$\nin 1/mm\n')
plt.tight_layout()
plt.savefig("Plots/PDF_definition_cbar" + "." + graph_type, format=graph_type)