mixle.analysis.biodiversity module

Biodiversity impact, habitat connectivity, and reclamation offsets (workstream N; N4 + N6).

Two related pieces share this module because they share the same underlying object – an N1 HabitatModel’s fitted suitability field:

  • N4 – connectivity (resistance_raster(), least_cost_corridor(), effective_conductance() / habitat_connectivity(), max_flow_connectivity(), fragmentation_impact()): graph resistance on a habitat-cost raster. Suitability is inverted into a per-cell movement cost; mixle.relations.ShortestPath gives the cheapest corridor between two patches, and the equivalent conductance graph feeds two genuinely different connectivity metrics: effective_conductance() (aliased as habitat_connectivity()) solves the Dirichlet/Laplacian boundary-value problem that circuit theory actually means by “effective conductance” / “effective current” (McRae et al. 2008’s Circuitscape method – fix the source patch cells at potential 1 and the sink patch cells at potential 0, solve for the harmonic potential everywhere else, and read off the current the network carries at that unit potential difference), while max_flow_connectivity() is the separate, bottleneck-only maximum-flow quantity over the same conductance graph. An earlier revision used max-flow as the circuit-theory metric (MXR-080-0071): that is wrong, because series conductances add reciprocally while max-flow only sees the tightest cut – the two agree only on a single-edge network. Use habitat_connectivity() for the population-viability proxy; max_flow_connectivity() is kept only because bottleneck throughput is occasionally independently useful on its own terms, never as a stand-in for the other. fragmentation_impact() scores a candidate development footprint (a mine footprint, in this module’s H3 worked instantiation) by how much it raises corridor resistance / drops habitat_connectivity() relative to the undisturbed baseline – the population-viability proxy N6 and J’s objective consume. Nothing about the connectivity graph is mining-specific: the same resistance-raster/least-cost-corridor/effective-conductance construction is standard landscape-ecology practice for any development footprint (a road, a subdivision, a mine) cutting through habitat.

  • N6 – reclamation ecology & biodiversity offsets (habitat_offset_liability(), no_net_loss_constraint()): prices the habitat impact of a development footprint as a liability the same shape as J6’s other priced terms (reclamation/remediation, health, carbon) and emits the companion no-net-loss hard constraint, so biodiversity offsets trade off against grade/cost/carbon inside ONE risk-adjusted objective instead of being a separate side calculation. habitat_offset_liability/no_net_loss_constraint work off the same “lost habitat-hectare- equivalents” quantity: the fitted suitability field (HabitatModel.mean, i.e. lambda_c) times per-cell area, summed over whatever footprint of cells a candidate disturbance plan (a mine plan, in this module’s worked instantiation) disturbs. Suitability, area, and every economic rate/cost feeding these two functions are validated finite and sign-correct before they reach a solver or accounting payload (MXR-080-0073) – see _lost_equivalents() and the _require_* helpers below.

Every function here reads only duck-typed attributes off habitat (.mean, optionally .cell_area) or takes plain arrays, so anything satisfying the IC-1 Posterior surface over a suitability field – in particular N1’s HabitatModel – works; the HabitatModel type hint is a forward reference (evaluated only under TYPE_CHECKING), so this module has no hard runtime dependency on mixle.analysis.sdm. effective_conductance() / habitat_connectivity() use SciPy’s sparse linear solver; max_flow_connectivity() and fragmentation_impact() use this module’s sparse float-capacity push-relabel implementation. Only least_cost_corridor() calls the frozen ShortestPath surface, imported lazily because mixle.relations transitively imports mixle.stats and this module can be loaded while the statistics/analysis import graph is still initializing.

resistance_raster(habitat, *, floor=1e-3)[source]

Per-cell movement cost from a fitted suitability field: cost_c = 1 / max(suitability_c, floor).

Reads only habitat.mean (N1’s fitted intensity field lambda_c); floor keeps near-zero suitability cells at a large-but-finite cost (near-impermeable) rather than blowing up to inf, so the resulting raster is always safe to feed straight into least_cost_corridor() / habitat_connectivity() without any further cleanup.

Raises:

ValueError – if floor is not finite and strictly positive, or if habitat.mean contains a NaN/Inf/negative entry. An earlier revision silently clamped negative suitability into the same near-impermeable cost as zero suitability (MXR-080-0073), hiding invalid upstream model output; a negative suitability now always raises instead of being clamped. Zero suitability remains a legitimate input (clamped to floor, as documented above).

Parameters:

  • habitat (HabitatModel)

  • floor (float)

Return type:

np.ndarray

least_cost_corridor(resistance, patch_a, patch_b)[source]

Cheapest movement path between two cells of a resistance raster (the least-cost corridor).

