Multi-sample shortest paths#

Here we will see how the multi-sample capabilities of CORNETO allow to exactly solve more complex problems on networks.

Imagine you have a directed graph G with a defined set of vertices and edges. Each edge in this graph has two different costs associated with it, one for each sample (let’s call these Sample 1 and Sample 2). These costs represent two possible scenarios or configurations. For each sample, we want to find the shortest path using the costs defined for that sample. We could run a standard shortest path method per sample. However, given that there can be multiple shortest paths with the same total cost, following that strategy, we cannot identify interesting solutions, for example, if there is a path that is the shortest for both samples. Without a multi-sample method, we need to enumerate all shortest paths per sample, and then compare them to find the answer. This is not only computationally expensive but also not feasible for large networks.

Using CORNETO, it is possible to create a multi shortest path problem to identify a pair of shortest paths—one for each sample—that minimize the number of “non-selected” edges across both samples. This objective can be interpreted in two ways:

  1. If there are multiple shortest paths with the same total cost for a given sample, prefer the one with the fewest edges (i.e., the most direct route).

  2. If the shortest paths differ between the two samples, prioritize the selection of paths that share the greatest number of common edges. This ensures that, despite different costs, the paths are as similar as possible.

We will see how to solve this problem using CORNETO’s multi-sample capabilities.

import corneto as cn

cn.info()
Installed version:v1.0.0.dev3 (latest stable: v1.0.0-alpha)
Available backends:CVXPY v1.6.0, PICOS v2.5.1
Default backend (corneto.opt):CVXPY
Installed solvers:CVXOPT, GLPK, GLPK_MI, HIGHS, SCIP, SCIPY
Graphviz version:v0.20.3
Installed path:/home/runner/work/corneto/corneto/corneto
Repository:https://github.com/saezlab/corneto
G = cn.Graph()
G.add_edge("A", "B")
G.add_edge("A", "C")
G.add_edge("A", "F")
G.add_edge("B", "D")
G.add_edge("B", "E")
G.add_edge("C", "F")
G.add_edge("D", "A")
G.add_edge("D", "C")
G.add_edge("D", "G")
G.add_edge("E", "D")
G.add_edge("E", "H")
G.add_edge("F", "G")
G.add_edge("G", "E")
G.add_edge("G", "H")
G.add_edge("H", "F")
G.plot()
../../_images/64708dcf8e21fbf201e5e2ced254ecf76f437aa9e12028dc0896e8c7f9c0eafd.svg
for i, e in enumerate(G.E):
    print(i, e)
0 (frozenset({'A'}), frozenset({'B'}))
1 (frozenset({'A'}), frozenset({'C'}))
2 (frozenset({'A'}), frozenset({'F'}))
3 (frozenset({'B'}), frozenset({'D'}))
4 (frozenset({'B'}), frozenset({'E'}))
5 (frozenset({'C'}), frozenset({'F'}))
6 (frozenset({'D'}), frozenset({'A'}))
7 (frozenset({'D'}), frozenset({'C'}))
8 (frozenset({'D'}), frozenset({'G'}))
9 (frozenset({'E'}), frozenset({'D'}))
10 (frozenset({'E'}), frozenset({'H'}))
11 (frozenset({'F'}), frozenset({'G'}))
12 (frozenset({'G'}), frozenset({'E'}))
13 (frozenset({'G'}), frozenset({'H'}))
14 (frozenset({'H'}), frozenset({'F'}))
# Weights for the 14 edges for the two samples
c1_weights = [3, 6, 7, 1, 7, 1, 1, 2, 6, 2, 5, 2, 1, 6, 1]
c2_weights = [5, 15, 17, 6, 14, 1, 2, 3, 6, 1, 3, 2, 2, 5, 2]
from corneto.methods import solve_shortest_path

# We solve independently each shortest path problem
edges, P, Gc = solve_shortest_path(G, "A", "H", edge_weights=c1_weights, solver="SCIPY")
Gc.edge_subgraph(edges).plot(orphan_edges=False)
../../_images/845e503ef5eeff51de2e3ac0d2a9ff73f6c61fec08a071d70b0082842d9e5cb7.svg
P.objectives[0].value
np.float64(15.0)
edges, P, Gc = solve_shortest_path(G, "A", "H", edge_weights=c2_weights, solver="SCIPY")
Gc.edge_subgraph(edges).plot(orphan_edges=False)
../../_images/0e35e72763cfaa3e9b6a43d268c4454b08cba574de910d183980a692b5a9b9bf.svg
P.objectives[0].value
np.float64(22.0)
from corneto.methods.shortest_path import create_multisample_shortest_path

