diff --git a/README.md b/README.md
index e5a9d4b..1fc2a20 100644
--- a/README.md
+++ b/README.md
@@ -75,13 +75,25 @@ noll_zk = params['noll_zk']
 
 ## Minimal usage (NeuralActiveOpticsSys)
 
+Production and RA-style runs should pass **checkpoint paths** so WaveNet and AlignNet weights come from disk only (no Hugging Face Hub or internet):
+
 ```python
 from tarts.NeuralActiveOpticsSys import NeuralActiveOpticsSys
 
-model = NeuralActiveOpticsSys(dataset_params='TARTS/python/tarts/dataset_params.yaml')
-# Use forward/forward_align/forward_shifts, or deploy_run for LSST exposures
+model = NeuralActiveOpticsSys(
+    dataset_params="TARTS/python/tarts/dataset_params.yaml",
+    wavenet_path="/path/to/wavenet.ckpt",
+    alignet_path="/path/to/alignnet.ckpt",
+    aggregatornet_path="/path/to/aggregator.ckpt",  # optional if aggregator_on=False
+)
 ```
 
+For local experiments without checkpoints, `NeuralActiveOpticsSys` defaults `pretrained=False` so timm/torchvision **does not** download ImageNet backbones at import time. Set `pretrained=True` only when you want those downloads (needs network and a writable Hugging Face cache).
+
+### Hugging Face cache / offline hosts
+
+Backbone downloads (when `pretrained=True`) use the Hugging Face Hub cache. Set **`HF_HOME`** to a writable directory (or **`TARTS_HF_HOME`**, which sets `HF_HOME` if it is not already set) **before** `import tarts` if `HOME` is not writable or points at a non-existent path. For air-gapped use, pre-populate that cache (or bake it into the image) and set **`HF_HUB_OFFLINE=1`**.
+
 ## Notes
 
 - When integrating with Rubin/LSST data, ensure the LSST science stack is available; this package itself does not enforce LSST dependencies by default.
diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
index 7216e98..2002f58 100644
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@@ -1,5 +1,5 @@
 {% set name = "tarts" %}
-{% set version = "0.0.0a7" %}  # Update this to your current version
+{% set version = "4.2" %}
 {% set python_min = "3.10" %}
 
 package:
diff --git a/python/tarts/NeuralActiveOpticsSys.py b/python/tarts/NeuralActiveOpticsSys.py
index 0250629..4365476 100644
--- a/python/tarts/NeuralActiveOpticsSys.py
+++ b/python/tarts/NeuralActiveOpticsSys.py
@@ -4,10 +4,14 @@
 import copy
 import logging
 import os
+import time
+import warnings
+from pathlib import Path
 from typing import Any, Dict, List
 
 # Third-party imports
 import joblib
+import numpy as np
 import pytorch_lightning as pl
 import torch
 import torch.nn.functional as F_loss
@@ -36,9 +40,50 @@
 
 logger = logging.getLogger(__name__)
 
+try:
+    from sklearn.exceptions import InconsistentVersionWarning as _SklearnInconsistentVersionWarning
+except ImportError:
+    _SklearnInconsistentVersionWarning = None
+
 
 class NeuralActiveOpticsSys(pl.LightningModule):