resistance is an n-dimensional cost raster (e.g. from resistance_raster()); patch_a/ patch_b are flat (row-major/np.ravel) cell indices into it. Delegates to mixle.relations.ShortestPath over the rook-adjacency grid graph with each edge weighted by _edge_cost() (mean resistance of its two endpoint cells); a cell resistance of inf (a mined footprint, fragmentation_impact()) makes every edge touching it non-finite and so unusable, effectively removing that cell from the graph. Lower resistance = better connected.

ShortestPath’s underlying engine (mixle.relations.best_first_paths) is a pure best-first search with no closed/visited set – intentional there, since it is also used to enumerate k best (not just the single best) paths, where revisiting a state via a costlier path is legitimate output. On a cyclic grid graph asked for only its single best path (k=1), that same lack of a closed set is combinatorially explosive (every back-and-forth detour is a distinct, never-pruned path). Since successors(state) is called exactly once per pop and pops come out in non-decreasing cost order, the first call for a given state is provably its cheapest arrival (Dijkstra’s invariant, valid here because every edge cost is non-negative) – so a small closed set local to this call, never touching relations.py, turns the search back into ordinary Dijkstra: each state’s successors are computed once and [] (already settled) on every later call.

Returns:

the total path resistance and the visited flat cell indices from patch_a to patch_b inclusive; (inf, []) if no finite-cost path exists.

Return type:

(corridor_resistance, path)

Parameters:
effective_conductance(resistance, sources, sinks)[source]

Circuit-theory effective conductance between sources and sinks over a resistance raster.

Solves the standard Dirichlet/Laplacian boundary-value problem for landscape connectivity (McRae et al. 2008; the “Circuitscape” formulation): every sources cell is held at potential 1 and every sinks cell at potential 0 – equivalent to shorting each terminal’s cells together with an ideal conductor, with no arbitrary coupling constant needed (unlike a max-flow super-source/-sink). The remaining free cells’ potentials solve the reduced linear system L_free @ v_free = -L_free,fixed @ v_fixed for the graph Laplacian L (edge weight = conductance = 1 / edge_cost, _edge_cost()); the total current flowing out of the source terminal at this unit potential difference is the effective conductance (I = G @ V, V = 1), read off directly from the Laplacian: current out of node i is (L @ v)[i], which is exactly 0 at every free node (that is what the linear solve enforces) and the true injected/withdrawn current at every fixed node.

This is a genuinely different quantity from maximum flow (max_flow_connectivity()): conductances in series add reciprocally (two unit conductances in series give an effective conductance of 0.5, not 1), so a max-flow solver – which only sees bottleneck capacity and ignores series resistance entirely – systematically overstates connectivity whenever current has to cross more than one edge to reach the sink. On a unit-conductance square between opposite corners this function returns the textbook value 1; a max-flow solver over the same graph returns 2 (the sum of the two parallel corner-to-corner paths’ raw capacities, MXR-080-0071).

Cells with no path at all to either terminal (an isolated island cut off by +inf-resistance neighbours on every side) carry no current and are excluded from the linear solve rather than left to make it singular; a raster whose sources have no path to its sinks at all correctly yields 0.0 (the linear system still solves – each side simply settles to its own terminal’s potential everywhere reachable, so no potential gradient, hence no current, exists anywhere).

Raises:

ValueError – on a NaN/negative/zero resistance entry (+inf, the impassable-cell sentinel, is fine), an empty raster, an out-of-bounds source/sink index, or overlapping source/sink sets.

Parameters:
Return type:

float

habitat_connectivity(resistance, sources, sinks)[source]

Domain-friendly alias for effective_conductance(), the circuit-theory connectivity metric.

Kept under this name because it is the pre-existing public entry point N4/N6/J callers and tests already use; see effective_conductance() for the full Dirichlet/Laplacian derivation and its McRae et al. (2008) citation. Prior to MXR-080-0071 this function actually computed maximum flow (now max_flow_connectivity()) despite the circuit-theory claim in its docstring – a different, systematically-inflated quantity for any network with more than one edge between source and sink. It now computes the real thing.

Parameters:
Return type:

float

max_flow_connectivity(resistance, sources, sinks)[source]

Maximum-flow connectivity between sources and sinks over a resistance raster’s conductance graph.

Builds a conductance graph (arc capacity 1 / edge_cost between rook-adjacent cells) with a super-source over sources and a super-sink over sinks, then returns mixle.relations.max_flow()’s value between them – the maximum flow (bottleneck-limited throughput) the habitat network can carry from the source patch(es) to the sink patch(es).

This is NOT the circuit-theory effective-conductance quantity (effective_conductance() / habitat_connectivity()): max flow only sees the tightest cut and entirely ignores series resistance, so it systematically overstates connectivity whenever current must cross more than one edge to reach the sink (MXR-080-0071). This function exists as a distinctly-named, independently useful bottleneck-throughput metric in its own right – never use it as a stand-in for the other.

Parameters:
Return type:

float

