MeanField

 71class RateSystem:
 72    """General mean-field solver with helper utilities for ERF and fixpoints.
 73
 74    Subclasses implement `_build_dynamics(...)` and then inherit the fixed-point
 75    solver, ERF sweep, and fixpoint analysis helpers.
 76    """
 77    # Minimum allowed variance to prevent numerical instabilities.
 78    # Chosen small enough to avoid affecting dynamics, only avoids divide-by-zero.
 79    VAR_EPS = 1e-12
 80    # Accept near-roots even when the underlying SciPy solver does not set success=True.
 81    RESIDUAL_ACCEPT_TOL = 5e-4
 82    # Smallest continuation step introduced when adaptively subdividing the ERF sweep.
 83    ADAPTIVE_CONTINUATION_MIN_STEP = 1e-3
 84
 85    def __init__(
 86        self,
 87        parameter: Dict,
 88        v_focus: float,
 89        *,
 90        focus_population: Optional[Union[int, Sequence[int]]] = None,
 91        prefer_jax: bool = True,
 92        root_tol: float = 1e-9,
 93        max_function_evals: int = 4000,
 94        max_newton_steps: Optional[int] = 1000,
 95        **network_kwargs,
 96    ) -> None:
 97        self.parameter = dict(parameter)
 98        self.v_focus = float(v_focus)
 99        self.prefer_jax = bool(prefer_jax)
