INRETS and bpr2 VD functions (#273)

* French INRETS VD function

First implementation of the French INRETS Volume Delay function.
Beta is not used but kept for better consistancy across VD functions.

* BPR2 implementation

First implementation of BPR2 Volume Delay function.
This doubles beta above capacity so that fewer vehicles are affected when capacity is exceeded.
Double Beta insteed of using a Beta' allow to use only 2 parameters as for other VD functions. Integration into QGIS GUI is also easier

aequilibrae/paths/AoN.pyx

@@ -4,7 +4,7 @@
  Purpose:    Implement shortest path and network loading routines
  Original Author:  Pedro Camargo (c@margo.co)
  Contributors:
- Last edited by: Pedro Camrgo
+ Last edited by: Arthur E.
  Website:    www.AequilibraE.com
  Repository:  https://github.com/AequilibraE/AequilibraE
  Created:    15/09/2013
@@ -21,7 +21,9 @@ from libcpp cimport bool
 # include 'parameters.pxi'
 include 'basic_path_finding.pyx'
 include 'bpr.pyx'
+include 'bpr2.pyx'
 include 'conical.pyx'
+include 'inrets.pyx'
 include 'parallel_numpy.pyx'
 include 'path_file_saving.pyx'
 

aequilibrae/paths/bpr.pyx

@@ -29,36 +29,38 @@ def delta_bpr(dbpr, link_flows, capacity, fftime, alpha, beta, cores):
 @cython.embedsignature(True)
 @cython.boundscheck(False)
 cpdef void bpr_cython(double[:] congested_time,
-                      double[:] link_flows,
-                      double [:] capacity,
-                      double [:] fftime,
-                      double[:] alpha,
-                      double [:] beta,
-                      int cores):
-  cdef long long i
-  cdef long long l = congested_time.shape[0]
+                    double[:] link_flows,
+                    double [:] capacity,
+                    double [:] fftime,
+                    double[:] alpha,
+                    double [:] beta,
+                    int cores):
+    cdef long long i
+    cdef long long l = congested_time.shape[0]
 
-  for i in prange(l, nogil=True, num_threads=cores):
-      if link_flows[i] > 0:
-        congested_time[i] = fftime[i] * (1 + alpha[i] * (pow(link_flows[i] / capacity[i], beta[i])))
-      else:
-        congested_time[i] = fftime[i]
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            congested_time[i] = fftime[i] * (1 + alpha[i] * (
+            pow(link_flows[i] / capacity[i], beta[i])))
+        else:
+            congested_time[i] = fftime[i]
 
 @cython.wraparound(False)
 @cython.embedsignature(True)
 @cython.boundscheck(False)
 cpdef void dbpr_cython(double[:] deltaresult,
-                       double[:] link_flows,
-                       double [:] capacity,
-                       double [:] fftime,
-                       double[:] alpha,
-                       double [:] beta,
-                       int cores):
-  cdef long long i
-  cdef long long l = deltaresult.shape[0]
+                    double[:] link_flows,
+                    double [:] capacity,
+                    double [:] fftime,
+                    double[:] alpha,
+                    double [:] beta,
+                    int cores):
+    cdef long long i
+    cdef long long l = deltaresult.shape[0]
 
-  for i in prange(l, nogil=True, num_threads=cores):
-    if link_flows[i] > 0:
-        deltaresult[i] = fftime[i] * (alpha[i] * beta[i] * (pow(link_flows[i] / capacity[i], beta[i]-1)))/ capacity[i]
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            deltaresult[i] = fftime[i] * (alpha[i] * beta[i] * (
+            pow(link_flows[i] / capacity[i], beta[i]-1))) / capacity[i]
     else:
-        deltaresult[i] = fftime[i]
+        deltaresult[i] = fftime[i]

aequilibrae/paths/bpr2.pyx

