345 lines
16 KiB
Python
Executable File
345 lines
16 KiB
Python
Executable File
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, ticker
|
|
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()
|
|
mortalities = list()
|
|
production = list()
|
|
initial_pop_distrib = population_distribution.copy()
|
|
|
|
for i in range(len(years)):
|
|
fraction_of_new_cars = correction_factor[i]
|
|
|
|
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")
|
|
|
|
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.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, mortality = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
|
|
else:
|
|
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_with_covid.xlsx")
|
|
|
|
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.copy()
|
|
|
|
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, mortality = iterate(covid_mortality, population_distribution, fraction_of_new_cars)
|
|
elif i == 1:
|
|
population_distribution, mortality = iterate(subsidy_mortality, population_distribution, fraction_of_new_cars)
|
|
else:
|
|
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/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()
|
|
#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()
|
|
|
|
total_amount = list()
|
|
production = list()
|
|
year = list()
|
|
title = list()
|
|
|
|
total_production = list()
|
|
total_title = list()
|
|
total_mortalities = list()
|
|
|
|
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, 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 = 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))
|
|
|
|
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, 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))
|
|
|
|
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, 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))
|
|
|
|
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, 3.5))
|
|
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.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)
|
|
# 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, 3.5))
|
|
plt.subplot(121)
|
|
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)
|
|
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()):
|
|
# 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",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("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
|
|
# 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()
|
|
|
|
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:])}") |