Final for Paper
This commit is contained in:
parent
e0df453875
commit
eb0b99e969
|
|
@ -237,5 +237,5 @@ dmypy.json
|
|||
# Cython debug symbols
|
||||
cython_debug/
|
||||
|
||||
out.xlsx
|
||||
animation.mp4
|
||||
output/*
|
||||
!output/.keep
|
||||
|
|
|
|||
|
|
@ -1,11 +1,6 @@
|
|||
import itertools
|
||||
import logging
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
from celluloid import Camera
|
||||
from scipy import stats
|
||||
import seaborn as sns
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
DEFAULT_AGES = np.arange(0, 50)
|
||||
DEFAULT_MORTALITY_DISTRIBUTION = 5 * stats.exponweib(a=1, c=3, loc=0, scale=30).pdf(DEFAULT_AGES)
|
||||
|
|
@ -53,6 +48,9 @@ def steady_state_population(mortality_distribution: np.array = None, total_popul
|
|||
new_cars = ((old_population) + 0) - population_before_creation
|
||||
population_distribution[0] = new_cars
|
||||
|
||||
if i > 1000:
|
||||
raise Exception("Failed to converge after 1000 iterations!")
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
logger.info(f"Steady-state population achieved after {i} iterations.")
|
||||
|
||||
|
|
@ -103,52 +101,6 @@ def iterate(mortality_distribution: np.array, population_distribution: np.array,
|
|||
|
||||
return population_distribution
|
||||
|
||||
def plot_mortality(mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION, subsidy_distribution = None,
|
||||
covid_distribution = None):
|
||||
logging.basicConfig(level=logging.DEBUG)
|
||||
|
||||
plt.figure()
|
||||
if subsidy_distribution is None:
|
||||
plot_ = sns.barplot(np.arange(len(mortality_distribution)), mortality_distribution, color='cornflowerblue')
|
||||
elif covid_distribution is None:
|
||||
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)), np.arange(len(mortality_distribution)))),
|
||||
"label": list(itertools.chain(["Normal"]*len(mortality_distribution), ["With Subsidy"]*len(mortality_distribution))),
|
||||
"Legend": np.hstack((mortality_distribution, subsidy_distribution))})
|
||||
else:
|
||||
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"label": list(itertools.chain(["Normal"] * len(mortality_distribution),
|
||||
["With Subsidy"] * len(mortality_distribution),
|
||||
["With Covid"] * len(mortality_distribution))),
|
||||
"Scenario": np.hstack((mortality_distribution,
|
||||
subsidy_distribution,
|
||||
covid_distribution))})
|
||||
plot_ = sns.barplot(data=df, x="Ages", y="Scenario", hue="label")
|
||||
plt.title("Mortality Distribution")
|
||||
plt.ylabel("Recycling probability")
|
||||
plt.xlabel("Age (years)")
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
def plot_population(population_distribution: np.ndarray):
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(np.arange(len(population_distribution)), population_distribution / 1e6, color='cornflowerblue')
|
||||
plt.title("Population Distribution")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Age (years)")
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
def sanity_check_for_subsidy():
|
||||
population_distribution = steady_state_population()
|
||||
|
|
@ -163,60 +115,3 @@ def sanity_check_for_subsidy():
|
|||
|
||||
assert np.abs((subsidy_new_cars - normal_year_new_cars) - 1e6) < 10000
|
||||
|
||||
def plot_population_development(years, total_pops: np.array):
|
||||
# Plot total vehicle pop
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(years, np.array(total_pops) / 1e6, color='cornflowerblue')
|
||||
plt.title("Total Vehicle Population")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Simulation year")
|
||||
plt.ylim(min(np.array(total_pops) / 1e6), max(np.array(total_pops) / 1e6))
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
def plot_production(years, production: np.array):
|
||||
df = pd.DataFrame({"Year": np.hstack((years, years)),
|
||||
"Production": np.hstack((production, [production[0]]*len(years))),
|
||||
"Scenario": itertools.chain([f"Subsidy ({np.sum(production)/1e6})"]*len(years), [f"Business as Usual ({(production[0]/1e6)*len(years)})"]*len(years))})
|
||||
df.Production /= 1e6
|
||||
|
||||
# Plot total vehicle pop
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(data=df, x="Year", y="Production", hue="Scenario")
|
||||
plt.title("Total Vehicle Production")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Simulation year")
|
||||
plt.ylim(0.8*min(np.array(production) / 1e6), 1.2*max(np.array(production) / 1e6))
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
def animate(years, population_distributions):
|
||||
ages = np.arange(len(population_distributions[0]))
|
||||
fig = plt.figure()
|
||||
camera = Camera(fig)
|
||||
# animation draws one data point at a time
|
||||
for i in range(len(population_distributions)):
|
||||
plot_ = sns.barplot(ages, population_distributions[i] / 1e6, color='cornflowerblue')
|
||||
plt.plot(ages, np.array(population_distributions[0]) / 1e6, color='blue')
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.title(f"Total Vehicle Population")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Vehicle Age")
|
||||
camera.snap()
|
||||
|
||||
anim = camera.animate(blit=False)
|
||||
anim.save('animation.mp4')
|
||||
7
demo.py
7
demo.py
|
|
@ -1,4 +1,9 @@
|
|||
from util import *
|
||||
from calculations import DEFAULT_MORTALITY_DISTRIBUTION, calculate_covid_mortality, calculate_subsidy_mortality, \
|
||||
iterate, sanity_check_for_subsidy, \
|
||||
steady_state_population
|
||||
from visualization import animate, plot_mortality, plot_population, plot_population_development, plot_production
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
if __name__ == "__main__":
|
||||
sanity_check_for_subsidy()
|
||||
|
|
|
|||
|
|
@ -0,0 +1,75 @@
|
|||
from calculations import iterate
|
||||
from calculations import steady_state_population
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
import seaborn as sns
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
|
||||
def real_mortality_distribution() -> np.ndarray:
|
||||
df = pd.read_excel("input/mortality.xlsx")
|
||||
mortality = np.array(df.Mortality)
|
||||
mortality[np.where(mortality>0)[0]] = 0
|
||||
mortality=np.hstack((mortality, mortality[-1]))
|
||||
return np.array(-1*mortality)
|
||||
|
||||
|
||||
def real_population_distribution(fraction_that_is_cars = 0.8413) -> np.ndarray:
|
||||
df = pd.read_excel("input/age distribution.xlsx")
|
||||
return np.array(df["Absolute Anzahl "])*fraction_that_is_cars
|
||||
|
||||
|
||||
def population_distribution_deltas(fraction_that_is_cars = 0.8413) -> np.ndarray:
|
||||
population_distribution = real_population_distribution()
|
||||
mortality_distribution = real_mortality_distribution()
|
||||
|
||||
df = pd.read_excel("input/new registration numbers.xlsx")
|
||||
years = np.arange(2020, 2031)
|
||||
new_registrations = np.array(df["scenario without Covid"][1:])*fraction_that_is_cars
|
||||
steady_state_registrationss = list()
|
||||
correction_factors = list()
|
||||
totals = list()
|
||||
|
||||
for i in range(len(years)):
|
||||
steady_state_registrations = iterate(real_mortality_distribution(), population_distribution.copy())[0]
|
||||
correction_factor = new_registrations[i]/steady_state_registrations
|
||||
correction_factors.append(correction_factor)
|
||||
population_distribution = iterate(real_mortality_distribution(), population_distribution, fraction_of_new_cars=correction_factor)
|
||||
totals.append(sum(population_distribution))
|
||||
steady_state_registrationss.append(steady_state_registrations)
|
||||
|
||||
df = pd.DataFrame({"Year": years,
|
||||
"Predicted Regsitrations": new_registrations,
|
||||
"Steady State Registrations": steady_state_registrationss,
|
||||
"Correction Factor": correction_factors,
|
||||
"Total": totals})
|
||||
|
||||
return df
|
||||
|
||||
if __name__ == "__main__":
|
||||
mortality_distribution = real_mortality_distribution()
|
||||
population_distribution = real_population_distribution()
|
||||
|
||||
steady_state_population_distribution = steady_state_population(mortality_distribution,
|
||||
total_population=sum(population_distribution))
|
||||
# Plot Age dtribution
|
||||
if True:
|
||||
plt.figure()
|
||||
df = pd.DataFrame({"Age": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"Data": np.hstack((population_distribution/1e6,
|
||||
steady_state_population_distribution/1e6)),
|
||||
"Scenario": np.hstack((["KBA Data"]*len(mortality_distribution),
|
||||
["Simulated Steady State"]*len(mortality_distribution)))})
|
||||
plot_ = sns.barplot(data=df, x="Age", y="Data", hue="Scenario")
|
||||
plt.title("Population Distribution")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Age (years)")
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
|
|
@ -0,0 +1,31 @@
|
|||
import pickle
|
||||
|
||||
import numpy as np
|
||||
from pyswarm import pso
|
||||
|
||||
from calculations import steady_state_population
|
||||
from load_real_data import real_mortality_distribution, real_population_distribution
|
||||
|
||||
|
||||
def optimized_mortality_distribution() -> np.ndarray:
|
||||
population_distribution = real_population_distribution()
|
||||
initial_mortality = real_mortality_distribution()
|
||||
upper_bound = 0.20 * np.ones(initial_mortality.shape)
|
||||
lower_bound = 0.01 * np.ones(initial_mortality.shape)
|
||||
|
||||
def _inner(x):
|
||||
mortality = np.array(x)
|
||||
try:
|
||||
steady_state_population_distribution = steady_state_population(mortality,
|
||||
total_population=sum(population_distribution))
|
||||
except Exception:
|
||||
return 1e100
|
||||
population_error = np.sum(np.abs(steady_state_population_distribution-population_distribution))/sum(population_distribution)
|
||||
return population_error
|
||||
|
||||
xopt, fopt = pso(_inner, lower_bound, upper_bound, swarmsize=1000, maxiter=100, debug=True)
|
||||
with open("output/swarm-mortality.pickle", "wb") as fp:
|
||||
pickle.dump(xopt, fp)
|
||||
|
||||
if __name__ == "__main__":
|
||||
optimized_mortality_distribution()
|
||||
|
|
@ -0,0 +1,253 @@
|
|||
import pickle
|
||||
|
||||
import matplotlib
|
||||
|
||||
matplotlib.use("pgf")
|
||||
matplotlib.rcParams.update({
|
||||
"pgf.texsystem": "pdflatex",
|
||||
'font.size' : 6.5,
|
||||
'font.family' : 'lmodern',
|
||||
'text.usetex': True,
|
||||
'pgf.rcfonts': False,
|
||||
})
|
||||
|
||||
import seaborn as sns
|
||||
|
||||
from calculations import calculate_covid_mortality, calculate_subsidy_mortality, iterate, steady_state_population
|
||||
import pandas as pd
|
||||
from matplotlib import pyplot as plt
|
||||
import numpy as np
|
||||
from load_real_data import population_distribution_deltas, real_mortality_distribution, real_population_distribution
|
||||
from visualization import plot_mortality
|
||||
|
||||
|
||||
def no_premium_no_covid(population_distribution, mortality_distribution):
|
||||
total_pops = list()
|
||||
distribs = list()
|
||||
production = list()
|
||||
initial_pop_distrib = population_distribution
|
||||
|
||||
for i in range(len(years)):
|
||||
fraction_of_new_cars = correction_factor[i]
|
||||
|
||||
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
|
||||
distribs.append(population_distribution.copy())
|
||||
total_pops.append(sum(population_distribution))
|
||||
production.append(population_distribution[0])
|
||||
|
||||
df = pd.DataFrame(index=np.hstack((years[0]-1,years)), data=np.vstack((initial_pop_distrib, distribs)))
|
||||
df.to_excel("output/no_premium_no_covid.xlsx")
|
||||
return df
|
||||
|
||||
def no_premium_with_covid(population_distribution, mortality_distribution):
|
||||
total_pops = list()
|
||||
distribs = list()
|
||||
production = list()
|
||||
initial_pop_distrib = population_distribution
|
||||
|
||||
for i in range(len(years)):
|
||||
fraction_of_new_cars = correction_factor[i]
|
||||
|
||||
if i == 0:
|
||||
covid_mortality = calculate_covid_mortality(population_distribution, mortality_distribution)
|
||||
population_distribution = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
|
||||
else:
|
||||
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
|
||||
|
||||
|
||||
distribs.append(population_distribution.copy())
|
||||
total_pops.append(sum(population_distribution))
|
||||
production.append(population_distribution[0])
|
||||
|
||||
df = pd.DataFrame(index=np.hstack((years[0]-1,years)), data=np.vstack((initial_pop_distrib, distribs)))
|
||||
df.to_excel("output/no_premium_with_covid.xlsx")
|
||||
return df
|
||||
|
||||
def premium_with_covid(population_distribution, mortality_distribution, subsidised_vehicle_count=1e6):
|
||||
total_pops = list()
|
||||
distribs = list()
|
||||
production = list()
|
||||
initial_pop_distrib = population_distribution
|
||||
|
||||
covid_mortality = calculate_covid_mortality(population_distribution, mortality_distribution)
|
||||
subsidy_mortality = calculate_subsidy_mortality(population_distribution, mortality_distribution, total_cars=subsidised_vehicle_count)
|
||||
|
||||
for i in range(len(years)):
|
||||
fraction_of_new_cars = correction_factor[i]
|
||||
|
||||
if i == 0:
|
||||
population_distribution = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
|
||||
elif i == 1:
|
||||
population_distribution = iterate(subsidy_mortality, population_distribution, fraction_of_new_cars)
|
||||
else:
|
||||
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
|
||||
|
||||
|
||||
distribs.append(population_distribution.copy())
|
||||
total_pops.append(sum(population_distribution))
|
||||
production.append(population_distribution[0])
|
||||
|
||||
df = pd.DataFrame(index=np.hstack((years[0]-1,years)), data=np.vstack((initial_pop_distrib, distribs)))
|
||||
df.to_excel(f"output/with_premium_with_covid_{subsidised_vehicle_count}.xlsx")
|
||||
return df
|
||||
|
||||
if __name__ == "__main__":
|
||||
mortality_distribution = real_mortality_distribution()
|
||||
#mortality_distribution = np.convolve(mortality_distribution, np.ones((3,))/3, mode='valid') # TODO: Remove
|
||||
|
||||
population_distribution = real_population_distribution()
|
||||
|
||||
steady_state_population_distribution = steady_state_population(mortality_distribution,
|
||||
total_population=sum(population_distribution))
|
||||
|
||||
# population_distribution = steady_state_population_distribution #TODO: Remove
|
||||
|
||||
years = population_distribution_deltas()["Year"]
|
||||
correction_factor = population_distribution_deltas()["Correction Factor"]
|
||||
# correction_factor = np.ones(correction_factor.shape) # TODO: Remove
|
||||
|
||||
if True:
|
||||
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
|
||||
plt.figure()
|
||||
df = pd.DataFrame({"Age": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"Data": np.hstack((population_distribution / 1e6,
|
||||
steady_state_population_distribution / 1e6)),
|
||||
"Scenario": np.hstack((["Actual Data 2019"] * len(mortality_distribution),
|
||||
["Using Normal Deregistration"] * len(mortality_distribution)))})
|
||||
plot_ = sns.lineplot(data=df, x="Age", y="Data", hue="Scenario")
|
||||
plt.title("In-Scope Population Distribution")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Age (years)")
|
||||
#for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
# if ind % 10 == 0: # every 10th label is kept
|
||||
# label.set_visible(True)
|
||||
# else:
|
||||
# label.set_visible(False)
|
||||
plt.savefig("output/population.png")
|
||||
plt.savefig("output/population.pgf")
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
production = list()
|
||||
year = list()
|
||||
title = list()
|
||||
|
||||
total_production = list()
|
||||
total_title = list()
|
||||
|
||||
df = no_premium_no_covid(population_distribution.copy(), mortality_distribution)
|
||||
production.extend(df[0]/1e6)
|
||||
year.extend(df.index)
|
||||
title.extend(["No Covid, No Premium"]*len(df.index))
|
||||
|
||||
total_production.append(np.sum(df[0][1:]/1e6))
|
||||
total_title.append("No Covid,\nNo Premium")
|
||||
|
||||
df = no_premium_with_covid(population_distribution.copy(), mortality_distribution)
|
||||
production.extend(df[0]/1e6)
|
||||
year.extend(df.index)
|
||||
title.extend(["With Covid, No Premium"]*len(df.index))
|
||||
|
||||
total_production.append(np.sum(df[0][1:]/1e6))
|
||||
total_title.append("With Covid,\nNo Premium")
|
||||
|
||||
subsidy_count = production[1]*1e6-df[0][2020]
|
||||
df = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
|
||||
production.extend(df[0]/1e6)
|
||||
year.extend(df.index)
|
||||
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
|
||||
|
||||
total_production.append(np.sum(df[0][1:]/1e6))
|
||||
total_title.append(f"Covid,\nPremium\n({subsidy_count/1e6:.2f}M\nvehicles)")
|
||||
|
||||
subsidy_count*=2
|
||||
df = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
|
||||
production.extend(df[0]/1e6)
|
||||
year.extend(df.index)
|
||||
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
|
||||
|
||||
total_production.append(np.sum(df[0][1:]/1e6))
|
||||
total_title.append(f"Covid,\nPremium\n({subsidy_count/1e6:.2f}M\nvehicles)")
|
||||
|
||||
subsidy_count*=2
|
||||
df = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
|
||||
production.extend(df[0]/1e6)
|
||||
year.extend(df.index)
|
||||
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
|
||||
|
||||
total_production.append(np.sum(df[0][1:]/1e6))
|
||||
total_title.append(f"Covid,\nPremium\n({subsidy_count/1e6:.2f}M\nvehicles)")
|
||||
|
||||
if True:
|
||||
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
|
||||
plt.figure(figsize=(6.5, 4))
|
||||
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"Case": np.hstack((["No Covid, No Premium"]*len(mortality_distribution),
|
||||
["With Covid, No Premium"]*len(mortality_distribution),
|
||||
[f"Covid, Premium ({subsidy_count/4e6:.2f}M vehicles)"]*len(mortality_distribution),
|
||||
[f"Covid, Premium ({subsidy_count/2e6:.2f}M vehicles)"]*len(mortality_distribution),
|
||||
[f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(mortality_distribution))),
|
||||
"Scenario": np.hstack((mortality_distribution,
|
||||
calculate_covid_mortality(population_distribution, mortality_distribution),
|
||||
calculate_subsidy_mortality(population_distribution, mortality_distribution, total_cars=subsidy_count/4),
|
||||
calculate_subsidy_mortality(population_distribution, mortality_distribution, total_cars=subsidy_count/2),
|
||||
calculate_subsidy_mortality(population_distribution, mortality_distribution, total_cars=subsidy_count)))})
|
||||
plot_ = sns.lineplot(data=df, x="Ages", y="Scenario", hue="Case")
|
||||
plot_.grid(b=True)
|
||||
#plt.title("Deregistration Distribution")
|
||||
plt.ylabel("Recycling probability")
|
||||
plt.xlabel("Age (years)")
|
||||
# for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
# if ind % 10 == 0: # every 10th label is kept
|
||||
# label.set_visible(True)
|
||||
# else:
|
||||
# label.set_visible(False)
|
||||
plt.savefig("output/mortality.png")
|
||||
plt.savefig("output/mortality.pgf")
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
if True:
|
||||
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
|
||||
plt.figure(figsize=(6.5, 2.5))
|
||||
plt.subplot(121)
|
||||
df = pd.DataFrame({"Year": year,
|
||||
"Data": production,
|
||||
"Scenario": title})
|
||||
plot_ = sns.lineplot(data=df, x="Year", y="Data", hue="Scenario")
|
||||
plot_.legend().set_visible(False)
|
||||
plot_.grid(b=True)
|
||||
plt.title("In-Scope Vehice Registrations per Year")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Year")
|
||||
#for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
# if ind % 10 == 0: # every 10th label is kept
|
||||
# label.set_visible(True)
|
||||
# else:
|
||||
# label.set_visible(False)
|
||||
|
||||
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
|
||||
plt.subplot(122)
|
||||
df = pd.DataFrame({"Data": total_production,
|
||||
"Scenario": total_title})
|
||||
plot_ = sns.barplot(data=df, x="Scenario", y="Data")
|
||||
plot_.grid(b=True, axis='y')
|
||||
plt.ylim(min(np.array(total_production))*0.99, max(np.array(total_production))*1.01)
|
||||
plt.title("In-Scope Vehicle Registrations 2020-2030")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
#plt.xlabel("Sce")
|
||||
#for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
# if ind % 10 == 0: # every 10th label is kept
|
||||
# label.set_visible(True)
|
||||
# else:
|
||||
# label.set_visible(False)
|
||||
plt.tight_layout()
|
||||
plt.savefig("output/productions.png")
|
||||
plt.savefig("output/productions.pgf")
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
|
@ -1,71 +1 @@
|
|||
from util import *
|
||||
|
||||
def real_mortality_distribution() -> np.ndarray:
|
||||
df = pd.read_excel("mortality.xlsx")
|
||||
mortality = df.Mortality
|
||||
mortality[0] = 0
|
||||
mortality[1] = 0
|
||||
mortality[2] = 0
|
||||
mortality[32:38] = mortality[31]
|
||||
mnp.hstack((mortality, mortality[31]))
|
||||
return np.array(-1*mortality)
|
||||
|
||||
def real_population_distribution() -> np.ndarray:
|
||||
df = pd.read_excel("age distribution.xlsx")
|
||||
return np.array(df["Absolute Anzahl "])
|
||||
|
||||
def population_distribution_deltas() -> np.ndarray:
|
||||
population_distribution = real_population_distribution()
|
||||
mortality_distribution = real_mortality_distribution()
|
||||
|
||||
df = pd.read_excel("new registration numbers.xlsx")
|
||||
years = np.arange(2020, 2031)
|
||||
new_registrations = np.array(df["scenario without Covid"][1:])
|
||||
steady_state_registrationss = list()
|
||||
correction_factors = list()
|
||||
totals = list()
|
||||
|
||||
for i in range(len(years)):
|
||||
steady_state_registrations = iterate(real_mortality_distribution(), population_distribution.copy())[0]
|
||||
correction_factor = new_registrations[i]/steady_state_registrations
|
||||
correction_factors.append(correction_factor)
|
||||
population_distribution = iterate(real_mortality_distribution(), population_distribution, fraction_of_new_cars=correction_factor)
|
||||
totals.append(sum(population_distribution))
|
||||
steady_state_registrationss.append(steady_state_registrations)
|
||||
|
||||
df = pd.DataFrame({"Year": years,
|
||||
"Predicted Regsitrations": new_registrations,
|
||||
"Steady State Registrations": steady_state_registrationss,
|
||||
"Correction Factor": correction_factors,
|
||||
"Total": totals})
|
||||
|
||||
return correction_factors
|
||||
|
||||
if __name__ == "__main__":
|
||||
mortality_distribution = real_mortality_distribution()
|
||||
population_distribution = real_population_distribution()
|
||||
|
||||
steady_state_population_distribution = steady_state_population(mortality_distribution,
|
||||
total_population=sum(population_distribution))
|
||||
|
||||
if True:
|
||||
plt.figure()
|
||||
df = pd.DataFrame({"Age": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"Data": np.hstack((population_distribution/1e6,
|
||||
steady_state_population_distribution/1e6)),
|
||||
"Scenario": np.hstack((["KBA Data"]*len(mortality_distribution),
|
||||
["Simulated Steady State"]*len(mortality_distribution)))})
|
||||
plot_ = sns.barplot(data=df, x="Age", y="Data", hue="Scenario")
|
||||
plt.title("Population Distribution")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Age (years)")
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
population_distribution_deltas()
|
||||
|
|
@ -0,0 +1,122 @@
|
|||
import itertools
|
||||
import logging
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import seaborn as sns
|
||||
from celluloid import Camera
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
from calculations import DEFAULT_MORTALITY_DISTRIBUTION
|
||||
|
||||
|
||||
def plot_mortality(mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION, subsidy_distribution = None,
|
||||
covid_distribution = None):
|
||||
logging.basicConfig(level=logging.DEBUG)
|
||||
|
||||
plt.figure()
|
||||
if subsidy_distribution is None:
|
||||
plot_ = sns.barplot(np.arange(len(mortality_distribution)), mortality_distribution, color='cornflowerblue')
|
||||
elif covid_distribution is None:
|
||||
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)), np.arange(len(mortality_distribution)))),
|
||||
"label": list(itertools.chain(["Normal"]*len(mortality_distribution), ["With Subsidy"]*len(mortality_distribution))),
|
||||
"Legend": np.hstack((mortality_distribution, subsidy_distribution))})
|
||||
else:
|
||||
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)),
|
||||
np.arange(len(mortality_distribution)))),
|
||||
"Case": list(itertools.chain(["Normal"] * len(mortality_distribution),
|
||||
["With Subsidy"] * len(mortality_distribution),
|
||||
["With Covid"] * len(mortality_distribution))),
|
||||
"Scenario": np.hstack((mortality_distribution,
|
||||
subsidy_distribution,
|
||||
covid_distribution))})
|
||||
plot_ = sns.lineplot(data=df, x="Ages", y="Scenario", hue="Case")
|
||||
plt.title("Deregistration Distribution")
|
||||
plt.ylabel("Recycling probability")
|
||||
plt.xlabel("Age (years)")
|
||||
#for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
# if ind % 10 == 0: # every 10th label is kept
|
||||
# label.set_visible(True)
|
||||
# else:
|
||||
# label.set_visible(False)
|
||||
plt.savefig("output/mortality.png")
|
||||
plt.savefig("output/mortality.pgf")
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
|
||||
def plot_population(population_distribution: np.ndarray):
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(np.arange(len(population_distribution)), population_distribution / 1e6, color='cornflowerblue')
|
||||
plt.title("Population Distribution")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Age (years)")
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
|
||||
def plot_population_development(years, total_pops: np.array):
|
||||
# Plot total vehicle pop
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(years, np.array(total_pops) / 1e6, color='cornflowerblue')
|
||||
plt.title("Total Vehicle Population")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Simulation year")
|
||||
plt.ylim(min(np.array(total_pops) / 1e6), max(np.array(total_pops) / 1e6))
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
|
||||
def plot_production(years, production: np.array):
|
||||
df = pd.DataFrame({"Year": np.hstack((years, years)),
|
||||
"Production": np.hstack((production, [production[0]]*len(years))),
|
||||
"Scenario": itertools.chain([f"Subsidy ({np.sum(production)/1e6})"]*len(years), [f"Business as Usual ({(production[0]/1e6)*len(years)})"]*len(years))})
|
||||
df.Production /= 1e6
|
||||
|
||||
# Plot total vehicle pop
|
||||
plt.figure()
|
||||
plot_ = sns.barplot(data=df, x="Year", y="Production", hue="Scenario")
|
||||
plt.title("Total Vehicle Production")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Simulation year")
|
||||
plt.ylim(0.8*min(np.array(production) / 1e6), 1.2*max(np.array(production) / 1e6))
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.show()
|
||||
plt.close()
|
||||
|
||||
|
||||
def animate(years, population_distributions):
|
||||
ages = np.arange(len(population_distributions[0]))
|
||||
fig = plt.figure()
|
||||
camera = Camera(fig)
|
||||
# animation draws one data point at a time
|
||||
for i in range(len(population_distributions)):
|
||||
plot_ = sns.barplot(ages, population_distributions[i] / 1e6, color='cornflowerblue')
|
||||
plt.plot(ages, np.array(population_distributions[0]) / 1e6, color='blue')
|
||||
for ind, label in enumerate(plot_.get_xticklabels()):
|
||||
if ind % 10 == 0: # every 10th label is kept
|
||||
label.set_visible(True)
|
||||
else:
|
||||
label.set_visible(False)
|
||||
plt.title(f"Total Vehicle Population")
|
||||
plt.ylabel("Vehicles (Millions)")
|
||||
plt.xlabel("Vehicle Age")
|
||||
camera.snap()
|
||||
|
||||
anim = camera.animate(blit=False)
|
||||
anim.save('animation.mp4')
|
||||
Loading…
Reference in New Issue
Block a user