100        self.root_tol = root_tol
101        self.max_function_evals = max_function_evals
102        default_steps = 256
103        if max_newton_steps is None:
104            self.max_steps = default_steps
105        else:
106            steps = int(max_newton_steps)
107            if steps <= 0:
108                raise ValueError("max_newton_steps must be positive.")
109            self.max_steps = steps
110        self.network_kwargs = dict(network_kwargs)
111        dynamics = self._build_dynamics(self.parameter, **self.network_kwargs)
112        if len(dynamics) == 4:
113            self.A, self.B, self.bias, self.tau = dynamics
114            self.C = np.zeros_like(self.B)
115        elif len(dynamics) == 5:
116            self.A, self.B, self.C, self.bias, self.tau = dynamics
117            if self.C is None:
118                self.C = np.zeros_like(self.B)
119        else:
120            raise ValueError("Expected _build_dynamics to return (A, B, bias, tau) or (A, B, C, bias, tau).")
121        self.population_count = int(self.A.shape[0])
122        if self.A.shape != self.B.shape or self.A.shape != self.C.shape:
123            raise ValueError("Connectivity mean and variance matrices must match in shape.")
124        if self.bias.shape[0] != self.population_count or self.tau.shape[0] != self.population_count:
125            raise ValueError("Bias and tau vectors must match the matrix dimensions.")
126        if focus_population is not None:
127            focus_config = focus_population
128        elif isinstance(self.parameter.get("focus_population"), (list, tuple, range)):
129            focus_config = self.parameter.get("focus_population")
130        else:
131            count = int(self.parameter.get("focus_count", 1) or 1)
132            focus_config = list(range(max(count, 1)))
133        self.focus_indices = self._resolve_focus_indices(focus_config)
134        self._initialize_group_constraints()
135        self.use_jax = self.prefer_jax and HAS_JAX and optx is not None
136        self._jax_args = None
137
138    # --- abstract hooks -------------------------------------------------
139    def _build_dynamics(
140        self, parameter: Dict, **network_kwargs
141    ) -> Union[
142        Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
143        Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray],
144    ]:
145        raise NotImplementedError("Derived classes must implement '_build_dynamics'.")
146
147    def _build_population_groups(self, focus: np.ndarray) -> List[np.ndarray]:
148        focus_set = set(int(idx) for idx in focus.tolist())
149        if not focus_set:
150            focus_set = {0}
151        groups: List[np.ndarray] = [np.array(sorted(focus_set), dtype=int)]
152        remaining = [idx for idx in range(self.population_count) if idx not in focus_set]
153        for idx in remaining:
154            groups.append(np.array([idx], dtype=int))
155        return groups
156
157    # --- solver core ----------------------------------------------------
158    def _initialize_group_constraints(self) -> None:
159        groups = self._build_population_groups(self.focus_indices)
160        covered = np.concatenate(groups)
161        if covered.size != self.population_count or not np.all(np.sort(covered) == np.arange(self.population_count)):
162            raise ValueError("Population groups must cover each population exactly once.")
163        self.groups: List[np.ndarray] = groups
164        self.group_count = len(groups)
165        self.focus_group_index = 0
166        self.solve_groups = [idx for idx in range(self.group_count) if idx != self.focus_group_index]
167        self.dim = len(self.solve_groups)
168        self.focus_population_mask = np.zeros(self.population_count, dtype=bool)
169        self.focus_population_mask[self.focus_indices] = True
170        self.group_membership = np.zeros((self.population_count, self.group_count), dtype=float)
171        for idx, members in enumerate(self.groups):
172            self.group_membership[members, idx] = 1.0
173        self.group_sizes = np.array([len(members) for members in self.groups], dtype=float)
174        self.group_inverse_sizes = np.reciprocal(self.group_sizes)
175        self.focus_vector = np.zeros(self.group_count, dtype=float)
176        self.focus_vector[self.focus_group_index] = 1.0
177        self.selector_matrix = np.zeros((self.group_count, self.dim), dtype=float)
178        for col, group_idx in enumerate(self.solve_groups):
179            self.selector_matrix[group_idx, col] = 1.0
180        self.residual_matrix = np.zeros((self.dim, self.population_count), dtype=float)
181        for row, group_idx in enumerate(self.solve_groups):
182            members = self.groups[group_idx]
183            self.residual_matrix[row, members] = 1.0 / len(members)
184
185    def _resolve_focus_indices(self, focus_config) -> np.ndarray:
186        if focus_config is None:
187            entries = [0]
188        elif isinstance(focus_config, (int, np.integer)):
189            entries = [int(focus_config)]
190        else:
191            entries = []
192            for value in focus_config:
193                if isinstance(value, slice):
194                    start = 0 if value.start is None else value.start
195                    stop = self.population_count if value.stop is None else value.stop
196                    step = 1 if value.step is None else value.step
197                    entries.extend(range(start, stop, step))
198                else:
199                    entries.append(int(value))
200        focus = sorted(set(entries))
201        if not focus:
202            focus = [0]
203        for idx in focus:
204            if idx < 0 or idx >= self.population_count:
205                raise ValueError(f"Focus index {idx} out of bounds for population size {self.population_count}.")
206        return np.asarray(focus, dtype=int)
207
208    def _full_rates_numpy(self, x: np.ndarray) -> np.ndarray:
209        arr = np.asarray(x, dtype=float).ravel()
210        if arr.size == self.population_count:
211            return arr
212        if arr.size == self.group_count:
213            group_values = arr.copy()
214            group_values[self.focus_group_index] = self.v_focus
215        else:
216            if arr.size != self.dim:
217                raise ValueError(
218                    f"Expected vector of length {self.dim}, {self.group_count}, or {self.population_count}, got {arr.size}."
219                )
220            group_values = self.selector_matrix @ arr
221            group_values += self.focus_vector * self.v_focus
222        return self.group_membership @ group_values
223
224    def _phi_numpy(self, full_rates: np.ndarray) -> np.ndarray:
225        mean = self.A.dot(full_rates) + self.bias
226        rates_sq = full_rates * full_rates
227        var = np.maximum(self.B.dot(full_rates) + self.C.dot(rates_sq), self.VAR_EPS)
228        return 0.5 * (1 - special.erf(-mean / np.sqrt(2.0 * var)))
229
230    def _phi_jacobian_numpy(self, full_rates: np.ndarray) -> np.ndarray:
231        mean = self.A.dot(full_rates) + self.bias
232        rates_sq = full_rates * full_rates
233        var = np.maximum(self.B.dot(full_rates) + self.C.dot(rates_sq), self.VAR_EPS)
234        inv_sqrt = 1.0 / np.sqrt(2.0 * var)
235        exp_term = np.exp(-(mean ** 2) / (2.0 * var))
236        coeff = (1.0 / np.sqrt(np.pi)) * exp_term * inv_sqrt
237        correction = mean / (2.0 * var)
238        dvar_dm = self.B + self.C * (2.0 * full_rates[None, :])
239        return coeff[:, None] * (self.A - correction[:, None] * dvar_dm)
240
241    def residual_numpy(self, x: np.ndarray) -> np.ndarray:
242        rates = self._full_rates_numpy(x)
243        phi = self._phi_numpy(rates)
244        residual = (phi - rates) / self.tau
245        if self.dim == 0:
246            return np.zeros((0,), dtype=float)
247        return self.residual_matrix @ residual
248
249    @staticmethod
250    def _jax_residual(x, v_focus, args):
251        A, B, C, bias, tau, focus_vector, selector, membership, reduction, var_eps = args
252        group_values = selector @ x + focus_vector * v_focus
253        full_rates = membership @ group_values
254        mean = A @ full_rates + bias
255        var = jnp.maximum(B @ full_rates + C @ (full_rates * full_rates), var_eps)
256        phi = 0.5 * (1 - jspecial.erf(-mean / jnp.sqrt(2.0 * var)))
257        residual = (phi - full_rates) / tau
258        return reduction @ residual
259
260    def solve(self, initial_guess: Optional[np.ndarray] = None) -> Tuple[np.ndarray, np.ndarray, bool]:
261        """Solve for a fixed point of the reduced mean-field system.
262
263        Returns
264        -------
265        tuple[np.ndarray, np.ndarray, bool]
266            `(x, residual, success)` where `x` is the reduced solution vector.
267
268        Expected output
269        ---------------
270        `success` is `True` when the nonlinear solver converged and `residual`
271        is close to zero.
272        """
273        initial = self._coerce_initial(initial_guess)
274        if self.dim == 0:
275            residual = self.residual_numpy(initial)
276            return initial, residual, True
277        if self.use_jax:
278            try:
279                return self._solve_with_optimistix(initial)
280            except SolverConvergenceError:
281                raise
282            except (ValueError, RuntimeError, TypeError, AttributeError) as e:
283                logger.warning(
284                    "Optimistix solver failed with %s: %s. Falling back to scipy solver.",
285                    type(e).__name__,
286                    str(e),
287                )
288                self.use_jax = False
289        return self._solve_with_scipy(initial)
290
291    @staticmethod
292    def _residual_norm(residual: np.ndarray) -> float:
293        arr = np.asarray(residual, dtype=float).ravel()
294        if arr.size == 0:
295            return 0.0
296        return float(np.linalg.norm(arr))
297
298    def _normalize_solver_result(
299        self,
300        x: np.ndarray,
301        residual: np.ndarray,
302        success: bool,
303        *,
304        method: str,
305    ) -> Tuple[np.ndarray, np.ndarray, bool]:
306        x_arr = np.asarray(x, dtype=float).ravel()
307        residual_arr = np.asarray(residual, dtype=float).ravel()
308        finite = np.isfinite(x_arr).all() and np.isfinite(residual_arr).all()
309        accepted = bool(success) and finite
310        if (not accepted) and finite:
311            residual_norm = self._residual_norm(residual_arr)
312            if residual_norm <= self.RESIDUAL_ACCEPT_TOL:
313                logger.debug(
314                    "Accepting near-root from %s at v_focus %.6f with residual norm %.3e.",
315                    method,
316                    float(self.v_focus),
317                    residual_norm,
318                )
319                accepted = True
320        return x_arr, residual_arr, accepted
321
322    def _solve_with_scipy(self, initial: np.ndarray) -> Tuple[np.ndarray, np.ndarray, bool]:
323        best_x = np.asarray(initial, dtype=float).ravel()
324        best_residual = np.asarray(self.residual_numpy(best_x), dtype=float).ravel()
325        best_norm = self._residual_norm(best_residual)
326
327        def remember(x: np.ndarray, residual: np.ndarray) -> None:
328            nonlocal best_x, best_residual, best_norm
329            if not (np.isfinite(x).all() and np.isfinite(residual).all()):
330                return
331            norm = self._residual_norm(residual)
332            if norm < best_norm:
333                best_x = np.asarray(x, dtype=float).ravel()
334                best_residual = np.asarray(residual, dtype=float).ravel()
335                best_norm = norm
336
337        hybr = optimize.root(
338            self.residual_numpy,
339            initial,
340            method="hybr",
341            tol=self.root_tol,
342            options={"maxfev": self.max_function_evals},
343        )
344        x_hybr, residual_hybr, success_hybr = self._normalize_solver_result(
345            hybr.x,
346            hybr.fun,
347            bool(hybr.success),
348            method="scipy.root(hybr)",
349        )
350        if success_hybr:
351            return x_hybr, residual_hybr, True
352        remember(x_hybr, residual_hybr)
353
354        lm_initial = x_hybr if np.isfinite(x_hybr).all() else np.asarray(initial, dtype=float).ravel()
355        lm = optimize.root(
356            self.residual_numpy,
357            lm_initial,
358            method="lm",
359            tol=self.root_tol,
360            options={"maxiter": self.max_function_evals},
361        )
362        x_lm, residual_lm, success_lm = self._normalize_solver_result(
363            lm.x,
364            lm.fun,
365            bool(lm.success),
366            method="scipy.root(lm)",
367        )
368        if success_lm:
369            return x_lm, residual_lm, True
370        remember(x_lm, residual_lm)
371
372        ls_initial = x_lm if np.isfinite(x_lm).all() else lm_initial
373        ls = optimize.least_squares(
374            self.residual_numpy,
375            ls_initial,
376            method="trf",
377            xtol=self.root_tol,
378            ftol=self.root_tol,
379            gtol=self.root_tol,
380            max_nfev=self.max_function_evals,
381        )
382        x_ls, residual_ls, success_ls = self._normalize_solver_result(
383            ls.x,
384            ls.fun,
385            bool(ls.success),
386            method="scipy.least_squares(trf)",
387        )
388        if success_ls:
389            return x_ls, residual_ls, True
390        remember(x_ls, residual_ls)
391
392        return best_x, best_residual, False
393
394    def _solve_with_optimistix(self, initial: np.ndarray) -> Tuple[np.ndarray, np.ndarray, bool]:
395        args = self._prepare_jax_args()
396        solver_entry = self._get_jax_solver()
397        if solver_entry is None:
398            raise RuntimeError("Optimistix solver is unavailable.")
399        x0 = jnp.asarray(initial, dtype=jnp.float64)
400        v_focus = jnp.asarray(float(self.v_focus), dtype=jnp.float64)
401        value, status, _ = solver_entry["single"](x0, v_focus, args)
402        success = bool(np.asarray(status, dtype=bool))
403        if not success:
404            raise SolverConvergenceError(self.v_focus, self.max_steps)
405        value_np = np.asarray(value, dtype=float)
406        return value_np, self.residual_numpy(value_np), success
407
408    def _prepare_jax_args(self):
409        if self._jax_args is None:
410            self._jax_args = (
411                jnp.asarray(self.A, dtype=jnp.float64),
412                jnp.asarray(self.B, dtype=jnp.float64),
413                jnp.asarray(self.C, dtype=jnp.float64),
414                jnp.asarray(self.bias, dtype=jnp.float64),
415                jnp.asarray(self.tau, dtype=jnp.float64),
416                jnp.asarray(self.focus_vector, dtype=jnp.float64),
417                jnp.asarray(self.selector_matrix, dtype=jnp.float64),
418                jnp.asarray(self.group_membership, dtype=jnp.float64),
419                jnp.asarray(self.residual_matrix, dtype=jnp.float64),
420                float(self.VAR_EPS),
421            )
422        return self._jax_args
423
424    def _get_jax_solver(self):
425        if not (self.use_jax and HAS_JAX and optx is not None):
426            return None
427        if self.dim == 0:
428            return None
429        key = (self.dim, float(self.root_tol), int(self.max_steps))
430        entry = JAX_SOLVER_CACHE.get(key)
431        if entry is None:
432            entry = _build_jax_solver_entry(self.dim, float(self.root_tol), int(self.max_steps))
433            JAX_SOLVER_CACHE[key] = entry
434        return entry
435
436    def phi_numpy(self, x: np.ndarray) -> np.ndarray:
437        """Evaluate the transfer function on the full-rate state implied by `x`.
438
439        Expected output
440        ---------------
441        Returns one activity value per population in the interval `[0, 1]`.
442        """
443        rates = self._full_rates_numpy(np.asarray(x, dtype=float))
444        return self._phi_numpy(rates)
445
446    def full_rates_numpy(self, x: np.ndarray) -> np.ndarray:
447        """Expand a reduced solver vector into one rate per population.
448
449        Expected output
450        ---------------
451        The returned array has length `population_count`.
452        """
453        return self._full_rates_numpy(np.asarray(x, dtype=float))
454
455    def jacobian_numpy(self, x: np.ndarray) -> np.ndarray:
456        """Return the Jacobian of the full mean-field residual at `x`.
457
458        Expected output
459        ---------------
460        The returned matrix has shape `(population_count, population_count)`.
461        """
462        rates = self._full_rates_numpy(np.asarray(x, dtype=float))
463        return (self._phi_jacobian_numpy(rates) - np.eye(self.population_count)) / self.tau[:, np.newaxis]
464
465    def focus_output(self, rates: np.ndarray) -> float:
466        """Average the rates over the configured focus populations."""
467        values = np.asarray(rates, dtype=float)[self.focus_indices]
468        return float(values.mean()) if values.size else float(rates[0])
469
470    def _reduce_full_rates(self, full_rates: np.ndarray) -> np.ndarray:
471        if self.dim == 0:
472            return np.zeros((0,), dtype=float)
473        full = np.asarray(full_rates, dtype=float).reshape((self.population_count,))
474        group_sums = self.group_membership.T @ full
475        group_means = group_sums * self.group_inverse_sizes
476        return group_means[self.solve_groups]
477
478    def _coerce_initial(self, initial_guess: Optional[np.ndarray]) -> np.ndarray:
479        if self.dim == 0:
480            return np.zeros((0,), dtype=float)
481        if initial_guess is None:
482            return np.full((self.dim,), 0.1, dtype=float)
483        arr = np.asarray(initial_guess, dtype=float).ravel()
484        if arr.size == self.dim:
485            return arr
486        if arr.size == self.group_count:
487            copy = arr.copy()
488            copy[self.focus_group_index] = self.v_focus
489            return copy[self.solve_groups]
490        if arr.size == self.population_count:
491            return self._reduce_full_rates(arr)
492        non_focus_count = self.population_count - len(self.focus_indices)
493        if arr.size == non_focus_count:
494            full_rates = np.empty(self.population_count, dtype=float)
495            full_rates[self.focus_population_mask] = self.v_focus
496            full_rates[~self.focus_population_mask] = arr
497            return self._reduce_full_rates(full_rates)
498        raise ValueError(
499            f"Initial guess must have length {self.dim}, {self.group_count}, {self.population_count}, or {non_focus_count}, got {arr.size}."
500        )
501
502    def solve_sequence(
503        self,
504        v_focus_values: Sequence[float],
505        initial_guess: Optional[np.ndarray] = None,
506    ) -> Optional[Tuple[np.ndarray, np.ndarray]]:
507        """Solve a sequence of focus inputs with the accelerated JAX path.
508
509        Expected output
510        ---------------
511        Returns `(solutions, success_flags)` when the JAX solver is available,
512        otherwise `None`.
513        """
514        values = np.asarray(list(v_focus_values), dtype=float)
515        if values.size == 0:
516            return np.zeros((0, self.dim)), np.zeros((0,), dtype=bool)
517        if self.dim == 0:
518            zeros = np.zeros((values.size, 0), dtype=float)
519            return zeros, np.ones((values.size,), dtype=bool)
520        if not (self.use_jax and HAS_JAX and optx is not None):
521            return None
522        solver_entry = self._get_jax_solver()
523        if solver_entry is None:
524            return None
525        args = self._prepare_jax_args()
526        initial = self._coerce_initial(initial_guess)
527        x0 = jnp.asarray(initial, dtype=jnp.float64)
528        v_seq = jnp.asarray(values, dtype=jnp.float64)
529        try:
530            solutions, statuses, _ = solver_entry["scan"](x0, v_seq, args)
531        except (ValueError, RuntimeError, TypeError, AttributeError) as exc:
532            logger.warning(
533                "Optimistix sweep failed with %s: %s. Falling back to sequential solver.",
534                type(exc).__name__,
535                str(exc),
536            )
537            self.use_jax = False
538            return None
539        solution_np = np.asarray(solutions, dtype=float)
540        success = np.asarray(statuses, dtype=bool)
541        return solution_np, success
542
543    # --- class helpers --------------------------------------------------
544    @classmethod
545    def generate_erf_curve(
546        cls,
547        parameter: Dict,
548        *,
549        start: float = 0.02,
550        end: float = 1.0,
551        step_number: int = 20,
552        retry_step: Optional[float] = None,
553        initial_guess: Optional[np.ndarray] = None,
554        fallback_initials: Optional[Sequence[float]] = None,
555        **network_kwargs,
556    ) -> ERFResult:
557        """Compute the ERF for ``parameter``.
558
559        When ``retry_step`` is provided, the solver retries inputs separated by
560        that value whenever convergence fails. If the ERF cannot be completed
561        the ``completed`` flag is ``False`` and no serialization should happen.
562        Parameters
563        ----------
564        start : float, optional
565            Lower bound of the input range for the ERF (default: 0.02).
566        end : float, optional
567            Upper bound of the input range for the ERF (default: 1.0).
568        step_number : int, optional
569            Number of steps between ``start`` and ``end`` (default: 20).
570        ...
571
572        Expected output
573        ---------------
574        Returns an :class:`ERFResult` whose `x_data` stores the driven focus
575        inputs and whose `y_data` stores the corresponding mean-field outputs.
576        """
577        ERF_EPS = 1e-12  # Tolerance for floating-point comparison of v_in against end
578        x_data: List[float] = []
579        y_data: List[float] = []
580        solves: List[np.ndarray] = []
581        step = (end - start) / max(step_number, 1)
582        adaptive_min_step = min(abs(step), cls.ADAPTIVE_CONTINUATION_MIN_STEP) if step != 0 else cls.ADAPTIVE_CONTINUATION_MIN_STEP
583        aborted = False
584        fallback_values = list(fallback_initials) if fallback_initials is not None else [0.02, 0.2, 0.5, 0.8, 0.98]
585
586        def solve_with_fallback_initials(
587            system: "RateSystem",
588            initial: Optional[np.ndarray],
589        ) -> Tuple[Optional[np.ndarray], Optional[np.ndarray], bool]:
590            best_solution: Optional[np.ndarray] = None
591            best_residual: Optional[np.ndarray] = None
592            best_norm = float("inf")
593            candidates: List[Optional[np.ndarray]] = [initial]
594            for seed in fallback_values:
595                if system.dim == 0:
596                    candidate = np.zeros((0,), dtype=float)
597                else:
598                    candidate = np.full((system.dim,), float(seed), dtype=float)
599                candidates.append(candidate)
600            for candidate in candidates:
601                try:
602                    solution, residual, success = system.solve(candidate)
603                except SolverConvergenceError as exc:
604                    logger.debug("Skipping v_in %.6f: %s", float(system.v_focus), str(exc))
605                    continue
606                residual_arr = np.asarray(residual, dtype=float).ravel()
607                if success:
608                    return np.asarray(solution, dtype=float).ravel(), residual_arr, True
609                if np.isfinite(residual_arr).all():
610                    residual_norm = cls._residual_norm(residual_arr)
611                    if residual_norm < best_norm:
612                        best_solution = np.asarray(solution, dtype=float).ravel()
613                        best_residual = residual_arr
614                        best_norm = residual_norm
615            return best_solution, best_residual, False
616
617        vector_values: Optional[List[float]] = None
618        if retry_step is None:
619            vector_values = []
620            if step == 0:
621                vector_values.append(float(start))
622            else:
623                cursor = float(start)
624                while cursor <= end + ERF_EPS:
625                    vector_values.append(cursor)
626                    cursor += step
627        initial_focus = float(start)
628        if vector_values:
629            initial_focus = float(vector_values[0])
630        system = cls(parameter, initial_focus, **network_kwargs)
631        current_initial = initial_guess
632        prefetched_solutions: Optional[np.ndarray] = None
633        prefix_limit = 0
634        if vector_values and system.use_jax:
635            seq_result = system.solve_sequence(vector_values, initial_guess=current_initial)
636            if seq_result is not None:
637                prefetched_solutions, success_flags = seq_result
638                failure = np.flatnonzero(~success_flags)
639                prefix_limit = int(failure[0]) if failure.size else len(vector_values)
640        if vector_values is not None:
641            idx = 0
642            while idx < len(vector_values):
643                v_in = vector_values[idx]
644                system.v_focus = float(v_in)
645                use_prefetched = prefetched_solutions is not None and idx < prefix_limit
646                if use_prefetched:
647                    solution = np.asarray(prefetched_solutions[idx], dtype=float)
648                    success = True
649                else:
650                    solution, residual, success = solve_with_fallback_initials(system, current_initial)
651                    if not success:
652                        prev_v = x_data[-1] if x_data else None
653                        if prev_v is not None:
654                            gap = float(v_in) - float(prev_v)
655                            if gap > adaptive_min_step + ERF_EPS:
656                                midpoint = float(prev_v) + 0.5 * gap
657                                vector_values.insert(idx, midpoint)
658                                prefetched_solutions = None
659                                prefix_limit = 0
660                                continue
661                        aborted = True
662                        break
663                phi_values = system.phi_numpy(solution)
664                x_data.append(system.v_focus)
665                y_data.append(system.focus_output(phi_values))
666                solves.append(solution)
667                current_initial = solution
668                idx += 1
669            completed = (not aborted) and (len(x_data) == len(vector_values))
670        else:
671            v_in = float(start)
672            next_value = v_in
673            while v_in <= end + ERF_EPS:
674                system.v_focus = float(v_in)
675                solution, residual, success = solve_with_fallback_initials(system, current_initial)
676                if not success:
677                    if retry_step is not None:
678                        v_in += retry_step
679                        next_value = v_in
680                        if solution is not None:
681                            current_initial = solution
682                        continue
683                    aborted = True
684                    break
685                phi_values = system.phi_numpy(solution)
686                x_data.append(system.v_focus)
687                y_data.append(system.focus_output(phi_values))
688                solves.append(solution)
689                current_initial = solution
690                if step == 0:
691                    next_value = float("inf")
692                    break
693                v_in = system.v_focus + step
694                next_value = v_in
695            completed = (not aborted) and (step == 0 or next_value > end + ERF_EPS)
696        return ERFResult(x_data=x_data, y_data=y_data, solves=solves, completed=completed)
697
698    @classmethod
699    def compute_fixpoints(
700        cls,
701        sweep_entry: Sequence,
702        *,
703        tol: float = 1e-3,
704        interpolation_steps: int = 10_000,
705        **network_kwargs,
706    ) -> Dict[float, Dict[str, Any]]:
707        """
708        Compute fixed points from an ERF sweep.
709        Fixed points with slope larger than 1 at the intersection with the identity line are considered unstable in the 1D map approximation.
710
711        Parameters
712        ----------
713        sweep_entry : sequence
714            Tuple (x_data, y_data, solves, parameter) as returned by generate_erf_curve.
715        tol : float, optional
716            Tolerance for detecting crossings of the identity line (x = y) and for
717            merging nearby crossings. Crossings where |x - y| <= tol are treated as
718            fixed points, and crossings within tol of each other are merged. Default
719            is 1e-3.
720        interpolation_steps : int, optional
721            Number of interpolation points used to refine the ERF before searching
722            for crossings. Larger values increase accuracy but also cost. Default
723            is 10_000.
724
725        Expected output
726        ---------------
727        Returns a dictionary keyed by fixed-point location. Each value contains
728        at least `stability`, `rates`, `residual_norm`, and `solver_success`.
729        """
730        SLOPE_STABILITY_THRESHOLD = 1.0 # |d(ERF)/dv| < 1 ⇒ stable in 1D; > 1 ⇒ unstable
731        x_data, y_data, solves, parameter = sweep_entry
732        x_interp, y_interp = interpolate_curve(x_data, y_data, steps=interpolation_steps)
733        if x_interp.size == 0 or y_interp.size == 0:
734            print("Skipping fixpoint analysis: empty ERF data.")
735            return {}
736        diff = x_interp - y_interp
737        crossings = []
738        prev_diff = diff[0]
739        for idx in range(1, len(diff)):
740            curr_diff = diff[idx]
741            cross_val = None
742            if np.abs(curr_diff) <= tol:
743                cross_val = y_interp[idx]
744            elif np.abs(prev_diff) <= tol:
745                cross_val = y_interp[idx - 1]
746            elif prev_diff * curr_diff < 0:
747                weight = prev_diff / (prev_diff - curr_diff)
748                cross_val = y_interp[idx - 1] + weight * (y_interp[idx] - y_interp[idx - 1])
749            if cross_val is not None:
750                if crossings and np.abs(cross_val - crossings[-1][0]) <= tol:
751                    crossings[-1] = (float(cross_val), idx)
752                else:
753                    crossings.append((float(cross_val), idx))
754            prev_diff = curr_diff
755        fixpoints: Dict[float, Dict[str, Any]] = {}
756        if not crossings:
757            return fixpoints
758        solves_array = [np.asarray(s, dtype=float) for s in solves]
759        v_out_old = np.asarray(y_data, dtype=float)
760        for cross_point, idx in crossings:
761            print(f"Cross-Point: {cross_point}")
762            if idx <= 0:
763                slope = np.inf
764            else:
765                slope = (y_interp[idx] - y_interp[idx - 1]) / (x_interp[idx] - x_interp[idx - 1])
766            entry: Dict[str, Any] = {
767                "stability": "unstable",
768                "rates": None,
769                "residual_norm": float("inf"),
770                "solver_success": False,
771                "slope": float(slope) if np.isfinite(slope) else float("inf"),
772                "included": False,
773            }
774            slope_unstable = not np.isfinite(slope) or slope > SLOPE_STABILITY_THRESHOLD
775            if len(solves_array) == 0:
776                entry["reason"] = "missing_erf_solution"
777                fixpoints[cross_point] = entry
778                continue
779            closest_idx = int(np.argmin(np.abs(v_out_old - cross_point)))
780            closest_idx = min(max(closest_idx, 0), len(solves_array) - 1)
781            initial = solves_array[closest_idx]
782            system = cls(parameter, cross_point, **network_kwargs)
783            try:
784                solve, residual, success = system.solve(initial)
785            except SolverConvergenceError as exc:
786                entry["solver_success"] = False
787                entry["residual_norm"] = float("inf")
788                entry["rates"] = None
789                entry["stability"] = "unstable"
790                entry["reason"] = "solver_failed"
791                entry["error"] = str(exc)
792                fixpoints[cross_point] = entry
793                continue
794            residual = np.asarray(residual, dtype=float)
795            residual_norm = float(np.linalg.norm(residual)) if residual.size else 0.0
796            if not np.isfinite(residual).all():
797                residual_norm = float("inf")
798            entry["solver_success"] = bool(success)
799            entry["residual_norm"] = residual_norm
800            if success and np.isfinite(residual).all():
801                entry["rates"] = system.full_rates_numpy(solve)
802            else:
803                entry["rates"] = None
804            if slope_unstable or not success or not np.isfinite(residual).all():
805                if not success or not np.isfinite(residual).all():
806                    print("Warning: convergence problems near cross-point")
807                stability = "unstable"
808            else:
809                jacobian = system.jacobian_numpy(solve)
810                if not np.isfinite(jacobian).all():
811                    stability = "unstable"
812                else:
813                    try:
814                        eigval = np.linalg.eigvals(jacobian)
815                    except np.linalg.LinAlgError:
816                        stability = "unstable"
817                    else:
818                        stability = "stable" if (eigval < 0).all() else "unstable"
819            entry["stability"] = stability
820            if "reason" not in entry:
821                entry["reason"] = "slope" if slope_unstable else None
822            fixpoints[cross_point] = entry
823        return fixpoints