From 68745efe9d39d4a1eef67e694beeb7bacffaabff Mon Sep 17 00:00:00 2001
From: Arthur Evrard <evrardarthur@orange.fr>
Date: Wed, 20 Oct 2021 15:51:10 +0200
Subject: [PATCH 1/5] 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.
---
 aequilibrae/paths/AoN.pyx               |  3 +-
 aequilibrae/paths/inrets.pyx            | 79 +++++++++++++++++++++++++
 aequilibrae/paths/traffic_assignment.py |  2 +-
 aequilibrae/paths/vdf.py                |  9 ++-
 docs/source/traffic_assignment.rst      | 12 +++-
 tests/aequilibrae/paths/test_vdf.py     |  2 +-
 6 files changed, 99 insertions(+), 8 deletions(-)
 create mode 100644 aequilibrae/paths/inrets.pyx

diff --git a/aequilibrae/paths/AoN.pyx b/aequilibrae/paths/AoN.pyx
index 918fbbf4c..334ae47ae 100644
--- a/aequilibrae/paths/AoN.pyx
+++ b/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
@@ -22,6 +22,7 @@ from libcpp cimport bool
 include 'basic_path_finding.pyx'
 include 'bpr.pyx'
 include 'conical.pyx'
+include 'inrets.pyx'
 include 'parallel_numpy.pyx'
 include 'path_file_saving.pyx'
 
diff --git a/aequilibrae/paths/inrets.pyx b/aequilibrae/paths/inrets.pyx
new file mode 100644
index 000000000..75e492b2d
--- /dev/null
+++ b/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]
diff --git a/aequilibrae/paths/traffic_assignment.py b/aequilibrae/paths/traffic_assignment.py
index 5f85d9466..0cad4bed3 100644
--- a/aequilibrae/paths/traffic_assignment.py
+++ b/aequilibrae/paths/traffic_assignment.py
@@ -257,7 +257,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", "CONICAL", "INRETS"]:
             for p1 in ["alpha", "beta"]:
                 if p1 not in par:
                     raise ValueError(f"{p1} should exist in the set of parameters provided")
diff --git a/aequilibrae/paths/vdf.py b/aequilibrae/paths/vdf.py
index 1cb48882d..3255b1b34 100644
--- a/aequilibrae/paths/vdf.py
+++ b/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, 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", "conical", "inrets"]
 
 
 class VDF:
@@ -17,7 +17,7 @@ class VDF:
 
         vdf = VDF()
         vdf.functions_available()
