Initial commit
This commit is contained in:
commit
e0df453875
|
|
@ -0,0 +1,241 @@
|
|||
# General
|
||||
.DS_Store
|
||||
.AppleDouble
|
||||
.LSOverride
|
||||
|
||||
# Icon must end with two \r
|
||||
Icon
|
||||
|
||||
|
||||
# Thumbnails
|
||||
._*
|
||||
|
||||
# Files that might appear in the root of a volume
|
||||
.DocumentRevisions-V100
|
||||
.fseventsd
|
||||
.Spotlight-V100
|
||||
.TemporaryItems
|
||||
.Trashes
|
||||
.VolumeIcon.icns
|
||||
.com.apple.timemachine.donotpresent
|
||||
|
||||
# Directories potentially created on remote AFP share
|
||||
.AppleDB
|
||||
.AppleDesktop
|
||||
Network Trash Folder
|
||||
Temporary Items
|
||||
.apdisk
|
||||
|
||||
# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio, WebStorm and Rider
|
||||
# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839
|
||||
|
||||
# User-specific stuff
|
||||
.idea/**/workspace.xml
|
||||
.idea/**/tasks.xml
|
||||
.idea/**/usage.statistics.xml
|
||||
.idea/**/dictionaries
|
||||
.idea/**/shelf
|
||||
|
||||
# Generated files
|
||||
.idea/**/contentModel.xml
|
||||
|
||||
# Sensitive or high-churn files
|
||||
.idea/**/dataSources/
|
||||
.idea/**/dataSources.ids
|
||||
.idea/**/dataSources.local.xml
|
||||
.idea/**/sqlDataSources.xml
|
||||
.idea/**/dynamic.xml
|
||||
.idea/**/uiDesigner.xml
|
||||
.idea/**/dbnavigator.xml
|
||||
|
||||
# Gradle
|
||||
.idea/**/gradle.xml
|
||||
.idea/**/libraries
|
||||
|
||||
# Gradle and Maven with auto-import
|
||||
# When using Gradle or Maven with auto-import, you should exclude module files,
|
||||
# since they will be recreated, and may cause churn. Uncomment if using
|
||||
# auto-import.
|
||||
# .idea/artifacts
|
||||
# .idea/compiler.xml
|
||||
# .idea/jarRepositories.xml
|
||||
# .idea/modules.xml
|
||||
# .idea/*.iml
|
||||
# .idea/modules
|
||||
# *.iml
|
||||
# *.ipr
|
||||
|
||||
# CMake
|
||||
cmake-build-*/
|
||||
|
||||
# Mongo Explorer plugin
|
||||
.idea/**/mongoSettings.xml
|
||||
|
||||
# File-based project format
|
||||
*.iws
|
||||
|
||||
# IntelliJ
|
||||
out/
|
||||
|
||||
# mpeltonen/sbt-idea plugin
|
||||
.idea_modules/
|
||||
|
||||
# JIRA plugin
|
||||
atlassian-ide-plugin.xml
|
||||
|
||||
# Cursive Clojure plugin
|
||||
.idea/replstate.xml
|
||||
|
||||
# Crashlytics plugin (for Android Studio and IntelliJ)
|
||||
com_crashlytics_export_strings.xml
|
||||
crashlytics.properties
|
||||
crashlytics-build.properties
|
||||
fabric.properties
|
||||
|
||||
# Editor-based Rest Client
|
||||
.idea/httpRequests
|
||||
|
||||
# Android studio 3.1+ serialized cache file
|
||||
.idea/caches/build_file_checksums.ser
|
||||
|
||||
# Byte-compiled / optimized / DLL files
|
||||
__pycache__/
|
||||
*.py[cod]
|
||||
*$py.class
|
||||
|
||||
# C extensions
|
||||
*.so
|
||||
|
||||
# Distribution / packaging
|
||||
.Python
|
||||
build/
|
||||
develop-eggs/
|
||||
dist/
|
||||
downloads/
|
||||
eggs/
|
||||
.eggs/
|
||||
lib/
|
||||
lib64/
|
||||
parts/
|
||||
sdist/
|
||||
var/
|
||||
wheels/
|
||||
share/python-wheels/
|
||||
*.egg-info/
|
||||
.installed.cfg
|
||||
*.egg
|
||||
MANIFEST
|
||||
|
||||
# PyInstaller
|
||||
# Usually these files are written by a python script from a template
|
||||
# before PyInstaller builds the exe, so as to inject date/other infos into it.
|
||||
*.manifest
|
||||
*.spec
|
||||
|
||||
# Installer logs
|
||||
pip-log.txt
|
||||
pip-delete-this-directory.txt
|
||||
|
||||
# Unit test / coverage reports
|
||||
htmlcov/
|
||||
.tox/
|
||||
.nox/
|
||||
.coverage
|
||||
.coverage.*
|
||||
.cache
|
||||
nosetests.xml
|
||||
coverage.xml
|
||||
*.cover
|
||||
*.py,cover
|
||||
.hypothesis/
|
||||
.pytest_cache/
|
||||
cover/
|
||||
|
||||
# Translations
|
||||
*.mo
|
||||
*.pot
|
||||
|
||||
# Django stuff:
|
||||
*.log
|
||||
local_settings.py
|
||||
db.sqlite3
|
||||
db.sqlite3-journal
|
||||
|
||||
# Flask stuff:
|
||||
instance/
|
||||
.webassets-cache
|
||||
|
||||
# Scrapy stuff:
|
||||
.scrapy
|
||||
|
||||
# Sphinx documentation
|
||||
docs/_build/
|
||||
|
||||
# PyBuilder
|
||||
.pybuilder/
|
||||
target/
|
||||
|
||||
# Jupyter Notebook
|
||||
.ipynb_checkpoints
|
||||
|
||||
# IPython
|
||||
profile_default/
|
||||
ipython_config.py
|
||||
|
||||
# pyenv
|
||||
# For a library or package, you might want to ignore these files since the code is
|
||||
# intended to run in multiple environments; otherwise, check them in:
|
||||
# .python-version
|
||||
|
||||
# pipenv
|
||||
# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
|
||||
# However, in case of collaboration, if having platform-specific dependencies or dependencies
|
||||
# having no cross-platform support, pipenv may install dependencies that don't work, or not
|
||||
# install all needed dependencies.
|
||||
#Pipfile.lock
|
||||
|
||||
# PEP 582; used by e.g. github.com/David-OConnor/pyflow
|
||||
__pypackages__/
|
||||
|
||||
# Celery stuff
|
||||
celerybeat-schedule
|
||||
celerybeat.pid
|
||||
|
||||
# SageMath parsed files
|
||||
*.sage.py
|
||||
|
||||
# Environments
|
||||
.env
|
||||
.venv
|
||||
env/
|
||||
venv/
|
||||
ENV/
|
||||
env.bak/
|
||||
venv.bak/
|
||||
|
||||
# Spyder project settings
|
||||
.spyderproject
|
||||
.spyproject
|
||||
|
||||
# Rope project settings
|
||||
.ropeproject
|
||||
|
||||
# mkdocs documentation
|
||||
/site
|
||||
|
||||
# mypy
|
||||
.mypy_cache/
|
||||
.dmypy.json
|
||||
dmypy.json
|
||||
|
||||
# Pyre type checker
|
||||
.pyre/
|
||||
|
||||
# pytype static type analyzer
|
||||
.pytype/
|
||||
|
||||
# Cython debug symbols
|
||||
cython_debug/
|
||||
|
||||
out.xlsx
|
||||
animation.mp4
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
# Default ignored files
|
||||
/shelf/
|
||||
/workspace.xml
|
||||
# Datasource local storage ignored files
|
||||
/dataSources/
|
||||
/dataSources.local.xml
|
||||
# Editor-based HTTP Client requests
|
||||
/httpRequests/
|
||||
|
|
@ -0,0 +1,10 @@
|
|||
<component name="InspectionProjectProfileManager">
|
||||
<profile version="1.0">
|
||||
<option name="myName" value="Project Default" />
|
||||
<inspection_tool class="SpellCheckingInspection" enabled="false" level="TYPO" enabled_by_default="false">
|
||||
<option name="processCode" value="true" />
|
||||
<option name="processLiterals" value="true" />
|
||||
<option name="processComments" value="true" />
|
||||
</inspection_tool>
|
||||
</profile>
|
||||
</component>
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
<component name="InspectionProjectProfileManager">
|
||||
<settings>
|
||||
<option name="USE_PROJECT_PROFILE" value="false" />
|
||||
<version value="1.0" />
|
||||
</settings>
|
||||
</component>
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project version="4">
|
||||
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (venv)" project-jdk-type="Python SDK" />
|
||||
</project>
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project version="4">
|
||||
<component name="ProjectModuleManager">
|
||||
<modules>
|
||||
<module fileurl="file://$PROJECT_DIR$/.idea/verschrottungswahrscheinlichkeit.iml" filepath="$PROJECT_DIR$/.idea/verschrottungswahrscheinlichkeit.iml" />
|
||||
</modules>
|
||||
</component>
|
||||
</project>
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project version="4">
|
||||
<component name="VcsDirectoryMappings">
|
||||
<mapping directory="$PROJECT_DIR$/../.." vcs="Git" />
|
||||
<mapping directory="$PROJECT_DIR$" vcs="Git" />
|
||||
</component>
|
||||
</project>
|
||||
|
|
@ -0,0 +1,15 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<module type="PYTHON_MODULE" version="4">
|
||||
<component name="NewModuleRootManager">
|
||||
<content url="file://$MODULE_DIR$" />
|
||||
<orderEntry type="inheritedJdk" />
|
||||
<orderEntry type="sourceFolder" forTests="false" />
|
||||
</component>
|
||||
<component name="PackageRequirementsSettings">
|
||||
<option name="requirementsPath" value="" />
|
||||
</component>
|
||||
<component name="PyDocumentationSettings">
|
||||
<option name="format" value="GOOGLE" />
|
||||
<option name="myDocStringFormat" value="Google" />
|
||||
</component>
|
||||
</module>
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
from util import *
|
||||
|
||||
if __name__ == "__main__":
|
||||
sanity_check_for_subsidy()
|
||||
|
||||
population_distribution = steady_state_population()
|
||||
plot_population(population_distribution)
|
||||
|
||||
subsidy_mortality = calculate_subsidy_mortality(population_distribution, DEFAULT_MORTALITY_DISTRIBUTION, total_cars=1e6)
|
||||
covid_mortality = calculate_covid_mortality(population_distribution, DEFAULT_MORTALITY_DISTRIBUTION)
|
||||
plot_mortality(DEFAULT_MORTALITY_DISTRIBUTION, subsidy_mortality, covid_mortality)
|
||||
|
||||
total_pops = list()
|
||||
distribs = list()
|
||||
production = list()
|
||||
duration = 50
|
||||
start = 2010
|
||||
years = np.arange(start, start+duration)
|
||||
|
||||
for i in range(duration):
|
||||
if i < 10:
|
||||
fraction_of_new_cars = 1.0
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION
|
||||
elif i == 10:
|
||||
fraction_of_new_cars = 1.00
|
||||
mortality_distribution = covid_mortality
|
||||
elif i == 11:
|
||||
fraction_of_new_cars = 1.00
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION# subsidy_mortality
|
||||
elif i > 11 and i <= 20:
|
||||
fraction_of_new_cars = 1.00
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION
|
||||
elif i > 20:
|
||||
fraction_of_new_cars = 1.0
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION
|
||||
else:
|
||||
fraction_of_new_cars = 1.0
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION
|
||||
|
||||
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=years, data=distribs)
|
||||
df.to_excel("out.xlsx")
|
||||
|
||||
plot_population_development(years, total_pops)
|
||||
plot_production(years, production)
|
||||
animate(years, distribs)
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,222 @@
|
|||
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)
|
||||
DEFAULT_TOTAL_POPULATION = 40e6
|
||||
|
||||
|
||||
def steady_state_population(mortality_distribution: np.array = None, total_population: int = DEFAULT_TOTAL_POPULATION):
|
||||
"""
|
||||
## Establishing the initial population
|
||||
|
||||
We use an array ranging from 0 to 50 years of age. We start with an equal distribution and run the simulation until
|
||||
it arrives at the steady-state solution. _In reality, we'd use the real data here_.
|
||||
|
||||
Args:
|
||||
mortality_distribution: The chance a vehicle that is _x_ years old leaves the population
|
||||
total_population: The total vehicle population
|
||||
|
||||
Returns:
|
||||
The steady-state population distribution, splitting the `total_population` vehicles across their ages.
|
||||
"""
|
||||
|
||||
# Initialize mortality
|
||||
if mortality_distribution is None:
|
||||
mortality_distribution = DEFAULT_MORTALITY_DISTRIBUTION
|
||||
|
||||
# Initialize population with uniform distribution
|
||||
population_distribution = np.ones(mortality_distribution.shape) * (total_population / len(mortality_distribution))
|
||||
old_population_distribution = np.zeros(mortality_distribution.shape)
|
||||
i = 0
|
||||
|
||||
while not np.all(np.isclose(population_distribution, old_population_distribution)):
|
||||
old_population = sum(population_distribution)
|
||||
old_population_distribution = population_distribution.copy()
|
||||
i += 1
|
||||
|
||||
# First, simulate the mortality and keep count of the new vehicles we need to produce to replace the existing ones
|
||||
new_cars = population_distribution[-1]
|
||||
for j in reversed(range(len(mortality_distribution) - 1)):
|
||||
population_distribution[j + 1] = (1 - mortality_distribution[j]) * population_distribution[j]
|
||||
new_cars += (mortality_distribution[j] * population_distribution[j])
|
||||
|
||||
# Last, add new vehicles in otder for the sum to make sense
|
||||
population_distribution[0] = 0
|
||||
population_before_creation = sum(population_distribution)
|
||||
new_cars = ((old_population) + 0) - population_before_creation
|
||||
population_distribution[0] = new_cars
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
logger.info(f"Steady-state population achieved after {i} iterations.")
|
||||
|
||||
return population_distribution
|
||||
|
||||
def calculate_subsidy_mortality(steady_state_population: np.array, normal_mortality_distribution: np.array,
|
||||
min_age: int = 10, total_cars: int = 1e6):
|
||||
"""
|
||||
Calculates a mortality distribution for the subsidy year so that `total_cars` with an age greater or equal then
|
||||
`min_age` leave the population.
|
||||
"""
|
||||
|
||||
# Count the total cars older than the minimum age:
|
||||
total_old_cars = np.sum(steady_state_population[min_age:])
|
||||
|
||||
assert total_old_cars > total_cars, "You cannot recycle more vehicles than exist."
|
||||
|
||||
# Calculate which fraction of them needs to be recycles
|
||||
recycling_fraction = total_cars/total_old_cars
|
||||
|
||||
subsidy_mortality_distribution = normal_mortality_distribution.copy()
|
||||
for i in np.arange(min_age, len(steady_state_population)):
|
||||
subsidy_mortality_distribution[i] += recycling_fraction
|
||||
assert subsidy_mortality_distribution[i] <= 1, "You cannot recycle more vehicles than exist."
|
||||
|
||||
return subsidy_mortality_distribution
|
||||
|
||||
def calculate_covid_mortality(steady_state_population: np.array, normal_mortality_distribution: np.array,
|
||||
reduction_in_production: int = 0.23):
|
||||
"""
|
||||
Calculates a mortality distribution for the subsidy year so that `total_cars` with an age greater or equal then
|
||||
`min_age` leave the population.
|
||||
"""
|
||||
return normal_mortality_distribution * (1-reduction_in_production)
|
||||
|
||||
def iterate(mortality_distribution: np.array, population_distribution: np.array, fraction_of_new_cars: float = 1.0):
|
||||
old_population = sum(population_distribution)
|
||||
|
||||
# First, simulate the mortality
|
||||
for j in reversed(range(len(mortality_distribution) - 1)):
|
||||
population_distribution[j + 1] = (1 - mortality_distribution[j]) * population_distribution[j]
|
||||
|
||||
# Last, add new vehicles in order for the sum to make sense
|
||||
population_distribution[0] = 0
|
||||
population_before_creation = sum(population_distribution)
|
||||
new_cars = old_population - population_before_creation
|
||||
population_distribution[0] = new_cars * fraction_of_new_cars
|
||||
|
||||
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()
|
||||
|
||||
subsidy_mortality = calculate_subsidy_mortality(population_distribution, DEFAULT_MORTALITY_DISTRIBUTION)
|
||||
|
||||
new_population_distribution = iterate(DEFAULT_MORTALITY_DISTRIBUTION, population_distribution)
|
||||
normal_year_new_cars = new_population_distribution[0]
|
||||
|
||||
subsidy_population_distribution = iterate(subsidy_mortality, population_distribution)
|
||||
subsidy_new_cars = subsidy_population_distribution[0]
|
||||
|
||||
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')
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
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()
|
||||
Loading…
Reference in New Issue
Block a user