Add inexact algorithm: possibly better for some cases.

Author: SebKrantz

src/models/solve_allocation_by_duality_cgc.jl

@@ -129,6 +129,7 @@ function hessian_duality_cgc(
         graph = auxdata.graph
         nodes = graph.nodes
         kappa = auxdata.kappa
+        inexact_algo = param.duality == 2
         beta = param.beta
         m1dbeta = -1 / beta
         sigma = param.sigma
@@ -147,9 +148,12 @@ function hessian_duality_cgc(
         # Constant term: first part of Bprime
         cons = m1dbeta * n1dnum1 * numbetam1
         ## Constant in T'
-        # cons2 = numbetam1 / (numbetam1+nu*beta)
-        # New: constant in T' -> More accurate!
-        cons3 = m1dbeta * n1dnum1 * (numbetam1+nu*beta) / (1+beta)
+        if inexact_algo
+            cons2 = numbetam1 / (numbetam1+nu*beta)
+        else
+            # New: constant in T' -> More accurate!
+            cons3 = m1dbeta * n1dnum1 * (numbetam1+nu*beta) / (1+beta)
+        end
 
         # Precompute elements
         res = recover_allocation_duality_cgc(x, auxdata)
@@ -203,13 +207,19 @@ function hessian_duality_cgc(
                     if jd == j
                         KPprimeAB1 = m1dbeta * Pjn[j, nd]^(-sigma) * PCj[j]^(sigma-1)
                         term += Qjkn[j, k, n] * (KPprimeAB1 - KPAprimeB1 + KPABprime1) # Derivative of Qjkn
-                        # term += T * (((1+beta) * KPprimeAB1 + cons2 * KPABprime1) * G + Gprime) # T'G + TG'
-                        term += T * ((1+beta) * KPprimeAB1 * G + Gprime) # T'G (first part) + TG'
-                        term += cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G
+                        if inexact_algo
+                            term += T * (((1+beta) * KPprimeAB1 + cons2 * KPABprime1) * G + Gprime) # T'G + TG'
+                        else
+                            term += T * ((1+beta) * KPprimeAB1 * G + Gprime) # T'G (first part) + TG'
+                            term += cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G
+                        end
                     elseif jd == k
                         term += Qjkn[j, k, n] * (KPAprimeB1 - KPABprime1) # Derivative of Qjkn
-                        # term -= T * cons2 * KPABprime1 * G # T'G: second part [B'(k) has opposite sign]
-                        term -= cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G 
+                        if inexact_algo
+                            term -= T * cons2 * KPABprime1 * G # T'G: second part [B'(k) has opposite sign]
+                        else
+                            term -= cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G 
+                        end
                     end
                 end
                 if Qjkn[k, j, n] > 0 # Flows in the direction of j