Trip Distribution#
In this example, we calibrate a Synthetic Gravity Model that same model plus IPF (Fratar/Furness).
[1]:
# Imports
from uuid import uuid4
from tempfile import gettempdir
from os.path import join
from aequilibrae.utils.create_example import create_example
import pandas as pd
import numpy as np
We create the example project inside our temp folder
[2]:
fldr = join(gettempdir(), uuid4().hex)
project = create_example(fldr)
[3]:
# We get the demand matrix directly from the project record
# so let's inspect what we have in the project
proj_matrices = project.matrices
print(proj_matrices.list())
name file_name cores procedure procedure_id \
0 demand_omx demand.omx 1 None None
1 demand_mc demand_mc.omx 3 None None
2 skims skims.omx 2 None None
3 demand_aem demand.aem 1 None None
timestamp description status
0 2020-11-24 08:47:18 Original data imported to OMX format
1 2021-02-24 00:51:35 None
2 None Example skim
3 2020-11-24 08:46:42 Original data imported to AEM format
[4]:
# We get the demand matrix
demand = proj_matrices.get_matrix("demand_omx")
demand.computational_view(["matrix"])
# And the impedance
impedance = proj_matrices.get_matrix("skims")
impedance.computational_view(["time_final"])
Let’s have a function to plot the Trip Length Frequency Distribution
[5]:
from math import log10, floor
import matplotlib.pyplot as plt
def plot_tlfd(demand, skim, name):
plt.clf()
b = floor(log10(skim.shape[0]) * 10)
n, bins, patches = plt.hist(
np.nan_to_num(skim.flatten(), 0),
bins=b,
weights=np.nan_to_num(demand.flatten()),
density=False,
facecolor="g",
alpha=0.75,
)
plt.xlabel("Trip length")
plt.ylabel("Probability")
plt.title("Trip-length frequency distribution")
plt.savefig(name, format="png")
return plt
[6]:
from aequilibrae.distribution import GravityCalibration
[7]:
for function in ["power", "expo"]:
gc = GravityCalibration(matrix=demand, impedance=impedance, function=function, nan_as_zero=True)
gc.calibrate()
model = gc.model
# We save the model
model.save(join(fldr, f"{function}_model.mod"))
# We can save an image for the resulting model
_ = plot_tlfd(gc.result_matrix.matrix_view, impedance.matrix_view, join(fldr, f"{function}_tfld.png"))
# We can save the result of applying the model as well
# We can also save the calibration report
with open(join(fldr, f"{function}_convergence.log"), "w") as otp:
for r in gc.report:
otp.write(r + "\n")
/opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/aequilibrae/distribution/gravity_application.py:322: RuntimeWarning: divide by zero encountered in power
self.output.matrix_view[i, :] = (np.power(self.impedance.matrix_view[i, :], -self.model.alpha) * p * a)[
/opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/aequilibrae/distribution/gravity_application.py:336: RuntimeWarning: invalid value encountered in multiply
self.output.matrix_view[:, :] = self.output.matrix_view[:, :] * non_inf
We save a trip length frequency distribution for the demand itself
[8]:
plt = plot_tlfd(demand.matrix_view, impedance.matrix_view, join(fldr, "demand_tfld.png"))
plt.show()
Forecast#
We create a set of ‘future’ vectors by applying some models and apply the model for both deterrence functions
[9]:
from aequilibrae.distribution import Ipf, GravityApplication, SyntheticGravityModel
from aequilibrae.matrix import AequilibraeData
import numpy as np
[10]:
zonal_data = pd.read_sql("Select zone_id, population, employment from zones order by zone_id", project.conn)
# We compute the vectors from our matrix
args = {
"file_path": join(fldr, "synthetic_future_vector.aed"),
"entries": demand.zones,
"field_names": ["origins", "destinations"],
"data_types": [np.float64, np.float64],
"memory_mode": True,
}
vectors = AequilibraeData()
vectors.create_empty(**args)
vectors.index[:] = zonal_data.zone_id[:]
# We apply a trivial regression-based model and balance the vectors
vectors.origins[:] = zonal_data.population[:] * 2.32
vectors.destinations[:] = zonal_data.employment[:] * 1.87
vectors.destinations *= vectors.origins.sum() / vectors.destinations.sum()
[11]:
# We simply apply the models to the same impedance matrix now
for function in ["power", "expo"]:
model = SyntheticGravityModel()
model.load(join(fldr, f"{function}_model.mod"))
outmatrix = join(proj_matrices.fldr, f"demand_{function}_model.aem")
args = {
"impedance": impedance,
"rows": vectors,
"row_field": "origins",
"model": model,
"columns": vectors,
"column_field": "destinations",
"nan_as_zero": True,
}
gravity = GravityApplication(**args)
gravity.apply()
# We get the output matrix and save it to OMX too,
gravity.save_to_project(name=f"demand_{function}_model_omx", file_name=f"demand_{function}_model.omx")
[12]:
# We update the matrices table/records and verify that the new matrices are indeed there
proj_matrices.update_database()
print(proj_matrices.list())
name file_name cores \
0 demand_omx demand.omx 1
1 demand_mc demand_mc.omx 3
2 skims skims.omx 2
3 demand_aem demand.aem 1
4 demand_power_model_omx demand_power_model.omx 1
5 demand_expo_model_omx demand_expo_model.omx 1
procedure procedure_id \
0 None None
1 None None
2 None None
3 None None
4 Synthetic gravity trip distribution 7d3b4204f7644e8f92c8176093613a2d
5 Synthetic gravity trip distribution 25631e86854749aea393b0ebf4a3239c
timestamp description \
0 2020-11-24 08:47:18 Original data imported to OMX format
1 2021-02-24 00:51:35 None
2 None Example skim
3 2020-11-24 08:46:42 Original data imported to AEM format
4 2023-10-09 07:32:49.833228 Synthetic gravity trip distribution. POWER
5 2023-10-09 07:32:50.024875 Synthetic gravity trip distribution. EXPO
status
0
1
2
3
4
5
We now run IPF for the future vectors
[13]:
args = {
"matrix": demand,
"rows": vectors,
"columns": vectors,
"column_field": "destinations",
"row_field": "origins",
"nan_as_zero": True,
}
ipf = Ipf(**args)
ipf.fit()
ipf.save_to_project(name="demand_ipf", file_name="demand_ipf.aem")
ipf.save_to_project(name="demand_ipf_omx", file_name="demand_ipf.omx")
[13]:
<aequilibrae.project.data.matrix_record.MatrixRecord at 0x7f55bc560cd0>
[14]:
print(proj_matrices.list())
name file_name cores \
0 demand_omx demand.omx 1
1 demand_mc demand_mc.omx 3
2 skims skims.omx 2
3 demand_aem demand.aem 1
4 demand_power_model_omx demand_power_model.omx 1
5 demand_expo_model_omx demand_expo_model.omx 1
6 demand_ipf demand_ipf.aem 1
7 demand_ipf_omx demand_ipf.omx 1
procedure procedure_id \
0 None None
1 None None
2 None None
3 None None
4 Synthetic gravity trip distribution 7d3b4204f7644e8f92c8176093613a2d
5 Synthetic gravity trip distribution 25631e86854749aea393b0ebf4a3239c
6 Iterative Proportional fitting ae41e127e96d45e088f6dd4e8a51a1e6
7 Iterative Proportional fitting ae41e127e96d45e088f6dd4e8a51a1e6
timestamp description \
0 2020-11-24 08:47:18 Original data imported to OMX format
1 2021-02-24 00:51:35 None
2 None Example skim
3 2020-11-24 08:46:42 Original data imported to AEM format
4 2023-10-09 07:32:49.833228 Synthetic gravity trip distribution. POWER
5 2023-10-09 07:32:50.024875 Synthetic gravity trip distribution. EXPO
6 2023-10-09 07:32:50.236655 None
7 2023-10-09 07:32:50.236655 None
status
0
1
2
3
4
5
6
7
[15]:
project.close()