fragmentation_impact(resistance, footprint_mask, sources, sinks)[source]

Habitat-connectivity impact of a development footprint: baseline vs. disturbed corridor/connectivity.

Sets every footprint_mask cell’s resistance to inf (impassable / zero conductance – an H3 mine plan, in this module’s worked instantiation) and recomputes both least_cost_corridor() (between the first sources/sinks cell, taken as the representative patch-to-patch corridor endpoints – a judgment call: the public API only takes cell sets here, not a single designated pair) and habitat_connectivity() (over the full sources/sinks sets) before and after. This is the population-viability proxy N6/J feed into the biodiversity-offset objective: a footprint that severs the only corridor raises corridor_resistance and drops connectivity sharply; one that misses every real corridor leaves both essentially unchanged.

Returns:

A dict with corridor_resistance_baseline/corridor_resistance_mined, connectivity_baseline/connectivity_mined, delta (connectivity_baseline - connectivity_mined, i.e. connectivity lost), and mincut_edges – the baseline network’s minimum-cut arcs (mixle.relations.min_cut()), i.e. the single weakest link a footprint would need to sever to maximally damage connectivity.

Parameters:
Return type:

dict

habitat_offset_liability(plan_footprint, habitat, *, offset_ratio, unit_offset_cost)[source]

Priced biodiversity-offset liability of disturbing plan_footprint (a J6 priced-objective term).

lost_equivalents = sum_{c in footprint} suitability_c * area_c (suitability = N1’s fitted HabitatModel.mean); the liability is offset_ratio * lost_equivalents * unit_offset_cost – an additive dollar term the same shape J6’s priced_liabilities already sums for carbon/health/ remediation (workstream-J.md J6). offset_ratio=0 or unit_offset_cost=0 reduces this to zero, i.e. no biodiversity-offset requirement – both remain valid, meaningful inputs and are not rejected.

Because lost_equivalents is linear in the boolean footprint, the per-cell rate offset_ratio * unit_offset_cost * suitability_c * area_c is itself a valid per-block deduction a MILP-based optimizer (H4/J6’s risk_adjusted_plan) can net directly out of expected per-block profit – this function is the scalar evaluator for a given (candidate or solved) footprint.

Raises:

ValueError – if offset_ratio/unit_offset_cost is negative or non-finite, or if _lost_equivalents rejects plan_footprint/habitat. A negative ratio or cost used to silently turn habitat damage into a “profit” (MXR-080-0073); both are now rejected. Also if the priced product itself is not representable in float64 – finite inputs near the float64 ceiling used to publish an infinite liability (MXR-080-1572).

Parameters:

  • plan_footprint (np.ndarray)

  • habitat (HabitatModel)

  • offset_ratio (float)

  • unit_offset_cost (float)

Return type:

float

no_net_loss_constraint(plan_footprint, habitat, *, offset_ratio)[source]

Hard no-net-loss constraint payload: created/purchased offsets >= offset_ratio * lost_equivalents.

Returns a dict carrying both the raw quantities (lost_equivalents, per_cell_lost_equivalents, required_offset) and a solver-agnostic linear-constraint row in this repo’s standard coeffs @ x <= bound convention (mixle.relations/mixle.stochastic_opt’s a_ub rows), expressed over a single offsets_created decision variable: coeffs=[-1.0], bound=-required_offset encodes -offsets_created <= -required_offset, i.e. offsets_created >= required_offset. Placing that row (and the offsets_created column it references) into the wider extraction/offset-purchase decision space is H4/J6’s job – this module never edits their MILP variable indexing, only hands them the row.

Because that row is emitted ALREADY normalized into the <= convention, its "sense" label is "<=" – the label describes the row that is actually here, not the inequality the row expresses about offsets_created. It used to read ">=" (MXR-080-1570), which made the payload internally contradictory: mixle.analysis.objective.hard_constraints() and mixle.stochastic_opt.risk_adjusted_plan() both trust the label and negate a ">=" row into <= form, so a second negation turned the no-net-loss FLOOR into the CEILING offsets_created <= required_offset – permitting zero offsets and forbidding any offset purchase above the bare minimum, i.e. exactly the opposite of a no-net-loss guarantee.

Raises:

ValueError – if offset_ratio is negative or non-finite, or if _lost_equivalents rejects plan_footprint/habitat (NaN/negative suitability or area). A NaN or negative lost_equivalents used to be able to propagate into required_offset, silently weakening or inverting a constraint that is supposed to be a hard floor (MXR-080-0073); both inputs are now validated before this payload is built. Also if required_offset itself is not representable in float64 – finite inputs near the float64 ceiling used to publish an infinite bound, which is not a hard constraint but an unsatisfiable row (MXR-080-1572).

Parameters:

  • plan_footprint (np.ndarray)

  • habitat (HabitatModel)

  • offset_ratio (float)

Return type:

dict