"""Black-box optimization of IIT quantities over substrate parameters.
Where :mod:`pyphi.landscape` analyzes the φ landscape along one axis, this
module searches it: :func:`optimize` runs a seeded population method over a
bounded box of connection weights, looking for a substrate that maximizes an
IIT quantity — by default the signed normalized system irreducibility φₛ,
which stays continuous across minimum-information-partition switches and so
gives a gradient-free search no discontinuities to trip on.
:func:`weight_axes` builds the search space for the common case: a map from a
parameter vector to a :func:`~pyphi.substrate_generator.build_substrate`
substrate, varying a chosen set of weight-matrix entries.
"""
from __future__ import annotations
import contextlib
import math
from collections.abc import Callable
from collections.abc import Sequence
from dataclasses import dataclass
from typing import Any
import numpy as np
import pandas as pd
from numpy.typing import NDArray
from pyphi.conf import config
from pyphi.conf import fallback
from pyphi.conf._helpers import check_parallel_flag
from pyphi.conf._helpers import check_progress_flag
from pyphi.display import Description
from pyphi.display import Displayable
from pyphi.display import Row
from pyphi.display import Section
from pyphi.display.numbers import format_value
from pyphi.landscape import _QUANTITIES
from pyphi.landscape import _UNREACHABLE
from pyphi.landscape import _optional_float
from pyphi.landscape import _part_id
from pyphi.serializable import Serializable
# Rank-based DE selection only needs unreachable to compare worse than any real
# candidate; a fixed large sentinel (internal minimization convention) suffices.
_UNREACHABLE_PENALTY = 1e18
[docs]
def weight_axes(
unit_functions: Any,
weights: NDArray[Any],
indices: Sequence[tuple[int, int]],
**kwargs: Any,
) -> Callable[[NDArray[Any]], Any]:
"""Return a parameter axis varying several weights of a generated substrate.
The vector analogue of :func:`pyphi.landscape.weight_axis`.
Parameters
----------
unit_functions : str or Callable or Iterable
Unit function(s), as accepted by
:func:`pyphi.substrate_generator.build_substrate`.
weights : ArrayLike
The base weight matrix; ``weights[i, j]`` is the connection from unit
``i`` to unit ``j``. Copied on every call — never mutated.
indices : Sequence[tuple[int, int]]
The ``(i, j)`` entries to vary, in the order the parameter vector
indexes them.
Returns
-------
_WeightAxis
A picklable callable mapping a length-``len(indices)`` vector θ to the
substrate built from the weight matrix with ``weights[indices[k]] =
θ[k]``. Additional keyword arguments forward to ``build_substrate`` on
every call.
Notes
-----
The returned axis is a picklable object, not a closure, so a population can
be evaluated across worker processes (``optimize(..., parallel=True)``).
Setting a weight to exactly 0 removes the connection from the derived
connectivity matrix (``cm = weights != 0``), a discrete topology change,
exactly as for :func:`pyphi.landscape.weight_axis`.
"""
from pyphi.substrate_generator import _check_unit_function_kwargs
_check_unit_function_kwargs(unit_functions, kwargs)
return _WeightAxis(
unit_functions=unit_functions,
base=np.array(weights, dtype=float),
entries=[(int(i), int(j)) for i, j in indices],
kwargs=dict(kwargs),
)
@dataclass(frozen=True)
class _WeightAxis:
"""A picklable vector → substrate map varying a fixed set of weight entries.
A module-level callable rather than a closure so the process backend can
pickle it when a population is evaluated in parallel.
"""
unit_functions: Any
base: NDArray[Any]
entries: list[tuple[int, int]]
kwargs: dict[str, Any]
def __call__(self, theta: NDArray[Any]) -> Any:
from pyphi.substrate_generator import build_substrate
varied = self.base.copy()
for (i, j), value in zip(self.entries, np.asarray(theta), strict=True):
varied[i, j] = value
return build_substrate(self.unit_functions, varied, **self.kwargs)
def _objective_value(sia: Any, objective: Any) -> float:
"""The natural (maximization-convention) objective scalar for one SIA."""
if callable(objective):
value: Any = objective(sia)
return float(value)
try:
value = getattr(sia, objective)
except AttributeError:
raise ValueError(
f"quantity {objective!r} is not available on "
f"{type(sia).__name__}; choose a quantity the requested "
f"formalism provides (e.g. quantity='phi')"
) from None
return _optional_float(value)
def _eval_one(
theta: NDArray[Any],
*,
builder: Callable[[NDArray[Any]], Any],
state: tuple[int, ...],
subset: Sequence[int] | None,
formalism: str | None,
quantity: Any,
) -> dict[str, Any]:
"""Evaluate one candidate; unreachable states are penalized, not raised.
``theta`` is first-positional and the rest keyword-only so a
``functools.partial`` over the fixed arguments is a picklable one-argument
worker for ``map_reduce``.
"""
from pyphi.analysis import analyze
from pyphi.direction import Direction
try:
sia: Any = analyze(
builder(np.asarray(theta)),
state,
subset=subset,
formalism=formalism,
compute="sia",
)
except _UNREACHABLE:
return {
"objective": math.nan,
"reachable": False,
"partition": None,
"cause_state": None,
"effect_state": None,
"partition_margin": math.nan,
"cause_state_margin": math.nan,
"effect_state_margin": math.nan,
"_sia": None,
}
# Selection margins and state specifications exist only on IIT 4.0
# SIAs; other formalisms' rows carry None/NaN in those columns.
system_state = getattr(sia, "system_state", None)
cause = system_state.cause if system_state is not None else None
effect = system_state.effect if system_state is not None else None
margins = getattr(sia, "state_margins", None)
return {
"objective": _objective_value(sia, quantity),
"reachable": True,
"partition": _part_id(sia.partition),
"cause_state": None if cause is None else tuple(int(x) for x in cause.state),
"effect_state": None if effect is None else tuple(int(x) for x in effect.state),
"partition_margin": _optional_float(getattr(sia, "partition_margin", None)),
"cause_state_margin": _optional_float(
margins[Direction.CAUSE] if margins is not None else None
),
"effect_state_margin": _optional_float(
margins[Direction.EFFECT] if margins is not None else None
),
"_sia": sia,
}
[docs]
@dataclass(frozen=True, eq=False, repr=False)
class OptimizationResult(Displayable, Serializable):
"""The outcome of an :func:`optimize` run.
Attributes
----------
best_params : NDArray
The winning parameter vector.
best_objective : float
Its objective value.
best_substrate : Substrate
``builder(best_params)``.
best_sia : SystemIrreducibilityAnalysis
The analysis at the best point — its selection identity and margins say
how robust the winning substrate is.
trajectory : pd.DataFrame
One row per evaluation: ``eval``, ``generation``, one column per
parameter dimension (``p0``, ``p1``, …), ``objective``,
``reachable``, the selection identity (``partition``,
``cause_state``, ``effect_state``) and the selection margins
(``partition_margin``, ``cause_state_margin``, ``effect_state_margin``).
bounds : list[tuple[float, float]]
The search box.
seed : int
The seed driving the run.
goal : str
``"maximize"`` or ``"minimize"``.
quantity_name : str
The optimized quantity's name, or ``"<callable>"``.
settings : dict
Backend and budget: ``backend``, ``popsize``, ``maxiter``, ``tol``.
config_snapshot : dict
PyPhi configuration at run time (``precision``, ``formalism``).
n_evaluations : int
Total candidate evaluations.
n_unreachable : int
How many candidates had a dynamically unreachable state.
"""
best_params: NDArray[Any]
best_objective: float
best_substrate: Any
best_sia: Any
trajectory: pd.DataFrame
bounds: list[tuple[float, float]]
seed: int
goal: str
quantity_name: str
settings: dict[str, Any]
config_snapshot: dict[str, Any]
n_evaluations: int
n_unreachable: int
def __eq__(self, other: object) -> bool:
# The generated dataclass __eq__ would compare ``best_params`` and
# ``trajectory`` with ``==``, which yields elementwise results whose
# truth values raise.
if not isinstance(other, OptimizationResult):
return NotImplemented
return (
np.array_equal(self.best_params, other.best_params)
and self.trajectory.equals(other.trajectory)
and self.best_objective == other.best_objective
and self.best_substrate == other.best_substrate
and self.best_sia == other.best_sia
and self.bounds == other.bounds
and self.seed == other.seed
and self.goal == other.goal
and self.quantity_name == other.quantity_name
and self.settings == other.settings
and self.config_snapshot == other.config_snapshot
and self.n_evaluations == other.n_evaluations
and self.n_unreachable == other.n_unreachable
)
def to_pandas(self) -> pd.DataFrame:
return self.trajectory
def _describe(self, verbosity: int) -> Description: # noqa: ARG002
params = ", ".join(format_value(float(p)) for p in self.best_params)
rows = (
Row("Goal", f"{self.goal} {self.quantity_name}"),
Row("Best objective", self.best_objective),
Row("Best parameters", f"({params})"),
Row("Evaluations", self.n_evaluations),
Row("Unreachable", self.n_unreachable),
Row("Seed", self.seed),
)
return Description(
title="OptimizationResult",
sections=(Section(rows=rows),),
compact=(
f"OptimizationResult({self.goal} {self.quantity_name}: "
f"{format_value(float(self.best_objective))})"
),
)
def _eval_batch(
builder: Callable[[NDArray[Any]], Any],
thetas: list[NDArray[Any]],
state: tuple[int, ...],
subset: Sequence[int] | None,
formalism: str | None,
quantity: Any,
parallel: bool,
progress: bool,
) -> list[dict[str, Any]]:
"""Evaluate a whole population, sequentially or through ``map_reduce``.
The parallel branch mirrors ``pyphi.sweeps._run_cells_parallel``: one level
of parallelism (``parallel=True`` on the outer ``map_reduce``, inner
``parallel=False`` in the worker config snapshot), one candidate per cell.
"""
if not parallel or len(thetas) <= 1:
return [
_eval_one(
theta,
builder=builder,
state=state,
subset=subset,
formalism=formalism,
quantity=quantity,
)
for theta in thetas
]
from functools import partial
from pyphi.conf import presets
from pyphi.parallel import map_reduce
preset = presets.by_name[formalism] if formalism is not None else {}
# partial binds the fixed args as keywords; theta stays the one positional,
# so cell_fn is a picklable one-argument worker.
cell_fn = partial(
_eval_one,
builder=builder,
state=state,
subset=subset,
formalism=formalism,
quantity=quantity,
)
results: list[Any] = []
with config.override(parallel=False, **preset):
results = map_reduce(
cell_fn,
thetas,
parallel=True,
ordered=True,
reduce_func=list,
progress=progress,
desc="optimize population",
chunksize=1,
)
if len(results) != len(thetas):
raise AssertionError("map_reduce flattened population results")
return results
def _check_inputs(
builder: Callable[[NDArray[Any]], Any],
box: list[tuple[float, float]],
state: tuple[int, ...],
subset: Sequence[int] | None,
formalism: str | None,
) -> None:
"""Build the substrate at the midpoint of ``box`` and the candidate system
on it, so a wrong bounds length, state or subset raises its own error
before the search starts rather than inside the optimizer.
An unreachable state at the midpoint is not an input error: the search
penalizes unreachable candidates.
"""
from pyphi.conf import presets
from pyphi.system import System
midpoint = np.array([(lo + hi) / 2 for lo, hi in box])
try:
substrate = builder(midpoint)
except (ValueError, TypeError, IndexError) as error:
raise ValueError(
f"builder could not build a substrate from the midpoint of bounds "
f"({len(box)} parameters): {error}"
) from error
with presets.use(formalism), contextlib.suppress(*_UNREACHABLE):
System(substrate, state, subset=None if subset is None else tuple(subset))
[docs]
def optimize(
builder: Callable[[NDArray[Any]], Any],
state: tuple[int, ...],
bounds: Sequence[tuple[float, float]],
*,
seed: int,
quantity: Any = "signed_normalized_phi",
goal: str = "maximize",
x0: NDArray[Any] | None = None,
popsize: int = 15,
maxiter: int = 100,
tol: float = 0.01,
subset: Sequence[int] | None = None,
formalism: str | None = None,
parallel: bool | None = None,
progress: bool | None = None,
) -> OptimizationResult:
"""Search a substrate's weights for a maximizer of an IIT quantity.
Runs a seeded, population-based black-box optimizer
(:func:`scipy.optimize.differential_evolution`) over a bounded box of
parameters, evaluating each candidate as one system irreducibility
analysis. The default quantity, signed normalized φₛ, is continuous across
minimum-information-partition switches, so the gradient-free search sees no
discontinuities.
Parameters
----------
builder : Callable[[NDArray], Substrate]
The parameter axis: a vector → substrate map. :func:`weight_axes` builds
one for the single-substrate weight case.
state : tuple[int, ...]
The substrate state analyzed at every candidate.
bounds : Sequence[tuple[float, float]]
``(low, high)`` per parameter dimension. ``len(bounds)`` fixes the
search dimension.
seed : int
Seeds an isolated ``np.random.default_rng`` driving the optimizer, and
is recorded on the result. Required.
quantity : str or Callable
The quantity to optimize: ``"signed_normalized_phi"`` (default),
``"phi"``, ``"signed_phi"``, or ``"normalized_phi"`` by name, as in
:func:`~pyphi.landscape.perturb`, or a callable ``SIA → float`` for
any other quantity.
goal : str
``"maximize"`` (default) or ``"minimize"``.
x0 : NDArray, optional
Seeds one member of the initial population, e.g. a known-good
substrate's weights.
popsize : int
Differential-evolution population multiplier.
maxiter : int
Maximum number of generations.
tol : float
Relative convergence tolerance.
subset : Sequence[int], optional
Candidate system node indices; ``None`` uses the whole substrate.
formalism : str, optional
Leave unset to compute IIT under the active configuration. Pass an
earlier version's name only to reproduce published results; it
applies to these evaluations only, as in :func:`pyphi.analyze`.
parallel : bool, optional
Evaluate each generation's population through ``map_reduce``; defaults to
``config.infrastructure.parallel``. Requires ``builder`` and a callable
``quantity`` to be picklable (:func:`weight_axes` and the named
quantities are).
progress : bool, optional
Show a progress bar; defaults to ``config.infrastructure.progress_bars``.
Returns
-------
OptimizationResult
Raises
------
ValueError
If ``quantity`` is a name that is not one of the four φₛ variants;
if ``goal`` is neither ``"maximize"`` nor ``"minimize"``; or if
``builder`` cannot build a substrate from ``len(bounds)`` parameters,
or ``state`` or ``subset`` does not fit that substrate. These are
checked once, at the midpoint of ``bounds``, before the search.
Notes
-----
The driver reads only the scalar objective to choose a move. It never
branches on a selection identity (the chosen partition or state), because
that identity flips between near-tied candidates from one grid point to the
next even where the objective is smooth. Selection identities and margins
are logged for inspection only.
With ``parallel=True`` the ``builder`` and ``quantity`` are sent to worker
processes, so both must be picklable. A closure builder or a lambda
quantity must be run with ``parallel=False``.
"""
from scipy.optimize import differential_evolution
if isinstance(quantity, str) and quantity not in _QUANTITIES:
raise ValueError(
f"unknown quantity {quantity!r}; expected one of "
f"{', '.join(sorted(_QUANTITIES))}, or a callable"
)
from pyphi.conf import presets
formalism = presets.canonical(formalism)
if goal not in ("maximize", "minimize"):
raise ValueError(f"goal must be 'maximize' or 'minimize', got {goal!r}")
box = [(float(lo), float(hi)) for lo, hi in bounds]
_check_inputs(builder, box, state, subset, formalism)
sign = -1.0 if goal == "maximize" else 1.0 # DE minimizes sign*objective
show = fallback(check_progress_flag(progress), config.infrastructure.progress_bars)
use_parallel = fallback(
check_parallel_flag(parallel), config.infrastructure.parallel
)
quantity_name = quantity if isinstance(quantity, str) else "<callable>"
rows: list[dict[str, Any]] = []
# Only the best candidate's SIA is retained; keeping every candidate's
# full SIA would grow memory linearly with the number of evaluations.
best_index: int | None = None
best_value = 0.0
best_sia_seen: Any = None
generation = 0
def batch_objective(population: NDArray[Any]) -> NDArray[Any]:
nonlocal generation, best_index, best_value, best_sia_seen
# scipy hands a (D, S) matrix under vectorized=True.
thetas = [population[:, i] for i in range(population.shape[1])]
evaluated = _eval_batch(
builder, thetas, state, subset, formalism, quantity, use_parallel, show
)
scores = np.empty(len(evaluated))
for member, (theta, row) in enumerate(zip(thetas, evaluated, strict=True)):
record = {
"eval": len(rows),
"generation": generation,
**{f"p{d}": float(theta[d]) for d in range(len(theta))},
"objective": row["objective"],
"reachable": row["reachable"],
"partition": row["partition"],
"cause_state": row["cause_state"],
"effect_state": row["effect_state"],
"partition_margin": row["partition_margin"],
"cause_state_margin": row["cause_state_margin"],
"effect_state_margin": row["effect_state_margin"],
}
rows.append(record)
value = record["objective"]
if value is not None and not math.isnan(float(value)):
improved = (
float(value) > best_value
if goal == "maximize"
else float(value) < best_value
)
if best_index is None or improved:
best_index = len(rows) - 1
best_value = float(value)
best_sia_seen = row["_sia"]
if row["reachable"]:
scores[member] = sign * row["objective"]
else:
scores[member] = _UNREACHABLE_PENALTY
generation += 1
return scores
# scipy replaces an exception raised by the objective with its own
# RuntimeError about the map-like callable; keep the original to re-raise.
failures: list[Exception] = []
def checked_objective(population: NDArray[Any]) -> NDArray[Any]:
try:
return batch_objective(population)
except Exception as error:
failures.append(error)
raise
try:
outcome = differential_evolution(
checked_objective,
box,
rng=np.random.default_rng(seed),
maxiter=maxiter,
popsize=popsize,
tol=tol,
x0=x0,
vectorized=True,
updating="deferred",
polish=False,
workers=1,
)
except Exception:
if failures:
raise failures[0] from None
raise
trajectory = pd.DataFrame(rows)
# Derive the best from the logged evaluations (not `outcome`) so
# best_params / best_sia / best_objective stay mutually consistent and
# match the table. The running best skips NaN (unreachable) rows and
# keeps the first occurrence on exact ties; the default RangeIndex makes
# `best_index` line up with the trajectory rows.
param_cols = [f"p{d}" for d in range(len(box))]
if best_index is not None:
best_row = best_index
best_objective = best_value
best_params = trajectory.iloc[best_row][param_cols].to_numpy(dtype=float)
best_sia = best_sia_seen
else:
# Every candidate was dynamically unreachable — no meaningful optimum.
best_objective = math.nan
best_params = np.asarray(outcome.x, dtype=float)
best_sia = None
best_substrate = builder(best_params)
return OptimizationResult(
best_params=best_params,
best_objective=float(best_objective),
best_substrate=best_substrate,
best_sia=best_sia,
trajectory=trajectory,
bounds=box,
seed=int(seed),
goal=goal,
quantity_name=quantity_name,
settings={
"backend": "scipy.differential_evolution",
"popsize": popsize,
"maxiter": maxiter,
"tol": tol,
},
config_snapshot={
"precision": config.numerics.precision,
"formalism": formalism,
},
n_evaluations=len(trajectory),
n_unreachable=int(np.count_nonzero(~trajectory["reachable"].to_numpy())),
)
