Skip to content
GitLab
Commit 4138f774 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

add simulation specific scenarios to zhang-scenarios.

parent 8b8e7875
Pipeline #87554 failed with stages
in 100 minutes and 17 seconds
Expand all Hide whitespace changes
Inline Side-by-side
Tools/VadereAnalysisTools/Plots/fundamentalDiagrams/OSM_calibration.ipynb
View file @ 4138f774
%% Cell type:markdown id: tags:
 
# Calibration of the Optimal Steps Model
 
This script is an attempt to recompute the results in silver-2016b page 51. The scenario [scenario](./../../../../VadereModelTests/TestOSM_calibration/rimea_04_calibration_osm.scenario) is based on the RiMEA-Test 4. We use the `Teleporter` to model a circular scenario and the parameter `useFreeSpaceOnly = false` to generate high densities. The following code plots all the necessary diagrams.
 
%% Cell type:code id: tags:
 
``` python
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from vadere_analysis_tool import ScenarioOutput, VadereProject
from scipy.optimize import curve_fit
 
sns.set(style="whitegrid")
sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.0})
 
def plotEvolution(vproject, ending, yval, ylabel):
plots = []
i = 1
plt.figure(1, figsize=(30, 10))
for outStr in project.output_dirs :
out = project.output_dirs[outStr]
dataFrames = []
for fileStr in out.files :
if fileStr.endswith(ending) :
ndf = pd.DataFrame();
df = out.files[fileStr]()
#df = df[df.velocity > 0]
ndf['density'] = df.density.astype(float)
ndf['velocity'] = df.velocity.astype(float)
ndf['timeStep'] = df.timeStep.astype(int)
ndf['scenario'] = [out.scenario['name']] * len(df.density)
dataFrames.append(ndf)
#concatFrames = pd.concat(dataFrames)
plt.subplot(2, 1, i)
plt.title(out.scenario['name'])
plt.xlabel('timeStep')
plt.ylabel(ylabel)
sns.lineplot(x='timeStep', y=yval, data=pd.concat(dataFrames))
i = i + 1
 
def plotFundamentalDiagram(vproject, ending):
dataFrames = []
for outStr in project.output_dirs :
out = project.output_dirs[outStr]
for fileStr in out.files :
if fileStr.endswith(ending) :
ndf = pd.DataFrame();
df = out.files[fileStr]()
ndf['density'] = df.density.astype(float)
ndf['velocity'] = df.velocity.astype(float)
ndf['scenario'] = [out.scenario['name']] * len(df.velocity)
#ndf = ndf[ndf.density < 7]
dataFrames.append(ndf)
concatFrames = pd.concat(dataFrames)
g = sns.relplot(x="density", y="velocity", hue="scenario", data=concatFrames,
height=10, aspect=2)
def plotFundamentalDiagramScatter(vproject, ending, sep=False, width = 10, height = 5):
dataFrames = []
index = 0
cols = 3
rows = len(vproject.output_dirs) / cols + 1
if not sep :
fig, axs = plt.subplots(int(rows), int(cols), figsize=(height*rows, width*cols), sharex=False, sharey=True)
for outStr in vproject.output_dirs :
out = vproject.output_dirs[outStr]
for fileStr in out.files :
if fileStr.endswith(ending) :
if not sep :
axes = axs[int(index / cols), int(index % cols)]
else :
fig = plt.figure(1, figsize=(width, height))
axes = plt.axes()
fig.add_axes(axes)
ndf = pd.DataFrame();
df = out.files[fileStr]()
ndf['density'] = df.density.astype(float)
ndf['velocity'] = df.velocity.astype(float)
ndf['scenario'] = [out.scenario['name']] * len(df.velocity)
#ndf = ndf[ndf.density < 7]
#plt.scatter()
#axes.set_title(out.scenario['name'])
axes.set_xlabel('density')
axes.set_ylabel('velocity')
axes.set_xticks([0,1,2,3,4,5,6])
axes.set_yticks([0,0.5,1,1.5,2,2.5])
axes.set_xlim(0,6)
axes.set_ylim(0,2.5)
axes.scatter(ndf['density'], ndf['velocity'], s=0.7, marker='*', color='#555555')
wm = plotWeidmann(axes)
popt, pcov = curve_fit(kladek, ndf['density'], ndf['velocity'], p0=(1.34, 1.913, 5.4))
print(str(popt[0]) + "," + str(popt[1]) + "," + str(popt[2]))
xx = np.linspace(0.1, 6, 1000)
yy = kladek(xx, *popt)
axes.plot(xx, yy, '--', c=sns.color_palette().as_hex()[1])
axes.legend(['Weidmann', 'regression', 'Simulated data'])
index = index + 1;
if sep :
fig.savefig("./"+out.scenario['name']+"_fundamental_diagram"+".png", bbox_inches='tight')
plt.show()
if not sep :
fig.savefig("./"+vproject.project_name+"_fundamental_diagrams"+".png", bbox_inches='tight')
 
```
 
%% Cell type:code id: tags:
 