-    """Transfer learning driven WaveNet."""
+    """Neural Active Optics System: AlignNet + WaveNet + AggregatorNet pipeline.
+
+    For inference speed you can:
+    - Compile submodules only: pass compile_models=True (each submodel is compiled).
+    - Compile the entire system: use create_compiled() or wrap with torch.compile(model)
+      after construction to fuse across submodules (may have graph breaks due to dynamic
+      control flow in forward).
+    """
+
+    @classmethod
+    def create_compiled(cls, *, backend=None, mode=None, **kwargs):
+        """Build a NeuralActiveOpticsSys and wrap it with torch.compile for full-system optimization.
+
+        Use this for deployment when you want one compiled graph across AlignNet, WaveNet,
+        and AggregatorNet. The first run may be slow (compilation); later runs are faster.
+
+        Parameters
+        ----------
+        backend : str, optional
+            torch.compile backend (e.g. "inductor" for CPU, None for default GPU).
+        mode : str, optional
+            torch.compile mode ("default", "reduce-overhead", "max-autotune").
+        **kwargs
+            Passed to NeuralActiveOpticsSys.__init__().
+
+        Returns
+        -------
+        compiled_model
+            torch.compile(NeuralActiveOpticsSys(**kwargs)), with optional backend and mode.
+        """
+        model = cls(**kwargs)
+        compile_kw = {}
+        if backend is not None:
+            compile_kw["backend"] = backend
+        if mode is not None:
+            compile_kw["mode"] = mode
+        return torch.compile(model, **compile_kw)
 
     def __init__(
         self,
@@ -49,7 +94,7 @@ def __init__(
         lr=1e-3,
         final_layer=None,
         aggregator_on=True,
-        pretrained=True,
+        pretrained=False,
         compile_models=False,
         ood_model_path=None,
         kmin=None,
@@ -76,16 +121,21 @@ def __init__(
             Whether to use aggregator network for final prediction. Defaults to True.
         pretrained : bool, optional
             Whether to use pre-trained CNN weights when creating new models (when checkpoint paths are None).
-            Set to False for deployment mode to avoid downloading weights. Defaults to True.
-            Note: This parameter is ignored when loading from checkpoint files.
+            Defaults to False so Summit/USDF/RA starts do not contact the Hugging Face Hub or require a
+            writable ImageNet-weight cache. Set True only when you intend to download backbone weights
+            (requires network and a writable ``HF_HOME`` / default cache). Ignored when loading from
+            checkpoint files (checkpoints carry full weights).
         compile_models : bool, optional
             Whether to apply torch.compile to the submodels (WaveNet, AlignNet, AggregatorNet).
             This can significantly speed up inference but may increase compilation time on first run.
             Automatically selects backend: "inductor" for CPU, default for GPU.
-            Defaults to False.
+            For full-system compilation (one fused graph), use create_compiled() or
+            torch.compile(model) after construction. Defaults to False.
         ood_model_path : str, optional
             Path to OOD detection model (joblib file). If provided, OOD detection will be performed
             during inference and scores will be stored in internal metadata. Defaults to None.
+            Prefer artifacts saved with a ``precision`` array (see ``training/fit_OOD_model.py``) so
+            inference does not unpickle sklearn estimators across version mismatches.
         kmin : float, optional
             Minimum frequency magnitude for frequency filtering. If None, frequency filtering is disabled.
             If set, applies frequency filter to each cropped image before passing to WaveNet.
@@ -121,12 +171,41 @@ def __init__(
         # Load OOD detection model if path is provided
         self.ood_model = None
         self.ood_mean = None
+        self.ood_precision: torch.Tensor | None = None
         if ood_model_path is not None and os.path.exists(ood_model_path):
             logger.info(f"Loading OOD detection model from {ood_model_path}...")
-            ood_data = joblib.load(ood_model_path)
-            self.ood_model = ood_data["cov_model"]
+            with warnings.catch_warnings():
+                # Suppress sklearn's "unpickle estimator from different version" warning.
+                # This warning can appear when `ood_model.joblib` was created by an older
+                # sklearn version (e.g. LedoitWolf) and is being unpickled at runtime.
+                #
+                # We suppress by category (if available) and also by message text as a fallback
+                # to make this robust across sklearn versions / import-path changes.
+                if _SklearnInconsistentVersionWarning is not None:
+                    warnings.simplefilter("ignore", _SklearnInconsistentVersionWarning)
+                warnings.filterwarnings(
+                    "ignore",
+                    message=r"Trying to unpickle estimator.*",
+                    category=Warning,
+                )
+                ood_data = joblib.load(ood_model_path)
             self.ood_mean = torch.tensor(ood_data["mean"], device=self.device_val, dtype=torch.float32)
-            logger.info("OOD detection model loaded successfully")
+            if ood_data.get("precision") is not None:
+                self.ood_precision = torch.as_tensor(
+                    ood_data["precision"], dtype=torch.float32, device=self.device_val
+                )
+                logger.info("OOD detector: using saved precision matrix (sklearn-free inference).")
+            elif ood_data.get("cov_model") is not None:
+                self.ood_model = ood_data["cov_model"]
+                logger.warning(
+                    "OOD artifact uses legacy pickled sklearn estimator; re-run training/fit_OOD_model.py "
+                    "to save a precision-based artifact and avoid sklearn version warnings."
+                )
+            else:
+                self.ood_mean = None
+                logger.warning("OOD file has no 'precision' or 'cov_model'; OOD detection disabled.")
+            if self.ood_mean is not None:
+                logger.info("OOD detection model loaded successfully")
         elif ood_model_path is not None:
             logger.warning(f"OOD model path provided but file not found: {ood_model_path}")
             logger.info("Continuing without OOD detection...")
@@ -135,22 +214,32 @@ def __init__(
         with open(dataset_params, "r") as yaml_file:
             params = yaml.safe_load(yaml_file)
 
+        cnn_model = params.get("cnn_model", "resnet34")
+        if (wavenet_path is None or alignet_path is None) and pretrained:
+            logger.warning(
+                "NeuralActiveOpticsSys: building WaveNet/AlignNet without checkpoint(s) and with "
+                "pretrained=True will download backbone weights (timm/Hugging Face Hub) and write "
+                "under HF_HOME or HOME. For production use checkpoint paths and leave pretrained=False."
+            )
+
         if wavenet_path is None:
-            self.wavenet_model = WaveNetSystem(pretrained=pretrained).to(self.device_val)
+            self.wavenet_model = WaveNetSystem(cnn_model=cnn_model, pretrained=pretrained).to(self.device_val)
         else:
-            # Always use checkpoint loading - the pretrained parameter doesn't matter
-            # when loading from checkpoint
+            # Override any legacy checkpoint hyperparameter so Lightning does
+            # not initialize a pretrained backbone before loading weights.
             self.wavenet_model = WaveNetSystem.load_from_checkpoint(
-                wavenet_path, map_location=str(self.device_val), strict=False
+                wavenet_path, map_location=str(self.device_val), strict=False, pretrained=False
             ).to(self.device_val)
 
         if alignet_path is None:
-            self.alignnet_model = AlignNetSystem(pretrained=pretrained).to(self.device_val)
+            self.alignnet_model = AlignNetSystem(cnn_model=cnn_model, pretrained=pretrained).to(
+                self.device_val
+            )
         else:
-            # Always use checkpoint loading - the pretrained parameter doesn't matter
-            # when loading from checkpoint
+            # Override any legacy checkpoint hyperparameter so Lightning does
+            # not initialize a pretrained backbone before loading weights.
             self.alignnet_model = AlignNetSystem.load_from_checkpoint(
-                alignet_path, map_location=str(self.device_val), strict=False
+                alignet_path, map_location=str(self.device_val), strict=False, pretrained=False
             ).to(self.device_val)
 
         self.max_seq_length = params["max_seq_len"]
@@ -191,6 +280,24 @@ def __init__(
         self.SCALE = params["adjustment_AlignNet"]
         self.num_zernikes = len(params["noll_zk"])
 
+        # LSSTCam deployment: Z4/focus correction from branch_prediction_by_band_detector
+        # (see focus_offset.yaml).
+        focus_offset_path = Path(__file__).resolve().parent / "focus_offset.yaml"
+        self._lsstcam_focus_offset_table: Dict[str, Dict[str, float]] | None = None
+        if focus_offset_path.is_file():
+            try:
+                with open(focus_offset_path, "r") as f:
+                    focus_data = yaml.safe_load(f)
+                z4_cfg = focus_data.get("z4_peak_pattern") or {}
+                tbl = z4_cfg.get("branch_prediction_by_band_detector")
+                if isinstance(tbl, dict) and tbl:
+                    self._lsstcam_focus_offset_table = tbl
+                    logger.info("Loaded LSSTCam focus offset table from %s", focus_offset_path)
+            except (OSError, yaml.YAMLError, TypeError, ValueError) as e:
+                logger.warning("Could not load focus offset config from %s: %s", focus_offset_path, e)
+        else:
+            logger.debug("Focus offset file not found at %s", focus_offset_path)
+
         # Apply torch.compile to submodels if requested
         if compile_models:
             # Determine compilation backend based on device
@@ -289,7 +396,8 @@ def get_internal_data(self):
         - band: Band information
         - SNR: Signal-to-noise ratio
         - centers: Center coordinates
-        - zernikes: Estimated Zernike coefficients
+        - zk_deviations_CCS: Legacy alias for the estimated Zernike coefficients
+        - zernikes: Estimated Zernike coefficients in CCS convention
         - ood_score: Out-of-distribution score (if OOD detection is enabled)
 
         Returns:
@@ -304,8 +412,8 @@ def get_internal_data(self):
         if len(self.fx) > 0 and self.cropped_image is not None and self.total_zernikes is not None:
             num_donuts = len(self.fx)
 
-            # Ensure these are not None before indexing
             for i in range(num_donuts):
+                zernikes = self.total_zernikes[i].clone().detach()
                 data_dict = {
                     "cropped_image": self.cropped_image[i].clone().detach(),
                     "fx": self.fx[i].clone().detach(),
@@ -314,7 +422,8 @@ def get_internal_data(self):
                     "band": self.band[i].clone().detach(),
                     "SNR": self.SNR[i].clone().detach(),
                     "centers": self.centers[i].clone().detach(),
-                    "zernikes": self.total_zernikes[i].clone().detach(),
+                    "zk_deviations_CCS": zernikes,
+                    "zernikes": zernikes,
                 }
                 # Add OOD score if available
                 if self.ood_scores is not None and i < len(self.ood_scores):
@@ -335,6 +444,9 @@ def get_internal_data(self):
                 "band": torch.tensor(0).clone().detach(),
                 "SNR": torch.tensor(0).clone().detach(),
                 "centers": torch.tensor([0, 0]).clone().detach(),
+                "zk_deviations_CCS": torch.full((self.num_zernikes,), float("nan"), device=self.device_val)
+                .clone()
+                .detach(),
                 "zernikes": torch.full((self.num_zernikes,), float("nan"), device=self.device_val)
                 .clone()
                 .detach(),
@@ -432,6 +544,24 @@ def filter_batch_images(self, images):
 
         return torch.stack(filtered_images, dim=0)
 
+    def _compute_ood_mahalanobis(self, features_np: np.ndarray) -> torch.Tensor:
+        """Mahalanobis distance of rows of ``features_np`` from the OOD training mean.
+
+        Uses a stored precision matrix when available (no sklearn at inference). Legacy artifacts
+        fall back to ``sklearn`` ``LedoitWolf.mahalanobis``.
+        """
+        if self.ood_mean is None:
+            raise RuntimeError("OOD mean is not loaded; cannot compute Mahalanobis distance.")
+        centered = np.ascontiguousarray(features_np - self.ood_mean.cpu().numpy(), dtype=np.float32)
+        if self.ood_precision is not None:
+            fc = torch.from_numpy(centered).to(device=self.device_val, dtype=torch.float32)
+            d2 = torch.einsum("bi,ij,bj->b", fc, self.ood_precision, fc)
+            return torch.sqrt(torch.clamp(d2, min=0.0))
+        if self.ood_model is None:
+            raise RuntimeError("OOD model is not loaded; cannot compute Mahalanobis distance.")
+        mahalanobis_dist = self.ood_model.mahalanobis(centered)
+        return torch.as_tensor(mahalanobis_dist, device=self.device_val, dtype=torch.float32)
+
     def identity(self, x):
         """Return the input unchanged (identity function).
 
@@ -675,17 +805,14 @@ def forward(
 
         # Compute OOD scores if OOD model is available
         ood_scores = None
-        if self.ood_model is not None and self.wavenet_model.wavenet.predictor_features is not None:
-            # Get predictor penultimate features and detach/convert to CPU for numpy
-            penultimate = self.wavenet_model.wavenet.predictor_features  # Shape: (batch_size, n_features)
-
-            # Detach from computation graph and move to CPU for numpy conversion
+        if (
+            self.ood_mean is not None
+            and (self.ood_precision is not None or self.ood_model is not None)
+            and self.wavenet_model.wavenet.predictor_features is not None
+        ):
+            penultimate = self.wavenet_model.wavenet.predictor_features
             features_np = penultimate.detach().cpu().numpy()
-            if self.ood_mean is not None:
-                features_centered = features_np - self.ood_mean.cpu().numpy()
-                mahalanobis_dist = self.ood_model.mahalanobis(features_centered)
-                # Store as tensor
-                ood_scores = torch.tensor(mahalanobis_dist, device=self.device_val, dtype=torch.float32)
+            ood_scores = self._compute_ood_mahalanobis(features_np)
 
         # Ensure all tensors are on the same device before concatenation
         device = total_zernikes.device
@@ -729,18 +856,17 @@ def forward(
 
         embedded_features = embedded_features[None, ...]
 
-        # Match training data mean computation:
-        # Training data: zk_mean1 = torch.mean(zk_pred1, dim=0) / 1000
-        # Since total_zernikes is already divided by 1000, we don't divide again
-        mean_zernike = torch.mean(total_zernikes, dim=0)
-
+        # Use robust median reduction across crop-level predictions.
+        # Since total_zernikes is already divided by 1000, we don't divide again.
+        median_zernike = torch.median(total_zernikes, dim=0).values
+        self.mean_zernike = median_zernike
         # (OPTIONAL) Check the types
         if self.aggregator_on:
             # Match training data: mean_zernike should NOT have convert_zernikes_deploy applied
             # Training data stores mean directly without conversion
-            final_zernike_prediction = self.aggregatornet_model((embedded_features, mean_zernike))
+            final_zernike_prediction = self.aggregatornet_model((embedded_features, median_zernike))
         else:
-            final_zernike_prediction = mean_zernike
+            final_zernike_prediction = median_zernike
         final_zernike_prediction = self.final_layer(final_zernike_prediction)
 
         return final_zernike_prediction
@@ -843,18 +969,14 @@ def forward_shifts(
 
         # Compute OOD scores if OOD model is available
         ood_scores = None
-        if self.ood_model is not None and self.wavenet_model.wavenet.predictor_features is not None:
-            # Get predictor penultimate features and detach/convert to CPU for numpy
-            penultimate = self.wavenet_model.wavenet.predictor_features  # Shape: (batch_size, n_features)
-
-            # Detach from computation graph and move to CPU for numpy conversion
+        if (
+            self.ood_mean is not None
+            and (self.ood_precision is not None or self.ood_model is not None)
+            and self.wavenet_model.wavenet.predictor_features is not None
+        ):
+            penultimate = self.wavenet_model.wavenet.predictor_features
             features_np = penultimate.detach().cpu().numpy()
-            if self.ood_mean is not None:
-                features_centered = features_np - self.ood_mean.cpu().numpy()
-                mahalanobis_dist = self.ood_model.mahalanobis(features_centered)
-
-                # Store as tensor
-                ood_scores = torch.tensor(mahalanobis_dist, device=self.device_val, dtype=torch.float32)
+            ood_scores = self._compute_ood_mahalanobis(features_np)
 
         # Ensure all tensors are on the same device before concatenation
         device = total_zernikes.device
@@ -889,18 +1011,17 @@ def forward_shifts(
 
         embedded_features = embedded_features[None, ...]
 
-        # Match training data mean computation:
-        # Training data: zk_mean1 = torch.mean(zk_pred1, dim=0) / 1000
-        # Since total_zernikes is already divided by 1000, we don't divide again
-        mean_zernike = torch.mean(total_zernikes, dim=0)
+        # Use robust median reduction across crop-level predictions.
+        # Since total_zernikes is already divided by 1000, we don't divide again.
+        median_zernike = torch.median(total_zernikes, dim=0).values
 
         # (OPTIONAL) Check the types
         if self.aggregator_on:
             # Match training data: mean_zernike should NOT have convert_zernikes_deploy applied
             # Training data stores mean directly without conversion
-            final_zernike_prediction = self.aggregatornet_model((embedded_features, mean_zernike))
+            final_zernike_prediction = self.aggregatornet_model((embedded_features, median_zernike))
         else:
-            final_zernike_prediction = mean_zernike
+            final_zernike_prediction = median_zernike
         final_zernike_prediction = self.final_layer(final_zernike_prediction)
 
         return final_zernike_prediction
@@ -979,6 +1100,61 @@ def configure_optimizers(self) -> Any:
             },
         }
 
+    @staticmethod
+    def _filter_str_to_band_letter(filter_str: str) -> str | None:
+        """Map FILTER keyword string to a single band letter (u–y), same heuristics as deploy_run."""
+        if "u" in filter_str:
+            return "u"
+        if "g" in filter_str:
+            return "g"
+        if "r" in filter_str:
+            return "r"
+        if "i" in filter_str:
+            return "i"
+        if "z" in filter_str:
+            return "z"
+        if "y" in filter_str:
+            return "y"
+        return None
+
+    def _apply_lsstcam_focus_offset_if_applicable(
+        self,
+        pred: torch.Tensor,
+        metadata,
+        detector_id: int | str,
+        band_letter: str | None,
+    ) -> torch.Tensor:
+        """If INSTRUME matches LSSTCam (case-insensitive), add focus_offset.yaml correction to index 0."""
+        if self._lsstcam_focus_offset_table is None or band_letter is None:
+            return pred
+        try:
+            instrume = metadata["INSTRUME"]
+        except (KeyError, TypeError):
+            return pred
+        instrume_str = (instrume if isinstance(instrume, str) else str(instrume)).strip()
+        if instrume_str.casefold() != "lsstcam":
+            return pred
+        by_det = self._lsstcam_focus_offset_table.get(band_letter)
+        if not by_det:
+            return pred
+        det_key = str(int(detector_id)) if isinstance(detector_id, (int, float)) else str(detector_id)
+        offset = by_det.get(det_key)
+        if offset is None:
+            return pred
+        off_f = float(offset)
+        logger.info(
+            "LSSTCam global Z4 offset (focus_offset.yaml): Z4 += %s (band=%s detector=%s; "
+            "applied to aggregated output index 0 only; INSTRUME=%r)",
+            off_f,
+            band_letter,
+            det_key,
+            instrume_str,
+        )
+        out = pred.clone()
+        off = torch.tensor(off_f, device=out.device, dtype=out.dtype)
+        out[..., 0] = out[..., 0] + off
+        return out
+
     def deploy_run(self, exposure, detectorName=None):
         """Run inference on a real LSST exposure for deployment.
 
@@ -1001,22 +1177,22 @@ def deploy_run(self, exposure, detectorName=None):
         2. Subtracts background
         3. Extracts metadata (filter, focal plane position)
         4. Computes field coordinates for all detected donuts
-        5. Runs forward pass to predict Zernike coefficients
         """
+        start_t = time.perf_counter()
         camera = LsstCam().getCamera()
         try:
             assembleCcdTask = AssembleCcdTask()
             new = assembleCcdTask.assembleCcd(exposure)
             SubtractBackground = subtractBackground.SubtractBackgroundTask()
             SubtractBackground.run(new)
-        except (AttributeError, RuntimeError, ValueError, LengthError) as e:
-            logger.warning(f"Switching to no CCD assembly: {e}")
+        except (AttributeError, RuntimeError, ValueError, LengthError):
             new = exposure
             SubtractBackground = subtractBackground.SubtractBackgroundTask()
             SubtractBackground.run(new)
         image = new.getImage().array
         header = exposure.metadata
         filter_name = header["FILTER"]
+        band_letter = self._filter_str_to_band_letter(filter_name)
         if detectorName is None:
             full_detectorName = header["RAFTBAY"] + "_" + header["CCDSLOT"]
             detectorName = MAP_DETECTOR_TO_NUMBER[full_detectorName]
@@ -1067,7 +1243,16 @@ def deploy_run(self, exposure, detectorName=None):
 
         image_tensor = F.to_tensor(image)[None, ...]
         with torch.no_grad():
-            pred = self.forward(image_tensor, field_x, field_y, focal_val, band_val)
+            pred = self.forward(
+                image_tensor,
+                field_x,
+                field_y,
+                focal_val,
+                band_val,
+            )
+        pred = self._apply_lsstcam_focus_offset_if_applicable(pred, header, detectorName, band_letter)
+        elapsed_s = time.perf_counter() - start_t
+        logger.info(f"NAOS.deploy_run: output computed in {elapsed_s:.3f}s")
         return pred
 
     def deploy_run_shifts(self, exposure, detectorName=None, shift_amount=5):
@@ -1107,6 +1292,7 @@ def deploy_run_shifts(self, exposure, detectorName=None, shift_amount=5):
         image = new.getImage().array
         header = exposure.metadata
         filter_name = header["FILTER"]
+        band_letter = self._filter_str_to_band_letter(filter_name)
         if detectorName is None:
             detectorName = header["CHIPID"]
         #  U G R I Z Y
@@ -1156,8 +1342,15 @@ def deploy_run_shifts(self, exposure, detectorName=None, shift_amount=5):
 
         image_tensor = F.to_tensor(image)[None, ...]
         with torch.no_grad():
-            pred = self.forward_shifts(image_tensor, field_x, field_y, focal_val, band_val, shift_amount)
-        return pred
+            pred = self.forward_shifts(
+                image_tensor,
+                field_x,
+                field_y,
+                focal_val,
+                band_val,
+                shift_amount,
+            )
+        return self._apply_lsstcam_focus_offset_if_applicable(pred, header, detectorName, band_letter)
 
     def deploy_detect(self, exposure, detectorName=None):
         """Run inference on a real LSST exposure for deployment.
diff --git a/python/tarts/__init__.py b/python/tarts/__init__.py
index 02e4af2..7e52551 100644
--- a/python/tarts/__init__.py
+++ b/python/tarts/__init__.py
@@ -10,6 +10,15 @@
 - WaveNet: For Zernike coefficient prediction from donut images
 - AggregatorNet: For combining predictions from multiple donuts
 - Utilities for data processing, model training, and deployment
+
+Deployment (Summit / USDF / RA)
+-------------------------------
+- Prefer loading ``NeuralActiveOpticsSys`` with WaveNet/AlignNet checkpoint paths so no
+  backbone weights are fetched at runtime.
+- By default ``pretrained=False`` on ``NeuralActiveOpticsSys`` to avoid Hugging Face Hub
+  downloads and cache writes under ``HOME``.
+- To redirect the Hub cache, set ``HF_HOME`` or ``TARTS_HF_HOME`` (alias used if ``HF_HOME``
+  is unset) before ``import tarts``.
 """
 
 import logging
@@ -29,6 +38,8 @@
     handler.setFormatter(formatter)
     _logger.addHandler(handler)
 
+from . import _runtime_env  # noqa: F401  # configures HF cache before timm/hub imports
+
 from .NeuralActiveOpticsSys import *  # noqa: F403, F401
 from .utils import *  # noqa: F403, F401
 from .dataloader import *  # noqa: F403, F401
diff --git a/python/tarts/_runtime_env.py b/python/tarts/_runtime_env.py
new file mode 100644
index 0000000..e230209
--- /dev/null
+++ b/python/tarts/_runtime_env.py
@@ -0,0 +1,29 @@
+"""Process-wide defaults applied before timm / huggingface_hub are imported.
+
+Summit and USDF often run with no internet and with ``HOME`` pointing at paths
+that are not writable or do not exist. timm backbones with ``pretrained=True``
+download via the Hugging Face Hub into ``HF_HOME`` (or defaults derived from
+``HOME``). Setting ``HF_HOME`` (or this package's alias) early avoids failed
+``makedirs`` under e.g. ``/home/saluser/.cache/...`` on RA hosts.
+"""
+
+from __future__ import annotations
+
+import os
+from pathlib import Path
+
+
+def apply_tarts_runtime_env() -> None:
+    """Apply environment tweaks that must run before Hub-backed models load.
+
+    ``TARTS_HF_HOME``
+        If set and ``HF_HOME`` is unset, ``HF_HOME`` is set to this path
+        (``~`` expanded). Operators may instead set ``HF_HOME`` directly; see
+        Hugging Face Hub documentation.
+    """
+    tarts_hf = os.environ.get("TARTS_HF_HOME")
+    if tarts_hf and not os.environ.get("HF_HOME"):
+        os.environ["HF_HOME"] = os.fspath(Path(tarts_hf).expanduser())
+
+
+apply_tarts_runtime_env()
diff --git a/python/tarts/aggregatornet_coral.py b/python/tarts/aggregatornet_coral.py
index 539abd2..9dd18cf 100644
--- a/python/tarts/aggregatornet_coral.py
+++ b/python/tarts/aggregatornet_coral.py
@@ -10,6 +10,7 @@
 from typing import Any, Tuple
 
 # Third-party imports
+import numpy as np
 import pytorch_lightning as pl
 import torch
 import torch.nn as nn
@@ -69,6 +70,11 @@ def __init__(
         tradeoff_scale: float = 0.001,
         threshold: float = 0.9,
         dare_gram_weight: float = 1.0,
+        exp_rise_x1: float = 0.15,
+        exp_rise_x2: float = 0.16,
+        weight_decay: float = 0.0,
+        coral_mode: bool = False,
+        rotate: bool = False,
     ):
         """Initialize the AggregatorNet_Coral model.
 
@@ -97,6 +103,22 @@ def __init__(
         dare_gram_weight : float, optional, default=1.0
             Overall weight to scale the DARE-GRAM loss relative to regression loss.
             Higher values prioritize domain adaptation over regression accuracy.
+        exp_rise_x1 : float, optional, default=0.15
+            Lower threshold for exponential rise function in DARE-GRAM scaling.
+            When mRSSE <= x1, scale_loss = 1.0 (full DARE-GRAM influence).
+        exp_rise_x2 : float, optional, default=0.16
+            Upper threshold for exponential rise function in DARE-GRAM scaling.
+            When mRSSE >= x2, scale_loss = 0.0 (no DARE-GRAM influence).
+            Exponential transition occurs between x1 and x2.
+        weight_decay : float, optional, default=0.0
+            Weight decay (L2 regularization) coefficient for the optimizer.
+            Helps prevent overfitting by penalizing large weights.
+        coral_mode : bool, optional, default=False
+            Whether coral/real data sampling is enabled for domain adaptation.
+            When True, the dataset samples real observational data alongside simulation data.
+        rotate : bool, optional, default=False
+            Whether to use rotated Zernikes (zk_true_camera) in camera frame.
+            When True, uses camera frame coordinates instead of telescope frame.
 
         Notes
         -----
@@ -116,6 +138,7 @@ def __init__(
             d_model=d_model,
             nhead=nhead,
             dim_feedforward=dim_feedforward,
+            dropout=0.2,  # 50% dropout for regularization
             batch_first=True,
         )
         self.transformer_encoder = nn.TransformerEncoder(encoder_layer, num_layers=num_layers)
@@ -125,6 +148,9 @@ def __init__(
         # Cache for transformer features (for domain adaptation)
         self.transformer_features = None
         self.val_mRSSE: torch.Tensor | None = None
+        # Accumulators for epoch-average validation mRSSE
+        self._val_mRSSE_sum: torch.Tensor | None = None
+        self._val_batch_count: int = 0
 
     def forward(self, x: tuple[torch.Tensor, torch.Tensor], cache_features: bool = False) -> torch.Tensor:
         """Forward pass of the AggregatorNet_Coral model.
@@ -162,7 +188,8 @@ def forward(self, x: tuple[torch.Tensor, torch.Tensor], cache_features: bool = F
         if cache_features:
             # Use the last token's features before FC layer
             # Convert to float32 for SVD operations in DARE-GRAM loss
-            self.transformer_features = x_tensor[:, -1, :].detach().float()
+            # IMPORTANT: Do NOT detach here - we need gradients for DARE-GRAM loss training
+            self.transformer_features = x_tensor[:, -1, :].float()
 
         x_tensor = x_tensor[:, -1, :]  # Take the last token's output
         x_tensor = self.fc(x_tensor)  # Predict the next token
@@ -221,6 +248,7 @@ def dare_gram_loss(self, features_source: torch.Tensor, features_target: torch.T
                 torch.use_deterministic_algorithms(True)
 
         # Determine rank k based on threshold
+        # Use at least eigenvalue at index 1, but prefer threshold if reached
         T = self.hparams.threshold
 
         # Find index where cumulative variance reaches threshold
@@ -258,6 +286,7 @@ def dare_gram_loss(self, features_source: torch.Tensor, features_target: torch.T
             return torch.tensor(0.0, device=self.device)
 
         # Compute pseudo-inverse with low-rank regularization
+        # Add minimum rtol to prevent numerical instability
         rtol_A = max((L_A[k] / L_A[0]).item(), 1e-6)
         rtol_B = max((L_B[k] / L_B[0]).item(), 1e-6)
         A_pinv = torch.linalg.pinv(cov_A, rtol=rtol_A)
@@ -285,36 +314,38 @@ def dare_gram_loss(self, features_source: torch.Tensor, features_target: torch.T
         return dare_gram_loss
 
     def exp_rise_flipped(self, loss, a=6.0):
-        """Exponentially rises from 0 at loss=0.10 to 1 at loss=0.09.
+        """Exponential rise function for DARE-GRAM scaling, flipped to give higher scale at lower loss.
 
-        Then flattens at 0 above 0.10 and 1 below 0.09.
+        Uses exp_rise_x1 and exp_rise_x2 hyperparameters to control transition region.
+        When loss <= x1: scale = 1.0 (full DARE-GRAM)
+        When loss >= x2: scale = 0.0 (no DARE-GRAM)
+        Between x1 and x2: exponential transition
         """
-        loss = torch.as_tensor(loss, dtype=torch.float32, device=self.device)
-
-        # CHANGE HERE:
-        x1 = torch.tensor(0.09, dtype=torch.float32, device=loss.device)  # Peak (1.0)
-        x2 = torch.tensor(0.10, dtype=torch.float32, device=loss.device)  # Start (0.0)
-
-        a = torch.tensor(a, dtype=torch.float32, device=loss.device)
-
-        f = torch.zeros_like(loss, dtype=torch.float32)
-
-        # Region 1: loss <= 0.09 (Before flip = 0)
-        f[loss <= x1] = 0.0
-
-        # Region 2: 0.09 < loss < 0.10
-        mask = (loss > x1) & (loss < x2)
-        t = (loss[mask] - x1) / (x2 - x1)
-        f[mask] = (1 - torch.exp(-a * t)) / (1 - torch.exp(-a))
-
-        # Region 3: loss >= 0.10 (Before flip = 1)
-        f[loss >= x2] = 1.0
+        # Convert loss to scalar
+        loss_val = loss.item() if isinstance(loss, torch.Tensor) else float(loss)
+
+        # Get hyperparameters (must be provided, will error if not)
+        x1 = float(self.hparams.exp_rise_x1)
+        x2 = float(self.hparams.exp_rise_x2)
+        a = float(a)
+
+        # Region 1: loss <= x1 → f = 0
+        if loss_val <= x1:
+            f = 0.0
+        # Region 3: loss >= x2 → f = 1
+        elif loss_val >= x2:
+            f = 1.0
+        # Region 2: x1 < loss < x2 → exponential rise
+        else:
+            t = (loss_val - x1) / (x2 - x1)
+            f = (1 - np.exp(-a * t)) / (1 - np.exp(-a))
 
-        # FLIP:
-        # 0.10 (which was 1) -> becomes 0
-        # 0.09 (which was 0) -> becomes 1
+        # Flip: high f (when loss is high) becomes low scale (when loss is high)
+        # We want: low loss (good) → high scale (more DARE-GRAM), high loss (bad) → low scale (less DARE-GRAM)
         f = -f + 1
-        return f
+
+        # Convert back to tensor on correct device
+        return torch.tensor(f, dtype=torch.float32, device=self.device)
 
     def calc_losses(
         self,
@@ -322,7 +353,7 @@ def calc_losses(
         batch_idx: int,
         use_coral: bool = False,
         add_dare_gram_to_loss: bool = True,
-    ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
+    ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
         """Calculate losses with optional DARE-GRAM domain adaptation.
 
         Parameters
@@ -343,14 +374,14 @@ def calc_losses(
         Returns
         -------
         tuple
-            (total_loss, mRSSE, dare_gram_loss)
+            (total_loss, mRSSE, dare_gram_loss, dare_gram_scale)
             dare_gram_loss is 0 if use_coral=False or coral data not available.
+            dare_gram_scale is the scaling factor from exp_rise_flipped function.
         """
         x, y = batch  # y is the target token
 
         # Check if coral data is present (8 elements in x vs 4 elements)
         has_coral = len(x) == 8
-
         if has_coral:
             x_input, x_mean, filter_name, chipid, coral_x_total, coral_x_mean, coral_filter, coral_chipid = x
         else:
@@ -370,11 +401,13 @@ def calc_losses(
         mRSSE = torch.sqrt(sse).mean()
 
         # DARE-GRAM loss if coral data is available
+        # Skip computation entirely if dare_gram_weight is 0.0 for faster training
         dare_gram_loss = torch.tensor(0.0, device=self.device)
-        if use_coral and has_coral:
+        if use_coral and has_coral and self.hparams.dare_gram_weight > 0.0:
             try:
                 # Forward pass on target/coral data
                 # Use eval mode to prevent BN updates (keep statistics source-domain only)
+                # NOTE: eval() doesn't disable gradients, only affects batch norm/dropout
                 was_training = self.training
                 try:
                     self.eval()
@@ -384,20 +417,27 @@ def calc_losses(
                     if was_training:
                         self.train()
 
-                # Compute DARE-GRAM loss
+                # Compute DARE-GRAM loss (now has gradients since we removed .detach() from features)
                 dare_gram_loss = self.dare_gram_loss(source_features, target_features)
             except (RuntimeError, ValueError, IndexError) as e:
                 logger.warning(f"DARE-GRAM loss computation failed: {e}")
+                # Use zero tensor (no gradients needed if loss computation failed)
                 dare_gram_loss = torch.tensor(0.0, device=self.device)
 
         # Add DARE-GRAM to loss only if requested
         if add_dare_gram_to_loss:
-            scale_loss = self.exp_rise_flipped(self.val_mRSSE if self.val_mRSSE is not None else mRSSE)
+            # Use epoch-averaged val_mRSSE if available, otherwise use current batch mRSSE
+            # (for first training batch before first validation completes)
+            mRSSE_for_scaling = self.val_mRSSE if self.val_mRSSE is not None else mRSSE
+            scale_loss = self.exp_rise_flipped(mRSSE_for_scaling)
             total_loss = regression_loss + self.hparams.dare_gram_weight * scale_loss * dare_gram_loss
         else:
             total_loss = regression_loss
+            # Still compute scale_loss for logging even if not added to loss
+            mRSSE_for_scaling = self.val_mRSSE if self.val_mRSSE is not None else mRSSE
+            scale_loss = self.exp_rise_flipped(mRSSE_for_scaling)
 
-        return total_loss, mRSSE, dare_gram_loss
+        return total_loss, mRSSE, dare_gram_loss, scale_loss
 
     def training_step(self, batch: tuple, batch_idx: int):
         """Perform a single training step with domain adaptation.
@@ -420,10 +460,11 @@ def training_step(self, batch: tuple, batch_idx: int):
         - DARE-GRAM loss is added to regression loss if coral data is available.
         - Training loss, mRSSE, and DARE-GRAM loss are logged for monitoring.
         """
-        loss, mRSSE, dare_gram_loss = self.calc_losses(batch, batch_idx, use_coral=True)
+        loss, mRSSE, dare_gram_loss, dare_gram_scale = self.calc_losses(batch, batch_idx, use_coral=True)
         self.log("train_loss", loss, prog_bar=True, sync_dist=True)
         self.log("train_mRSSE", mRSSE, sync_dist=True)
         self.log("train_dare_gram_loss", dare_gram_loss, sync_dist=True)
+        self.log("train_dare_gram_scale", dare_gram_scale, sync_dist=True)
         return loss
 
     def validation_step(self, batch, batch_idx):
@@ -446,17 +487,40 @@ def validation_step(self, batch, batch_idx):
         - DARE-GRAM is computed for logging but not added to validation loss.
         - This allows monitoring domain adaptation without affecting validation metrics.
         - Validation loss, mRSSE, and DARE-GRAM loss are logged for monitoring.
+        - mRSSE is accumulated across batches to compute epoch average.
         """
         # Compute DARE-GRAM for logging but don't add it to validation loss
-        loss, mRSSE, dare_gram_loss = self.calc_losses(
+        loss, mRSSE, dare_gram_loss, dare_gram_scale = self.calc_losses(
             batch, batch_idx, use_coral=True, add_dare_gram_to_loss=False
         )
         self.log("val_loss", loss, prog_bar=True, sync_dist=True)
         self.log("val_mRSSE", mRSSE, prog_bar=True, sync_dist=True)
         self.log("val_dare_gram_loss", dare_gram_loss, sync_dist=True)
-        self.val_mRSSE = mRSSE.clone().detach()  # Store a copy to avoid tensor reference issues
+        self.log("val_dare_gram_scale", dare_gram_scale, sync_dist=True)
+
+        # Accumulate mRSSE for epoch-average calculation
+        mRSSE_detached = mRSSE.detach()
+        if self._val_mRSSE_sum is None:
+            self._val_mRSSE_sum = torch.zeros_like(mRSSE_detached)
+            self._val_batch_count = 0
+        self._val_mRSSE_sum += mRSSE_detached
+        self._val_batch_count += 1
+
         return loss
 
+    def on_validation_epoch_end(self):
+        """Called at the end of validation epoch to compute epoch-average mRSSE.
+
+        Computes the average validation mRSSE across all batches in the epoch
+        and stores it for use in the next training epoch's DARE-GRAM scaling.
+        """
+        if self._val_mRSSE_sum is not None and self._val_batch_count > 0:
+            # Store epoch-averaged mRSSE for use in next training epoch
+            self.val_mRSSE = (self._val_mRSSE_sum / self._val_batch_count).clone().detach()
+            # Reset accumulators for next epoch
+            self._val_mRSSE_sum = None
+            self._val_batch_count = 0
+
     def loss_fn(self, x, y):
         """Compute the loss using the Root Sum of Squared Errors (mRSSE).
 
@@ -495,7 +559,16 @@ def configure_optimizers(self) -> Any:
         return {
             "optimizer": optimizer,
             "lr_scheduler": {
-                "scheduler": ReduceLROnPlateau(optimizer),
+                "scheduler": ReduceLROnPlateau(
+                    optimizer,
+                    mode="min",
+                    factor=0.5,  # Reduce LR by half when plateau
+                    patience=5,  # Wait 5 validation checks (epochs) before reducing
+                    threshold=1e-3,  # Minimum change to qualify as an improvement
+                    threshold_mode="abs",  # Use absolute threshold
+                    min_lr=1e-7,  # Minimum learning rate
+                    verbose=True,
+                ),
                 "monitor": "val_loss",
                 "frequency": 1,
             },
diff --git a/python/tarts/alignnet.py b/python/tarts/alignnet.py
index a62f3c2..71e942f 100644
--- a/python/tarts/alignnet.py
+++ b/python/tarts/alignnet.py
@@ -1,5 +1,8 @@
 """Neural network to predict donut placement coefficients from donut images and positions."""
 
+# Configure Hugging Face cache (``HF_HOME`` / ``TARTS_HF_HOME``) before timm pulls huggingface_hub.
+from . import _runtime_env  # noqa: F401
+
 # Third-party imports
 import timm
 import torch
@@ -45,7 +48,7 @@ def __init__(
         freeze_cnn: bool = False,
         n_predictor_layers: tuple = (256,),
         device="cuda",
-        pretrained: bool = True,
+        pretrained: bool = False,
     ) -> None:
         """Initialize the NeuralAlignment model.
 
@@ -63,8 +66,8 @@ def __init__(
             The output layer is fixed to have 2 nodes.
         device : str, optional, default='cuda'
             The device to use for computations ('cuda' or 'cpu').
-        pretrained : bool, optional, default=True
-            Whether to use pre-trained CNN weights. Set to False to avoid downloading weights.
+        pretrained : bool, optional, default=False
+            Whether to use pre-trained CNN weights. Set to True only when downloads are acceptable.
 
         Notes
         -----
@@ -213,8 +216,14 @@ def _reshape_image(self, image: torch.Tensor) -> torch.Tensor:
 
         # Get the number of input channels required by the CNN
         if self.is_timm_model:
-            # timm models (MobileNet, etc.)
-            n_channels = self.cnn.conv_stem.in_channels
+            # timm: MobileNet/ConvNeXt use conv_stem, ResNet uses conv1
+            first_conv = getattr(self.cnn, "conv_stem", None) or getattr(self.cnn, "conv1", None)
+            if first_conv is None:
+                msg = (
+                    f"AlignNet CNN backbone {type(self.cnn).__name__} has no " "'conv_stem' or 'conv1' layer."
+                )
+                raise AttributeError(msg)
+            n_channels = first_conv.in_channels
         else:
             # torchvision models (ResNet, etc.)
             n_channels = self.cnn.conv1.in_channels
diff --git a/python/tarts/dataloader.py b/python/tarts/dataloader.py
index 320f69d..6a0f9b8 100644
--- a/python/tarts/dataloader.py
+++ b/python/tarts/dataloader.py
@@ -2,6 +2,7 @@
 
 # Standard library imports
 import glob
+import hashlib
 import logging
 import os
 import pickle
@@ -14,13 +15,52 @@
 from astropy.table import Table
 from torch.utils.data import Dataset
 
+try:
+    from tqdm import tqdm
+except ImportError:
+
+    def tqdm(iterable, **kwargs):
+        """Fallback when tqdm is not installed: return iterable unchanged."""
+        return iterable
+
+
 # Local/application imports
 from .constants import DEFAULT_NOLL_ZK, DEFAULT_TRAIN_FRACTION
-from .utils import shift_offcenter, transform_inputs, augment_data_torch, add_random_hot_pixel
+from .utils import shift_offcenter, transform_inputs
 
 logger = logging.getLogger(__name__)
 
 
+def _plot_seqnum_histogram(selected_seqnums: np.ndarray, save_path: str) -> None:
+    """Plot histogram of SEQNUM counts in the selected training set.
+
+    Used to verify uniform coverage across perturbations (one bin per SEQNUM).
+    """
+    try:
+        import matplotlib
+
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+    except ImportError:
+        logger.warning("matplotlib not available; skipping SEQNUM histogram.")
+        return
+    fig, ax = plt.subplots(figsize=(10, 4))
+    # Exclude sentinel -1 if present
+    vals = selected_seqnums[selected_seqnums >= 0]
+    if vals.size == 0:
+        logger.warning("No valid SEQNUMs for histogram.")
+        return
+    ax.hist(vals, bins=min(100, max(50, int(np.ptp(vals)) + 1)), color="steelblue", edgecolor="white")
+    ax.set_xlabel("SEQNUM")
+    ax.set_ylabel("Sample count")
+    ax.set_title("SEQNUM distribution in training set (should be roughly uniform per bin)")
+    plt.tight_layout()
+    os.makedirs(os.path.dirname(save_path) or ".", exist_ok=True)
+    plt.savefig(save_path, dpi=150)
+    plt.close()
+    logger.info(f"SEQNUM histogram saved to {save_path}")
+
+
 class Donuts(Dataset):
     """PRETRAINING DATASET Batoid Sims.
 
@@ -291,12 +331,10 @@ def __init__(
         coral_filepath: str = "/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/coral/",
         coral_mode: bool = False,
         mask_mode: bool = False,
-        augment: bool = False,
-        augment_scale: float = 0.1,
-        augment_kmin: float = 0.1,
+        rotate: bool = False,
         **kwargs: Any,
     ) -> None:
-        """Load the simulated ImSim donuts and zernikes in a Pytorch Dataset.
+        r"""Load the simulated ImSim donuts and zernikes in a Pytorch Dataset.
 
         Parameters
         ----------
@@ -318,17 +356,29 @@ def __init__(
             Whether to enable coral mode for domain adaptation.
         mask_mode: bool, default=False
             Whether to use mask mode for zernike extraction.
-        augment: bool, default=False
-            Whether to apply frequency-based augmentation. Requires coral_mode=True.
-            When enabled, randomly augments 50% of samples by injecting high-frequency
-            structure from coral images into simulation images.
-        augment_scale: float, default=0.1
-            (Deprecated) Scaling factor parameter. The scale is now randomly sampled
-            uniformly between 0.1 and 1.0 for each augmentation. This parameter is
-            kept for backward compatibility but is not used.
-        augment_kmin: float, default=0.1
-            Minimum frequency magnitude to extract from coral images for augmentation.
+        rotate: bool, default=False
+            Whether to use rotated Zernike targets from the dataset when available.
+        train_fraction: float, optional
+            Fraction of stamps per SEQNUM (0.0-1.0). Always uses full coverage of all
+            unique SEQNUMs (perturbations); this controls how many stamps per SEQNUM.
+            For train, if None, uses DEFAULT_TRAIN_FRACTION. For val, if None, uses full data.
+        min_seqnums: int, default=10001
+            Unused; kept for API compatibility. Full coverage uses all unique SEQNUMs
+            found in the directory (number is logged).
+        rotation_direction: Optional[str], default=None
+            Required when rotate=True. Must be \"pos\" or \"neg\" to select
+            `zk_true_camera_pos` or `zk_true_camera_neg` from the dataset.
         """
+        train_fraction = kwargs.pop("train_fraction", None)
+        rotation_direction = kwargs.pop("rotation_direction", None)
+        _ = kwargs.pop("min_seqnums", 10001)  # Unused; kept for API compatibility
+        plot_seqnum_hist = kwargs.pop("plot_seqnum_hist", None)
+        # Optional: disable SEQNUM-balanced sampling and just use all images.
+        seqnum_balanced = kwargs.pop("seqnum_balanced", True)
+        # Ignore deprecated augmentation params for backward compatibility
+        kwargs.pop("augment", None)
+        kwargs.pop("augment_scale", None)
+        kwargs.pop("augment_kmin", None)
         self.settings = {
             "mode": mode,
             "transform": transform,
@@ -348,18 +398,17 @@ def __init__(
         logger.debug(f"Image directory: {self.image_dir}")
         self.coral_filepath = coral_filepath
         self.coral_mode = coral_mode
-        self.augment = augment
-        self.augment_scale = augment_scale
-        self.augment_kmin = augment_kmin
-
-        # Augmentation requires coral_mode to be enabled
-        if self.augment and not self.coral_mode:
-            logger.warning("augment=True requires coral_mode=True. Disabling augmentation.")
-            self.augment = False
+        self.rotate = rotate
+        # Enforce explicit rotation direction when rotate=True to avoid silent mismatches
+        if self.rotate and rotation_direction not in ("pos", "neg"):
+            raise ValueError(
+                f"rotation_direction must be 'pos' or 'neg' when rotate=True "
+                f"(got {rotation_direction!r})."
+            )
+        self.rotation_direction = rotation_direction
 
         if coral_mode:
             self.coral_image_files = []
-            # Use a temporary variable to avoid modifying the original path
             coral_data_path = coral_filepath
             if self.settings["mode"] == "train":
                 coral_data_path += "/train"
@@ -370,19 +419,103 @@ def __init__(
                     file_path = os.path.join(root, file)
                     self.coral_image_files.append(file_path)
 
-        if self.settings["mode"] == "train":
-            # self.image_files = self.image_files[: int(0.1 * len(self.image_files))]
-
-            self.image_files = self.image_files[: int(DEFAULT_TRAIN_FRACTION * len(self.image_files))]
+        self.used_indices = list(range(len(self.image_files)))
+        if not seqnum_balanced:
+            # Unbalanced path: optional subsample by train_fraction (no SEQNUM balancing).
+            frac = train_fraction if train_fraction is not None else DEFAULT_TRAIN_FRACTION
+            n_total = len(self.image_files)
+            n_keep = max(1, int(frac * n_total))
+            if n_keep < n_total:
+                chosen = np.random.choice(n_total, size=n_keep, replace=False)
+                self.used_indices = np.random.permutation(chosen).tolist()
+                logger.info(
+                    f"Subsampled to {len(self.used_indices)} / {n_total} samples (train_fraction={frac})."
+                )
+        else:
+            # SEQNUM-based balanced sampling: always full coverage of all perturbations.
+            # train_fraction controls stamps per SEQNUM (at least 1 per SEQNUM).
+            frac = train_fraction if train_fraction is not None else DEFAULT_TRAIN_FRACTION
+            # Store SEQNUM cache outside the data directory so it is never read as data.
+            _cache_dir = os.path.join(os.path.expanduser("~"), ".cache", "tarts", "seqnum")
+            os.makedirs(_cache_dir, exist_ok=True)
+            _path_hash = hashlib.sha256(os.path.abspath(self.image_dir).encode()).hexdigest()[:16]
+            cache_path = os.path.join(_cache_dir, f"{_path_hash}.pkl")
+            seqnum_per_index: Optional[List[int]] = None
+            if os.path.exists(cache_path):
+                try:
+                    with open(cache_path, "rb") as f:
+                        cache = pickle.load(f)
+                    path_to_seqnum = cache.get("path_to_seqnum")
+                    if path_to_seqnum is not None and set(path_to_seqnum.keys()) == set(self.image_files):
+                        seqnum_per_index = [path_to_seqnum[fp] for fp in self.image_files]
+                        logger.info(
+                            f"Loaded SEQNUM cache from {cache_path} ({len(seqnum_per_index)} entries)."
+                        )
+                    elif cache.get("image_files") == self.image_files:
+                        # Legacy format: same order required
+                        seqnum_per_index = cache.get("seqnum_per_index")
+                        if seqnum_per_index is not None:
+                            logger.info(
+                                f"Loaded SEQNUM cache from {cache_path} ({len(seqnum_per_index)} entries)."
+                            )
+                except (pickle.PickleError, OSError) as e:
+                    logger.warning(f"Could not load SEQNUM cache: {e}")
+            if seqnum_per_index is None:
+                seqnum_per_index = []
+                for i, fp in enumerate(tqdm(self.image_files, desc="Building SEQNUM index")):
+                    try:
+                        data = np.load(fp, allow_pickle=True)
+                        try:
+                            seqnum = int(data["SEQNUM"].item()) if "SEQNUM" in data else -1
+                        finally:
+                            if hasattr(data, "close"):
+                                data.close()
+                    except (KeyError, ValueError, OSError) as e:
+                        logger.warning(f"Failed to load SEQNUM from {fp}: {e}")
+                        seqnum = -1
+                    seqnum_per_index.append(seqnum)
+                try:
+                    path_to_seqnum = {fp: seqnum_per_index[i] for i, fp in enumerate(self.image_files)}
+                    with open(cache_path, "wb") as f:
+                        pickle.dump({"path_to_seqnum": path_to_seqnum}, f)
+                    logger.info(f"Saved SEQNUM cache to {cache_path}.")
+                except OSError as e:
+                    logger.warning(f"Could not save SEQNUM cache: {e}")
+            seqnum_to_indices: Dict[int, List[int]] = {}
+            for i, seqnum in enumerate(seqnum_per_index):
+                seqnum_to_indices.setdefault(seqnum, []).append(i)
+            num_unique = len(seqnum_to_indices) - (1 if -1 in seqnum_to_indices else 0)
+            logger.info(f"Found {num_unique} unique SEQNUMs in dataset.")
+            selected: List[Tuple[int, int]] = []  # (file_index, seqnum)
+            for seqnum, indices in seqnum_to_indices.items():
+                n_take = max(1, min(len(indices), int(frac * len(indices))))
+                chosen = np.random.choice(len(indices), size=n_take, replace=False)
+                for j in chosen:
+                    selected.append((indices[int(j)], seqnum))
+            np.random.shuffle(selected)
+            self.used_indices = [s[0] for s in selected]
+            selected_seqnums = np.array([s[1] for s in selected])
+            if len(self.used_indices) < num_unique:
+                logger.warning(
+                    f"Selected {len(self.used_indices)} samples < {num_unique} unique SEQNUMs; "
+                    "some perturbation bins may be unrepresented."
+                )
+            # Optional: plot SEQNUM distribution to verify uniform coverage
+            if plot_seqnum_hist and self.settings["mode"] == "train":
+                _plot_seqnum_histogram(
+                    selected_seqnums,
+                    plot_seqnum_hist if isinstance(plot_seqnum_hist, str) else "seqnum_histogram.png",
+                )
 
         if noll_zk is None:
             noll_zk = DEFAULT_NOLL_ZK
         self.noll_zk = np.array(noll_zk) - 4
         self.adjustment_factor = adjustment_factor
+        self.rotate = rotate
 
     def __len__(self) -> int:
         """Return length of this Dataset."""
-        return len(self.image_files)
+        return len(self.used_indices)
 
     def sample_coral(self) -> Dict[str, Any]:
         """Sample a coral image randomly."""
@@ -407,7 +540,13 @@ def sample_coral(self) -> Dict[str, Any]:
                 if img_file in corrupted_files:
                     continue
 
-                state = np.load(img_file, allow_pickle=True)
+                # Load coral file (avoid 'with' for NumPy 2.x / dict return)
+                loaded_state = np.load(img_file, allow_pickle=True)
+                try:
+                    state = {key: loaded_state[key].copy() for key in loaded_state.keys()}
+                finally:
+                    if hasattr(loaded_state, "close"):
+                        loaded_state.close()
                 break  # Successfully loaded, exit retry loop
             except (EOFError, IOError, OSError) as e:
                 # File is corrupted, truncated, or missing - delete it
@@ -514,30 +653,71 @@ def __getitem__(self, idx: int) -> Dict[str, torch.Tensor]:
                 obsID: the observation ID
                 objID: the object ID
         """
-        # get the image file
-        img_file = self.image_files[idx]
+        # get the image file (use index mapping for SEQNUM-balanced sampling)
+        img_file = self.image_files[self.used_indices[idx]]
 
+        # Load npz; avoid 'with' so we support np.load returning dict (e.g. NumPy 2.x / .npy)
         state = np.load(img_file, allow_pickle=True)
+        try:
+            # get the donut locations
+            fx, fy = (
+                torch.tensor(state["field_x"]) * np.pi / 180,
+                torch.tensor(state["field_y"]) * np.pi / 180,
+            )
 
-        # get the donut locations
-        fx, fy = (
-            torch.tensor(state["field_x"]) * np.pi / 180,
-            torch.tensor(state["field_y"]) * np.pi / 180,
-        )
-
-        # get the intra/extra flag
-        intra = torch.tensor(state["intra"]).int()
-
-        if self.mask_mode:
-            zernikes = torch.tensor(state["zk_true"])
-            zernikes = zernikes[:, 0]
-            zernikes = zernikes[None, :]
-        else:
-            zernikes = torch.tensor(state["zk_true"])[:, self.noll_zk]
-
-        band_tensor = torch.tensor(state["band"]).int()
-        band = band_tensor.item()
-        img = torch.tensor(state["image_aligned"])
+            # get the intra/extra flag
+            intra = torch.tensor(state["intra"]).int()
+
+            # Choose which Zernike target to use based on rotation settings
+            zk_true_key = None
+            is_rotated = False
+            if self.rotate:
+                # Prefer new explicit positive/negative rotation fields when available
+                if self.rotation_direction == "pos" and "zk_true_camera_pos" in state:
+                    zk_true_key = "zk_true_camera_pos"
+                    is_rotated = True
+                elif self.rotation_direction == "neg" and "zk_true_camera_neg" in state:
+                    zk_true_key = "zk_true_camera_neg"
+                    is_rotated = True
+                # Legacy single-rotation field
+                elif "zk_true_camera" in state:
+                    zk_true_key = "zk_true_camera"
+                    is_rotated = True
+                # Fall back to unrotated if nothing else is present
+                elif "zk_true" in state:
+                    zk_true_key = "zk_true"
+                    is_rotated = False
+                else:
+                    raise KeyError(
+                        f"No suitable rotated Zernikes found in npz file when rotate=True. File: {img_file}"
+                    )
+            else:
+                if "zk_true" in state:
+                    zk_true_key = "zk_true"
+                    is_rotated = False
+                else:
+                    raise KeyError(f"'zk_true' not found in npz file. File: {img_file}")
+
+            if self.mask_mode:
+                zernikes = torch.tensor(state[zk_true_key])
+                # zk_true_camera has 1 less dimension, so no need to index [:, 0]
+                if not is_rotated:
+                    zernikes = zernikes[:, 0]
+                zernikes = zernikes[None, :]
+            else:
+                zernikes = torch.tensor(state[zk_true_key])
+                # zk_true_camera has 1 less dimension, so indexing is different
+                if is_rotated:
+                    zernikes = zernikes[self.noll_zk]
+                else:
+                    zernikes = zernikes[:, self.noll_zk]
+
+            band_tensor = torch.tensor(state["band"]).int()
+            band = band_tensor.item()
+            img = torch.tensor(state["image_aligned"])
+        finally:
+            if hasattr(state, "close"):
+                state.close()
 
         # standardize all the inputs for the neural net
         if self.settings["transform"]:
@@ -582,7 +762,11 @@ def __getitem__(self, idx: int) -> Dict[str, torch.Tensor]:
             offset_amount = [0, 0]
             offset_vec = np.array([0.0, 0.0])
             offset_r = np.array(0.0)
-        zernikes = zernikes.float()[0, :]
+        # zk_true_camera is already 1D, so no need to index [0, :]
+        if zernikes.dim() > 1:
+            zernikes = zernikes.float()[0, :]
+        else:
+            zernikes = zernikes.float()
 
         output = {
             "image": img,
@@ -599,50 +783,6 @@ def __getitem__(self, idx: int) -> Dict[str, torch.Tensor]:
             coral_output = self.sample_coral()
             output.update(coral_output)
 
-            # Apply augmentations if enabled (only in training mode)
-            if self.augment and self.settings["mode"] == "train":
-                # Apply frequency-based augmentation 50% of the time
-                if torch.rand(1).item() > 0.5:
-                    try:
-                        # Randomly sample scale between 0.1 and 1.0 uniformly
-                        random_scale = torch.empty(1).uniform_(0.5, 1.0).item()
-                        # Apply augmentation (function handles device, shape, and dtype internally)
-                        augmented_img = augment_data_torch(
-                            img, coral_output["coral_image"], scale=random_scale, kmin=self.augment_kmin
-                        )
-
-                        # Ensure augmented image matches original image's dtype and device
-                        augmented_img = augmented_img.to(dtype=img.dtype, device=img.device)
-
-                        # Restore original shape if img had extra dimensions
-                        original_shape = img.shape
-                        if augmented_img.shape != original_shape:
-                            # Add back channel/batch dimensions if needed
-                            while len(augmented_img.shape) < len(original_shape):
-                                augmented_img = augmented_img.unsqueeze(0)
-
-                        output["image"] = augmented_img
-                    except Exception as e:
-                        # If augmentation fails, use original image and log warning
-                        logger.warning(f"Frequency augmentation failed, using original image: {e}")
-                        # output["image"] remains as img
-
-                # Apply hot pixel augmentation 10% of the time (independent of frequency augmentation)
-                # This can be applied to the original image or the frequency-augmented image
-                try:
-                    current_img = output.get("image", img)
-                    hot_pixel_img = add_random_hot_pixel(
-                        current_img, sigma=0.05, prob=0.1, min_scale=5.0, max_scale=20
-                    )
-                    # Ensure dtype and device match (function preserves device, but ensure dtype consistency)
-                    if hot_pixel_img.dtype != img.dtype:
-                        hot_pixel_img = hot_pixel_img.to(dtype=img.dtype)
-                    output["image"] = hot_pixel_img
-                except Exception as e:
-                    # If hot pixel augmentation fails, keep current image and log warning
-                    logger.warning(f"Hot pixel augmentation failed: {e}")
-                    # output["image"] remains as current_img
-
         return output
 
 
@@ -711,11 +851,13 @@ def __init__(
         train=True,
         data_dir="/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator/",
         alpha=1e-3,
+        rotate=False,
+        rotation_direction: str | None = None,
         return_true=False,
         coral_mode=False,
         coral_filepath="/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator_real/",
     ):
-        """Initialize the zernikeDataset.
+        r"""Initialize the zernikeDataset.
 
         Parameters
         ----------
@@ -729,6 +871,12 @@ def __init__(
             Parameter used for adjusting Zernike coefficients during processing.
         return_true : bool, optional, default=False
             Whether to return the true Zernike coefficients or estimated coefficients.
+        rotate : bool, optional, default=False
+            Whether to use rotated Zernikes in the camera frame when available.
+        rotation_direction : {"pos", "neg"}, optional
+            When rotate=True, selects which camera-frame Zernikes to use:
+            - \"pos\" → zk_true_camera_pos (positive RTP rotation)
+            - \"neg\" → zk_true_camera_neg (negative RTP rotation)
         coral_mode : bool, optional, default=False
             Whether to enable coral mode for real data sampling.
         coral_filepath : str, optional
@@ -742,7 +890,7 @@ def __init__(
         if train:
             self.image_dir = data_dir + "/train"
         else:
-            self.image_dir = data_dir + "/train"
+            self.image_dir = data_dir + "/val"
 
         self.filename = []
         for root, _, files in os.walk(self.image_dir):
@@ -750,11 +898,8 @@ def __init__(
                 file_path = os.path.join(root, file)
                 self.filename.append(file_path)
 
-        if train:
-            self.filename = self.filename[: int(0.8 * len(self.filename))]
-        else:
-            self.filename = self.filename[int(0.8 * len(self.filename)) :]
-
+        # Use all files in the directory (no additional train/val split)
+        # The split is already handled by using different directories (train/val)
         self.num_samples = len(self.filename)
         self.alpha = alpha
         self.return_true = return_true
@@ -778,6 +923,14 @@ def __init__(
         # Always use CPU in dataset - PyTorch Lightning handles GPU transfer
         # This avoids CUDA reinitialization issues with num_workers > 0
         self.device = torch.device("cpu")
+        self.rotate = rotate
+        # Enforce explicit rotation direction when using rotated Zernikes
+        if self.rotate and rotation_direction not in ("pos", "neg"):
+            raise ValueError(
+                f"rotation_direction must be 'pos' or 'neg' when rotate=True "
+                f"(got {rotation_direction!r})."
+            )
+        self.rotation_direction = rotation_direction
 
     def __len__(self) -> int:
         """Return the number of samples in the dataset."""
@@ -958,8 +1111,28 @@ def __getitem__(self, idx: int) -> Union[
                 "header": json.loads(str(npz_data["header_json"])),
             }
             # Add conditional fields if they exist
-            if "zk_true" in npz_data:
-                loaded_data["zk_true"] = torch.from_numpy(npz_data["zk_true"])
+            # When rotate=True, use camera-frame Zernikes; enforce that the correct
+            # explicit positive/negative field is present for a hard guarantee.
+            if self.rotate:
+                if self.rotation_direction == "neg":
+                    if "zk_true_camera_neg" not in npz_data:
+                        raise KeyError(
+                            "Expected 'zk_true_camera_neg' in npz file when "
+                            "rotate=True and rotation_direction='neg'."
+                        )
+                    loaded_data["zk_true"] = torch.from_numpy(npz_data["zk_true_camera_neg"])
+                elif self.rotation_direction == "pos":
+                    if "zk_true_camera_pos" not in npz_data:
+                        raise KeyError(
+                            "Expected 'zk_true_camera_pos' in npz file when "
+                            "rotate=True and rotation_direction='pos'."
+                        )
+                    loaded_data["zk_true"] = torch.from_numpy(npz_data["zk_true_camera_pos"])
+            else:
+                if "zk_true" in npz_data:
+                    loaded_data["zk_true"] = torch.from_numpy(npz_data["zk_true"])
+                else:
+                    raise KeyError("'zk_true' not found in npz file when rotate=False")
         except (IOError, OSError, RuntimeError, KeyError) as e:
             logger.error(
                 f"Error loading file {self.filename[idx] if idx < len(self.filename) else 'unknown'}: {e}"
diff --git a/python/tarts/dataset_params.yaml b/python/tarts/dataset_params.yaml
index be35512..9dd7cc2 100644
--- a/python/tarts/dataset_params.yaml
+++ b/python/tarts/dataset_params.yaml
@@ -1,20 +1,23 @@
 version: 0
 seed: 42  # Random seed for reproducibility across all training scripts
+# Must be > 10k to ensure full coverage of all perturbation bins (SEQNUMs) when sampling.
+# Each SEQNUM = unique Zernike perturbation; dataset has 10k such perturbations.
+min_seqnums: 10001
 adjustment_WaveNet: 10
 adjustment_AlignNet: 120
 refinements: 1
 CROP_SIZE : 160
 deg_per_pix : 0.00005555555
 mm_pix : 0.01
-alpha: 1
+alpha: 0.5
 max_seq_len: 200
-butler_repo_path: 'butler'
+butler_repo_path: 'butler_training_fix'
 
 # Model architecture parameters
-cnn_model: 'mobilenetv4_conv_small'  # Options: mobilenetv4_conv_small, resnet18, resnet34, resnet50
-finetune_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain'
-finetune_filepath_train: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain/train/'
-finetune_filepath_val: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain/val/'
+cnn_model: 'mobilenetv4_conv_small'  # Options: mobilenetv4_conv_small, mobilenetv4_conv_medium, mobilenetv4_conv_large, mobilenetv4_hybrid_large, resnet18, resnet34, resnet50
+finetune_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain_42'
+finetune_filepath_train: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain_42/train/'
+finetune_filepath_val: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/simulation_pretrain_42/val/'
 coral_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/coral/'
 coral_filepath_train: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/coral/train/'
 coral_filepath_val: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/coral/val/'
@@ -23,19 +26,28 @@ coral_filepath_usdf_val: '~/val'
 aggregator_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator/'
 aggregator_train_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator/train/'
 aggregator_val_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator/val/'
-fullframe_train_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_plane_ml_opsim_0925/train/'
-fullframe_val_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_plane_ml_opsim_0925/val/'
-fullframe_test_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_plane_ml_opsim_0925/test/'
+fullframe_train_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_frame_final/train/'
+fullframe_val_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_frame_final/val/'
+fullframe_test_filepath: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/full_frame_final/test/'
 aggregator_real_filepath: "/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator_real"
 aggregator_real_train_filepath: "/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator_real/train"
 aggregator_real_val_filepath: "/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/aggregator_real/val"
+batoid_filepath_train: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/batoid/train/'
+batoid_filepath_val: '/media/peterma/mnt2/peterma/research/LSST_FULL_FRAME/batoid/val/'
+# Optional: path to DOF table (.npy shape (N, n_dof) or .csv). If set, each sample uses a row; if unset, all samples use zero DOF.
+# batoid_dof_table: '/path/to/dof_table.npy'
 
 # AggregatorNet model parameters
 aggregator_model:
-  d_model: 128
+  d_model: 32
   nhead: 4
   num_layers: 6
-  dim_feedforward: 512
+  dim_feedforward: 32
+# aggregator_model:
+#   d_model: 128
+#   nhead: 4
+#   num_layers: 6
+#   dim_feedforward: 512
 
 noll_zk:
 - 4
diff --git a/python/tarts/focus_offset.yaml b/python/tarts/focus_offset.yaml
new file mode 100644
index 0000000..f718758
--- /dev/null
+++ b/python/tarts/focus_offset.yaml
@@ -0,0 +1,89 @@
+z4_peak_pattern:
+  description: "Detector-polarity branch model for Z4 signed error (pred - truth)."
+  source:
+    by_band_peaks:
+      u: [0.028133]
+      g: [0.013092]
+      r: [-0.04707, 0.048187]
+      i: [-0.12729, 0.13342]
+      z: [-0.36794, 0.37407]
+      y: [-0.9934645787151222, -0.9602425377593652, 0.8337476738515093, 0.8171366533736306]
+    # Peaks keyed by raft half (intra/extra pairs share the same branch: 191/192, 195/196, …).
+    by_detector_peaks:
+      "191": [0.8337476738515093, 0.37407, 0.13342, 0.043173]
+      "192": [0.8337476738515093, 0.37407, 0.13342, 0.043173]
+      "195": [-0.9934645787151222, -0.37295, -0.12729, -0.042057]
+      "196": [-0.9934645787151222, -0.37295, -0.12729, -0.042057]
+      "199": [-0.9602425377593652, -0.36794, -0.12227, -0.016989]
+      "200": [-0.9602425377593652, -0.36794, -0.12227, -0.016989]
+      "203": [0.8171366533736306, 0.36905, 0.1284, 0.048187]
+      "204": [0.8171366533736306, 0.36905, 0.1284, 0.048187]
+
+  detector_polarity_groups:
+    positive_branch: [191, 192, 203, 204]
+    negative_branch: [195, 196, 199, 200]
+
+  branch_amplitude_by_band:
+    u: 0.028133
+    g: 0.013092
+    r: 0.048187
+    i: 0.13342
+    z: 0.37407
+    y: 0.8969951058054373
+
+  # Applied in NeuralActiveOpticsSys: same offset for intra/extra chip of each half (191↔192, etc.).
+  branch_prediction_by_band_detector:
+    u:
+      "191": 0.028133
+      "192": 0.028133
+      "195": -0.028133
+      "196": -0.028133
+      "199": -0.028133
+      "200": -0.028133
+      "203": 0.028133
+      "204": 0.028133
+    g:
+      "191": 0.013092
+      "192": 0.013092
+      "195": -0.013092
+      "196": -0.013092
+      "199": -0.013092
+      "200": -0.013092
+      "203": 0.013092
+      "204": 0.013092
+    r:
+      "191": 0.048187
+      "192": 0.048187
+      "195": -0.048187
+      "196": -0.048187
+      "199": -0.048187
+      "200": -0.048187
+      "203": 0.048187
+      "204": 0.048187
+    i:
+      "191": 0.13342
+      "192": 0.13342
+      "195": -0.13342
+      "196": -0.13342
+      "199": -0.13342
+      "200": -0.13342
+      "203": 0.13342
+      "204": 0.13342
+    z:
+      "191": 0.37407
+      "192": 0.37407
+      "195": -0.37407
+      "196": -0.37407
+      "199": -0.37407
+      "200": -0.37407
+      "203": 0.37407
+      "204": 0.37407
+    y:
+      "191": 0.8337476738515093
+      "192": 0.8337476738515093
+      "195": -0.9934645787151222
+      "196": -0.9934645787151222
+      "199": -0.9602425377593652
+      "200": -0.9602425377593652
+      "203": 0.8171366533736306
+      "204": 0.8171366533736306
diff --git a/python/tarts/focus_offset_danish.yaml b/python/tarts/focus_offset_danish.yaml
new file mode 100644
index 0000000..f718758
--- /dev/null
+++ b/python/tarts/focus_offset_danish.yaml
@@ -0,0 +1,89 @@
+z4_peak_pattern:
+  description: "Detector-polarity branch model for Z4 signed error (pred - truth)."
+  source:
+    by_band_peaks:
+      u: [0.028133]
+      g: [0.013092]
+      r: [-0.04707, 0.048187]
+      i: [-0.12729, 0.13342]
+      z: [-0.36794, 0.37407]
+      y: [-0.9934645787151222, -0.9602425377593652, 0.8337476738515093, 0.8171366533736306]
+    # Peaks keyed by raft half (intra/extra pairs share the same branch: 191/192, 195/196, …).
+    by_detector_peaks:
+      "191": [0.8337476738515093, 0.37407, 0.13342, 0.043173]
+      "192": [0.8337476738515093, 0.37407, 0.13342, 0.043173]
+      "195": [-0.9934645787151222, -0.37295, -0.12729, -0.042057]
+      "196": [-0.9934645787151222, -0.37295, -0.12729, -0.042057]
+      "199": [-0.9602425377593652, -0.36794, -0.12227, -0.016989]
+      "200": [-0.9602425377593652, -0.36794, -0.12227, -0.016989]
+      "203": [0.8171366533736306, 0.36905, 0.1284, 0.048187]
+      "204": [0.8171366533736306, 0.36905, 0.1284, 0.048187]
+
+  detector_polarity_groups:
+    positive_branch: [191, 192, 203, 204]
+    negative_branch: [195, 196, 199, 200]
+
+  branch_amplitude_by_band:
+    u: 0.028133
+    g: 0.013092
+    r: 0.048187
+    i: 0.13342
+    z: 0.37407
+    y: 0.8969951058054373
+
+  # Applied in NeuralActiveOpticsSys: same offset for intra/extra chip of each half (191↔192, etc.).
+  branch_prediction_by_band_detector:
+    u:
+      "191": 0.028133
+      "192": 0.028133
+      "195": -0.028133
+      "196": -0.028133
+      "199": -0.028133
+      "200": -0.028133
+      "203": 0.028133
+      "204": 0.028133
+    g:
+      "191": 0.013092
+      "192": 0.013092
+      "195": -0.013092
+      "196": -0.013092
+      "199": -0.013092
+      "200": -0.013092
+      "203": 0.013092
+      "204": 0.013092
+    r:
+      "191": 0.048187
+      "192": 0.048187
+      "195": -0.048187
+      "196": -0.048187
+      "199": -0.048187
+      "200": -0.048187
+      "203": 0.048187
+      "204": 0.048187
+    i:
+      "191": 0.13342
+      "192": 0.13342
+      "195": -0.13342
+      "196": -0.13342
+      "199": -0.13342
+      "200": -0.13342
+      "203": 0.13342
+      "204": 0.13342
+    z:
+      "191": 0.37407
+      "192": 0.37407
+      "195": -0.37407
+      "196": -0.37407
+      "199": -0.37407
+      "200": -0.37407
+      "203": 0.37407
+      "204": 0.37407
+    y:
+      "191": 0.8337476738515093
+      "192": 0.8337476738515093
+      "195": -0.9934645787151222
+      "196": -0.9934645787151222
+      "199": -0.9602425377593652
+      "200": -0.9602425377593652
+      "203": 0.8171366533736306
+      "204": 0.8171366533736306
diff --git a/python/tarts/intrinsics_corners.parquet b/python/tarts/intrinsics_corners.parquet
new file mode 100644
index 0000000..333c8b9
Binary files /dev/null and b/python/tarts/intrinsics_corners.parquet differ
diff --git a/python/tarts/lightning_wavenet.py b/python/tarts/lightning_wavenet.py
index 70b2997..3030386 100644
--- a/python/tarts/lightning_wavenet.py
+++ b/python/tarts/lightning_wavenet.py
@@ -134,6 +134,9 @@ def _build_loader(self, mode: str, shuffle: bool = False, drop_last: bool = True
             pin_memory=self.hparams.pin_memory,
             shuffle=shuffle,
             drop_last=drop_last,
+            prefetch_factor=(
+                4 if self.hparams.num_workers > 0 else None
+            ),  # Prefetch more batches to prevent GPU waiting
         )
 
     def train_dataloader(self) -> DataLoader:
diff --git a/python/tarts/lightning_wavenet_coral.py b/python/tarts/lightning_wavenet_coral.py
index 472ec20..c0967d9 100644
--- a/python/tarts/lightning_wavenet_coral.py
+++ b/python/tarts/lightning_wavenet_coral.py
@@ -6,7 +6,8 @@
 
 # Standard library imports
 import logging
-from typing import Any, Dict, Tuple
+import numpy as np
+from typing import Any, Dict, Optional, Tuple
 
 # Third-party imports
 import pytorch_lightning as pl
@@ -58,6 +59,9 @@ def __init__(
         tradeoff_scale: float = 0.001,
         threshold: float = 0.9,
         dare_gram_weight: float = 1.0,
+        exp_rise_x1: float = 0.15,
+        exp_rise_x2: float = 0.16,
+        backbone_lr: Optional[float] = None,
     ) -> None:
         """Create the WaveNet with DARE-GRAM domain adaptation.
 
@@ -92,6 +96,15 @@ def __init__(
         dare_gram_weight: float, default=1.0
             Overall weight to scale the DARE-GRAM loss relative to regression loss.
             Higher values prioritize domain adaptation over regression accuracy.
+        exp_rise_x1: float, default=0.15
+            Lower threshold for exponential rise function in DARE-GRAM scaling.
+            When mRSSE <= x1, scale_loss = 1.0 (full DARE-GRAM influence).
+        exp_rise_x2: float, default=0.16
+            Upper threshold for exponential rise function in DARE-GRAM scaling.
+            When mRSSE >= x2, scale_loss = 0.0 (no DARE-GRAM influence).
+            Exponential transition occurs between x1 and x2.
+        backbone_lr: float, optional
+            If set, use this LR for the CNN backbone and lr for the predictor (e.g. 1e-5 for ResNet101).
         """
         super().__init__()
         self.save_hyperparameters()
@@ -108,6 +121,16 @@ def __init__(
         self.inputShape = DEFAULT_INPUT_SHAPE
         self.val_mRSSE: torch.Tensor | None = None
 
+    def _dare_gram_enabled(self) -> bool:
+        """Return True if DARE-GRAM can contribute non-zero loss.
+
+        This is used to skip expensive DARE-GRAM computations (SVD/pinv, extra forward)
+        when the effective weight is zero.
+        """
+        return float(self.hparams.dare_gram_weight) != 0.0 and (
+            float(self.hparams.tradeoff_angle) != 0.0 or float(self.hparams.tradeoff_scale) != 0.0
+        )
+
     def dare_gram_loss(self, features_source: torch.Tensor, features_target: torch.Tensor) -> torch.Tensor:
         """Compute DARE-GRAM loss between source and target features.
 
@@ -235,9 +258,9 @@ def predict_step(self, batch: Dict[str, Any], batch_idx: int) -> Tuple[torch.Ten
         return zk_pred, zk_true
 
     def exp_rise_flipped(self, loss, a=6.0):
-        """Exponentially rises from 0 at loss=0.13 to 1 at loss=0.16.
+        """Exponentially rises from 0 at loss=x1 to 1 at loss=x2.
 
-        Then flattens at 1 above 0.16 and 0 below 0.13.
+        Then flattens at 1 above x2 and 0 below x1.
 
         Parameters
         ----------
@@ -246,28 +269,31 @@ def exp_rise_flipped(self, loss, a=6.0):
         a: float
             Controls steepness (larger = sharper rise).
         """
-        # Convert loss to tensor and ensure float32 dtype
-        loss = torch.as_tensor(loss, dtype=torch.float32, device=self.device_val)
-
-        # Ensure scalar values are float32 tensors
-        x1 = torch.tensor(0.13, dtype=torch.float32, device=loss.device)
-        x2 = torch.tensor(0.16, dtype=torch.float32, device=loss.device)
-        a = torch.tensor(a, dtype=torch.float32, device=loss.device)
-
-        f = torch.zeros_like(loss, dtype=torch.float32)
-
-        # Region 1: loss <= 0.13 → f = 0
-        f[loss <= x1] = 0.0
-
-        # Region 2: 0.13 < loss < 0.16 → exponential rise
-        mask = (loss > x1) & (loss < x2)
-        t = (loss[mask] - x1) / (x2 - x1)
-        f[mask] = (1 - torch.exp(-a * t)) / (1 - torch.exp(-a))
+        # Convert loss to scalar
+        loss_val = loss.item() if isinstance(loss, torch.Tensor) else float(loss)
+
+        # Get hyperparameters (must be provided, will error if not)
+        x1 = float(self.hparams.exp_rise_x1)
+        x2 = float(self.hparams.exp_rise_x2)
+        a = float(a)
+
+        # Region 1: loss <= x1 → f = 0
+        if loss_val <= x1:
+            f = 0.0
+        # Region 3: loss >= x2 → f = 1
+        elif loss_val >= x2:
+            f = 1.0
+        # Region 2: x1 < loss < x2 → exponential rise
+        else:
+            t = (loss_val - x1) / (x2 - x1)
+            f = (1 - np.exp(-a * t)) / (1 - np.exp(-a))
 
-        # Region 3: loss >= 0.16 → f = 1
-        f[loss >= x2] = 1.0
+        # Flip: high f (when loss is high) becomes low scale (when loss is high)
+        # We want: low loss (good) → high scale (more DARE-GRAM), high loss (bad) → low scale (less DARE-GRAM)
         f = -f + 1
-        return f
+
+        # Convert back to tensor on correct device
+        return torch.tensor(f, dtype=torch.float32, device=self.device_val)
 
     def calc_losses(
         self,
@@ -275,7 +301,7 @@ def calc_losses(
         batch_idx: int,
         use_coral: bool = False,
         add_dare_gram_to_loss: bool = True,
-    ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
+    ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
         """Predict Zernikes and calculate losses with optional DARE-GRAM.
 
         Parameters
@@ -293,8 +319,9 @@ def calc_losses(
         Returns
         -------
         tuple
-            (total_loss, mRSSE, dare_gram_loss)
+            (total_loss, mRSSE, dare_gram_loss, dare_gram_scale)
             dare_gram_loss is 0 if use_coral=False.
+            dare_gram_scale is the scaling factor from exp_rise_flipped function.
         """
         zk_pred, zk_true = self.predict_step(batch, batch_idx)
 
@@ -310,14 +337,18 @@ def calc_losses(
         regression_loss = sse.mean() + self.hparams.alpha * A.square().sum()
         mRSSE = torch.sqrt(sse).mean()
 
+        # Compute DARE-GRAM scaling first. If it is 0, DARE-GRAM is effectively disabled for
+        # this step and we can skip its expensive computation entirely (extra forward + SVD/pinv).
+        scale_loss = self.exp_rise_flipped(self.val_mRSSE if self.val_mRSSE is not None else mRSSE)
+        scale_is_zero = float(scale_loss.item()) == 0.0
+
         # DARE-GRAM loss if coral data is available
+        # Skip computation entirely if the effective DARE-GRAM weight is 0.0 for faster training
         dare_gram_loss = torch.tensor(0.0, device=self.device_val)
-        if use_coral and "coral_image" in batch:
+        if use_coral and "coral_image" in batch and self._dare_gram_enabled() and not scale_is_zero:
             try:
-                # Forward pass on source data to get features
-                zk_pred_s, _ = self.predict_step(batch, batch_idx)
-
-                # Forward pass on target/coral data
+                # Features from the source forward were already produced by predict_step() above
+                # (WaveNet stores them on self.wavenet.predictor_features).
                 coral_img = batch["coral_image"]
                 coral_fx = batch["coral_field_x"]
                 coral_fy = batch["coral_field_y"]
@@ -345,19 +376,19 @@ def calc_losses(
 
         # Add DARE-GRAM to loss only if requested
         if add_dare_gram_to_loss:
-            scale_loss = self.exp_rise_flipped(self.val_mRSSE if self.val_mRSSE is not None else mRSSE)
             total_loss = regression_loss + self.hparams.dare_gram_weight * scale_loss * dare_gram_loss
         else:
             total_loss = regression_loss
 
-        return total_loss, mRSSE, dare_gram_loss
+        return total_loss, mRSSE, dare_gram_loss, scale_loss
 
     def training_step(self, batch: Dict[str, Any], batch_idx: int) -> torch.Tensor:
         """Execute training step on a batch."""
-        loss, mRSSE, dare_gram_loss = self.calc_losses(batch, batch_idx, use_coral=True)
+        loss, mRSSE, dare_gram_loss, dare_gram_scale = self.calc_losses(batch, batch_idx, use_coral=True)
         self.log("train_loss", loss, sync_dist=True, prog_bar=True)
         self.log("train_mRSSE", mRSSE, sync_dist=True)
         self.log("train_dare_gram_loss", dare_gram_loss, sync_dist=True)
+        self.log("train_dare_gram_scale", dare_gram_scale, sync_dist=True)
 
         return loss
 
@@ -365,27 +396,51 @@ def validation_step(self, batch: Dict[str, Any], batch_idx: int) -> torch.Tensor
         """Execute validation step on a batch."""
         # Compute DARE-GRAM for logging but don't add it to validation loss
         # This allows monitoring domain adaptation without affecting validation metrics
-        loss, mRSSE, dare_gram_loss = self.calc_losses(
+        loss, mRSSE, dare_gram_loss, dare_gram_scale = self.calc_losses(
             batch, batch_idx, use_coral=True, add_dare_gram_to_loss=False
         )
         self.log("val_loss", loss, sync_dist=True, prog_bar=True)
         self.log("val_mRSSE", mRSSE, sync_dist=True)
         self.log("val_dare_gram_loss", dare_gram_loss, sync_dist=True)
-        self.val_mRSSE = mRSSE.clone().detach()  # Store a copy to avoid tensor reference issues
+        self.log("val_dare_gram_scale", dare_gram_scale, sync_dist=True)
         return loss
 
+    def on_validation_epoch_end(self):
+        """Update val_mRSSE with epoch-averaged value after validation completes."""
+        # Get the epoch-averaged val_mRSSE from logged metrics
+        if "val_mRSSE" in self.trainer.callback_metrics:
+            epoch_avg_val_mRSSE = self.trainer.callback_metrics["val_mRSSE"]
+            self.val_mRSSE = epoch_avg_val_mRSSE.clone().detach()
+
     def configure_optimizers(self) -> Any:
-        """Configure the optimizer."""
-        optimizer = torch.optim.AdamW(
-            self.parameters(), lr=self.hparams.lr, weight_decay=self.hparams.weight_decay
-        )
+        """Configure the optimizer. Uses separate LRs for backbone vs predictor if backbone_lr is set."""
+        backbone_lr = getattr(self.hparams, "backbone_lr", None)
+        if backbone_lr is not None:
+            param_groups = [
+                {"params": self.wavenet.cnn.parameters(), "lr": float(backbone_lr)},
+                {"params": self.wavenet.predictor.parameters(), "lr": self.hparams.lr},
+            ]
+            optimizer = torch.optim.AdamW(param_groups, weight_decay=self.hparams.weight_decay)
+        else:
+            optimizer = torch.optim.AdamW(
+                self.parameters(), lr=self.hparams.lr, weight_decay=self.hparams.weight_decay
+            )
 
         if self.hparams.lr_schedule:
             return {
                 "optimizer": optimizer,
                 "lr_scheduler": {
-                    "scheduler": ReduceLROnPlateau(optimizer),
-                    "monitor": "val_loss",
+                    "scheduler": ReduceLROnPlateau(
+                        optimizer,
+                        mode="min",
+                        factor=0.5,  # Reduce LR by half when plateau
+                        patience=5,  # Wait 5 validation checks (epochs) before reducing
+                        threshold=1e-3,  # Minimum change to qualify as an improvement
+                        threshold_mode="abs",  # Use absolute threshold
+                        min_lr=1e-7,  # Minimum learning rate
+                        verbose=True,
+                    ),
+                    "monitor": "val_mRSSE",  # Monitor val_mRSSE instead of val_loss for better stability
                     "frequency": 1,
                 },
             }
diff --git a/python/tarts/wavenet.py b/python/tarts/wavenet.py
index 787b48c..59d954c 100644
--- a/python/tarts/wavenet.py
+++ b/python/tarts/wavenet.py
@@ -1,5 +1,8 @@
 """Neural network to predict zernike coefficients from donut images and positions."""
 
+# Configure Hugging Face cache (``HF_HOME`` / ``TARTS_HF_HOME``) before timm pulls huggingface_hub.
+from . import _runtime_env  # noqa: F401
+
 # Third-party imports
 import timm
 import torch
@@ -20,7 +23,7 @@ def __init__(
         n_predictor_layers: tuple = (256,),
         n_zernikes: int = 25,
         device="cuda",
-        pretrained: bool = True,
+        pretrained: bool = False,
     ) -> None:
         """Create the WaveNet.
 
@@ -39,8 +42,8 @@ def __init__(
             Number of Zernike coefficients to predict.
         device: str, default='cuda'
             Device to run the model on ('cuda' or 'cpu').
-        pretrained: bool, default=True
-            Whether to use pre-trained CNN weights. Set to False to avoid downloading weights.
+        pretrained: bool, default=False
+            Whether to use pre-trained CNN weights. Set to True only when downloads are acceptable.
         """
         super().__init__()
         # Check CUDA availability and convert device string to torch.device
@@ -171,8 +174,14 @@ def _reshape_image(self, image: torch.Tensor) -> torch.Tensor:
 
         # Get the number of input channels required by the CNN
         if self.is_timm_model:
-            # timm models (MobileNet, etc.)
-            n_channels = self.cnn.conv_stem.in_channels
+            # timm: MobileNet/ConvNeXt use conv_stem, ResNet uses conv1
+            first_conv = getattr(self.cnn, "conv_stem", None) or getattr(self.cnn, "conv1", None)
+            if first_conv is None:
+                msg = (
+                    f"WaveNet CNN backbone {type(self.cnn).__name__} has no " "'conv_stem' or 'conv1' layer."
+                )
+                raise AttributeError(msg)
+            n_channels = first_conv.in_channels
         else:
             # torchvision models (ResNet, etc.)
             n_channels = self.cnn.conv1.in_channels
diff --git a/setup.py b/setup.py
index dc28735..f247bf3 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@
 
 setup(
     name="tarts",
-    version="0.0.0a7",
+    version="4.2",
     packages=find_packages(),
     install_requires=[],
     author="Peter Xiangyuan Ma",