Feature / save long-running sim outputs to cache. Manager preempts sending identical work in favor of pulling from cache

.claude/skills/generate-scripts/references/results_metadata.md

@@ -39,4 +39,3 @@ If the minimum objective value is exactly 0.0, check whether those rows have
 `sim_ended == True`. Unevaluated rows often have fields initialized to zero.
 This is common for the last few rows when the simulation budget is exhausted — they were
 allocated by the generator but never evaluated.
-

.github/workflows/extra.yml

@@ -76,6 +76,7 @@ jobs:
           rm ./libensemble/tests/unit_tests/test_ufunc_runners.py # needs globus-compute
           rm ./libensemble/tests/regression_tests/test_gpCAM.py # needs gpcam, which doesn't build on 3.13
           rm ./libensemble/tests/regression_tests/test_asktell_gpCAM.py # needs gpcam, which doesn't build on 3.13
+          rm ./libensemble/tests/regression_tests/test_persistent_gp_multitask_ax.py # needs ax-platform, which doesn't yet support 3.14
           rm ./libensemble/tests/regression_tests/test_optimas_ax_mf.py # needs ax-platform, which doesn't yet support 3.14
           rm ./libensemble/tests/regression_tests/test_optimas_ax_sf.py # needs ax-platform, which doesn't yet support 3.14
 

.pre-commit-config.yaml

@@ -37,4 +37,6 @@ repos:
   rev: v1.19.1
   hooks:
     - id: mypy
+      pass_filenames: false
+      args: [--package=libensemble.utils]
       exclude: ^docs/conf\.py$|libensemble/utils/(launcher|loc_stack|runners|pydantic|output_directory)\.py$|libensemble/tests/(regression_tests|functionality_tests|unit_tests|scaling_tests)/.*

docs/data_structures/libE_specs/libE_specs_general.rst

@@ -38,3 +38,20 @@ General
 **gen_workers** [list of ints]:
     List of workers that should run only generators. All other workers will run
     only simulator functions.
+
+**cache_long_sims** [bool] = ``False``:
+    Cache simulation results with runtimes >1s to disk. Subsequent runs of the same
+    base script with the same command-line arguments will access this cache.
+
+    Upon the generator creating points already in the cache, those points will be skipped from
+    being sent for evaluation. Instead the corresponding cached results are retrieved and returned
+    to the generator.
+
+    The cache is saved in ``cache_dir``, and by default is named after the joined command-line arguments.
+
+**cache_dir** [str] = ``"~/.cache/libensemble"``:
+    The directory to store the cache file. Defaults to ``~/.cache/libensemble``.
+
+**cache_name** [str] = ``"_".join(sys.argv)``:
+    The name of the cache file. Stored in ``cache_dir``, and by default is named after the
+    joined command-line arguments.

libensemble/alloc_funcs/start_only_persistent.py