@@ -0,0 +1,75 @@
+from libc.math cimport pow
+from cython.parallel import prange
+
+def bpr2(congested_times, link_flows, capacity, fftime, alpha, beta, cores):
+    cdef int c = cores
+
+    cdef double [:] congested_view = congested_times
+    cdef double [:] link_flows_view = link_flows
+    cdef double [:] capacity_view = capacity
+    cdef double [:] fftime_view = fftime
+    cdef double [:] alpha_view = alpha
+    cdef double [:] beta_view = beta
+
+    bpr2_cython(congested_view, link_flows_view, capacity_view, fftime_view, alpha_view, beta_view, c)
+
+def delta_bpr2(dbpr2, link_flows, capacity, fftime, alpha, beta, cores):
+    cdef int c = cores
+
+    cdef double [:] dbpr2_view = dbpr2
+    cdef double [:] link_flows_view = link_flows
+    cdef double [:] capacity_view = capacity
+    cdef double [:] fftime_view = fftime
+    cdef double [:] alpha_view = alpha
+    cdef double [:] beta_view = beta
+
+    dbpr2_cython(dbpr2_view, link_flows_view, capacity_view, fftime_view, alpha_view, beta_view, c)
+
+@cython.wraparound(False)
+@cython.embedsignature(True)
+@cython.boundscheck(False)
+cpdef void bpr2_cython(double[:] congested_time,
+                    double[:] link_flows,
+                    double [:] capacity,
+                    double [:] fftime,
+                    double[:] alpha,
+                    double [:] beta,
+                    int cores):
+    cdef long long i
+    cdef long long l = congested_time.shape[0]
+
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            if link_flows[i] > capacity[i]:
+                congested_time[i] = fftime[i] * (1 + alpha[i] * (
+                    pow(link_flows[i] / capacity[i], 2*beta[i])))
+            else:
+                congested_time[i] = fftime[i] * (1 + alpha[i] * (
+                    pow(link_flows[i] / capacity[i], beta[i])))
+        else:
+            congested_time[i] = fftime[i]
+
+@cython.wraparound(False)
+@cython.embedsignature(True)
+@cython.boundscheck(False)
+cpdef void dbpr2_cython(double[:] deltaresult,
+                    double[:] link_flows,
+                    double [:] capacity,
+                    double [:] fftime,
+                    double[:] alpha,
+                    double [:] beta,
+                    int cores):
+    cdef long long i
+    cdef long long l = deltaresult.shape[0]
+
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            if link_flows[i] > capacity[i]:
+                deltaresult[i] = fftime[i] * (alpha[i] * 2 * beta[i] * (
+                    pow(link_flows[i] / capacity[i], (2*beta[i])-1))) / ( 
+                    capacity[i])
+            else:
+                deltaresult[i] = fftime[i] * (alpha[i] * beta[i] * (
+                    pow(link_flows[i] / capacity[i], beta[i]-1))) / capacity[i]
+        else:
+            deltaresult[i] = fftime[i]

aequilibrae/paths/conical.pyx

@@ -29,23 +29,24 @@ def delta_conical(dbpr, link_flows, capacity, fftime, alpha, beta, cores):
 @cython.embedsignature(True)
 @cython.boundscheck(False)
 cpdef void conical_cython(double[:] congested_time,
-                      double[:] link_flows,
-                      double [:] capacity,
-                      double [:] fftime,
-                      double[:] alpha,
-                      double [:] beta,
-                      int cores):
-  cdef long long i
-  cdef long long l = congested_time.shape[0]
+                        double[:] link_flows,
+                        double [:] capacity,
+                        double [:] fftime,
+                        double[:] alpha,
+                        double [:] beta,
+                        int cores):
+    cdef long long i
+    cdef long long l = congested_time.shape[0]
 
-  for i in prange(l, nogil=True, num_threads=cores):
-      if link_flows[i] > 0:
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
 
-          congested_time[i] = fftime[i] * (
-                  sqrt(pow(alpha[i], 2) * pow(1 - link_flows[i] / capacity[i], 2) + pow(beta[i], 2)) - alpha[i] * (
-                  1 - link_flows[i] / capacity[i]) - beta[i] + 2)
-      else:
-        congested_time[i] = fftime[i]
+            congested_time[i] = fftime[i] * (
+                sqrt(pow(alpha[i], 2) * pow(1 - link_flows[i] / capacity[i], 2)\
+                + pow(beta[i], 2)) - alpha[i] * (
+                1 - link_flows[i] / capacity[i]) - beta[i] + 2)
+        else:
+            congested_time[i] = fftime[i]
 
 @cython.wraparound(False)
 @cython.embedsignature(True)
@@ -57,14 +58,15 @@ cpdef void dconical_cython(double[:] deltaresult,
                        double[:] alpha,
                        double [:] beta,
                        int cores):
