Getting Started#
The package stochastic_matching allows you to:
Build a stochastic model (from scratch, from one of the provided models, or by some combination).
Attach arrival rates on nodes to it.
Study the theoretical stability of the system.
Explore the feasible matching rates.
Use a simulator to study the behavior of the system under pre-defined or custom-made policies.
Graph building#
Two types of graphs can be used:
Simple graphs, where edges are pairs of distinct nodes.
Hypergraphs, where edges are arbitrary non-empty subsets of nodes, possibly with multiplicity.
Simple graphs can be created from an adjacency matrix.
[1]:
import numpy as np
import stochastic_matching as sm
adja = np.array([[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]])
diamond = sm.Model(adjacency=adja)
The graph can be displayed.
[2]:
diamond.show_graph()
Hypergraphs can be build from their incidence matrix. For example, the following hypergraph is made of a line of three nodes with mono-edges at both ends.
[3]:
incidence = np.array([[1, 1, 0, 0], [0, 1, 1, 0], [0, 0, 1, 1]])
graph = sm.Model(incidence=incidence, names="alpha")
On display, edges are represented by small black nodes, and the edges displayed alway link one node to one edge.
[4]:
graph.show_graph()
The package provides several graph generators to create graphs without manually describing them.
[5]:
paw = sm.Tadpole()
paw.show_graph()
[6]:
candy = sm.HyperPaddle(names="alpha")
candy.show_graph()
See the Graphs notebook for details.
Stability and achievable matching rates#
To study a matching queues system, you need to associate a graph to arrival rates on nodes.
[7]:
diamond = sm.CycleChain(names="alpha", rates=[2, 3, 2, 1])
From then you can see if a stable solution exists.
[8]:
diamond.stabilizable
[8]:
True
You can access the kernel space of solutions.
[9]:
diamond.kernel
[9]:
Kernels of a graph with 4 nodes and 5 edges.
Node dimension is 0.
Edge dimension is 1
Type: Surjective-only
Node kernel:
[]
Edge kernel:
[[ 1 -1 0 -1 1]]
You can have the default solution, which is the product of the arrival rates by the pseudo-inverse of the incidence matrix. Note that the default solution may not be positive even if a positive solution exists.
[10]:
diamond.base_flow
[10]:
array([1.25, 0.75, 1. , 0.75, 0.25])
You can compute the maximin solution, which maximizes the minimum of the rates over all edges.
[11]:
diamond.maximin
[11]:
array([1.5, 0.5, 1. , 0.5, 0.5])
You can find a flow that maximizes a given edge. The following maximizes the last edge.
[12]:
diamond.optimize_edge(4, 1)
[12]:
array([2., 0., 1., 0., 1.])
We can also minimize it.
[13]:
diamond.optimize_edge(4, -1)
[13]:
array([1., 1., 1., 1., 0.])
Another possibility is to ask to maximize some reward based on weights on the edges.
[14]:
diamond.optimize_rates([1, 2, 2, 2, 1])
[14]:
array([1., 1., 1., 1., 0.])
[15]:
diamond.optimize_rates([1, 2, 2, 2, 4])
[15]:
array([2., 0., 1., 0., 1.])
We can also provide the incompressible flow, i.e. for each edge the minimal rate that goes through in any solution.
[16]:
diamond.incompressible_flow()
[16]:
array([1., 0., 1., 0., 0.])
Flows can be displayed using the show_flow method. If no flow is provided, the base flow is displayed.
[17]:
diamond.show_flow()
[18]:
diamond.show_flow(flow=diamond.maximin)
A color code shows the problematic nodes with flow conservation issues and non-positive edges.
[19]:
diamond.show_flow(flow=diamond.optimize_edge(4, 1))
If you want to have a closer look to the polytope of solutions, you can display the kernel:
[20]:
diamond.show_kernel(disp_flow=True, flow=diamond.maximin)
You can access the vertices of the polytope.
[21]:
diamond.vertices
[21]:
[{'kernel_coordinates': array([-0.25]),
'edge_coordinates': array([1., 1., 1., 1., 0.]),
'null_edges': [4],
'bijective': True},
{'kernel_coordinates': array([0.75]),
'edge_coordinates': array([2., 0., 1., 0., 1.]),
'null_edges': [1, 3],
'bijective': False}]
You can display a vertex.
[22]:
diamond.show_vertex(0)
[23]:
diamond.show_vertex(1)
See the Algebra notebook for more details.
Simulations#
To run a simulation, you just need to specify:
A policy (using the name of an implemented policy or providing a custom one)
Number of arrival rates to simulate (default to 1,000,000)
The maximal queue size (default to 1,000)
A random seed (optional)
The simulator returns an estimate of the stability of the policy: False if a queue was maxed before the end, True otherwise.
[24]:
diamond.run("longest", n_steps=10000000)
[24]:
True
The results of simulation are stored in a logs attribute that contains:
The number of matchings for each edge
For each node and state, the number of steps spent by that node in that state
Number of steps performed
A printable version, plogs, is available.
[25]:
diamond.simulator.plogs
Arrivals: [2497426 3748647 2502895 1251032]
Traffic: [1664184 833241 1251543 832920 418111]
Queues: [[7418931 1425915 652538 ... 0 0 0]
[6668530 1336514 800062 ... 0 0 0]
[8886833 742688 246867 ... 0 0 0]
[7672998 1455916 557020 ... 0 0 0]]
Steps done: 10000000
The flow can be displayed.
[26]:
diamond.show_flow()
The CCDFs can be computed and shown.
[27]:
diamond.simulator.ccdf[:, :20]
[27]:
array([[1.000000e+00, 2.581069e-01, 1.155154e-01, 5.026160e-02,
2.141220e-02, 8.995300e-03, 3.730000e-03, 1.524100e-03,
6.317000e-04, 2.628000e-04, 1.046000e-04, 3.570000e-05,
1.180000e-05, 3.900000e-06, 1.100000e-06, 2.000000e-07,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00],
[1.000000e+00, 3.331470e-01, 1.994956e-01, 1.194894e-01,
7.159350e-02, 4.286670e-02, 2.560520e-02, 1.531590e-02,
9.173200e-03, 5.459200e-03, 3.252300e-03, 1.954700e-03,
1.170700e-03, 7.005000e-04, 4.128000e-04, 2.319000e-04,
1.277000e-04, 7.000000e-05, 4.020000e-05, 2.450000e-05],
[1.000000e+00, 1.113167e-01, 3.704790e-02, 1.236120e-02,
4.097700e-03, 1.349100e-03, 4.402000e-04, 1.369000e-04,
4.160000e-05, 1.150000e-05, 2.500000e-06, 6.000000e-07,
1.000000e-07, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00],
[1.000000e+00, 2.327002e-01, 8.710860e-02, 3.140660e-02,
1.115720e-02, 3.904000e-03, 1.312800e-03, 4.458000e-04,
1.518000e-04, 4.920000e-05, 1.160000e-05, 1.500000e-06,
2.000000e-07, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00]])
[28]:
fig = diamond.simulator.show_ccdf(sort=True)
The queues occupancies are given by:
[29]:
diamond.simulator.avg_queues
[29]:
array([0.4605973, 0.8301682, 0.166806 , 0.3682495])
They can be displayed (optionally sorted).
[30]:
fig = diamond.simulator.show_average_queues(sort=True)
It is possible to use the arrival rates to express the results as waiting times.
[31]:
fig = diamond.simulator.show_average_queues(sort=True, as_time=True)
See the Simulation notebook for more details.