-      ['bpr', 'conical']
+      ['bpr', 'conical', 'inrets']
 
     """
 
@@ -36,6 +36,9 @@ def __setattr__(self, instance, value) -> None:
             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:
diff --git a/docs/source/traffic_assignment.rst b/docs/source/traffic_assignment.rst
index 642a84e6f..b9c993aa2 100644
--- a/docs/source/traffic_assignment.rst
+++ b/docs/source/traffic_assignment.rst
@@ -73,14 +73,22 @@ 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
+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
diff --git a/tests/aequilibrae/paths/test_vdf.py b/tests/aequilibrae/paths/test_vdf.py
index c1e1085d8..5c6a70f9a 100644
--- a/tests/aequilibrae/paths/test_vdf.py
+++ b/tests/aequilibrae/paths/test_vdf.py
@@ -5,7 +5,7 @@
 class TestVDF(TestCase):
     def test_functions_available(self):
         v = VDF()
-        self.assertEqual(v.functions_available(), ["bpr", "conical"], "VDF class returning wrong availability")
+        self.assertEqual(v.functions_available(), ["bpr", "conical", "inrets"], "VDF class returning wrong availability")
         self.assertEqual(v.apply_vdf, None, "VDF is missing term")
         self.assertEqual(v.apply_derivative, None, "VDF is missing term")
 

From 5010e7d4681eb0d59b4da55faa53e4f5911efa73 Mon Sep 17 00:00:00 2001
From: Arthur Evrard <evrardarthur@orange.fr>
Date: Wed, 20 Oct 2021 16:09:42 +0200
Subject: [PATCH 2/5] 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 :stuck_out_tongue_winking_eye:
---
 aequilibrae/paths/AoN.pyx               |  1 +
 aequilibrae/paths/bpr2.pyx              | 70 +++++++++++++++++++++++++
 aequilibrae/paths/traffic_assignment.py |  2 +-
 aequilibrae/paths/vdf.py                |  9 ++--
 docs/source/traffic_assignment.rst      |  4 +-
 tests/aequilibrae/paths/test_vdf.py     |  2 +-
 6 files changed, 82 insertions(+), 6 deletions(-)
 create mode 100644 aequilibrae/paths/bpr2.pyx

diff --git a/aequilibrae/paths/AoN.pyx b/aequilibrae/paths/AoN.pyx
index 334ae47ae..9ad9906b2 100644
--- a/aequilibrae/paths/AoN.pyx
+++ b/aequilibrae/paths/AoN.pyx
@@ -21,6 +21,7 @@ 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'
diff --git a/aequilibrae/paths/bpr2.pyx b/aequilibrae/paths/bpr2.pyx
new file mode 100644
index 000000000..27a30ad19
--- /dev/null
+++ b/aequilibrae/paths/bpr2.pyx
@@ -0,0 +1,70 @@
+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]
\ No newline at end of file
diff --git a/aequilibrae/paths/traffic_assignment.py b/aequilibrae/paths/traffic_assignment.py
index 0cad4bed3..bb90cd3c4 100644
--- a/aequilibrae/paths/traffic_assignment.py
+++ b/aequilibrae/paths/traffic_assignment.py
@@ -257,7 +257,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", "INRETS"]:
+        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")
diff --git a/aequilibrae/paths/vdf.py b/aequilibrae/paths/vdf.py
index 3255b1b34..bb846a48e 100644
--- a/aequilibrae/paths/vdf.py
+++ b/aequilibrae/paths/vdf.py
@@ -1,11 +1,11 @@
 from aequilibrae import logger
 
 try:
-    from aequilibrae.paths.AoN import bpr, delta_bpr, conical, delta_conical, inrets, delta_inrets
+    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", "inrets"]
+all_vdf_functions = ["bpr", "bpr2", "conical", "inrets"]
 
 
 class VDF:
@@ -17,7 +17,7 @@ class VDF:
 
         vdf = VDF()
         vdf.functions_available()
-      ['bpr', 'conical', 'inrets']
+      ['bpr', 'bpr2', 'conical', 'inrets']
 
     """
 
@@ -33,6 +33,9 @@ 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
diff --git a/docs/source/traffic_assignment.rst b/docs/source/traffic_assignment.rst
index b9c993aa2..0b0e2ea3a 100644
--- a/docs/source/traffic_assignment.rst
+++ b/docs/source/traffic_assignment.rst
@@ -77,6 +77,8 @@ 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)`
 
+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)`
@@ -96,7 +98,7 @@ 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
diff --git a/tests/aequilibrae/paths/test_vdf.py b/tests/aequilibrae/paths/test_vdf.py
index 5c6a70f9a..161a19bc4 100644
--- a/tests/aequilibrae/paths/test_vdf.py
+++ b/tests/aequilibrae/paths/test_vdf.py
@@ -5,7 +5,7 @@
 class TestVDF(TestCase):
     def test_functions_available(self):
         v = VDF()
-        self.assertEqual(v.functions_available(), ["bpr", "conical", "inrets"], "VDF class returning wrong availability")
+        self.assertEqual(v.functions_available(), ["bpr","bpr2" , "conical", "inrets"], "VDF class returning wrong availability")
         self.assertEqual(v.apply_vdf, None, "VDF is missing term")
         self.assertEqual(v.apply_derivative, None, "VDF is missing term")
 

From 7df6a8cf4e408024c3c97a291aba21881438393d Mon Sep 17 00:00:00 2001
From: Arthur Evrard <evrardarthur@orange.fr>
Date: Thu, 21 Oct 2021 13:09:02 +0200
Subject: [PATCH 3/5] Style issues 1

Correction of style issues base on flake8 for .py files (E501, line too long , is ignore).
Do we need to correct all E501 ?
Anything to check for pyx and rst files ?
---
 aequilibrae/paths/traffic_assignment.py | 1 -
 tests/aequilibrae/paths/test_vdf.py     | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/aequilibrae/paths/traffic_assignment.py b/aequilibrae/paths/traffic_assignment.py