-  cdef long long i
-  cdef long long l = deltaresult.shape[0]
+    cdef long long i
+    cdef long long l = deltaresult.shape[0]
 
-  for i in prange(l, nogil=True, num_threads=cores):
-    if link_flows[i] > 0:
-        deltaresult[i] = fftime[i] * ((alpha[i] / capacity[i]) - (
-                (pow(alpha[i], 2) * (1 - link_flows[i] / capacity[i])) / (capacity[i] * sqrt(
-            pow(alpha[i], 2) * pow(1 - link_flows[i] / capacity[i], 2) + pow(beta[i], 2)))))
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            deltaresult[i] = fftime[i] * ((alpha[i] / capacity[i]) - (
+                    (pow(alpha[i], 2) * (1 - link_flows[i] / capacity[i])) / (
+                    capacity[i] * sqrt(pow(alpha[i], 2) * pow(
+                    1 - link_flows[i] / capacity[i], 2) + pow(beta[i], 2)))))
 
-    else:
-        deltaresult[i] = fftime[i]
+        else:
+            deltaresult[i] = fftime[i]

aequilibrae/paths/inrets.pyx

@@ -0,0 +1,79 @@
+from libc.math cimport pow
+from cython.parallel import prange
+
+def inrets(congested_times, link_flows, capacity, fftime, alpha, beta, cores):
+    cdef int c = cores
+
+    cdef double [:] congested_view = congested_times
+    cdef double [:] link_flows_view = link_flows
+    cdef double [:] capacity_view = capacity
+    cdef double [:] fftime_view = fftime
+    cdef double [:] alpha_view = alpha
+    cdef double [:] beta_view = beta
+
+    inrets_cython(congested_view, link_flows_view, capacity_view, fftime_view, alpha_view, beta_view, c)
+
+def delta_inrets(dbpr, link_flows, capacity, fftime, alpha, beta, cores):
+    cdef int c = cores
+
+    cdef double [:] dbpr_view = dbpr
+    cdef double [:] link_flows_view = link_flows
+    cdef double [:] capacity_view = capacity
+    cdef double [:] fftime_view = fftime
+    cdef double [:] alpha_view = alpha
+    cdef double [:] beta_view = beta
+
+    dinrets_cython(dbpr_view, link_flows_view, capacity_view, fftime_view, alpha_view, beta_view, c)
+
+@cython.wraparound(False)
+@cython.embedsignature(True)
+@cython.boundscheck(False)
+cpdef void inrets_cython(double[:] congested_time,
+                    double[:] link_flows,
+                    double [:] capacity,
+                    double [:] fftime,
+                    double[:] alpha,
+                    double [:] beta,
+                    int cores):
+    cdef long long i
+    cdef long long l = congested_time.shape[0]
+
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            if link_flows[i] > capacity[i]:
+                congested_time[i] = fftime[i] * (
+                    (1.1 - alpha[i])/0.1) * (
+                    pow(link_flows[i] / capacity[i],2) )
+            else:
+                congested_time[i] = fftime[i] * (
+                    1.1 - (alpha[i]*(link_flows[i] / capacity[i]))) / (
+                    1.1 - (link_flows[i] / capacity[i]) )
+        else:
+            congested_time[i] = fftime[i]
+
+@cython.wraparound(False)
+@cython.embedsignature(True)
+@cython.boundscheck(False)
+cpdef void dinrets_cython(double[:] deltaresult,
+                        double[:] link_flows,
+                        double [:] capacity,
+                        double [:] fftime,
+                        double[:] alpha,
+                        double [:] beta,
+                        int cores):
+    cdef long long i
+    cdef long long l = deltaresult.shape[0]
+
+    for i in prange(l, nogil=True, num_threads=cores):
+        if link_flows[i] > 0:
+            if link_flows[i] > capacity[i]:
+                deltaresult[i] = fftime[i] * (
+                    (-20)*(alpha[i]-1.1)*link_flows[i]) / (
+                    pow(capacity[i],2))
+            else:
+                deltaresult[i] = fftime[i] * (
+                    (-110)*(alpha[i]-1)*capacity[i]) / (
+                    pow((11*capacity[i])-(10*link_flows[i]),2))
+
+        else:
+            deltaresult[i] = fftime[i]

aequilibrae/paths/traffic_assignment.py

@@ -8,7 +8,6 @@
 import numpy as np
 import pandas as pd
 from aequilibrae.project.database_connection import ENVIRON_VAR
