PERF pre-allocate arrays in NARX _get_hc_ids/_get_jc_ids (#242)

* PERF replace quadratic np.vstack/np.append with pre-allocated arrays in NARX ID construction

_get_hc_ids and _get_jc_ids grew arrays with np.vstack inside loops,
causing O(n²) allocation. Two-pass approach: count rows first,
pre-allocate exact size, fill in-place. 1.5x faster, 50% less RAM.

* vectorize and comments

* fix np astype

* PERF avoid result copies in NARX _get_jc_ids/_get_hc_ids

Use astype(copy=False) and free the boolean mask before allocating
result arrays, so peak RAM no longer exceeds the previous incremental
implementation at meaningful sizes (-15% to -27% peak for large term
counts). Outputs are byte-identical; jac/hess tests pass.

* CI re-trigger (OpenML 504 timeouts)

* readability

* scikit-learn<1.9

* pin pixi 0.70.0

* fix readthedocs

---------

Co-authored-by: sikai zhang <matthew.szhang91@gmail.com>

Author: sebastianbreguel

.github/workflows/asv.yml

@@ -36,6 +36,7 @@ jobs:
 
       - uses: prefix-dev/setup-pixi@v0.9.6
         with:
+          pixi-version: v0.70.0
           environments: dev
           cache: true
 
@@ -102,6 +103,7 @@ jobs:
 
       - uses: prefix-dev/setup-pixi@v0.9.6
         with:
+          pixi-version: v0.70.0
           environments: dev
           cache: true
 

.github/workflows/static.yml

@@ -11,6 +11,7 @@ jobs:
       - uses: actions/checkout@v6
       - uses: prefix-dev/setup-pixi@v0.9.6
         with:
+          pixi-version: v0.70.0
           environments: static
           cache: true
 

.github/workflows/test.yml

@@ -25,6 +25,7 @@ jobs:
       - uses: actions/checkout@v6
       - uses: prefix-dev/setup-pixi@v0.9.6
         with:
+          pixi-version: v0.70.0
           environments: >-
             dev
             docs

.github/workflows/wheel.yml

@@ -10,6 +10,7 @@ jobs:
       - uses: actions/checkout@v6
       - uses: prefix-dev/setup-pixi@v0.9.6
         with:
+          pixi-version: v0.70.0
           environments: dev
           cache: true
       - name: Re-install local
@@ -52,8 +53,6 @@ jobs:
           CIBW_ARCHS_LINUX: auto
           CIBW_ARCHS_MACOS: x86_64 arm64
           CIBW_ARCHS_WINDOWS: auto64
-          # Include free-threaded support
-          CIBW_ENABLE: cpython-freethreading
           CIBW_BUILD_FRONTEND: "build[uv]"
       - name: Upload package
         uses: actions/upload-artifact@v7

.readthedocs.yml

@@ -6,8 +6,8 @@ build:
    jobs:
       create_environment:
          - asdf plugin add pixi
-         - asdf install pixi latest
-         - asdf global pixi latest
+         - asdf install pixi 0.70.0
+         - asdf global pixi 0.70.0
       install:
          - pixi install -e docs
       build:

README.rst

@@ -108,12 +108,6 @@ fastcan can be used for system identification.
 In particular, we provide a submodule `fastcan.narx` to build Nonlinear AutoRegressive eXogenous (NARX) models.
 For more information, check this `NARX model example <https://fastcan.readthedocs.io/en/latest/auto_examples/plot_narx.html>`_.
 
-
-Support Free-Threaded Wheels
-----------------------------
-fastcan has support for free-threaded (also known as nogil) CPython 3.13.
-For more information about free-threaded CPython, check `how to install a free-threaded CPython <https://py-free-threading.github.io/installing_cpython/>`_.
-
 Support WASM Wheels
 -------------------
 fastcan is compiled to WebAssembly (WASM) wheels using `pyodide <https://github.com/pyodide/pyodide>`_.

fastcan/narx/_base.py

@@ -718,59 +718,54 @@ def _get_hc_ids(
         d2ydx2[k, i] += Constant terms
         """
         n_degree = jac_feat_ids.shape[1]
-        hess_yyd_ids = np.zeros((0, 3), dtype=np.int32)
-        hess_yd_ids = np.zeros((0, 2), dtype=np.int32)
-        hess_coef_ids = np.zeros(0, dtype=np.int32)
-        hess_feat_ids = np.zeros((0, n_degree), dtype=np.int32)
-        hess_delay_ids = np.zeros((0, n_degree), dtype=np.int32)
 
         if mode < 2:
             return (
-                hess_yyd_ids,
-                hess_yd_ids,
-                hess_coef_ids,
-                hess_feat_ids,
-                hess_delay_ids,
+                np.zeros((0, 3), dtype=np.int32),  # hess_yyd_ids
+                np.zeros((0, 2), dtype=np.int32),  # hess_yd_ids
+                np.zeros(0, dtype=np.int32),  # hess_coef_ids
+                np.zeros((0, n_degree), dtype=np.int32),  # hess_feat_ids
+                np.zeros((0, n_degree), dtype=np.int32),  # hess_delay_ids
             )
 
-        for yyd_id, coef_id, feat_ids, delay_ids in zip(
-            jac_yyd_ids, jac_coef_ids, jac_feat_ids, jac_delay_ids
-        ):
-            # In jac, x * term will generate a term with coef 1.
-            # In hess, it will be
-            # d term / dx = 1 * d y_in * yi * yj
-            hess_yyd_ids = np.vstack([hess_yyd_ids, yyd_id])
-            hess_yd_ids = np.vstack([hess_yd_ids, [-1, -1]])  # empty
-            # constant 1 handled in _update_hc
-            hess_coef_ids = np.append(hess_coef_ids, coef_id)
-            hess_feat_ids = np.vstack([hess_feat_ids, feat_ids])
-            hess_delay_ids = np.vstack([hess_delay_ids, delay_ids])
-            for var_id, (feat_id, delay_id) in enumerate(zip(feat_ids, delay_ids)):
-                #  d JC / dx = coef * d y_in * d yi * yj
-                # hess_yyd_ids: y_out and y_in
-                # hess_coef_ids: coef
-                # hess_feat_ids: yj ..
-                # hess_delay_ids: yj ..
-                # feat_ids and delay_ids contain spaceholder -1 to keep poly_degree size
-                # Skip input x and spaceholder -1
-                if feat_id >= n_features_in and delay_id > 0:
-                    # when feat_id is output y, drop it from hess_feat_ids
-                    hess_yd_ids = np.vstack(
-                        [hess_yd_ids, [feat_id - n_features_in, delay_id]]
-                    )
-                    hess_yyd_ids = np.vstack([hess_yyd_ids, yyd_id])
-                    hess_coef_ids = np.append(hess_coef_ids, coef_id)
-                    hess_feat_ids = np.vstack([hess_feat_ids, feat_ids])
-                    hess_delay_ids = np.vstack([hess_delay_ids, delay_ids])
-                    hess_feat_ids[-1][var_id] = -1
-                    hess_delay_ids[-1][var_id] = -1
+        # In jac, x * term will generate two parts: 1 * term + x * d term / dx
+        # In hess, the part of a term with coef 1 will be
+        # d term / dx = 1 * d y_in * yi * yj
+        # constant 1 in hess_coef_ids handled in _update_hc
+        # The second part will be
+        #  d JC / dx = coef * d y_in * d yi * yj
+        # hess_yyd_ids: y_out and y_in
+        # hess_coef_ids: coef
+        # hess_feat_ids: yj ..
+        # hess_delay_ids: yj ..
+        # feat_ids and delay_ids contain spaceholder -1 to keep poly_degree size
+        # Skip input x and spaceholder -1
+        n_d_coefs = len(jac_yyd_ids)
+        mask = (jac_feat_ids >= n_features_in) & (jac_delay_ids > 0)
+        d_terms_ids, var_ids = np.nonzero(mask)
+
+        hess_yyd_ids = np.vstack([jac_yyd_ids, jac_yyd_ids[d_terms_ids]])
+        yd_d_coefs = np.full((n_d_coefs, 2), -1, dtype=np.int32)
+        yd_d_terms = np.column_stack(
+            [jac_feat_ids[mask] - n_features_in, jac_delay_ids[mask]]
+        )
+        hess_yd_ids = np.vstack([yd_d_coefs, yd_d_terms])
+        hess_coef_ids = np.concatenate([jac_coef_ids, jac_coef_ids[d_terms_ids]])
+
+        hess_feat_ids = np.vstack([jac_feat_ids, jac_feat_ids[d_terms_ids]])
+        hess_delay_ids = np.vstack([jac_delay_ids, jac_delay_ids[d_terms_ids]])
+
+        row_indices = np.arange(n_d_coefs, len(hess_yyd_ids))
+        # when feat_id is output y, drop it from hess_feat_ids
+        hess_feat_ids[row_indices, var_ids] = -1
+        hess_delay_ids[row_indices, var_ids] = -1
 
         return (
-            hess_yyd_ids.astype(np.int32),
-            hess_yd_ids.astype(np.int32),
-            hess_coef_ids.astype(np.int32),
-            hess_feat_ids.astype(np.int32),
-            hess_delay_ids.astype(np.int32),
+            hess_yyd_ids.astype(np.int32, copy=False),
+            hess_yd_ids.astype(np.int32, copy=False),
+            hess_coef_ids.astype(np.int32, copy=False),
+            hess_feat_ids.astype(np.int32, copy=False),
+            hess_delay_ids.astype(np.int32, copy=False),
         )
 
     @staticmethod
@@ -802,35 +797,26 @@ def _get_jc_ids(feat_ids, delay_ids, output_ids, n_features_in):
         axis-2 (j) input y: dy0(k-d)/dx, dy1(k-d)/dx, ..., dyn(k-d)/dx
         """
 
-        n_degree = feat_ids.shape[1]
-        jac_yyd_ids = np.zeros((0, 3), dtype=np.int32)
-        jac_coef_ids = np.zeros(0, dtype=int)
-        jac_feat_ids = np.zeros((0, n_degree), dtype=np.int32)
-        jac_delay_ids = np.zeros((0, n_degree), dtype=np.int32)
+        mask = (feat_ids >= n_features_in) & (delay_ids > 0)
+        jac_coef_ids, var_ids = np.nonzero(mask)
+        n_rows = jac_coef_ids.shape[0]
 
-        for coef_id, (term_feat_ids, term_delay_ids) in enumerate(
-            zip(feat_ids, delay_ids)
-        ):
-            out_y_id = output_ids[coef_id]  # y(k, id), output
-            for var_id, (feat_id, delay_id) in enumerate(
-                zip(term_feat_ids, term_delay_ids)
-            ):
-                if feat_id >= n_features_in and delay_id > 0:
-                    in_y_id = feat_id - n_features_in  # y(k-d, id), input
-                    jac_yyd_ids = np.vstack(
-                        [jac_yyd_ids, [out_y_id, in_y_id, delay_id]]
-                    )
-                    jac_coef_ids = np.append(jac_coef_ids, coef_id)
-                    jac_feat_ids = np.vstack([jac_feat_ids, term_feat_ids])
-                    jac_delay_ids = np.vstack([jac_delay_ids, term_delay_ids])
-                    jac_feat_ids[-1][var_id] = -1
-                    jac_delay_ids[-1][var_id] = -1
+        jac_yyd_ids = np.empty((n_rows, 3), dtype=np.int32)
+        jac_yyd_ids[:, 0] = output_ids[jac_coef_ids]  # y(k, id), output
+        jac_yyd_ids[:, 1] = feat_ids[mask] - n_features_in  # y(k-d, id), input
+        jac_yyd_ids[:, 2] = delay_ids[mask]
+
+        jac_feat_ids = feat_ids[jac_coef_ids]
+        jac_delay_ids = delay_ids[jac_coef_ids]
+        row_indices = np.arange(n_rows)
+        jac_feat_ids[row_indices, var_ids] = -1
+        jac_delay_ids[row_indices, var_ids] = -1
 
         return (
-            jac_yyd_ids.astype(np.int32),
-            jac_coef_ids.astype(np.int32),
-            jac_feat_ids.astype(np.int32),
-            jac_delay_ids.astype(np.int32),
+            jac_yyd_ids.astype(np.int32, copy=False),
+            jac_coef_ids.astype(np.int32, copy=False),
+            jac_feat_ids.astype(np.int32, copy=False),
+            jac_delay_ids.astype(np.int32, copy=False)
         )
 
     @staticmethod