Finaler for Paper

This commit is contained in:
Ludger Heide 2020-10-15 16:16:28 +02:00
parent eb0b99e969
commit ed9ed964b5
5 changed files with 133 additions and 37 deletions

View File

@ -88,10 +88,12 @@ def calculate_covid_mortality(steady_state_population: np.array, normal_mortalit
def iterate(mortality_distribution: np.array, population_distribution: np.array, fraction_of_new_cars: float = 1.0):
old_population = sum(population_distribution)
mortality_totals = mortality_distribution * population_distribution
population_distribution -= mortality_totals
# First, simulate the mortality
for j in reversed(range(len(mortality_distribution) - 1)):
population_distribution[j + 1] = (1 - mortality_distribution[j]) * population_distribution[j]
population_distribution[j + 1] = population_distribution[j]
# Last, add new vehicles in order for the sum to make sense
population_distribution[0] = 0
@ -99,7 +101,9 @@ def iterate(mortality_distribution: np.array, population_distribution: np.array,
new_cars = old_population - population_before_creation
population_distribution[0] = new_cars * fraction_of_new_cars
return population_distribution
# assert new_cars == sum(mortality_totals)
return population_distribution, mortality_totals
def sanity_check_for_subsidy():

Binary file not shown.

Binary file not shown.

View File

@ -6,12 +6,12 @@ import seaborn as sns
from matplotlib import pyplot as plt
def real_mortality_distribution() -> np.ndarray:
def real_mortality_distribution(scaling_factor = 1.51) -> 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)
return np.array(-1*mortality)*scaling_factor
def real_population_distribution(fraction_that_is_cars = 0.8413) -> np.ndarray:
@ -31,10 +31,10 @@ def population_distribution_deltas(fraction_that_is_cars = 0.8413) -> np.ndarray
totals = list()
for i in range(len(years)):
steady_state_registrations = iterate(real_mortality_distribution(), population_distribution.copy())[0]
steady_state_registrations = iterate(real_mortality_distribution(), population_distribution.copy())[0][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)
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)

View File