-from aequilibrae.paths.all_or_nothing import allOrNothing
 from aequilibrae.paths.linear_approximation import LinearApproximation
 from aequilibrae.paths.vdf import VDF, all_vdf_functions
 from aequilibrae.paths.traffic_class import TrafficClass
@@ -257,7 +256,7 @@ def set_vdf_parameters(self, par: dict) -> None:
             raise Exception("Before setting vdf parameters, you need to set traffic classes and choose a VDF function")
         self.__dict__["vdf_parameters"] = par
         pars = []
-        if self.vdf.function in ["BPR", "CONICAL"]:
+        if self.vdf.function in ["BPR", "BPR2", "CONICAL", "INRETS"]:
             for p1 in ["alpha", "beta"]:
                 if p1 not in par:
                     raise ValueError(f"{p1} should exist in the set of parameters provided")

aequilibrae/paths/vdf.py

@@ -1,11 +1,11 @@
 from aequilibrae import logger
 
 try:
-    from aequilibrae.paths.AoN import bpr, delta_bpr, conical, delta_conical
+    from aequilibrae.paths.AoN import bpr, delta_bpr, bpr2, delta_bpr2, conical, delta_conical, inrets, delta_inrets
 except ImportError as ie:
     logger.warning(f"Could not import procedures from the binary. {ie.args}")
 
-all_vdf_functions = ["bpr", "conical"]
+all_vdf_functions = ["bpr", "bpr2", "conical", "inrets"]
 
 
 class VDF:
@@ -17,7 +17,7 @@ class VDF:
 
         vdf = VDF()
         vdf.functions_available()
-      ['bpr', 'conical']
+      ['bpr', 'bpr2', 'conical', 'inrets']
 
     """
 
@@ -33,9 +33,15 @@ def __setattr__(self, instance, value) -> None:
             if value == "BPR":
                 self.__dict__["apply_vdf"] = bpr
                 self.__dict__["apply_derivative"] = delta_bpr
+            elif value == "BPR2":
+                self.__dict__["apply_vdf"] = bpr2
+                self.__dict__["apply_derivative"] = delta_bpr2
             elif value == "CONICAL":
                 self.__dict__["apply_vdf"] = conical
                 self.__dict__["apply_derivative"] = delta_conical
+            elif value == "INRETS":
+                self.__dict__["apply_vdf"] = inrets
+                self.__dict__["apply_derivative"] = delta_inrets
             else:
                 raise ValueError("VDF function not available")
         else:

docs/source/traffic_assignment.rst

@@ -73,22 +73,32 @@ To begin building the assignment it is easy:
 Volume Delay Function
 +++++++++++++++++++++
 
-For now, the only VDF functions available in AequilibraE are the BPR
+For now, the only VDF functions available in AequilibraE are the BPR,
 
 :math:`CongestedTime_{i} = FreeFlowTime_{i} * (1 + \alpha * (\frac{Volume_{i}}{Capacity_{i}})^\beta)`
 
-and Spiess' conical
+BPR2 which double beta when traffic flow is over the link capacity,
+
+Spiess' conical,
 
 :math:`CongestedTime_{i} = FreeFlowTime_{i} * (2 + \sqrt[2][\alpha^2*(1- \frac{Volume_{i}}{Capacity_{i}})^2 + \beta^2] - \alpha *(1-\frac{Volume_{i}}{Capacity_{i}})-\beta)`
 
+and French INRETS (alpha < 1)
+
+Before capacity
+:math:`CongestedTime_{i} = FreeFlowTime_{i} * \frac{1.1- (\alpha *\frac{Volume_{i}}{Capacity_{i}})}{1.1-\frac{Volume_{i}}{Capacity_{i}}}`
+
+and after capacity
+:math:`CongestedTime_{i} = FreeFlowTime_{i} * \frac{1.1- \alpha}{0.1} * (\frac{Volume_{i}}{Capacity_{i}})^2`
+
 More functions will be added as needed/requested/possible.
 
 Setting the volume delay function is one of the first things you should do after
 instantiating an assignment problem in AequilibraE, and it is as simple as:
 
 ::
 
-    assig.set_vdf('CONICAL')
+    assig.set_vdf('BPR')
 
 The implementation of the VDF functions in AequilibraE is written in Cython and
 fully multi-threaded, and therefore descent methods that may evaluate such