fix(analysis): fixed-point classification, nullcline selection, GD batches, plot kwargs (#849)

fix(analysis): fixed-point classification, nullcline selection, GD batch count, plot kwargs

- 2D star vs degenerate-node stability classification was inverted and the STAR
  branch was unreachable dead code; correctly distinguish star (scalar*I) from
  defective node (Critical)
- phase-plane select_candidates='nullclines' tested the fy condition twice,
  silently dropping all fx-nullcline candidate points; test fx OR fy (High)
- slow_points num_opt_loops = int(num_opt/num_batch) was 0 when num_opt<num_batch
  -> empty losses -> concatenate crash; use max(1, ceil(...)) (Medium)
- streamplot used plot_style.get('linewidth') (not pop), passing linewidth twice
  -> TypeError; use pop (Medium)
- set_markersize wrote to a typo local, leaving the module global stale (Medium)

Findings recorded in docs/issues-found-20260619-analysis.md

Author: chaoming0625

brainpy/analysis/highdim/slow_points.py

@@ -376,7 +376,10 @@ def batch_train(start_i, n_batch):
             print(f"Optimizing with {optimizer} to find fixed points:")
         opt_losses = []
         do_stop = False
-        num_opt_loops = int(num_opt / num_batch)
+        # Always run at least one batch, even when ``num_opt < num_batch`` (where
+        # ``int(num_opt / num_batch)`` would otherwise be 0 and leave ``opt_losses``
+        # empty, crashing the subsequent ``jnp.concatenate``).
+        num_opt_loops = max(1, int(np.ceil(num_opt / num_batch)))
         for oidx in range(num_opt_loops):
             if do_stop:
                 break

brainpy/analysis/highdim/slow_points_coverage_test.py

@@ -155,6 +155,19 @@ def step(x, scale):
     assert f.opt_losses.ndim == 1
 
 
+def test_gd_small_num_opt_runs_at_least_one_batch():
+    """Regression for P13-M1: ``num_opt < num_batch`` previously gave
+    ``num_opt_loops = int(num_opt / num_batch) == 0``, leaving ``opt_losses``
+    empty and crashing on ``jnp.concatenate``.  At least one batch must run."""
+    f = bp.analysis.SlowPointFinder(f_cell=_linear_step, f_type=CONTINUOUS, verbose=False)
+    rng = bm.random.RandomState(2)
+    # num_opt (50) < num_batch (100) -> would crash before the fix.
+    f.find_fps_with_gd_method(rng.random((4, 2)) * 0.2, num_opt=50, num_batch=100)
+    assert f.num_fps == 4
+    assert f.opt_losses.ndim == 1
+    assert f.opt_losses.shape[0] == 100
+
+
 # --------------------------------------------------------------------------- #
 # full pipeline for a callable system
 # --------------------------------------------------------------------------- #

brainpy/analysis/lowdim/lowdim_phase_plane.py

@@ -226,7 +226,10 @@ def plot_vector_field(self, with_plot=True, with_return=False,
             elif plot_method == 'streamplot':
                 if plot_style is None:
                     plot_style = dict(arrowsize=1.2, density=1, color='thistle')
-                linewidth = plot_style.get('linewidth', None)
+                # Pop (not get) so a user-supplied ``linewidth`` is consumed and
+                # only passed once -- otherwise it is forwarded both explicitly
+                # and via ``**plot_style`` -> TypeError (multiple values).
+                linewidth = plot_style.pop('linewidth', None)
                 if linewidth is None:
                     if (not np.isnan(dx).any()) and (not np.isnan(dy).any()):
                         min_width, max_width = 0.5, 5.5
@@ -326,7 +329,7 @@ def plot_fixed_point(self, with_plot=True, with_return=False, show=False,
                 candidates = jnp.vstack(candidates)
             elif select_candidates == 'nullclines':
                 candidates = [self.analyzed_results[key][0] for key in self.analyzed_results.keys()
-                              if key.startswith(C.fy_nullcline_points) or key.startswith(C.fy_nullcline_points)]
+                              if key.startswith(C.fx_nullcline_points) or key.startswith(C.fy_nullcline_points)]
                 if len(candidates) == 0:
                     raise errors.AnalyzerError(f'No nullcline points are found, please call '
                                                f'".{self.plot_nullcline.__name__}()" first.')

brainpy/analysis/lowdim/lowdim_phase_plane_coverage_test.py

@@ -111,17 +111,18 @@ def test_pp2d_quiver_and_return():
     assert dx.shape == dy.shape
 
 
-# NOTE (defect): passing ``linewidth`` inside ``plot_style`` to the
-# ``streamplot`` branch crashes.  The code reads ``plot_style.get('linewidth')``
-# but never pops it, so ``pyplot.streamplot(..., linewidth=linewidth,
-# **plot_style)`` passes ``linewidth`` twice -> ``TypeError: got multiple values
-# for keyword argument 'linewidth'`` (lowdim_phase_plane.py:235).
-def test_pp2d_streamplot_custom_linewidth_is_buggy():
+# Regression for P13-M2: passing ``linewidth`` inside ``plot_style`` to the
+# ``streamplot`` branch used to crash because the code read
+# ``plot_style.get('linewidth')`` without popping it, so ``linewidth`` was
+# forwarded twice (explicitly and via ``**plot_style``) ->
+# ``TypeError: got multiple values for keyword argument 'linewidth'``.
+# The fix pops the key, so a user-supplied ``linewidth`` is honoured exactly once.
+def test_pp2d_streamplot_custom_linewidth():
     pp = bp.analysis.PhasePlane2D(model=_fhn_2d(),
                                   target_vars={'V': [-3., 3.], 'w': [-1., 3.]},
                                   resolutions=0.2)
-    with pytest.raises(TypeError):
-        pp.plot_vector_field(plot_method='streamplot', plot_style=dict(linewidth=1.0))
+    # should no longer raise
+    pp.plot_vector_field(plot_method='streamplot', plot_style=dict(linewidth=1.0))
 
 
 def test_pp2d_unknown_plot_method_raises():
@@ -170,6 +171,37 @@ def test_pp2d_fixed_point_after_nullcline():
     assert fps is not None
 
 
+def test_pp2d_nullclines_selector_uses_both_nullclines():
+    """Regression for P13-H2: the ``select_candidates='nullclines'`` branch must
+    gather candidates from *both* the fx- and fy-nullcline point sets.  The old
+    code tested ``startswith(fy_...)`` twice, so fx-nullcline points were dropped.
+    """
+    import jax.numpy as jnp
+    from brainpy.analysis import constants as C
+
+    pp = bp.analysis.PhasePlane2D(model=_fhn_2d(),
+                                  target_vars={'V': [-3., 3.], 'w': [-1., 3.]},
+                                  resolutions=0.1)
+    pp.plot_nullcline()
+
+    fx_keys = [k for k in pp.analyzed_results if k.startswith(C.fx_nullcline_points)]
+    fy_keys = [k for k in pp.analyzed_results if k.startswith(C.fy_nullcline_points)]
+    assert len(fx_keys) > 0 and len(fy_keys) > 0
+    n_fx = sum(pp.analyzed_results[k][0].shape[0] for k in fx_keys)
+    n_fy = sum(pp.analyzed_results[k][0].shape[0] for k in fy_keys)
+
+    # Reproduce the candidate gathering used by ``select_candidates='nullclines'``.
+    candidates = [pp.analyzed_results[k][0] for k in pp.analyzed_results.keys()
+                  if k.startswith(C.fx_nullcline_points) or k.startswith(C.fy_nullcline_points)]
+    candidates = jnp.vstack(candidates)
+    # union of both nullcline candidate sets, not just one of them
+    assert candidates.shape[0] == n_fx + n_fy
+    assert candidates.shape[0] > n_fy
+
+    fps = pp.plot_fixed_point(select_candidates='nullclines', with_return=True)
+    assert fps is not None
+
+
 # --------------------------------------------------------------------------- #
 # PhasePlane2D trajectory + limit cycle
 # --------------------------------------------------------------------------- #

brainpy/analysis/plotstyle.py

@@ -76,6 +76,6 @@ def set_markersize(markersize):
     if not isinstance(markersize, int):
         raise TypeError(f"Must be an integer, but got {type(markersize)}: {markersize}")
     global _markersize
-    __markersize = markersize
+    _markersize = markersize
     for key in tuple(plot_schema.keys()):
         plot_schema[key]['markersize'] = markersize

brainpy/analysis/plotstyle_test.py

@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+# Copyright 2025 BrainX Ecosystem Limited. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+# ==============================================================================
+import pytest
+
+from brainpy.analysis import plotstyle
+
+
+def test_set_markersize_updates_global():
+    """Regression for P13-H1: ``set_markersize`` must update the module global
+    ``_markersize`` (it previously assigned to a typo local ``__markersize``)."""
+    original = plotstyle._markersize
+    try:
+        plotstyle.set_markersize(33)
+        # per-key schema entries updated
+        assert plotstyle.plot_schema[plotstyle.SADDLE_NODE]['markersize'] == 33
+        # module global updated, not just the per-key dicts
+        assert plotstyle._markersize == 33
+    finally:
+        plotstyle.set_markersize(original)
+
+
+def test_set_markersize_type_check():
+    with pytest.raises(TypeError):
+        plotstyle.set_markersize(1.5)
+
+
+def test_set_plot_schema_validations():
+    with pytest.raises(TypeError):
+        plotstyle.set_plot_schema(123)
+    with pytest.raises(KeyError):
+        plotstyle.set_plot_schema('not-a-real-fixed-point-type')
+    # valid update
+    plotstyle.set_plot_schema(plotstyle.SADDLE_NODE, color='black')
+    assert plotstyle.plot_schema[plotstyle.SADDLE_NODE]['color'] == 'black'

brainpy/analysis/stability.py

@@ -145,22 +145,25 @@ def stability_analysis(derivatives):
                 elif e > 0:
                     return UNSTABLE_NODE_2D
                 else:
-                    w = np.linalg.eigvals(derivatives)
-                    if w[0] == w[1]:
-                        return UNSTABLE_DEGENERATE_2D
-                    else:
+                    # Repeated eigenvalue. A star (proper) node has a full
+                    # eigenspace, which happens iff the matrix is a scalar
+                    # multiple of the identity (b == c == 0 and a == d).
+                    # Otherwise the matrix is defective -> degenerate (improper) node.
+                    if b == 0 and c == 0 and a == d:
                         return UNSTABLE_STAR_2D
+                    else:
+                        return UNSTABLE_DEGENERATE_2D
             else:
                 if e < 0:
                     return STABLE_FOCUS_2D
                 elif e > 0:
                     return STABLE_NODE_2D
                 else:
-                    w = np.linalg.eigvals(derivatives)
-                    if w[0] == w[1]:
-                        return STABLE_DEGENERATE_2D
-                    else:
+                    # Repeated eigenvalue. See the unstable branch above.
+                    if b == 0 and c == 0 and a == d:
                         return STABLE_STAR_2D
+                    else:
+                        return STABLE_DEGENERATE_2D
 
     elif np.size(derivatives) == 9:  # 3D dynamical system
         eigenvalues = np.linalg.eigvals(np.array(derivatives))

brainpy/analysis/stability_coverage_test.py

@@ -88,14 +88,13 @@ def test_2d_unstable_focus_and_node():
 
 
 def test_2d_unstable_degenerate_and_star():
-    # p > 0, e == 0. Distinct construction so eigenvalues equal -> degenerate.
+    # p > 0, e == 0. A scalar multiple of the identity is a proper (star) node:
+    # full 2-D eigenspace.  (P13-C1: previously mislabelled as degenerate.)
     J = [[2., 0.], [0., 2.]]  # trace 4, det 4, e = 16 - 16 = 0; eigvals both 2
-    assert stability_analysis(J) == st.UNSTABLE_DEGENERATE_2D
-    # p > 0, e == 0 but eigvals not literally equal (numerical) -> star branch.
-    # Use a Jordan-like block so np.linalg.eigvals returns slightly different.
+    assert stability_analysis(J) == st.UNSTABLE_STAR_2D
+    # p > 0, e == 0 with a defective (Jordan) block -> degenerate (improper) node.
     J = [[2., 1.], [0., 2.]]  # trace 4, det 4, e = 0; defective matrix
-    res = stability_analysis(J)
-    assert res in (st.UNSTABLE_DEGENERATE_2D, st.UNSTABLE_STAR_2D)
+    assert stability_analysis(J) == st.UNSTABLE_DEGENERATE_2D
 
 
 def test_2d_stable_focus_and_node():
@@ -108,13 +107,13 @@ def test_2d_stable_focus_and_node():
 
 
 def test_2d_stable_degenerate_and_star():
-    # p < 0, e == 0 -> stable degenerate (eigvals equal)
+    # p < 0, e == 0. Scalar multiple of identity -> proper (star) node.
+    # (P13-C1: previously mislabelled as degenerate.)
     J = [[-2., 0.], [0., -2.]]  # trace -4, det 4, e = 0
-    assert stability_analysis(J) == st.STABLE_DEGENERATE_2D
-    # defective -> degenerate or star
+    assert stability_analysis(J) == st.STABLE_STAR_2D
+    # defective (Jordan) block -> degenerate (improper) node.
     J = [[-2., 1.], [0., -2.]]
-    res = stability_analysis(J)
-    assert res in (st.STABLE_DEGENERATE_2D, st.STABLE_STAR_2D)
+    assert stability_analysis(J) == st.STABLE_DEGENERATE_2D
 
 
 # --------------------------------------------------------------------------- #

brainpy/analysis/stability_test.py

@@ -13,10 +13,44 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 # ==============================================================================
+import numpy as np
+
 from brainpy.analysis.stability import *
 
 
 def test_d1():
     assert stability_analysis(1.) == UNSTABLE_POINT_1D
     assert stability_analysis(-1.) == STABLE_POINT_1D
     assert stability_analysis(0.) == SADDLE_NODE
+
+
+# --------------------------------------------------------------------------- #
+# 2D star (proper) vs degenerate (improper) node — regression for P13-C1.
+#
+# A *star* node has a full 2-D eigenspace, i.e. the Jacobian is a scalar
+# multiple of the identity (b == c == 0 and a == d). A *degenerate*
+# (improper) node has a repeated eigenvalue but is defective (single
+# eigenvector), e.g. a Jordan block.
+# --------------------------------------------------------------------------- #
+def test_2d_stable_star_proper_node():
+    # -I : repeated eigenvalue -1, full eigenspace -> stable STAR.
+    J = np.array([[-1., 0.], [0., -1.]])
+    assert stability_analysis(J) == STABLE_STAR_2D
+
+
+def test_2d_stable_degenerate_defective():
+    # Jordan block with repeated eigenvalue -1, defective -> stable DEGENERATE.
+    J = np.array([[-1., 1.], [0., -1.]])
+    assert stability_analysis(J) == STABLE_DEGENERATE_2D
+
+
+def test_2d_unstable_star_proper_node():
+    # 2*I : repeated eigenvalue +2, full eigenspace -> unstable STAR.
+    J = np.array([[2., 0.], [0., 2.]])
+    assert stability_analysis(J) == UNSTABLE_STAR_2D
+
+
+def test_2d_unstable_degenerate_defective():
+    # Jordan block with repeated eigenvalue +2, defective -> unstable DEGENERATE.
+    J = np.array([[2., 1.], [0., 2.]])
+    assert stability_analysis(J) == UNSTABLE_DEGENERATE_2D