index bb90cd3c4..ef6906bf5 100644
--- a/aequilibrae/paths/traffic_assignment.py
+++ b/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
diff --git a/tests/aequilibrae/paths/test_vdf.py b/tests/aequilibrae/paths/test_vdf.py
index 161a19bc4..fea84c4b9 100644
--- a/tests/aequilibrae/paths/test_vdf.py
+++ b/tests/aequilibrae/paths/test_vdf.py
@@ -5,7 +5,7 @@
 class TestVDF(TestCase):
     def test_functions_available(self):
         v = VDF()
-        self.assertEqual(v.functions_available(), ["bpr","bpr2" , "conical", "inrets"], "VDF class returning wrong availability")
+        self.assertEqual(v.functions_available(), ["bpr", "bpr2", "conical", "inrets"], "VDF class returning wrong availability")
         self.assertEqual(v.apply_vdf, None, "VDF is missing term")
         self.assertEqual(v.apply_derivative, None, "VDF is missing term")
 

From cc1e90fdee9b45a84d95d7fdf6d6fa726f929ba8 Mon Sep 17 00:00:00 2001
From: Arthur Evrard <evrardarthur@orange.fr>
Date: Thu, 21 Oct 2021 17:25:27 +0200
Subject: [PATCH 4/5] Indentation correction and inrets test

---
 aequilibrae/paths/bpr.pyx              | 50 ++++++++---------
 aequilibrae/paths/bpr2.pyx             | 67 ++++++++++++----------
 aequilibrae/paths/conical.pyx          | 48 ++++++++--------
 aequilibrae/paths/inrets.pyx           | 78 +++++++++++++-------------
 tests/aequilibrae/paths/test_inrets.py | 60 ++++++++++++++++++++
 5 files changed, 184 insertions(+), 119 deletions(-)
 create mode 100644 tests/aequilibrae/paths/test_inrets.py

diff --git a/aequilibrae/paths/bpr.pyx b/aequilibrae/paths/bpr.pyx
index b8a3d63c9..f52bc0ea0 100644
--- a/aequilibrae/paths/bpr.pyx
+++ b/aequilibrae/paths/bpr.pyx
@@ -29,36 +29,36 @@ 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]
\ No newline at end of file
+        deltaresult[i] = fftime[i]
diff --git a/aequilibrae/paths/bpr2.pyx b/aequilibrae/paths/bpr2.pyx
index 27a30ad19..ec1246d36 100644
--- a/aequilibrae/paths/bpr2.pyx
+++ b/aequilibrae/paths/bpr2.pyx
@@ -29,42 +29,47 @@ def delta_bpr2(dbpr2, link_flows, capacity, fftime, alpha, beta, cores):
 @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]
+                    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])))
+    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] * (1 + alpha[i] * (pow(link_flows[i] / capacity[i], beta[i])))
-      else:
-        congested_time[i] = fftime[i]
+            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]
+                    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]
\ No newline at end of file
+    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]
diff --git a/aequilibrae/paths/conical.pyx b/aequilibrae/paths/conical.pyx
index 1124201b4..fd98a5a38 100644
--- a/aequilibrae/paths/conical.pyx
+++ b/aequilibrae/paths/conical.pyx
@@ -29,23 +29,23 @@ 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 +57,14 @@ 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]
diff --git a/aequilibrae/paths/inrets.pyx b/aequilibrae/paths/inrets.pyx
index 75e492b2d..152abae90 100644
--- a/aequilibrae/paths/inrets.pyx
+++ b/aequilibrae/paths/inrets.pyx
@@ -29,51 +29,51 @@ def delta_inrets(dbpr, link_flows, capacity, fftime, alpha, beta, cores):
 @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]