@ -5,7 +5,7 @@ import matplotlib
matplotlib.use("pgf")
matplotlib.rcParams.update({
"pgf.texsystem": "pdflatex",
'font.size' : 6.5,
#'font.size' : 6.5,
'font.family' : 'lmodern',
'text.usetex': True,
'pgf.rcfonts': False,
@ -15,7 +15,7 @@ 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
from matplotlib import pyplot as plt, ticker
import numpy as np
from load_real_data import population_distribution_deltas, real_mortality_distribution, real_population_distribution
from visualization import plot_mortality
@ -24,50 +24,63 @@ from visualization import plot_mortality
def no_premium_no_covid(population_distribution, mortality_distribution):
total_pops = list()
distribs = list()
mortalities = list()
production = list()
initial_pop_distrib = population_distribution
initial_pop_distrib = population_distribution.copy()
for i in range(len(years)):
fraction_of_new_cars = correction_factor[i]
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
population_distribution, mortality = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
mortalities.append(mortality.copy())
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
df2 = pd.DataFrame(index=np.hstack((years, years.iloc[-1] + 1)), data=np.vstack((mortalities, mortality_distribution*population_distribution)))
df2.to_excel("output/no_premium_no_covid_mortalities.xlsx")
return df, df2
def no_premium_with_covid(population_distribution, mortality_distribution):
total_pops = list()
distribs = list()
mortalities = list()
production = list()
initial_pop_distrib = population_distribution
initial_pop_distrib = population_distribution.copy()
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)
population_distribution, mortality = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
else:
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
population_distribution, mortality = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
mortalities.append(mortality.copy())
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 = 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
df2 = pd.DataFrame(index=np.hstack((years, years.iloc[-1] + 1)), data=np.vstack((mortalities, mortality_distribution*population_distribution)))
df2.to_excel("output/no_premium_with_covid_mortalities.xlsx")
return df, df2
def premium_with_covid(population_distribution, mortality_distribution, subsidised_vehicle_count=1e6):
total_pops = list()
distribs = list()
mortalities = list()
production = list()
initial_pop_distrib = population_distribution
initial_pop_distrib = population_distribution.copy()
covid_mortality = calculate_covid_mortality(population_distribution, mortality_distribution)
subsidy_mortality = calculate_subsidy_mortality(population_distribution, mortality_distribution, total_cars=subsidised_vehicle_count)
@ -76,20 +89,25 @@ def premium_with_covid(population_distribution, mortality_distribution, subsidis
fraction_of_new_cars = correction_factor[i]
if i == 0:
population_distribution = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
population_distribution, mortality = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
elif i == 1:
population_distribution = iterate(subsidy_mortality, population_distribution, fraction_of_new_cars)
population_distribution, mortality = iterate(subsidy_mortality, population_distribution, fraction_of_new_cars)
else:
population_distribution = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
population_distribution, mortality = iterate(mortality_distribution, population_distribution, fraction_of_new_cars)
mortalities.append(mortality.copy())
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
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/1e6:.2f}.xlsx")
df2 = pd.DataFrame(index=np.hstack((years, years.iloc[-1] + 1)), data=np.vstack((mortalities, mortality_distribution*population_distribution)))
df2.to_excel(f"output/with_premium_with_covid_mortalities_{subsidised_vehicle_count/1e6:.2f}.xlsx")
return df, df2
if __name__ == "__main__":
mortality_distribution = real_mortality_distribution()
@ -104,7 +122,7 @@ if __name__ == "__main__":
years = population_distribution_deltas()["Year"]
correction_factor = population_distribution_deltas()["Correction Factor"]
# correction_factor = np.ones(correction_factor.shape) # TODO: Remove
correction_factor = np.ones(correction_factor.shape) # TODO: Remove
if True:
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
@ -129,32 +147,42 @@ if __name__ == "__main__":
plt.show()
plt.close()
total_amount = list()
production = list()
year = list()
title = list()
total_production = list()
total_title = list()
total_mortalities = list()
df = no_premium_no_covid(population_distribution.copy(), mortality_distribution)
df, df2 = no_premium_no_covid(population_distribution.copy(), mortality_distribution)
no_covid_no_premium_production = df[0][2020]
normal_production = df[0]
production.extend(df[0]/1e6)
total_amount.extend(np.sum(df.values, axis=1))
total_mortalities.extend(np.sum(df2.values, axis=1))
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)
df, df2 = no_premium_with_covid(population_distribution.copy(), mortality_distribution)
production.extend(df[0]/1e6)
total_amount.extend(np.sum(df.values, axis=1))
total_mortalities.extend(np.sum(df2.values, axis=1))
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)
subsidy_count = no_covid_no_premium_production - df[0][2020]
df, df2 = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
production.extend(df[0]/1e6)
total_amount.extend(np.sum(df.values, axis=1))
total_mortalities.extend(np.sum(df2.values, axis=1))
year.extend(df.index)
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
@ -162,8 +190,10 @@ if __name__ == "__main__":
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)
df, df2 = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
production.extend(df[0]/1e6)
total_amount.extend(np.sum(df.values, axis=1))
total_mortalities.extend(np.sum(df2.values, axis=1))
year.extend(df.index)
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
@ -171,8 +201,10 @@ if __name__ == "__main__":
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)
df, df2 = premium_with_covid(population_distribution.copy(), mortality_distribution, subsidised_vehicle_count=subsidy_count)
production.extend(df[0]/1e6)
total_amount.extend(np.sum(df.values, axis=1))
total_mortalities.extend(np.sum(df2.values, axis=1))
year.extend(df.index)
title.extend([f"Covid, Premium ({subsidy_count/1e6:.2f}M vehicles)"]*len(df.index))
@ -181,7 +213,7 @@ if __name__ == "__main__":
if True:
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
plt.figure(figsize=(6.5, 4))
plt.figure(figsize=(6.5, 3.5))
df = pd.DataFrame({"Ages": np.hstack((np.arange(len(mortality_distribution)),
np.arange(len(mortality_distribution)),
np.arange(len(mortality_distribution)),
@ -199,9 +231,12 @@ if __name__ == "__main__":
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.legend(fontsize=6.5)
plt.gcf().axes[0].set_axisbelow(True)
#plt.title("Deregistration Distribution")
plt.ylabel("Recycling probability")
plt.xlabel("Age (years)")
plt.tight_layout()
# for ind, label in enumerate(plot_.get_xticklabels()):
# if ind % 10 == 0: # every 10th label is kept
# label.set_visible(True)
@ -214,15 +249,20 @@ if __name__ == "__main__":
if True:
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
plt.figure(figsize=(6.5, 2.5))
plt.figure(figsize=(6, 3.5))
plt.subplot(121)
df = pd.DataFrame({"Year": year,
df = pd.DataFrame({"Year": np.array(year),
"Data": production,
"Scenario": title})
df = df.drop(np.where(df.index % 12 == 0)[0])
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")
plot_.xaxis.set_major_locator(ticker.MultipleLocator(2))
plot_.xaxis.set_major_formatter(ticker.ScalarFormatter())
plt.legend(fontsize=6.5)
plt.gcf().axes[0].set_axisbelow(True)
plt.title("Per Year")
plt.ylabel("Vehicles (Millions)")
plt.xlabel("Year")
#for ind, label in enumerate(plot_.get_xticklabels()):
@ -235,11 +275,16 @@ if __name__ == "__main__":
plt.subplot(122)
df = pd.DataFrame({"Data": total_production,
"Scenario": total_title})
plot_ = sns.barplot(data=df, x="Scenario", y="Data")
plot_ = sns.barplot(data=df, x="Scenario", y="Data",edgecolor='black',linewidth=0.8)
plot_.grid(b=True, axis='y')
#plot_.set_xticklabels(plt.xticks()[1], size=6.5)
plt.xticks([], [])
plt.gcf().axes[1].set_axisbelow(True)
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.title("Total 2020-2030")
plt.ylabel("Vehicles (Millions)")
#plt.legend(fontsize=6.5)
plt.suptitle("In-Scope Vehicle Registrations")
#plt.xlabel("Sce")
#for ind, label in enumerate(plot_.get_xticklabels()):
# if ind % 10 == 0: # every 10th label is kept
@ -250,4 +295,51 @@ if __name__ == "__main__":
plt.savefig("output/productions.png")
plt.savefig("output/productions.pgf")
plt.show()
plt.close()
plt.close()
if True:
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
df = pd.DataFrame({"Year": year,
"Data": total_amount,
"Scenario": title})
plot_ = sns.lineplot(data=df, x="Year", y="Data", hue="Scenario")
plot_.legend().set_visible(False)
plot_.grid(b=True)
plt.ylim(min(np.array(total_amount))*0.99, max(np.array(total_amount))*1.01)
plt.title("Total vehicle count")
plt.ylabel("Vehicles (Millions)")
plt.xlabel("Year")
plt.tight_layout()
plt.show()
plt.savefig("output/amount.png")
plt.savefig("output/amount.pgf")
plt.close()
# 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)
if True:
sns.set_palette(sns.color_palette(("grey", "firebrick", "peru", "darkorange", "gold")))
df = pd.DataFrame({"Year": np.hstack((df2.index, df2.index,df2.index,df2.index,df2.index)),
"Data": np.array(total_mortalities)/1e6,
"Scenario": title})
plot_ = sns.lineplot(data=df, x="Year", y="Data", hue="Scenario")
plot_.legend().set_visible(False)
plot_.grid(b=True)
plt.title("Total recycled vehicles")
plt.ylabel("Vehicles (Millions)")
plt.xlabel("Year")
plt.tight_layout()
plt.show()
plt.savefig("output/recycled.png")
plt.savefig("output/amount.pgf")
plt.close()
# 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)
print(f"Production scaling factor: {(3455546*0.8413)/np.mean(normal_production[1:])}")