diff --git a/.cspell/custom-dictionary-workspace.txt b/.cspell/custom-dictionary-workspace.txt
index 1de20418b..b12724784 100644
--- a/.cspell/custom-dictionary-workspace.txt
+++ b/.cspell/custom-dictionary-workspace.txt
@@ -420,6 +420,7 @@ Roboto
 rowspan
 rstart
 rstrip
+rtol
 rtype
 ruamel
 Rvrt
diff --git a/apps/predbat/load_predictor.py b/apps/predbat/load_predictor.py
index 25bf6cc85..da0d75f9a 100644
--- a/apps/predbat/load_predictor.py
+++ b/apps/predbat/load_predictor.py
@@ -611,14 +611,23 @@ def _create_dataset(self, load_minutes, now_utc, pv_minutes=None, temp_minutes=N
         # Validation window: most recent chunks (chunk_idx 0 … validation_end_chunk-1)
         validation_end_chunk = validation_holdout_hours * 60 // CHUNK_MINUTES
 
-        X_train_list = []
+        # Preallocate the feature matrices and write each sample straight into its row. Holding
+        # a Python list of per-sample arrays and then copying it into np.array() keeps both
+        # alive at once, which for a three week window is an extra 70MB on top of the 66MB
+        # matrix. Rows are sized for the maximum possible sample count and the filled prefix is
+        # returned; samples are only skipped where the history has gaps, so the slack is small.
+        max_train_rows = max(max_chunk_idx - LOOKBACK_STEPS, 0)
+        max_val_rows = max(validation_end_chunk, 0)
+        X_train_all = np.empty((max_train_rows, TOTAL_FEATURES), dtype=np.float32)
+        X_val_all = np.empty((max_val_rows, TOTAL_FEATURES), dtype=np.float32)
+        train_rows = 0
+        val_rows = 0
         y_train_list = []
         weight_list = []
-        X_val_list = []
         y_val_list = []
 
-        def _build_sample(target_chunk_idx):
-            """Build one (features, target) sample centred on target_chunk_idx."""
+        def _build_sample(target_chunk_idx, out_row):
+            """Write one sample's features into out_row and return its target, or None on a gap."""
             lookback_start = target_chunk_idx + 1
             lookback_values = []
             pv_lookback_values = []
@@ -635,12 +644,12 @@ def _build_sample(target_chunk_idx):
                     import_rate_lookback.append(chunked_import.get(lb_idx, 0.0))
                     export_rate_lookback.append(chunked_export.get(lb_idx, 0.0))
                 else:
-                    return None, None  # Gap in data - skip
+                    return None  # Gap in data - skip
 
             if len(lookback_values) != LOOKBACK_STEPS:
-                return None, None
+                return None
             if target_chunk_idx not in chunked_energy:
-                return None, None
+                return None
 
             target_value = chunked_energy[target_chunk_idx]
 
@@ -651,25 +660,22 @@ def _build_sample(target_chunk_idx):
             day_of_year = target_time.timetuple().tm_yday
             time_features = self._create_time_features(minute_of_day, day_of_week, day_of_year)
 
-            features = np.concatenate(
-                [
-                    np.array(lookback_values, dtype=np.float32),
-                    np.array(pv_lookback_values, dtype=np.float32),
-                    np.array(temp_lookback_values, dtype=np.float32),
-                    np.array(import_rate_lookback, dtype=np.float32),
-                    np.array(export_rate_lookback, dtype=np.float32),
-                    time_features,
-                ]
-            )
-            return features, np.array([target_value], dtype=np.float32)
+            # Assigning the lists straight into row slices converts them in place, avoiding the
+            # five temporary arrays and the concatenated result this used to build per sample
+            offset = 0
+            for values in (lookback_values, pv_lookback_values, temp_lookback_values, import_rate_lookback, export_rate_lookback):
+                out_row[offset : offset + LOOKBACK_STEPS] = values
+                offset += LOOKBACK_STEPS
+            out_row[offset:] = time_features
+            return np.array([target_value], dtype=np.float32)
 
         # Training samples: all available chunks
         for target_chunk_idx in range(0, max_chunk_idx - LOOKBACK_STEPS):
-            features, target = _build_sample(target_chunk_idx)
-            if features is None:
+            target = _build_sample(target_chunk_idx, X_train_all[train_rows])
+            if target is None:
                 continue
 
-            X_train_list.append(features)
+            train_rows += 1
             y_train_list.append(target)
 
             # Time-decay weighting (older samples get lower weight)
@@ -678,23 +684,23 @@ def _build_sample(target_chunk_idx):
 
         # Validation samples: most recent validation_end_chunk chunks
         for target_chunk_idx in range(0, validation_end_chunk):
-            features, target = _build_sample(target_chunk_idx)
-            if features is None:
+            target = _build_sample(target_chunk_idx, X_val_all[val_rows])
+            if target is None:
                 continue
-            X_val_list.append(features)
+            val_rows += 1
             y_val_list.append(target)
 
-        if not X_train_list:
+        if not train_rows:
             return None, None, None, None, None
 
-        X_train = np.array(X_train_list, dtype=np.float32)
+        X_train = X_train_all[:train_rows]
         y_train = np.array(y_train_list, dtype=np.float32)
         train_weights = np.array(weight_list, dtype=np.float32)
 
         # Normalize weights to sum to number of samples
         train_weights = train_weights * len(train_weights) / np.sum(train_weights)
 
-        X_val = np.array(X_val_list, dtype=np.float32) if X_val_list else None
+        X_val = X_val_all[:val_rows] if val_rows else None
         y_val = np.array(y_val_list, dtype=np.float32) if y_val_list else None
 
         return X_train, y_train, train_weights, X_val, y_val
@@ -808,6 +814,43 @@ def _ar_rollout_diagnostic(self, load_minutes, now_utc, pv_minutes=None, temp_mi
 
         return float(np.mean(errors)), float(np.mean(biases))
 
+    @staticmethod
+    def _feature_mean_std(X, block=512):
+        """
+        Compute per-feature mean and std over row blocks, accumulating in float64.
+
+        np.mean/np.std on a float32 array accumulate in float32, which loses low bits on the
+        feature groups that sit far from zero, and np.std materialises the deviations array
+        internally - another full copy of a 66MB training matrix to produce 1,446 numbers.
+        Two passes over row blocks avoids both: no temporary larger than one block, and
+        float64 accumulators that match a float64 reference.
+
+        Args:
+            X: Feature array, shape (samples, features)
+            block: Rows accumulated per pass. The float64 deviation buffer is
+                   block * features * 8 bytes, so this bounds the working set: 512 rows of
+                   1,446 features is 5.9MB against 47.4MB at 4096.
+
+        Returns:
+            Tuple of (mean, std) as float32 arrays
+        """
+        rows = len(X)
+        total = np.zeros(X.shape[1], dtype=np.float64)
+        for start in range(0, rows, block):
+            total += X[start : start + block].sum(axis=0, dtype=np.float64)
+        mean = total / rows
+
+        squares = np.zeros(X.shape[1], dtype=np.float64)
+        for start in range(0, rows, block):
+            deviation = X[start : start + block].astype(np.float64)
+            deviation -= mean
+            # Square into the deviation buffer rather than allocating a second one of the
+            # same size; the values are identical either way
+            squares += np.square(deviation, out=deviation).sum(axis=0)
+        std = np.sqrt(squares / rows)
+
+        return mean.astype(np.float32), std.astype(np.float32)
+
     def _get_min_std_array(self, n_features):
         """
         Return the per-feature minimum std array used to prevent extreme normalization.
@@ -818,7 +861,10 @@ def _get_min_std_array(self, n_features):
         Returns:
             numpy array of minimum std values, shape (n_features,)
         """
-        min_std = np.ones(n_features) * 1e-8  # Default fallback
+        # float32 to match the dataset and weights: a float64 minimum promotes the whole
+        # normalised feature matrix to float64 via np.maximum, doubling its size and running
+        # training in double precision against float32 weights
+        min_std = np.ones(n_features, dtype=np.float32) * 1e-8  # Default fallback
         if n_features == TOTAL_FEATURES:
             min_std[0:LOOKBACK_STEPS] = 0.01  # Load energy (kWh)
             min_std[LOOKBACK_STEPS : 2 * LOOKBACK_STEPS] = 0.01  # PV energy (kWh)
@@ -858,7 +904,7 @@ def _log_normalization_stats(self, label=""):
 
         self.log("ML Predictor: Normalization stats [{}] target(mean={:.4f} std={:.4f}) {}".format(label, self.target_mean if self.target_mean is not None else 0, self.target_std if self.target_std is not None else 0, " ".join(parts)))
 
-    def _normalize_features(self, X, fit=False, ema_alpha=0.0):
+    def _normalize_features(self, X, fit=False, ema_alpha=0.0, in_place=False):
         """
         Normalize features using z-score normalization with feature-specific minimum stds.
 
@@ -868,13 +914,17 @@ def _normalize_features(self, X, fit=False, ema_alpha=0.0):
             ema_alpha: If > 0 and existing params exist, blend new stats with old via EMA
                        (new = alpha * new_stats + (1-alpha) * old_stats). Used during
                        fine-tuning to track feature distribution drift without sudden jumps.
+            in_place: If True, write the normalised values back into X instead of building a
+                      new array. A training feature matrix is around 66MB, and the caller
+                      keeps no use for the un-normalised values, so copying doubles the
+                      resident cost of a training pass for nothing. Only pass this when the
+                      caller is finished with the array it hands in.
 
         Returns:
             Normalized feature array
         """
         if fit:
-            self.feature_mean = np.mean(X, axis=0)
-            self.feature_std = np.std(X, axis=0)
+            self.feature_mean, self.feature_std = self._feature_mean_std(X)
 
             # Clamp std to per-feature minimums to prevent extreme normalization
             self.feature_std = np.maximum(self.feature_std, self._get_min_std_array(len(self.feature_std)))
@@ -882,8 +932,7 @@ def _normalize_features(self, X, fit=False, ema_alpha=0.0):
 
         elif ema_alpha > 0 and self.feature_mean is not None and self.feature_std is not None:
             # EMA update: blend new statistics with existing to track distribution drift
-            new_mean = np.mean(X, axis=0)
-            new_std = np.std(X, axis=0)
+            new_mean, new_std = self._feature_mean_std(X)
 
             # Apply same min-std clamping to new stats before blending
             new_std = np.maximum(new_std, self._get_min_std_array(len(new_std)))
@@ -896,6 +945,11 @@ def _normalize_features(self, X, fit=False, ema_alpha=0.0):
         if self.feature_mean is None or self.feature_std is None:
             return X
 
+        if in_place:
+            X -= self.feature_mean
+            X /= self.feature_std
+            return X
+
         return (X - self.feature_mean) / self.feature_std
 
     def _normalize_targets(self, y, fit=False):
@@ -1039,13 +1093,13 @@ def train(
         # On initial train: fit normalization from scratch
         # On fine-tune: apply EMA update to track distribution drift gradually
         if is_initial or not self.model_initialized:
-            X_train_norm = self._normalize_features(X_train, fit=True)
+            X_train_norm = self._normalize_features(X_train, fit=True, in_place=True)
             y_train_norm = self._normalize_targets(y_train, fit=True)
         else:
-            X_train_norm = self._normalize_features(X_train, fit=False, ema_alpha=norm_ema_alpha)
+            X_train_norm = self._normalize_features(X_train, fit=False, ema_alpha=norm_ema_alpha, in_place=True)
             y_train_norm = self._normalize_targets(y_train, fit=False)
             self.log("ML Predictor: Applied EMA normalization update (alpha={}) to track feature drift".format(norm_ema_alpha))
-        X_val_norm = self._normalize_features(X_val, fit=False)
+        X_val_norm = self._normalize_features(X_val, fit=False, in_place=True)
         y_val_norm = self._normalize_targets(y_val, fit=False)
 
         # Initialise weights if needed
diff --git a/apps/predbat/tests/test_ml_memory.py b/apps/predbat/tests/test_ml_memory.py
new file mode 100644
index 000000000..711308fdb
--- /dev/null
+++ b/apps/predbat/tests/test_ml_memory.py
@@ -0,0 +1,329 @@
+# fmt: off
+"""
+ML Training Memory Tests
+
+Training a curriculum pass over three weeks of history builds a feature matrix of
+around 11,500 samples by 1,446 float32 features - 66.7MB per copy. These tests pin the
+numerics of dataset construction and normalisation so the memory work cannot change
+them, and assert the normalisation reuses its input buffer rather than copying it.
+"""
+
+import math
+from datetime import datetime, timezone
+
+import numpy as np
+
+from load_predictor import LoadPredictor, TOTAL_FEATURES
+
+
+def _predictor_with_stats(samples=64, seed=7):
+    """Build a predictor with normalisation stats fitted to a fixed synthetic matrix."""
+    np.random.seed(seed)
+    features = (np.random.randn(samples, TOTAL_FEATURES) * 3.0 + 1.5).astype(np.float32)
+    predictor = LoadPredictor(learning_rate=0.001)
+    return predictor, features
+
+
+def test_normalize_features_values_are_stable(my_predbat=None):
+    """Normalisation produces the documented z-score, whatever the implementation"""
+    print("\n=== Testing _normalize_features() numerics ===")
+    failed = 0
+
+    predictor, features = _predictor_with_stats()
+    reference = features.copy()
+
+    result = predictor._normalize_features(features.copy(), fit=True)
+
+    expected = (reference - predictor.feature_mean) / predictor.feature_std
+    if not np.allclose(result, expected, rtol=1e-5, atol=1e-6):
+        worst = float(np.max(np.abs(result - expected)))
+        print("ERROR: normalised values drifted from the z-score definition, max diff {}".format(worst))
+        failed += 1
+    else:
+        print("✓ normalised output matches (X - mean) / std")
+
+    # The fitted statistics themselves must stay accurate against a float64 reference
+    ref_mean = reference.astype(np.float64).mean(axis=0)
+    if not np.allclose(predictor.feature_mean, ref_mean, rtol=1e-4, atol=1e-5):
+        worst = float(np.max(np.abs(predictor.feature_mean - ref_mean)))
+        print("ERROR: fitted mean drifted from the float64 reference, max diff {}".format(worst))
+        failed += 1
+    else:
+        print("✓ fitted mean matches a float64 reference")
+
+    ref_std = reference.astype(np.float64).std(axis=0)
+    if not np.allclose(predictor.feature_std, ref_std, rtol=1e-3, atol=1e-4):
+        worst = float(np.max(np.abs(predictor.feature_std - ref_std)))
+        print("ERROR: fitted std drifted from the float64 reference, max diff {}".format(worst))
+        failed += 1
+    else:
+        print("✓ fitted std matches a float64 reference")
+
+    return failed
+
+
+def test_normalize_features_reuses_the_input_buffer(my_predbat=None):
+    """in_place normalisation must not allocate a second copy of the feature matrix"""
+    print("\n=== Testing _normalize_features() in-place behaviour ===")
+    failed = 0
+
+    predictor, features = _predictor_with_stats()
+    predictor._normalize_features(features.copy(), fit=True)
+
+    subject = features.copy()
+    result = predictor._normalize_features(subject, fit=False, in_place=True)
+
+    if result is not subject:
+        print("ERROR: in_place normalisation returned a new array; at 11,500 samples that is a 66.7MB copy")
+        return failed + 1
+    print("✓ in_place normalisation writes into the array it was given")
+
+    expected = (features - predictor.feature_mean) / predictor.feature_std
+    if not np.allclose(result, expected, rtol=1e-5, atol=1e-6):
+        print("ERROR: in_place normalisation produced different values to the copying path")
+        failed += 1
+    else:
+        print("✓ in_place result matches the copying path")
+
+    # The default must keep the caller's array untouched
+    untouched = features.copy()
+    original = untouched.copy()
+    predictor._normalize_features(untouched, fit=False)
+    if not np.array_equal(untouched, original):
+        print("ERROR: default call mutated the caller's array")
+        failed += 1
+    else:
+        print("✓ default call leaves the caller's array untouched")
+
+    return failed
+
+
+def test_normalisation_stays_float32(my_predbat=None):
+    """Normalising must not promote the feature matrix to float64
+
+    The network's weights, the dataset and the activations are all float32. A float64
+    minimum-std array promotes the whole normalised matrix to float64, doubling a 66.7MB
+    training matrix and running every matmul in double precision against float32 weights.
+    """
+    print("\n=== Testing normalisation dtype ===")
+    failed = 0
+
+    predictor, features = _predictor_with_stats()
+    predictor._normalize_features(features.copy(), fit=True)
+
+    if predictor._get_min_std_array(TOTAL_FEATURES).dtype != np.float32:
+        print("ERROR: minimum-std array is {}, expected float32".format(predictor._get_min_std_array(TOTAL_FEATURES).dtype))
+        failed += 1
+    else:
+        print("✓ minimum-std array is float32")
+
+    if predictor.feature_std.dtype != np.float32:
+        print("ERROR: fitted std is {}, expected float32".format(predictor.feature_std.dtype))
+        failed += 1
+    else:
+        print("✓ fitted std stays float32")
+
+    result = predictor._normalize_features(features.copy(), fit=False)
+    if result.dtype != np.float32:
+        print("ERROR: normalised matrix is {}, expected float32 - a 66.7MB matrix becomes 133.4MB".format(result.dtype))
+        failed += 1
+    else:
+        print("✓ normalised matrix stays float32")
+
+    return failed
+
+
+def _synthetic_history(days=5):
+    """Build deterministic load/pv/temp/rate dicts spanning the given number of days."""
+    minutes = days * 24 * 60
+    load, pv, temp, imp, exp = {}, {}, {}, {}, {}
+    for minute in range(0, minutes, 5):
+        phase = minute / 180.0
+        load[minute] = 0.05 + 0.03 * math.sin(phase)
+        pv[minute] = max(0.0, 0.09 * math.sin(phase / 4.0))
+        temp[minute] = 12.0 + 5.0 * math.cos(phase / 6.0)
+        imp[minute] = 20.0 + 4.0 * math.sin(phase / 3.0)
+        exp[minute] = 8.0 + 2.0 * math.cos(phase / 5.0)
+    return load, pv, temp, imp, exp
+
+
+def test_create_dataset_output_is_stable(my_predbat=None):
+    """Building the feature matrix must produce the same data whatever the construction
+
+    This pins the dataset so the memory work - preallocating the matrix instead of
+    building a Python list of per-sample rows - cannot silently change what is trained on.
+    """
+    print("\n=== Testing _create_dataset() output ===")
+    failed = 0
+
+    load, pv, temp, imp, exp = _synthetic_history(days=5)
+    now = datetime(2026, 8, 17, 12, 0, 0, tzinfo=timezone.utc)
+    predictor = LoadPredictor(learning_rate=0.001)
+
+    X_train, y_train, weights, X_val, y_val = predictor._create_dataset(load, now, pv_minutes=pv, temp_minutes=temp, import_rates=imp, export_rates=exp, validation_holdout_hours=24)
+
+    if X_train is None:
+        print("ERROR: no dataset produced from 5 days of synthetic history")
+        return failed + 1
+
+    for name, array, dtype in (("X_train", X_train, np.float32), ("y_train", y_train, np.float32), ("weights", weights, np.float32), ("X_val", X_val, np.float32)):
+        if array.dtype != dtype:
+            print("ERROR: {} dtype is {}, expected {}".format(name, array.dtype, dtype))
+            failed += 1
+
+    if X_train.shape[1] != TOTAL_FEATURES:
+        print("ERROR: X_train has {} features, expected {}".format(X_train.shape[1], TOTAL_FEATURES))
+        failed += 1
+    elif X_train.shape[0] != y_train.shape[0] or X_train.shape[0] != weights.shape[0]:
+        print("ERROR: ragged dataset: X {} y {} weights {}".format(X_train.shape, y_train.shape, weights.shape))
+        failed += 1
+    else:
+        print("✓ {} samples x {} features, aligned with targets and weights".format(X_train.shape[0], X_train.shape[1]))
+
+    if not np.all(np.isfinite(X_train)):
+        print("ERROR: X_train contains non-finite values")
+        failed += 1
+    else:
+        print("✓ all feature values finite")
+
+    # Time-decay weights are normalised to sum to the sample count
+    if abs(float(np.sum(weights)) - len(weights)) > 0.01 * len(weights):
+        print("ERROR: weights sum to {}, expected about {}".format(float(np.sum(weights)), len(weights)))
+        failed += 1
+    else:
+        print("✓ time-decay weights normalised to the sample count")
+
+    # A checksum over the feature matrix catches any change in row content or ordering.
+    # Recorded from the list-of-rows construction this replaced; if numpy ever changes its
+    # summation this may need re-recording, but it must not move for a refactor.
+    checksum = float(np.sum(X_train.astype(np.float64)))
+    expected_checksum = 12977609.094096
+    if abs(checksum - expected_checksum) > abs(expected_checksum) * 1e-9:
+        print("ERROR: feature matrix changed - checksum {:.6f}, expected {:.6f}".format(checksum, expected_checksum))
+        failed += 1
+    else:
+        print("✓ feature matrix checksum unchanged ({:.6f})".format(checksum))
+
+    if X_train.shape[0] != 1151:
+        print("ERROR: sample count changed: {} rows, expected 1151".format(X_train.shape[0]))
+        failed += 1
+    else:
+        print("✓ sample count unchanged (1151)")
+
+    return failed
+
+
+def test_fitted_statistics_are_accurate_on_large_datasets(my_predbat=None):
+    """Fitting normalisation stats must not lose precision to float32 accumulation
+
+    np.mean/np.std accumulate in float32 for a float32 input, which drifts on the feature
+    groups that sit far from zero - temperature around 12C, import rates around 20p. It also
+    materialises the deviations array internally, another full copy of a 66.7MB matrix just
+    to produce 1,446 numbers. Accumulating over row blocks in float64 fixes both.
+    """
+    print("\n=== Testing fitted statistics accuracy ===")
+    failed = 0
+
+    np.random.seed(3)
+    # Large mean, small spread: the case where a float32 accumulator loses its low bits
+    features = (np.random.randn(20000, TOTAL_FEATURES // 32).astype(np.float32) * 0.01 + 500.0).astype(np.float32)
+
+    predictor = LoadPredictor(learning_rate=0.001)
+    predictor._normalize_features(features.copy(), fit=True)
+
+    reference = features.astype(np.float64)
+    mean_error = float(np.max(np.abs(predictor.feature_mean - reference.mean(axis=0))))
+    std_error = float(np.max(np.abs(predictor.feature_std - reference.std(axis=0))))
+
+    # A float32 accumulator lands around 2.2e-04 on this data. float64 accumulation reduces
+    # that to the float32 representation floor of the stored result - eps * 500 is about 6e-05
+    # - so the threshold sits between the two rather than at zero.
+    if mean_error > 5e-5:
+        print("ERROR: fitted mean drifted {:.3e} from the float64 reference (float32 accumulation)".format(mean_error))
+        failed += 1
+    else:
+        print("✓ fitted mean accurate to {:.3e}".format(mean_error))
+
+    if std_error > 1e-6:
+        print("ERROR: fitted std drifted {:.3e} from the float64 reference".format(std_error))
+        failed += 1
+    else:
+        print("✓ fitted std accurate to {:.3e}".format(std_error))
+
+    return failed
+
+
+def test_create_dataset_skips_gaps_in_the_history(my_predbat=None):
+    """Samples whose lookback window crosses a gap must be skipped, not half-written
+
+    Real history has gaps where the sensor stopped reporting. A sample spanning one cannot
+    be built, so the row is skipped - the dataset must stay dense and the targets must line
+    up with it.
+    """
+    print("\n=== Testing _create_dataset() with gaps ===")
+    failed = 0
+
+    load, pv, temp, imp, exp = _synthetic_history(days=6)
+    # Punch two holes in the load data, as a sensor dropout would
+    for minute in range(2 * 24 * 60, 2 * 24 * 60 + 600, 5):
+        load.pop(minute, None)
+    for minute in range(4 * 24 * 60, 4 * 24 * 60 + 300, 5):
+        load.pop(minute, None)
+
+    now = datetime(2026, 8, 17, 12, 0, 0, tzinfo=timezone.utc)
+    predictor = LoadPredictor(learning_rate=0.001)
+
+    result = predictor._create_dataset(load, now, pv_minutes=pv, temp_minutes=temp, import_rates=imp, export_rates=exp, validation_holdout_hours=24)
+    X_train, y_train, weights, X_val, y_val = result
+
+    if X_train is None:
+        print("ERROR: gaps caused the whole dataset to be discarded")
+        return failed + 1
+    print("✓ dataset still built with gaps present")
+
+    if X_train.shape[0] != y_train.shape[0]:
+        print("ERROR: {} feature rows but {} targets - skipped rows were still counted".format(X_train.shape[0], y_train.shape[0]))
+        failed += 1
+    elif X_train.shape[0] != weights.shape[0]:
+        print("ERROR: {} feature rows but {} weights".format(X_train.shape[0], weights.shape[0]))
+        failed += 1
+    else:
+        print("✓ {} rows, targets and weights all aligned".format(X_train.shape[0]))
+
+    if y_train.dtype != np.float32 or not np.all(np.isfinite(y_train)):
+        print("ERROR: targets are {} and finite={}".format(y_train.dtype, bool(np.all(np.isfinite(y_train)))))
+        failed += 1
+    else:
+        print("✓ targets are finite float32")
+
+    if not np.all(np.isfinite(X_train)):
+        print("ERROR: feature matrix has non-finite values - a skipped row was left unwritten")
+        failed += 1
+    else:
+        print("✓ every feature row fully written")
+
+    return failed
+
+
+def run_ml_memory_tests(my_predbat=None):
+    """Run all ML training memory tests"""
+    print("\n" + "=" * 80)
+    print("ML Training Memory Tests")
+    print("=" * 80)
+
+    failed = 0
+    failed += test_normalize_features_values_are_stable(my_predbat)
+    failed += test_normalize_features_reuses_the_input_buffer(my_predbat)
+    failed += test_normalisation_stays_float32(my_predbat)
+    failed += test_create_dataset_output_is_stable(my_predbat)
+    failed += test_fitted_statistics_are_accurate_on_large_datasets(my_predbat)
+    failed += test_create_dataset_skips_gaps_in_the_history(my_predbat)
+
+    print("\n" + "=" * 80)
+    if failed == 0:
+        print("✅ All ML training memory tests passed!")
+    else:
+        print(f"❌ {failed} ML training memory test(s) failed")
+    print("=" * 80)
+
+    return failed
diff --git a/apps/predbat/unit_test.py b/apps/predbat/unit_test.py
index 7453dfa97..72f5ddfb8 100644
--- a/apps/predbat/unit_test.py
+++ b/apps/predbat/unit_test.py
@@ -193,6 +193,7 @@
 from tests.test_mock_base import test_mock_base_all
 from tests.test_solis import run_solis_tests
 from tests.test_load_ml import test_load_ml
+from tests.test_ml_memory import run_ml_memory_tests
 from tests.test_temperature import test_temperature
 from tests.test_oauth_mixin import run_oauth_mixin_tests
 from tests.test_fox_oauth import run_fox_oauth_tests
@@ -547,6 +548,8 @@ def main():
         ("tariff_catalogue", test_tariff_catalogue, "Tariff catalogue tests", False),
         ("annual_integration", run_annual_integration_isolated, "Annual prediction integration tests", True),
         ("load_ml", test_load_ml, "ML Load Forecaster tests (MLP, training, persistence, validation)", True),
+        # ML training memory: dataset construction, normalisation dtype and statistics accuracy
+        ("ml_memory", run_ml_memory_tests, "ML training memory tests", False),
         ("random", run_random_scenario_tests, "Random scenario plan regression against the committed baseline", False),
     ]