import random
from amplpy import AMPL
import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib import gridspec
from matplotlib.ticker import PercentFormatter
import numpy
import pandas as pd
random.seed(a=100)
n_bins = 150
dpi = 120
solver = "highs"
# Conditional value at risk (for losses) is the expected value of the worst alpha% tail of the realizations
# in the random variable (represented in this case as a list of values)
# alpha: a value in ]0,1[
# data: a list of values
def CVaR(data, alpha):
data.sort()
n = len(data)
m = int(alpha * n)
return sum(data[m:-1]) / len(data[m:-1])
# Prints the double histogram of a random variable (represented in this case as a list of values)
# First the empirical density distribution, and then the empirical accumulated probability
def doubleHist(data, xlabel=None):
# plt.figure(figsize=(14.5,6),dpi=100)
plt.figure(figsize=(10, 6), dpi=dpi)
gs = gridspec.GridSpec(2, 1, height_ratios=[5, 1], wspace=0.25)
ax = plt.subplot(gs[0])
ax3 = plt.subplot(gs[1])
ax.grid(True)
ax.set_title(xlabel)
ax.set_ylabel("Likelihood of occurrence", color="blue")
ax2 = ax.twinx()
ax2.set_ylabel("Cumulative probability", color="red")
# plot the density and cumulative histogram
n, bins, patches = ax.hist(
data,
n_bins,
density=True,
alpha=0.65,
cumulative=False,
label="Density",
color="blue",
)
n, bins, patches = ax2.hist(
data,
n_bins,
density=True,
histtype="step",
cumulative=True,
label="CDF",
color="red",
)
ax3.boxplot(data, vert=False, whis=[5, 95])
# ax3.axis("off")
ax3.tick_params(
left=False,
right=False,
labelleft=False,
labelbottom=False,
bottom=False,
top=True,
)
# tidy up the figure
plt.show()
# Print basic info on a given solution
def solutionStats(model, losses):
global cost
losses.sort()
samples = len(losses)
order = model.var["order"].value()
print(
"Model Value",
model.get_objective("obj").value(),
"order Qty",
order,
"investment",
order * cost,
"#samples",
samples,
)
print("Losses stats:")
print(
"CVaR(95%,90%,80%,70%)",
CVaR(losses, 0.95),
CVaR(losses, 0.90),
CVaR(losses, 0.80),
CVaR(losses, 0.70),
)
print(
"min",
min(losses),
"max",
max(losses),
"average",
sum(losses) / samples,
"VaR_75%",
losses[int(0.75 * samples)],
)
# Build a dataframe with expected values and CVaR values specified in the cvars list
def getTable(data, colnames, cvars):
info = []
temp = ["Expected Value"]
for sol in data:
temp.append(sum(sol) / len(sol))
info.append(temp)
for c in cvars:
temp = ["CVaR_" + str(c) + "%"]
for sol in data:
temp.append(CVaR(sol, c / 100))
info.append(temp)
return pd.DataFrame(info, columns=["Metric"] + colnames)
# Get the list of losses from the AMPL instance
def getLosses(ampl):
return [-v[1] for v in ampl.var["profit"].getValues()]