# SPDX-License-Identifier: MIT
# https://gitlab.windenergy.dtu.dk/TOPFARM/OptiWindNet/
import abc
import logging
import math
import os
import sys
from collections.abc import Iterator, Mapping
from dataclasses import asdict, dataclass, field
from enum import StrEnum, auto
from inspect import cleandoc
from itertools import chain
from pathlib import Path
from textwrap import indent
from typing import TYPE_CHECKING, Any
import networkx as nx
from makefun import with_signature
from ..interarraylib import (
G_from_S,
_ring_split_position,
bfs_subtree_loads,
directed_links,
)
from ..pathfinding import PathFinder
from ..types import Topology
_lggr = logging.getLogger(__name__)
error, info, warn = _lggr.error, _lggr.info, _lggr.warning
def physical_core_count() -> int:
"""Count physical cores available to this process (cross-platform).
On Linux, reads sysfs topology to count physical cores within the
process affinity set. On other platforms, falls back to psutil.
"""
if sys.platform == 'linux':
try:
affinity = os.sched_getaffinity(0)
except OSError:
affinity = None
if affinity is not None:
physical_cores = set()
for cpu_id in affinity:
topo = Path(f'/sys/devices/system/cpu/cpu{cpu_id}/topology')
try:
pkg = (topo / 'physical_package_id').read_text().strip()
core = (topo / 'core_id').read_text().strip()
physical_cores.add((pkg, core))
except OSError:
# sysfs unavailable (e.g. container), fall through
physical_cores = None
break
if physical_cores is not None:
return len(physical_cores)
# Windows, macOS, or Linux without sysfs
import psutil
return psutil.cpu_count(logical=False) or os.cpu_count() or 1
def _identifier_from_class_name(c: type) -> str:
"Convert a camel-case class name to a snake-case identifier"
s = c.__name__
return s[0].lower() + ''.join('_' + c.lower() if c.isupper() else c for c in s[1:])
[docs]
class OWNWarmupFailed(Exception):
pass
[docs]
class OWNSolutionNotFound(Exception):
pass
[docs]
class FeederRoute(StrEnum):
"If feeder routes must be ``'straight'`` or can be detoured (``'segmented'``)."
STRAIGHT = auto()
SEGMENTED = auto()
DEFAULT = SEGMENTED
[docs]
class FeederLimit(StrEnum):
"""Whether to limit the number of feeders. Both ``'specified'`` (an upper
bound) and ``'exactly'`` (an exact count) require the additional kwarg
``'max_feeders'``. Option ``'balanced'`` is only enforceable if the feeder
count is pinned to a single value, i.e. ``'minimum'`` or ``'exactly'``.
"""
UNLIMITED = auto()
EXACTLY = auto()
SPECIFIED = auto()
MINIMUM = auto()
MIN_PLUS1 = auto()
MIN_PLUS2 = auto()
MIN_PLUS3 = auto()
DEFAULT = UNLIMITED
def feeder_and_load_bounds(
T: int,
capacity: int,
feeder_limit: FeederLimit,
max_feeders: int,
balanced: bool,
feeders_per_subtree: int = 1,
) -> tuple[int, int | None, int | None, int | None]:
"""Derive the feeder-count and feeder-load bounds a model must enforce.
Bounds are returned in units of **subtrees** (the quantity the model's
feeder-sum constraint counts): one feeder per subtree for radial/branched
topologies, but a RINGED subtree is a cycle with two feeders. The
caller signals this with ``feeders_per_subtree`` (2 for RINGED), so that a
user-supplied ``max_feeders`` — always expressed as the number of physical
**substation connections** — is converted to the subtree count the model
constrains, and the returned bounds are multiplied back by
``feeders_per_subtree`` for reporting.
The subtree count is bounded below by ``min_feeders = ceil(T/capacity)``
regardless of ``feeder_limit`` (a valid inequality). ``feeders_ub`` is
``None`` when the count is unbounded above; when it equals ``feeders_lb``,
the count is pinned and callers should emit an equality constraint.
Balanced subtrees (loads differing at most by one unit) are only expressible
with a pinned feeder count ``F``, in which case the loads must lie in
``{T // F, ceil(T / F)}``. A load bound of ``None`` means "do not emit":
either ``balanced`` is off, or it is not enforceable (a warning is issued),
or the bound is already implied by the flow variable's own bounds.
Returns:
``(feeders_lb, feeders_ub, load_lb, load_ub)`` in subtree-count units.
"""
if feeders_per_subtree != 1 and max_feeders:
if max_feeders % feeders_per_subtree:
raise ValueError(
f'max_feeders ({max_feeders}) must be a multiple of '
f'{feeders_per_subtree} for a RINGED topology (each ring uses '
f'{feeders_per_subtree} substation connections)'
)
max_feeders //= feeders_per_subtree
min_feeders = math.ceil(T / capacity)
if feeder_limit is FeederLimit.UNLIMITED:
feeders_lb, feeders_ub = min_feeders, None
elif feeder_limit is FeederLimit.MINIMUM:
feeders_lb = feeders_ub = min_feeders
elif feeder_limit is FeederLimit.EXACTLY:
if max_feeders < min_feeders:
raise ValueError('max_feeders is below the minimum necessary')
if max_feeders > T:
raise ValueError('max_feeders is above the number of terminals')
feeders_lb = feeders_ub = max_feeders
elif feeder_limit is FeederLimit.SPECIFIED:
if max_feeders < min_feeders:
raise ValueError('max_feeders is below the minimum necessary')
feeders_lb, feeders_ub = min_feeders, max_feeders
elif feeder_limit in (
FeederLimit.MIN_PLUS1,
FeederLimit.MIN_PLUS2,
FeederLimit.MIN_PLUS3,
):
plus = int(feeder_limit.value[-1])
feeders_lb, feeders_ub = min_feeders, min_feeders + plus
else:
raise NotImplementedError('Unknown value:', feeder_limit)
if not balanced:
return feeders_lb, feeders_ub, None, None
if feeders_lb != feeders_ub:
warn(
'Model option <balanced = True> requires a single possible feeder'
f' count, but <feeder_limit = {feeder_limit.value.upper()}> allows a'
' range: model will not enforce balanced subtrees.'
)
return feeders_lb, feeders_ub, None, None
F = feeders_lb
load_lb, load_ub = T // F, math.ceil(T / F)
# bounds at the extremes are already implied by the flow variable's bounds
return (
feeders_lb,
feeders_ub,
load_lb if load_lb > 1 else None,
load_ub if load_ub < capacity else None,
)
[docs]
class ModelOptions(dict):
"""Hold options for the modelling of the cable routing problem.
Use ModelOptions.help() to get the options and their permitted and default
values. Use ModelOptions() without any parameters to use the defaults.
"""
hints = {
_identifier_from_class_name(kind): kind
for kind in (Topology, FeederRoute, FeederLimit)
}
# this has to be kept in sync with make_min_length_model()
simple = dict(
balanced=(
bool,
False,
'Whether to enforce balanced subtrees (subtree loads differ at most '
'by one unit).',
),
max_feeders=(
int,
0,
'Number of feeders: the maximum if <feeder_limit = "specified">, '
'the exact count if <feeder_limit = "exactly">',
),
)
# `with_signature` rewrites `__init__` at import time so that `help()`,
# `inspect.signature` and IDE introspection show the real options. Type
# checkers cannot follow that, so the same signature is declared statically
# under TYPE_CHECKING -- keep the two in sync (they are both derived from
# `hints` and `simple`).
if TYPE_CHECKING:
def __init__(
self,
*,
topology: Topology | str = ...,
feeder_route: FeederRoute | str = ...,
feeder_limit: FeederLimit | str = ...,
balanced: bool = ...,
max_feeders: int = ...,
) -> None: ...
else:
@with_signature(
'__init__(self, *, '
+ ', '.join(
chain(
(
f'{k}: {v.__name__} = "{v.DEFAULT.value}"'
for k, v in hints.items()
),
(
f'{name}: {kind.__name__} = {default}'
for name, (kind, default, _) in simple.items()
),
)
)
+ ')'
)
def __init__(self, **kwargs):
# dispatch on the key, not on the value's type: a str passed for an
# int option is a type error, not an enum to look up
for k, v in kwargs.items():
kind = self.hints.get(k)
if kind is not None:
kwargs[k] = kind(v)
else:
expected = self.simple[k][0]
if not isinstance(v, expected):
raise TypeError(
f'{k} must be {expected.__name__}, got '
f'{type(v).__name__}: {v!r}'
)
super().__init__(kwargs)
# Options are a value object: every value is coerced and every key is
# present once __init__ returns, and the whole library reads them on that
# basis (`is Topology.RINGED` is false for a plain 'ringed' str). Mutating
# the mapping afterwards would bypass the coercion, so it is refused --
# build a new instance instead.
def _immutable(self, *args, **kwargs):
raise TypeError(
f'{type(self).__name__} is immutable: construct a new one instead of '
'modifying it in place'
)
__setitem__ = __delitem__ = __ior__ = _immutable
clear = pop = popitem = setdefault = update = _immutable
@classmethod
def _from_pickle(cls, values):
return cls(**values)
[docs]
def __reduce__(self):
return type(self)._from_pickle, (dict(self),)
[docs]
@classmethod
def help(cls):
for k, v in cls.hints.items():
doc = indent(cleandoc(v.__doc__ or ''), ' ')
print(
f'{k} in {{'
+ ', '.join(
f'"{m}"' for n, m in v.__members__.items() if n != 'DEFAULT'
)
+ f'}} default: {cls.hints[k].DEFAULT.value}\n'
f'{doc}\n'
)
for name, (kind, default, desc) in cls.simple.items():
print(f'{name} [{kind.__name__}] default: {default}\n {desc}\n')
_Link = tuple[int, int]
[docs]
@dataclass(slots=True)
class SolutionInfo:
runtime: float
bound: float
objective: float
relgap: float
termination: str
def warmstart_links(
metadata: ModelMetadata, S: nx.Graph
) -> Iterator[tuple[Any, Any | None, int]]:
"""Yield ``(link_var, flow_var, flow)`` for every link ``S`` activates.
Uses ``metadata``'s own link/flow maps as the index, orienting ``S`` the way
the flow formulation expects — radializing RINGED rings into directed chains
— via :func:`.interarraylib.directed_links`. ``flow_var`` is ``None`` for a
ring's closing feeder (a link variable with no flow one). Everything not
yielded is inactive (0); each backend seeds that however its warm-start API
requires, then applies these.
Raises:
OWNWarmupFailed: ``S`` uses a link the model lacks.
"""
for source, sink, flow in directed_links(S):
key = (source, sink)
if key not in metadata.link_:
raise OWNWarmupFailed(f'warmup_model() failed: model lacks S link {key}')
yield (
metadata.link_[key],
metadata.flow_[key] if key in metadata.flow_ else None,
flow,
)
def _finalize_ringed_mip_S(
S: nx.Graph,
closing_links: list[tuple[int, int]],
A: nx.Graph | None,
) -> int:
"""Close decoded MILP paths into canonical rings and recalculate their loads."""
rings: list[tuple[int, int, int, int, int, int]] = []
for close_root, tail in closing_links:
subtree = S.nodes[tail]['subtree']
ordered = [tail]
prev, curr = None, tail
while True:
nxts = [n for n in S[curr] if n != prev]
if len(nxts) != 1:
raise ValueError(
f'RINGED MILP path at terminal {curr} has {len(nxts)} continuations'
)
nxt = nxts[0]
if nxt < 0:
head_root = nxt
break
ordered.append(nxt)
prev, curr = curr, nxt
ordered.reverse()
head, n = ordered[0], len(ordered)
if n == 1:
if close_root != head_root:
S.add_edge(close_root, tail, load=0, reverse=False)
else:
m = _ring_split_position(ordered, A)
u, v = ordered[m - 1], ordered[m]
S[u][v].update(load=0, reverse=False)
# bfs_subtree_loads overwrites this placeholder with the arm load.
S.add_edge(close_root, tail, load=1, reverse=False)
rings.append((subtree, head_root, head, close_root, tail, n))
for nodeD in S.nodes.values():
nodeD.pop('load', None)
roots = set(range(-S.graph['R'], 0))
visited = roots.copy()
for root in roots:
S.nodes[root]['load'] = 0
for subtree, head_root, head, close_root, tail, n in rings:
bfs_subtree_loads(S, head_root, [head], subtree, visited)
if n > 1:
bfs_subtree_loads(S, close_root, [tail], subtree, visited)
max_load = max(
(edgeD['load'] for u, v, edgeD in S.edges(data=True) if u < 0 or v < 0),
default=0,
)
S.graph['max_load'] = max_load
S.graph['has_loads'] = True
return max_load
[docs]
class Solver(abc.ABC):
"Common interface to multiple MILP solvers"
name: str
metadata: ModelMetadata
solver: Any
options: dict[str, Any]
stopping: dict[str, Any]
model_options: ModelOptions
solution_info: SolutionInfo
applied_options: dict[str, Any]
@abc.abstractmethod
def _link_val(self, var: Any) -> int | bool:
"Get the value of a link variable from the current solution."
pass
@abc.abstractmethod
def _flow_val(self, var: Any) -> int:
"Get the value of a flow variable from the current solution."
pass
[docs]
@abc.abstractmethod
def set_problem(
self,
P: nx.PlanarEmbedding,
A: nx.Graph,
capacity: int,
model_options: Mapping[str, Any],
warmstart: nx.Graph | None = None,
):
"""Define the problem geometry, available edges and tree properties
Args:
P: planar embedding of the location
A: available edges for the location
capacity: maximum number of terminals in a subtree
model_options: tree properties - see ModelOptions.help()
warmstart: initial feasible solution to pass to solver
"""
pass
[docs]
@abc.abstractmethod
def solve(
self,
time_limit: float,
mip_gap: float,
options: dict[str, Any] = {},
verbose: bool = False,
) -> SolutionInfo:
"""Run the MILP solver search.
Args:
time_limit: maximum time (s) the solver is allowed to run.
mip_gap: relative difference from incumbent solution to lower bound
at which the search may be stopped before ``time_limit`` is reached.
options: additional options to pass to solver (see solver manual).
Returns:
General information about the solution search (use ``get_solution()`` for
the actual solution).
"""
pass
[docs]
@abc.abstractmethod
def get_incumbent_topology(self) -> nx.Graph:
"""Return the best model-objective incumbent as topology ``S``.
This method does not route or rank solution-pool entries by detoured
length. Use :meth:`get_solution` for the routed, post-processed result.
"""
pass
[docs]
@abc.abstractmethod
def get_solution(self, A: nx.Graph | None = None) -> tuple[nx.Graph, nx.Graph]:
"""Output solution topology A and routeset G.
Args:
A: optionally replace the A given via set_problem() (if normalized A)
Returns:
Topology graph S and routeset G.
"""
pass
def _make_graph_attributes(self) -> dict[str, Any]:
metadata, solution_info = self.metadata, self.solution_info
solver_details = self.applied_options.copy()
# the method_options dict is extracted by db utility function packmethod()
method_options = dict(
solver_name=self.name,
fun_fingerprint=metadata.fun_fingerprint,
**self.stopping,
**metadata.model_options,
)
# remaining graph attributes (key=value) are stored in db.RouteSet[].misc
attr = dict(
**asdict(solution_info),
method_options=method_options,
solver_details=solver_details,
)
if 'max_feeders' in method_options:
solver_details['max_feeders'] = method_options.pop('max_feeders')
if metadata.warmed_by:
attr['warmstart'] = metadata.warmed_by
return attr
def _topology_from_mip_sol(self):
"""Create a topology graph from the solution to the MILP model.
Returns:
Graph topology ``S`` from the solution.
"""
metadata = self.metadata
topology = metadata.model_options['topology']
R = metadata.R
S = nx.Graph(R=R, T=metadata.T)
# ensure roots are added, even if some are not connected
S.add_nodes_from(range(-R, 0))
# Get active links and if flow is reversed (i.e. from small to big)
rev_from_link = {
(u, v): u < v
for (u, v), var in metadata.link_.items()
if self._link_val(var)
}
flow_links = {
link: reverse
for link, reverse in rev_from_link.items()
if link in metadata.flow_
}
S.add_weighted_edges_from(
((u, v, self._flow_val(metadata.flow_[u, v])) for u, v in flow_links),
weight='load',
)
nx.set_edge_attributes(S, flow_links, name='reverse')
# Propagate the decoded flow-tree attributes for every topology. A RINGED
# model is a path forest until its flowless starsʹ links are restored.
subtree = -1
max_load = 0
for r in range(-R, 0):
for u, v in nx.edge_dfs(S, r):
S.nodes[v]['load'] = S[u][v]['load']
if u == r:
subtree += 1
S.nodes[v]['subtree'] = subtree
rootload = 0
for nbr in S.neighbors(r):
subtree_load = S.nodes[nbr]['load']
max_load = max(max_load, subtree_load)
rootload += subtree_load
S.nodes[r]['load'] = rootload
if topology is Topology.RINGED:
closing_links = [
link for link in rev_from_link if link not in metadata.flow_
]
max_load = _finalize_ringed_mip_S(
S, closing_links, getattr(self, 'A', None)
)
S.graph.update(
topology=topology,
capacity=metadata.capacity,
max_load=max_load,
has_loads=True,
creator='MILP.' + self.name,
solver_details={},
)
return S
class PoolHandler(abc.ABC):
name: str
num_solutions: int
model_options: ModelOptions
solution_info: SolutionInfo
@abc.abstractmethod
def _objective_at(self, index: int) -> float:
"Get objective value from solution pool at position ``index``"
pass
@abc.abstractmethod
def _topology_from_mip_pool(self) -> nx.Graph:
"Build topology from the pool solution at the last requested position"
pass
def _incumbent_topology_from_pool(self) -> nx.Graph:
"""Decode pool entry zero, which must hold the best model objective."""
try:
expected_objective = self.solution_info.objective
except AttributeError as exc:
exc.args += ('.solve() must be called before solution retrieval',)
raise
objective = self._objective_at(0)
if not math.isclose(objective, expected_objective, rel_tol=1e-6, abs_tol=1e-6):
raise ValueError(
'Best solution-pool objective does not agree with solution_info: '
f'{objective} != {expected_objective}'
)
return self._topology_from_mip_pool()
def _investigate_pool(
self, P: nx.PlanarEmbedding, A: nx.Graph
) -> tuple[nx.Graph, nx.Graph]:
"""Go through the solver's solutions checking which has the shortest length
after applying the detours with PathFinder."""
Λ = float('inf')
S = G = None
num_solutions = self.num_solutions
info(f'Solution pool has {num_solutions} solutions.')
for i in range(num_solutions):
λ = self._objective_at(i)
if λ > Λ:
info(
f"#{i} halted pool search: objective ({λ:.3f}) > incumbent's length"
)
break
Sʹ = self._topology_from_mip_pool()
Gʹ = PathFinder(G_from_S(Sʹ, A), planar=P, A=A).create_detours()
Λʹ = Gʹ.size(weight='length')
if Λʹ < Λ:
S, G, Λ = Sʹ, Gʹ, Λʹ
G.graph['pool_entry'] = i, λ
info(f'#{i} -> incumbent (objective: {λ:.3f}, length: {Λ:.3f})')
else:
info(f'#{i} discarded (objective: {λ:.3f}, length: {Λ:.3f})')
if S is None or G is None:
raise ValueError('Solution pool has no usable solution.')
G.graph['pool_count'] = num_solutions
return S, G