"""Utility functions."""
from collections.abc import Sequence
import numpy as np
from numpy.typing import ArrayLike, NDArray
from scipy.optimize import Bounds, LinearConstraint, milp
from tno.quantum.problems.mot import DetectionSequence
def capacity_constraints(
nd: list[int],
) -> tuple[NDArray[np.int32], NDArray[np.int32], NDArray[np.int32]]:
r"""Creates capacity constraints \sum x_ijt <= 1.
Returns a constraint matrix A, left-hand-side b_l and right-hand-side b_u
such that b_l <= A <= b_u
Args:
nd: the list of number of detections per frame
Returns:
A tuple (A, b_l, b_u)
"""
rows = []
T = len(nd) # noqa: N806
num_vars = np.sum([nd[j] * nd[j + 1] for j in range(T - 1)])
# capacity constraints intermediate nodes
ptr = nd[0] * nd[1]
for t in range(1, T - 1):
for i in range(nd[t]):
row = np.zeros(num_vars, dtype=np.int32)
start = ptr + i * nd[t + 1]
end = start + nd[t + 1]
row[start:end] = 1
rows.append(row)
ptr += nd[t] * nd[t + 1]
a = np.array(rows, dtype=np.int32)
b_u = np.ones(a.shape[0], dtype=np.int32) # \leq 1
b_l = np.full_like(b_u.astype(int), -999999)
return a, b_l, b_u
def flow_constraints(
nd: list[int],
) -> tuple[NDArray[np.int32], NDArray[np.int32], NDArray[np.int32]]:
r"""Creates flow constraints \sum x_ijt = 0.
Returns a constraint matrix A, left-hand-side b_l and right-hand-side b_u
such that b_l <= A <= b_u
Args:
nd: the list of number of detections per frame
Returns:
A tuple (A, b_l, b_u)
"""
rows = []
T = len(nd) # noqa: N806
for t in range(T - 2):
pad_left = np.zeros(
np.sum([nd[j] * nd[j + 1] for j in range(t)], dtype=np.int16)
)
pad_right = np.zeros(
np.sum([nd[j] * nd[j + 1] for j in range(t + 2, T - 1)], dtype=np.int16)
)
for i in range(nd[t + 1]):
mat1 = np.zeros((nd[t], nd[t + 1]))
mat1[:, i] = 1.0
mat2 = np.zeros((nd[t + 1], nd[t + 2]))
mat2[i, :] = -1.0
row = np.concatenate((pad_left, mat1.flatten(), mat2.flatten(), pad_right))
rows.append(row)
a = np.array(rows, dtype=np.int32)
b_u = np.zeros(a.shape[0], dtype=np.int32) # = 0
b_l = np.zeros(a.shape[0], dtype=np.int32)
return a, b_u, b_l
def boundary_constraints(
nd: list[int],
) -> tuple[NDArray[np.int32], NDArray[np.int32], NDArray[np.int32]]:
r"""Creates boundary constraints for t=0 and t=T-1.
Returns a constraint matrix A, left-hand-side b_l and right-hand-side b_u
such that b_l <= A <= b_u
Args:
nd: the list of number of detections per frame
Returns:
A tuple (A, b_l, b_u)
"""
T = len(nd) # noqa: N806
num_vars = np.sum([nd[j] * nd[j + 1] for j in range(T - 1)])
rows = []
for i in range(nd[0]):
row = np.zeros(num_vars)
start = i * nd[1]
end = start + nd[1]
row[start:end] = 1.0
rows.append(row)
a = np.array(rows, dtype=np.int32)
b_u = np.ones(nd[0], dtype=np.int32) # = 1
b_l = np.ones(nd[0], dtype=np.int32)
return a, b_u, b_l
[docs]
def calculate_cost(det_sequence: DetectionSequence) -> NDArray[np.float64]:
"""Calculates the cost array of a detection sequence.
The cost array is computed based on the Intersection over Union (IoU) between
bounding boxes in consecutive frames of a detection sequence.
Args:
det_sequence: A sequence of detections (FrameDets) for each frame,
where each FrameDets contains bounding boxes in the format xyxy.
Returns:
The cost values computed between all pairs of bounding boxes in consecutive
frames, flattened into a 1D array.
"""
cost_array = []
for t in range(len(det_sequence) - 1):
for i in range(det_sequence[t].size):
for j in range(det_sequence[t + 1].size):
box1 = det_sequence[t].xyxy[i]
box2 = det_sequence[t + 1].xyxy[j]
cost_array.append(compute_iou(box1, box2))
return 1 - np.array(cost_array)
[docs]
def compute_iou(box1: Sequence[float], box2: Sequence[float]) -> float:
"""Computes the Intersection over Union (IoU) between two bounding boxes.
The bounding boxes are represented in the format [x, y, x, y], with the first x, y
being the bottom-left corner coordinates and the second pair being the top-right
corner coordinates. The IoU value is a number between 0 and 1, where 0 means no
overlap and 1 means complete overlap.
Args:
box1 (Sequence[float]): The first bounding box in the format [x, y, x, y].
box2 (Sequence[float]): The second bounding box in the format [x, y, x, y].
Returns:
The IoU value between the two bounding boxes.
"""
if box1[1] > box1[3] or box2[1] > box2[3]:
msg = "Wrong box coordinate format. Must be bottom-left/top-right."
raise ValueError(msg)
# Compute intersection coordinates
inter_x1 = max(box1[0], box2[0])
inter_y1 = max(box1[1], box2[1])
inter_x2 = min(box1[2], box2[2])
inter_y2 = min(box1[3], box2[3])
inter_width = max(0, inter_x2 - inter_x1)
inter_height = max(0, inter_y2 - inter_y1)
inter_area = inter_width * inter_height
area1 = (box1[2] - box1[0]) * (box1[3] - box1[1])
area2 = (box2[2] - box2[0]) * (box2[3] - box2[1])
union_area = area1 + area2 - inter_area
return inter_area / union_area if union_area > 0 else 0.0
def solve_binary_program(
A: ArrayLike, # noqa: N803
b_l: ArrayLike,
b_u: ArrayLike,
c: ArrayLike,
) -> NDArray[np.int64]:
r"""Solves a binary linear programming problem using scipy's milp solver.
Solves binary problem of the form
.. math::
&\textrm{minimize } c\cdot x,
&\textrm{subject to } b_l \leq Ax \leq b_u,
&\textrm{with } x \textrm{ binary}.
Args:
A: constraint matrix
b_l: left hand side of constraints
b_u: right hand side of constraints
c: cost function
Returns:
Solution vector of x of the binary linear program.
"""
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)
bounds = Bounds(0, 1)
res = milp(
c=c,
constraints=constraints,
integrality=integrality,
bounds=bounds,
)
if not res.success:
msg = f"MILP not solved to optimality. {res.message}"
raise RuntimeError(msg)
return np.array(res.x, dtype=np.int64)