samples = [("A", "H"), ("A", "H")]
edge_weights = [c1_weights, c2_weights]
P, Gc = create_multisample_shortest_path(
    G, samples, edge_weights=edge_weights, lam=0.01
)

Gc.plot()
../../_images/ebf2fb31993bc89ed036c88fc10c8af4b950f549b2f8575420b2d57d7e650cd8.svg
P.solve(solver="SCIPY", verbosity=1)
===============================================================================
                                     CVXPY                                     
                                     v1.6.0                                    
===============================================================================
(CVXPY) Feb 08 11:47:11 AM: Your problem has 51 variables, 126 constraints, and 0 parameters.
(CVXPY) Feb 08 11:47:11 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Feb 08 11:47:11 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Feb 08 11:47:11 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Feb 08 11:47:11 AM: Your problem is compiled with the CPP canonicalization backend.
-------------------------------------------------------------------------------
                                  Compilation                                  
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:11 AM: Compiling problem (target solver=SCIPY).
(CVXPY) Feb 08 11:47:11 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCIPY
(CVXPY) Feb 08 11:47:11 AM: Applying reduction Dcp2Cone
(CVXPY) Feb 08 11:47:11 AM: Applying reduction CvxAttr2Constr
(CVXPY) Feb 08 11:47:11 AM: Applying reduction ConeMatrixStuffing
(CVXPY) Feb 08 11:47:11 AM: Applying reduction SCIPY
(CVXPY) Feb 08 11:47:11 AM: Finished problem compilation (took 1.152e-02 seconds).
-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:11 AM: Invoking solver SCIPY  to obtain a solution.
Solver terminated with message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
-------------------------------------------------------------------------------
                                    Summary                                    
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:11 AM: Problem status: optimal
(CVXPY) Feb 08 11:47:11 AM: Optimal value: 3.705e+01
(CVXPY) Feb 08 11:47:11 AM: Compilation took 1.152e-02 seconds
(CVXPY) Feb 08 11:47:11 AM: Solver (including time spent in interface) took 1.936e-03 seconds
Problem(Minimize(Expression(AFFINE, UNKNOWN, ())), [Inequality(Constant(CONSTANT, ZERO, (17, 2))), Inequality(Variable((17, 2), _flow)), Equality(Expression(AFFINE, UNKNOWN, (8, 2)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (2,)), Constant(CONSTANT, NONNEGATIVE, ())), Equality(Expression(AFFINE, UNKNOWN, (2,)), Constant(CONSTANT, NONNEGATIVE, ())), Equality(Expression(AFFINE, UNKNOWN, (2,)), Constant(CONSTANT, NONNEGATIVE, ())), Equality(Expression(AFFINE, UNKNOWN, (2,)), Constant(CONSTANT, NONNEGATIVE, ())), Inequality(Expression(AFFINE, UNKNOWN, (17,))), Inequality(Variable((17,), active_edge, boolean=True))])
# We obtain the same optimal cost for each shortest path, but now we identified a shortest path that is common to both samples
P.objectives[0].value, P.objectives[1].value, P.objectives[2].value
(np.float64(15.0), np.float64(22.0), np.float64(5.0))
# Even though flow is continuous, we have discrete values which indicate the shortest paths.
P.expr.flow.value
array([[ 1.,  1.],
       [-0., -0.],
       [-0., -0.],
       [-0., -0.],
       [ 1.,  1.],
       [ 0.,  0.],
       [-0., -0.],
       [-0., -0.],
       [-0., -0.],
       [-0., -0.],
       [ 1.,  1.],
       [ 0.,  0.],
       [-0., -0.],
       [-0., -0.],
       [-0., -0.],
       [ 1.,  1.],
       [ 1.,  1.]])
import numpy as np

sum(np.sum(P.expr.flow.value, axis=1) > 0)
np.int64(5)
Gc.edge_subgraph(np.flatnonzero(P.expr.flow.value[:, 0] > 0)).plot()
../../_images/d04abbfb2e9ea8b908ab9a6139f9121edde15eeee4600640c6c3095e9a169d2b.svg
Gc.edge_subgraph(np.flatnonzero(P.expr.flow.value[:, 1] > 0)).plot()
../../_images/d04abbfb2e9ea8b908ab9a6139f9121edde15eeee4600640c6c3095e9a169d2b.svg