@@ -68,17 +68,20 @@ def only_persistent_gens(W, H, sim_specs, gen_specs, alloc_specs, persis_info, l
     gen_count = support.count_persis_gens()
     Work = {}
 
-    # Asynchronous return to generator
-    async_return = user.get("async_return", False) and sum(H["sim_ended"]) >= initial_batch_size
+    # Asynchronous return to generator.
+    # Use the manager-maintained counter instead of re-scanning the full H array.
+    async_return = user.get("async_return", False) and libE_info["sim_ended_count"] >= initial_batch_size
 
     if gen_count < persis_info.get("num_gens_started", 0):
         # When a persistent worker is done, trigger a shutdown (returning exit condition of 1)
         return Work, persis_info, 1
 
-    # Give evaluated results back to a running persistent gen
+    # Give evaluated results back to a running persistent gen.
+    # Compute the sim_ended & ~gen_informed mask once; AND with per-worker gen_inds inside the loop.
+    pending_sim = H["sim_ended"] & ~H["gen_informed"]
     for wid in support.avail_worker_ids(persistent=EVAL_GEN_TAG, active_recv=active_recv_gen):
         gen_inds = H["gen_worker"] == wid
-        returned_but_not_given = np.logical_and.reduce((H["sim_ended"], ~H["gen_informed"], gen_inds))
+        returned_but_not_given = pending_sim & gen_inds
         if np.any(returned_but_not_given):
             if async_return or support.all_sim_ended(H, gen_inds):
                 point_ids = np.where(returned_but_not_given)[0]

libensemble/gen_classes/sampling.py

@@ -156,9 +156,7 @@ class UniformSampleWithVariableResources(LibensembleGenerator):
        path was tested with the default alloc.
     """
 
-    def __init__(
-        self, vocs: VOCS, max_resource_sets: int, random_seed: int = 1, *args, **kwargs
-    ):
+    def __init__(self, vocs: VOCS, max_resource_sets: int, random_seed: int = 1, *args, **kwargs):
         super().__init__(vocs, *args, **kwargs)
         self.rng = np.random.default_rng(random_seed)
         self.max_rsets = max_resource_sets

libensemble/history.py

@@ -1,5 +1,7 @@
+import json
 import logging
 import time
+from pathlib import Path
 
 import numpy as np
 import numpy.typing as npt
@@ -102,6 +104,8 @@ def __init__(
         self.index = len(H0)
         self.grow_count = 0
         self.safe_mode = False
+        self.use_cache = False
+        self.cache_set = False
 
         self.sim_started_count = np.sum(H["sim_started"])
         self.sim_ended_count = np.sum(H["sim_ended"])
@@ -115,6 +119,64 @@ def __init__(
         self.last_started = -1
         self.last_ended = -1
 
+    def init_cache(
+        self,
+        cache_name: str,
+        cache_dir: str | Path,
+        spec_hash: str | None = None,
+    ) -> None:
+        self.cache_dir = Path(cache_dir).expanduser()
+        self.cache_dir.mkdir(parents=True, exist_ok=True)
+        self.cache = self.cache_dir / Path(cache_name + ".npy")
+        self.cache_meta = self.cache_dir / Path(cache_name + ".meta.json")
+        self.spec_hash = spec_hash
+        self.use_cache = True
+        self.cache_set = False
+
+        # Precompute the sorted user-field names and their dtypes once, so
+        # _shelf_longrunning_sims doesn't recompute them on every sim return.
+        libE_field_names = {k[0] for k in libE_fields}
+        self.cache_keys = sorted([n for n in self.H.dtype.names if n not in libE_field_names])
+        self.cache_dtype = np.dtype(sorted([(n, self.H.dtype.fields[n][0]) for n in self.cache_keys]))
+
+        # Buffer for new entries collected during this run; deduplicated via bytes key.
+        self._cache_buffer: list = []
+        self._cache_seen: set = set()
+
+        # Validate any existing cache against the configuration hash.
+        cache_valid = False
+        if self.cache.exists():
+            if self.cache_meta.exists():
+                try:
+                    with open(self.cache_meta) as f:
+                        meta = json.load(f)
+                    if meta.get("spec_hash") == spec_hash:
+                        cache_valid = True
+                except (json.JSONDecodeError, KeyError):
+                    pass
+            if not cache_valid:
+                logger.debug(
+                    "Cache hash mismatch or missing metadata — starting fresh: %s",
+                    self.cache.name,
+                )
+                self.cache.unlink(missing_ok=True)
+
+        if not self.cache.exists():
+            self.cache.touch()
+
+        try:
+            self.in_cache = np.load(self.cache, allow_pickle=True)
+        except EOFError:
+            self.in_cache = None
+
+        # Pre-populate the seen-set from any on-disk entries so we don't re-add them.
+        # Also mark cache_set=True immediately when there is existing data — the manager
+        # uses this flag to decide whether to scan the cache when dispatching sim work.
+        if self.in_cache is not None and len(self.in_cache) > 0:
+            for row in self.in_cache:
+                self._cache_seen.add(row.tobytes())
+            self.cache_set = True
+
     def _append_new_fields(self, H_f: npt.NDArray) -> None:
         dtype_new = np.dtype(list(set(self.H.dtype.descr + np.lib.recfunctions.repack_fields(H_f).dtype.descr)))
         H_new = np.zeros(len(self.H), dtype=dtype_new)
@@ -123,6 +185,82 @@ def _append_new_fields(self, H_f: npt.NDArray) -> None:
             H_new[field][: len(self.H)] = self.H[field]
         self.H = H_new
 
+    def _shelf_longrunning_sims(self, index):
+        """Cache any f values that ran for more than a second.
+
+        Uses a bytes-keyed set for O(1) deduplication instead of np.unique on
+        every insertion, and accumulates new entries in a plain Python list that
+        is only materialised into a structured array at save_cache() time.
+        """
+        if self.H[index]["sim_ended_time"] - self.H[index]["sim_started_time"] <= 1:
+            return
+        entry = np.array([self.H[index][self.cache_keys]], dtype=self.cache_dtype)
+        key = entry[0].tobytes()
+        if key in self._cache_seen:
+            return
+        self._cache_seen.add(key)
+        self._cache_buffer.append(entry)
+        self.cache_set = True
+
+    def _materialize_cache(self) -> npt.NDArray | None:
+        """Combine the on-disk cache with any buffered new entries into one array."""
+        parts = []
+        if self.in_cache is not None:
+            parts.append(self.in_cache)
+        if self._cache_buffer:
+            parts.append(np.concatenate(self._cache_buffer))
+        if not parts:
+            return None
+        return np.concatenate(parts) if len(parts) > 1 else parts[0]
+
+    def save_cache(self) -> None:
+        if self.use_cache and self.cache_set:
+            combined = self._materialize_cache()
+            if combined is not None:
+                np.save(self.cache, combined, allow_pickle=True)
+                if self.spec_hash:
+                    with open(self.cache_meta, "w") as f:
+                        json.dump({"spec_hash": self.spec_hash}, f)
+
+    def get_shelved_sims(self) -> npt.NDArray:
+        combined = self._materialize_cache()
+        return combined if combined is not None else np.load(self.cache, allow_pickle=True)
+
+    @staticmethod
+    def _classify_fields(fields, returned_H, H):
+        """Partition returned fields into three buckets for update_history_f.
+
+        Returns
+        -------
+        scalar_fields : list[str]
+            Fields whose per-row value is a scalar or object (can be assigned
+            with a single fancy-indexed write across all rows).
+        uniform_fields : list[str]
+            Fixed-shape array fields whose shape exactly matches H's storage
+            shape (can also be assigned in one fancy-indexed write).
+        ragged_fields : list[str]
+            Fixed-shape array fields that are *smaller* than H's storage shape
+            (need per-row slice assignment).
+        """
+        scalar_fields = []
+        uniform_fields = []
+        ragged_fields = []
+        for field in fields:
+            if field in protected_libE_fields:
+                continue
+            dt = returned_H.dtype[field]
+            if dt.shape == () or dt.hasobject:
+                scalar_fields.append(field)
+            else:
+                # Compare element shape: returned vs H's allocated shape
+                h_shape = H.dtype[field].shape
+                r_shape = dt.shape
+                if r_shape == h_shape:
+                    uniform_fields.append(field)
+                else:
+                    ragged_fields.append(field)
+        return scalar_fields, uniform_fields, ragged_fields
+
     def update_history_f(self, D: dict, kill_canceled_sims: bool = False) -> None:
         """
         Updates the history after points have been evaluated
@@ -135,30 +273,42 @@ def update_history_f(self, D: dict, kill_canceled_sims: bool = False) -> None:
         if returned_H is not None and any([field not in self.H.dtype.names for field in returned_H.dtype.names]):
             self._append_new_fields(returned_H)
 
-        for j, ind in enumerate(new_inds):
+        if self.safe_mode:
             for field in fields:
-                if field in protected_libE_fields:
-                    if self.safe_mode:
-                        assert False, "The field '" + field + "' is protected"
-                    continue
+                assert field not in protected_libE_fields, "The field '" + field + "' is protected"
 
-                if np.isscalar(returned_H[field][j]) or returned_H.dtype[field].hasobject:
-                    self.H[field][ind] = returned_H[field][j]
-                else:
-                    # len or np.size
+        new_inds = np.asarray(new_inds)
+
+        if fields and returned_H is not None:
+            scalar_fields, uniform_fields, ragged_fields = self._classify_fields(fields, returned_H, self.H)
+
+            # Vectorized assignment for scalar and object fields (one op per field)
+            for field in scalar_fields:
+                self.H[field][new_inds] = returned_H[field]
+
+            # Vectorized assignment for fixed-shape array fields that exactly match H's shape
+            for field in uniform_fields:
+                self.H[field][new_inds] = returned_H[field]
+
+            # Per-row loop only for ragged (partial-fill) array fields
+            for j, ind in enumerate(new_inds):
+                for field in ragged_fields:
                     H0_size = len(returned_H[field][j])
                     assert H0_size <= len(self.H[field][ind]), (
                         "History update Error: Too many values received for " + field
                     )
                     assert H0_size, "History update Error: No values in this field " + field
-                    if H0_size == len(self.H[field][ind]):
-                        self.H[field][ind] = returned_H[field][j]  # ref
-                    else:
-                        self.H[field][ind][:H0_size] = returned_H[field][j]  # Slice View
-
-            self.H["sim_ended"][ind] = True
-            self.H["sim_ended_time"][ind] = time.time()
-            self.sim_ended_count += 1
+                    self.H[field][ind][:H0_size] = returned_H[field][j]
+
+        # Batch-update bookkeeping fields for all returned rows at once
+        t = time.time()
+        self.H["sim_ended"][new_inds] = True
+        self.H["sim_ended_time"][new_inds] = t
+        self.sim_ended_count += len(new_inds)
+
+        if self.use_cache:
+            for ind in new_inds:
+                self._shelf_longrunning_sims(ind)
 
         if kill_canceled_sims:
             for j in range(self.last_ended + 1, np.max(new_inds) + 1):

libensemble/libE.py

@@ -137,7 +137,7 @@
 from libensemble.tools.alloc_support import AllocSupport
 from libensemble.tools.tools import _USER_SIM_ID_WARNING
 from libensemble.utils import launcher
-from libensemble.utils.misc import specs_dump
+from libensemble.utils.misc import compute_config_hash, specs_dump
 from libensemble.utils.timer import Timer
 from libensemble.version import __version__
 from libensemble.worker import worker_main
@@ -235,12 +235,27 @@ def libE(
         exit_criteria=exit_criteria,
     )
 
+    # Compute a deterministic hash of the full configuration for cache integrity.
+    spec_hash = compute_config_hash(
+        sim_specs=ensemble.sim_specs,
+        gen_specs=ensemble.gen_specs,
+        alloc_specs=ensemble.alloc_specs,
+        libE_specs=ensemble.libE_specs,
+        exit_criteria=ensemble.exit_criteria,
+        H0=H0,
+    )
+
     (sim_specs, gen_specs, alloc_specs, libE_specs) = [
         specs_dump(spec, by_alias=True)
         for spec in [ensemble.sim_specs, ensemble.gen_specs, ensemble.alloc_specs, ensemble.libE_specs]
     ]
     exit_criteria = specs_dump(ensemble.exit_criteria, by_alias=True, exclude_none=True)
 
+    # Inject spec hash and auto-generate cache name when not explicitly provided.
+    libE_specs["_spec_hash"] = spec_hash
+    if libE_specs.get("cache_long_sims") and not libE_specs.get("cache_name"):
+        libE_specs["cache_name"] = f".libe_cache_{spec_hash[:16]}"
+
     # Restore objects that don't survive serialization via model_dump
     if hasattr(ensemble.sim_specs, "simulator") and ensemble.sim_specs.simulator is not None:
         sim_specs["simulator"] = ensemble.sim_specs.simulator