+                    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) )
+    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] * (
-                  1.1 - (alpha[i]*(link_flows[i] / capacity[i]))) / (
-                  1.1 - (link_flows[i] / capacity[i]) )
-      else:
-        congested_time[i] = fftime[i]
+            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]
+                        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))
+    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]
+        else:
+            deltaresult[i] = fftime[i]
diff --git a/tests/aequilibrae/paths/test_inrets.py b/tests/aequilibrae/paths/test_inrets.py
new file mode 100644
index 000000000..1a65d4612
--- /dev/null
+++ b/tests/aequilibrae/paths/test_inrets.py
@@ -0,0 +1,60 @@
+from unittest import TestCase
+from aequilibrae.paths.AoN import inrets, delta_inrets
+from multiprocessing import cpu_count
+import numpy as np
+
+class TestInrets(TestCase):
+    def test_inrets_funtion(self):
+        cores = cpu_count()
+
+        alpha = np.zeros(11)
+        beta = np.zeros(11)
+        fftime = np.zeros(11)
+        capacity = np.zeros(11)
+        congested_times = np.zeros(11)
+        delta = np.zeros(11)
+
+        alpha.fill(0.95)
+        beta.fill(0)
+        fftime.fill(1)
+        capacity.fill(1)
+        link_flows = np.arange(11).astype(float) * 0.2
+
+        inrets(congested_times, link_flows, capacity, fftime, alpha, beta, cores)
+
+        should_be = np.array(
+            [
+                1,
+                1.011111111,
+                1.028571429,
+                1.06,
+                1.133333333,
+                1.5,
+                2.16,
+                2.94,
+                3.84,
+                4.86,
+                6,
+            ]
+        )
+
+        for i in range(11):
+            self.assertAlmostEqual(should_be[i], congested_times[i], 5, "Inrets is wrong")
+
+        # Let's check the derivative for sections of the curve
+        dx = 0.00000001
+        for i in range(1, 11):
+            link_flows.fill(1 * 0.2 * i)
+            link_flows += np.arange(11) * dx
+
+            inrets(congested_times, link_flows, capacity, fftime, alpha, beta, cores)
+            delta_inrets(delta, link_flows, capacity, fftime, alpha, beta, cores)
+
+            # The derivative needs to be monotonically increasing.
+            self.assertGreater(delta[1], delta[0], "Delta is not increasing as it should")
+
+            # We check if the analytical solution matches the numerical differentiation
+            dydx = (congested_times[1] - congested_times[0]) / dx
+            self.assertAlmostEqual(dydx, delta[1], 6, "Problems with derivative for the inrets vdf")
+            print(dydx, delta[1])
+            
\ No newline at end of file

From 8eb6986e61335f985d497ca9caea3f6bd60530f5 Mon Sep 17 00:00:00 2001
From: Arthur Evrard <evrardarthur@orange.fr>
Date: Mon, 15 Nov 2021 09:49:38 +0100
Subject: [PATCH 5/5] Fix break line error

---
 aequilibrae/paths/bpr.pyx     | 6 ++++--
 aequilibrae/paths/bpr2.pyx    | 6 +++---
 aequilibrae/paths/conical.pyx | 8 +++++---
 3 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/aequilibrae/paths/bpr.pyx b/aequilibrae/paths/bpr.pyx
index f52bc0ea0..c3f26e17e 100644
--- a/aequilibrae/paths/bpr.pyx
+++ b/aequilibrae/paths/bpr.pyx
@@ -40,7 +40,8 @@ cpdef void bpr_cython(double[:] congested_time,
 
     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])))
+            congested_time[i] = fftime[i] * (1 + alpha[i] * (
+            pow(link_flows[i] / capacity[i], beta[i])))
         else:
             congested_time[i] = fftime[i]
 
@@ -59,6 +60,7 @@ cpdef void dbpr_cython(double[:] deltaresult,
 
     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]
+            deltaresult[i] = fftime[i] * (alpha[i] * beta[i] * (
+            pow(link_flows[i] / capacity[i], beta[i]-1))) / capacity[i]
     else:
         deltaresult[i] = fftime[i]
diff --git a/aequilibrae/paths/bpr2.pyx b/aequilibrae/paths/bpr2.pyx
index ec1246d36..160a5edda 100644
--- a/aequilibrae/paths/bpr2.pyx
+++ b/aequilibrae/paths/bpr2.pyx
@@ -66,10 +66,10 @@ cpdef void dbpr2_cython(double[:] deltaresult,
         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]
+                    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]
+                    pow(link_flows[i] / capacity[i], beta[i]-1))) / capacity[i]
         else:
             deltaresult[i] = fftime[i]
diff --git a/aequilibrae/paths/conical.pyx b/aequilibrae/paths/conical.pyx
index fd98a5a38..e6eaca3ca 100644
--- a/aequilibrae/paths/conical.pyx
+++ b/aequilibrae/paths/conical.pyx
@@ -42,7 +42,8 @@ cpdef void conical_cython(double[:] congested_time,
         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] * (
+                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]
@@ -63,8 +64,9 @@ cpdef void dconical_cython(double[:] deltaresult,
     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)))))
+                    (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]