``` python
def plotWeidmann(axes):
wmaxDensity = 5.4
wmeanVelocity = 1.34
wgamma = 1.913
wx = np.linspace(0.1, wmaxDensity, 100)
return plotKladek(wx, wmeanVelocity, wgamma, wmaxDensity, axes)
 
def plotKladek(x, v, gamma, pmax, axes):
result, = axes.plot(x, kladek(x, v, gamma, pmax), c=sns.color_palette().as_hex()[0])
return result
 
def kladek(x, v, gamma, pmax):
return v * (1 - np.exp(-gamma * (1/x - 1/pmax)))
```
 
%% Cell type:markdown id: tags:
 
First we load the project
 
%% Cell type:code id: tags:
 
``` python
projectFolder = "./../../../../VadereModelCalibration/TestOSM_calibration/"
project = VadereProject(projectFolder)
#out = project.named_output.C_050_180_180_2018_11_26_16_30_29_355()'
```
 
%% Output
 
loaded 5 out of 5 output directories.
 
%% Cell type:markdown id: tags:
 
# Meassurement methods
All methods are described in zhang-2011.
 
## Method A Plots
The computation of the velocity is slightly different i.e. we use the velocity computed by the current and last foot step.
 
%% Cell type:code id: tags:
 
``` python
# transform data frame
plotFundamentalDiagramScatter(project, "aTimeStep.fundamentalDiagram", True, 6, 6)
```
 
%% Cell type:markdown id: tags:
 
## Method B Plots
This method does not work for this scenario since agents run multiple times through the same measurement area.
 
%% Cell type:markdown id: tags:
 
## Method C Plots
The computation of the velocity is slightly different i.e. we use the velocity computed by the current and last foot step.
 
%% Cell type:code id: tags:
 
``` python
plotFundamentalDiagramScatter(project, "cTimeStep.fundamentalDiagram", True, 6, 6)
```
 
%% Output
 
1.09598977853,2.60307653064,5.7172600023
 
 
1.04365172849,2.77601030962,6.17771344067
 
 
1.00950690582,4.0231057037,5.42205556099
 
 
1.05243584192,2.9673646152,6.25808694348
 
 
1.06095434794,4.32504288205,5.47219183859
 
 
%% Cell type:markdown id: tags:
 
## Method D Plots
The computation of the velocity is slightly different i.e. we use the velocity computed by the current and last foot step.
 
%% Cell type:code id: tags:
 
``` python
myPlots = plotFundamentalDiagramScatter(project, "dTimeStep.fundamentalDiagram", True, 6, 6)
```
 
%% Output
 
1.31153677879,1.63591017016,6.31329791983
 
 
1.27929507782,1.61554754755,7.41274319187
 
 
1.23002056156,2.22943029934,6.31005306265
 
 
1.29467559573,1.6733716044,7.86471179297
 
 
1.23583074427,2.72867755264,5.83620246417
 
 
%% Cell type:markdown id: tags:
 
## Method E Plots
This method is similar to method D but the density is defined by $$\langle \rho \rangle = \frac{1}{N} \sum\limits_{i=1}^{N} A_i ,$$
and the velocity is defined by
$$\langle v \rangle = \sum\limits_{i=1}^{N} \frac{1}{A_i} \left( \sum\limits_{i=1}^{N} A_i v_i(t) \right) ,$$
 
where $N$ is the number of pedestrians inside the measurement area and $A_i$ is the area of the voronoi cell of agent $i$.
 
%% Cell type:code id: tags:
 
``` python
plotFundamentalDiagramScatter(project, "eTimeStep.fundamentalDiagram", False, 10, 10)
```
 
%% Output
 
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-827-2128e4d6f9ae> in <module>()
----> 1 plotFundamentalDiagramScatter(project, "eTimeStep.fundamentalDiagram", False, 10, 10)
 
<ipython-input-822-60976e09b448> in plotFundamentalDiagramScatter(vproject, ending, sep, width, height)
83 axes.scatter(ndf['density'], ndf['velocity'], s=0.7, marker='*', color='#555555')
84 wm = plotWeidmann(axes)
---> 85 popt, pcov = curve_fit(kladek, ndf['density'], ndf['velocity'], p0=(1.34, 1.913, 5.4))
86 print(str(popt[0]) + "," + str(popt[1]) + "," + str(popt[2]))
87 xx = np.linspace(0.1, 6, 1000)
~/anaconda3/lib/python3.6/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
699 # NaNs can not be handled
700 if check_finite:
--> 701 ydata = np.asarray_chkfinite(ydata)
702 else:
703 ydata = np.asarray(ydata)
~/anaconda3/lib/python3.6/site-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
1213 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
1214 raise ValueError(
-> 1215 "array must not contain infs or NaNs")
1216 return a
1217
ValueError: array must not contain infs or NaNs
 
 
%% Cell type:markdown id: tags:
 
## Density and velocity evolution
 
%% Cell type:code id: tags:
 
``` python
plotEvolution(project, "dTimeStep.fundamentalDiagram", "density", "density")
#plotEvolution(project, "aTimeStep.fundamentalDiagram", "density", "density")
```
 
%% Output