diff --git a/.gitignore b/.gitignore
index dd1a8c9..255aae8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,3 +32,4 @@ _skbuild/
 __pycache__/
 *.pyc
 logs/
+scripts/
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f074542..205cfcc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -147,7 +147,6 @@ if(DEEPITY_ENABLE_CUDA)
         enable_language(CUDA)
         set(CMAKE_CUDA_STANDARD 17)
         set(CMAKE_CUDA_STANDARD_REQUIRED ON)
-        set(CMAKE_CUDA_ARCHITECTURES 86)
         message(STATUS "CUDA support enabled.")
     else()
         message(STATUS "CUDA toolkit not found. Building CPU-only.")
@@ -328,14 +327,12 @@ set_target_properties(ActivationBenchmark PROPERTIES
     RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin
 )
 
-add_executable(CUDAIm2ColVerify tests/tCUDAIm2ColVerify.cpp)
-target_link_libraries(CUDAIm2ColVerify PRIVATE Deepity)
-set_target_properties(CUDAIm2ColVerify PROPERTIES
+add_executable(CUDALaunchErrorVerify tests/tCUDALaunchErrorVerify.cpp)
+target_link_libraries(CUDALaunchErrorVerify PRIVATE Deepity)
+set_target_properties(CUDALaunchErrorVerify PROPERTIES
     RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin
 )
 
-add_test(NAME CUDAIm2ColVerify COMMAND CUDAIm2ColVerify)
-
 add_test(NAME ActivationBenchmark COMMAND ActivationBenchmark)
 
 # --- Compiler flags -------------------------------------------------
diff --git a/CMakePresets.json b/CMakePresets.json
index c9f7e6c..5cedd29 100644
--- a/CMakePresets.json
+++ b/CMakePresets.json
@@ -51,7 +51,7 @@
                 "CMAKE_BUILD_TYPE": "Release",
                 "Python_EXECUTABLE": "${sourceDir}/.venv/bin/python3.12",
                 "DEEPITY_ENABLE_CUDA": "ON",
-                "CMAKE_CUDA_ARCHITECTURES": "86",
+                "CMAKE_CUDA_ARCHITECTURES": "80",
                 "CMAKE_CUDA_STANDARD": "17"
             }
         }
diff --git a/deepity_build/reporting/__pycache__/__init__.cpython-312.pyc b/deepity_build/reporting/__pycache__/__init__.cpython-312.pyc
index 705866f..640c1a8 100644
Binary files a/deepity_build/reporting/__pycache__/__init__.cpython-312.pyc and b/deepity_build/reporting/__pycache__/__init__.cpython-312.pyc differ
diff --git a/deepity_build/reporting/__pycache__/__init__.cpython-314.pyc b/deepity_build/reporting/__pycache__/__init__.cpython-314.pyc
deleted file mode 100644
index b70394a..0000000
Binary files a/deepity_build/reporting/__pycache__/__init__.cpython-314.pyc and /dev/null differ
diff --git a/deepity_build/reporting/__pycache__/base.cpython-312.pyc b/deepity_build/reporting/__pycache__/base.cpython-312.pyc
index 92f9220..77fa6ec 100644
Binary files a/deepity_build/reporting/__pycache__/base.cpython-312.pyc and b/deepity_build/reporting/__pycache__/base.cpython-312.pyc differ
diff --git a/deepity_build/reporting/__pycache__/base.cpython-314.pyc b/deepity_build/reporting/__pycache__/base.cpython-314.pyc
deleted file mode 100644
index 5859bf8..0000000
Binary files a/deepity_build/reporting/__pycache__/base.cpython-314.pyc and /dev/null differ
diff --git a/deepity_build/reporting/__pycache__/rich_reporter.cpython-312.pyc b/deepity_build/reporting/__pycache__/rich_reporter.cpython-312.pyc
index 3efd3f3..a16f640 100644
Binary files a/deepity_build/reporting/__pycache__/rich_reporter.cpython-312.pyc and b/deepity_build/reporting/__pycache__/rich_reporter.cpython-312.pyc differ
diff --git a/deepity_build/reporting/__pycache__/rich_reporter.cpython-314.pyc b/deepity_build/reporting/__pycache__/rich_reporter.cpython-314.pyc
deleted file mode 100644
index c965ce9..0000000
Binary files a/deepity_build/reporting/__pycache__/rich_reporter.cpython-314.pyc and /dev/null differ
diff --git a/include/deepity/backend/CUDABackend.h b/include/deepity/backend/CUDABackend.h
index 03e08a6..a8ad2a5 100644
--- a/include/deepity/backend/CUDABackend.h
+++ b/include/deepity/backend/CUDABackend.h
@@ -91,7 +91,8 @@ namespace Deep
         cudaGraphExec_t graphExec = nullptr;
 #endif
         bool hasGraph = false;
-        void *workspace = nullptr;
         float *onesVector = nullptr;
+        size_t onesCapacity = 0; // <--- Add this line
+        float *workspace = nullptr;
     };
 }
\ No newline at end of file
diff --git a/include/deepity/utils/ActivationType.h b/include/deepity/utils/ActivationType.h
new file mode 100644
index 0000000..6244e12
--- /dev/null
+++ b/include/deepity/utils/ActivationType.h
@@ -0,0 +1,35 @@
+#pragma once
+#include <cstdint>
+
+/**
+ * @file ActivationType.h
+ * @brief Just the ActivationType enum, deliberately split out of
+ * Activations.h -- zero dependency on <immintrin.h>/<sleef.h>/<omp.h>.
+ *
+ * Why this exists: IComputeBackend.h's methods take ActivationType as a
+ * plain enum parameter and never touch the actual CPU SIMD function
+ * implementations (ActivationFn/DerivativeFn/To_Fn/etc. all live in the
+ * full Activations.h) -- so IComputeBackend.h, and everything that
+ * includes it (crucially CUDABackend.h/.cu, compiled by nvcc), never
+ * needed the SIMD-intrinsics half of Activations.h at all.
+*/
+
+namespace Deep
+{
+    enum class ActivationType : uint8_t
+    {
+        RELU,
+        dRELU,
+        GELU,
+        dGELU,
+        SIGMOID,
+        dSIGMOID,
+        eSIGMOID,
+        d_eSIGMOID,
+        TANH,
+        dTANH,
+        LINEAR,
+        dLINEAR,
+        NONE
+    };
+}
\ No newline at end of file
diff --git a/include/deepity/utils/Activations.h b/include/deepity/utils/Activations.h
index 70760a1..e5eedca 100644
--- a/include/deepity/utils/Activations.h
+++ b/include/deepity/utils/Activations.h
@@ -6,6 +6,7 @@
 #include <sleef.h>
 #include <omp.h>
 #include <algorithm>
+#include <deepity/utils/ActivationType.h>
 
 #if defined(_MSC_VER)
 #define RESTRICT __restrict
@@ -46,23 +47,6 @@ namespace Deep
     // afford to mutate src in place -- src stays untouched throughout.
     using DerivativeFn2 = void (*)(float *RESTRICT, const float *RESTRICT, size_t);
 
-    enum class ActivationType : uint8_t
-    {
-        RELU,
-        dRELU,
-        GELU,
-        dGELU,
-        SIGMOID,
-        dSIGMOID,
-        eSIGMOID,
-        d_eSIGMOID,
-        TANH,
-        dTANH,
-        LINEAR,
-        dLINEAR,
-        NONE
-    };
-
     static inline void relu(float *, size_t) noexcept;
     static inline void gelu(float *, size_t) noexcept;
     static inline void sigmoid(float *, size_t) noexcept;
diff --git a/mnist.py b/mnist.py
index b7f47e4..4df1fa2 100644
--- a/mnist.py
+++ b/mnist.py
@@ -2,8 +2,10 @@
 import os
 import sys
 from pydeepity import DKPPCN
+from pydeepity.layer import Linear, Sigmoid
 from time import perf_counter
 
+
 def load_full_mnist():
     import gzip
     import urllib.request
@@ -39,6 +41,7 @@ def load_full_mnist():
     Y_train[np.arange(y_train_labels.shape[0]), y_train_labels] = 1.0 - eps
     return X_train, Y_train, X_test, y_test_labels
 
+
 def main() -> None:
     SEED = int(sys.argv[1]) if len(sys.argv) > 1 else 7
     EPOCHS = int(sys.argv[2]) if len(sys.argv) > 2 else 50
@@ -47,23 +50,29 @@ def main() -> None:
     X_train, Y_train, X_test, y_test_labels = load_full_mnist()
 
     BATCH_SIZE = 250
-    TERMINAL_SIZE = 10
     LR = 0.00373
     IR = 0.15
     FL = 1e-3
-    LMBDA = 1e-4  # The crucial Kolen-Pollack alignment decay
+    LMBDA = 1e-4
     DECAY_RATE = 0.94
 
-    print(f"\nBuilding network (784->512->512->10), seed={SEED}...")
-    net = DKPPCN(batch_size=BATCH_SIZE, device="gpu")
-    net.add_layer(784, 512, TERMINAL_SIZE, lr=LR, ir=IR, fl=FL, lmbda=LMBDA, act="linear")
-#     net.add_layer(512, 512, TERMINAL_SIZE, lr=LR, ir=IR, fl=FL, lmbda=LMBDA, act="sigmoid")
-    net.add_layer(512, TERMINAL_SIZE, TERMINAL_SIZE, lr=LR, ir=IR, fl=FL, lmbda=LMBDA, act="sigmoid")
-    net.add_layer(TERMINAL_SIZE, 0, TERMINAL_SIZE, lr=LR, ir=IR, fl=FL, lmbda=LMBDA, act="linear")
-    net.set_optimizer("ADAM")
-    net.set_psi_optimizer("ADAM")
-    net.compile()
-    net.randomize_weights()
+    print(f"\nBuilding network (784->512->10), seed={SEED}...")
+
+    net = DKPPCN(
+        Linear(784, 512),
+        Linear(512, 10),
+        Sigmoid(),
+        batch_size=BATCH_SIZE,
+        device="gpu",
+    )
+    net.configure(
+        learning_rate=LR,
+        inference_rate=IR,
+        feedback_rate=FL,
+        lmbda=LMBDA,
+        optimizer="ADAMW",
+        psi_optimizer="ADAMW",
+    )
 
     print(f"\n*** FULL DKP-PC RUN ***")
     print(f"Training DKPPCN: {EPOCHS} epochs, inference_steps={INFERENCE_STEPS}, ")
@@ -72,12 +81,11 @@ def main() -> None:
     rng = np.random.default_rng(SEED)
     n_batches = len(X_train) // BATCH_SIZE
     start_time = perf_counter()
-    epoch_accs = []
 
     for epoch in range(EPOCHS):
         current_lr = LR * (DECAY_RATE ** epoch)
         net.set_learning_rate(current_lr)
-        
+
         current_fl = FL * (DECAY_RATE ** epoch)
         net.set_feedback_rate(current_fl)
 
@@ -107,7 +115,6 @@ def main() -> None:
             total += BATCH_SIZE
 
         epoch_acc = 100.0 * correct / total
-        epoch_accs.append(epoch_acc)
         avg_energy = epoch_energy / n_batches
         elapsed = perf_counter() - start_time
         print(f"Epoch {epoch+1}/{EPOCHS} | Time: {elapsed:.1f}s | Acc: {epoch_acc:.2f}% | Avg energy: {avg_energy:.4f}")
diff --git a/pgo_workload.py b/pgo_workload.py
index 3943d0a..21e2a4d 100644
--- a/pgo_workload.py
+++ b/pgo_workload.py
@@ -1,18 +1,3 @@
-"""
-Short, dedicated workload for PGO profile collection -- run by
-deepity_build/cli.py's --pgo pass, not meant to be invoked directly for
-training. Deliberately NOT the full mnist.py run: PGO only needs to see
-which code paths are hot (branch outcomes, call frequency), and the
-settling loop's structure repeats identically on batch 1 and batch 234,
-so a few dozen batches already captures the same information a full
-15-epoch run would. This cost is paid on every --pgo build, so keeping
-it short matters.
-
-Matches the real network configuration from mnist.py exactly (same
-architecture, same activations) so the instrumented binary actually
-exercises the same code paths real training does -- a mismatched
-architecture would profile the wrong thing.
-"""
 import numpy as np
 import os
 from pydeepity import SimplePCN
diff --git a/pydeepity/DKPPCN.py b/pydeepity/DKPPCN.py
index 7e5f789..cf3a0dc 100644
--- a/pydeepity/DKPPCN.py
+++ b/pydeepity/DKPPCN.py
@@ -53,8 +53,10 @@ def _validate_architecture(self) -> None:
         if not isinstance(self.architecture[0], Linear):
             raise TypeError("Architecture must begin with a Linear layer.")
 
-        if not isinstance(self.architecture[-1], Linear):
-            raise TypeError("Architecture must end with a Linear layer.")
+        if not isinstance(self.architecture[-1], Linear) and not isinstance(self.architecture[-1], Activation):
+            raise TypeError(
+                "Network architecture must end with a Linear or Activation layer."
+            )
 
         previous_linear = None
 
@@ -95,10 +97,8 @@ def _build_backend(self) -> None:
             else:
                 activation = "linear"
 
-            # 1. Define the derivative string
             activation_deriv = "d" + activation
 
-            # 2. Pass them to the backend using the exact kwarg names it expects
             super().add_layer(
                 layer.in_n,
                 layer.out_n,
@@ -107,10 +107,22 @@ def _build_backend(self) -> None:
                 ir=self._inference_rate,
                 fl=self._feedback_rate,
                 lmbda=self._lambda,
-                activation=activation,             # Changed from act=
-                activation_deriv=activation_deriv, # Added derivative
+                activation=activation,
+                activation_deriv=activation_deriv,
             )
 
+        super().add_layer(
+            terminal_size,
+            0,
+            terminal_size,
+            lr=self._learning_rate,
+            ir=self._inference_rate,
+            fl=self._feedback_rate,
+            lmbda=self._lambda,
+            activation="linear",
+            activation_deriv="dlinear",
+        )
+
     def _terminal_size(self) -> int:
         for layer in reversed(self.architecture):
             if isinstance(layer, Linear):
diff --git a/run.slurm b/run.slurm
new file mode 100644
index 0000000..f6c6693
--- /dev/null
+++ b/run.slurm
@@ -0,0 +1,37 @@
+#!/bin/bash
+#SBATCH --job-name=deepity-mnist-dkppcn
+#SBATCH --account=pas0350
+#SBATCH --gpus=1
+#SBATCH --time=01:00:00
+#SBATCH --mem=16G
+#SBATCH --output=slurm-%j.out
+#SBATCH --error=slurm-%j.err
+
+# Slurm batch script for running Deepity's mnist.py (DKPPCN) on an A100
+# GPU on OSC's Ascend cluster. Matches everything confirmed working
+# tonight: cuda/12.9.1 (12.4.1 hit a real, documented nvcc bug parsing
+# GCC 11's amxtileintrin.h -- see include-patches/ or the earlier fix
+# discussion if this resurfaces on a different node/module set), the
+# project's .venv for python3.12, account pas0350.
+#
+# Submit with: sbatch run_mnist_gpu.slurm
+# Check status: squeue -u $USER
+# Watch output live: tail -f slurm-<jobid>.out
+
+set -euo pipefail
+
+module load cuda/12.9.1
+
+cd "$SLURM_SUBMIT_DIR"
+
+# Activate the project's venv (contains python3.12 + nanobind + the
+# built pydeepity extension module).
+source .venv/bin/activate
+
+echo "=== Job started on $(hostname) at $(date) ==="
+echo "=== GPU visible to this job: ==="
+nvidia-smi --query-gpu=name,memory.total --format=csv
+
+python mnist.py
+
+echo "=== Job finished at $(date) ==="
diff --git a/src/backend/CUDABackend.cu b/src/backend/CUDABackend.cu
index 4c9ed4d..f77b859 100644
--- a/src/backend/CUDABackend.cu
+++ b/src/backend/CUDABackend.cu
@@ -7,6 +7,15 @@
 #ifdef DEEPITY_USE_CUDA
 #include <curand_kernel.h>
 
+#define CHECK_CUDA_LAUNCH() \
+    do { \
+        cudaError_t err = cudaGetLastError(); \
+        if (err != cudaSuccess) { \
+            std::cerr << "CUDA error at " << __FILE__ << ":" << __LINE__ \
+                      << " -> " << cudaGetErrorString(err) << std::endl; \
+        } \
+    } while (0)
+
 namespace Deep
 {
     CUDABackend::CUDABackend()
@@ -14,10 +23,6 @@ namespace Deep
         cudaStreamCreate(&this->stream);
         cublasCreate(&this->handle);
         cublasSetStream(this->handle, this->stream);
-
-        //        constexpr size_t WORKSPACE_SIZE = 4 * 1024 * 1024;
-        //        cudaMalloc(&workspace, WORKSPACE_SIZE);
-        //        cublasSetWorkspace(handle, workspace, WORKSPACE_SIZE);
     }
 
     CUDABackend::~CUDABackend()
@@ -27,6 +32,8 @@ namespace Deep
             cudaGraphExecDestroy(graphExec);
             cudaGraphDestroy(graph);
         }
+        if (onesVector)
+            cudaFree(onesVector);
         if (workspace)
             cudaFree(workspace);
         cublasDestroy(this->handle);
@@ -96,28 +103,32 @@ namespace Deep
 
     void CUDABackend::Free(float *ptr) noexcept
     {
-        if (cudaFree(ptr) != cudaSuccess)
-            std::cerr << "Could not free memory from CUDA backend.\n";
+        if (ptr)
+            cudaFree(ptr);
     }
 
     void CUDABackend::Zero(float *ptr, size_t numFloats) noexcept
     {
-        cudaMemsetAsync(ptr, 0, numFloats * sizeof(float), stream);
+        if (ptr && numFloats > 0)
+            cudaMemsetAsync(ptr, 0, numFloats * sizeof(float), stream);
     }
 
     void CUDABackend::Copy(float *dst, const float *src, size_t numFloats) noexcept
     {
-        cudaMemcpyAsync(dst, src, numFloats * sizeof(float), cudaMemcpyDefault, stream);
+        if (dst && src && numFloats > 0)
+            cudaMemcpyAsync(dst, src, numFloats * sizeof(float), cudaMemcpyDefault, stream);
     }
 
     void CUDABackend::CopyFromHost(float *deviceDst, const float *hostSrc, size_t numFloats) noexcept
     {
-        cudaMemcpy(deviceDst, hostSrc, numFloats * sizeof(float), cudaMemcpyHostToDevice);
+        if (deviceDst && hostSrc && numFloats > 0)
+            cudaMemcpy(deviceDst, hostSrc, numFloats * sizeof(float), cudaMemcpyHostToDevice);
     }
 
     void CUDABackend::CopyToHost(float *hostDst, const float *deviceSrc, size_t numFloats) noexcept
     {
-        cudaMemcpy(hostDst, deviceSrc, numFloats * sizeof(float), cudaMemcpyDeviceToHost);
+        if (hostDst && deviceSrc && numFloats > 0)
+            cudaMemcpy(hostDst, deviceSrc, numFloats * sizeof(float), cudaMemcpyDeviceToHost);
     }
 
     __global__ void normal_generation(curandState *state, float *random_numbers,
@@ -143,8 +154,7 @@ namespace Deep
 
     void CUDABackend::RandomizeNormal(float *buf, size_t n, float mean, float stddev, uint32_t seed) noexcept
     {
-        if (n == 0)
-            return;
+        if (!buf || n == 0) return;
 
         curandState *state = nullptr;
         if (cudaMalloc(&state, n * sizeof(curandState)) != cudaSuccess)
@@ -153,38 +163,24 @@ namespace Deep
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
         normal_generation<<<blocks, BLOCK_SIZE, 0, stream>>>(state, buf, n, mean, stddev, seed);
+        CHECK_CUDA_LAUNCH();
         cudaStreamSynchronize(stream);
         cudaFree(state);
     }
 
     void CUDABackend::RandomizeUniform(float *buf, size_t n, float min, float max, uint32_t seed) noexcept
     {
-        if (n == 0)
-            return;
+        if (!buf || n == 0) return;
 
         curandState *state = nullptr;
-        cudaError_t mallocErr = cudaMalloc(&state, n * sizeof(curandState));
-        if (mallocErr != cudaSuccess)
-        {
-            std::cerr << "curandState cudaMalloc failed for n=" << n
-                      << " (" << n * sizeof(curandState) << " bytes): "
-                      << cudaGetErrorString(mallocErr) << "\n";
+        if (cudaMalloc(&state, n * sizeof(curandState)) != cudaSuccess)
             return;
-        }
 
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
         uniform_generation<<<blocks, BLOCK_SIZE, 0, stream>>>(state, buf, n, min, max - min, seed);
-
-        cudaError_t launchErr = cudaGetLastError();
-        if (launchErr != cudaSuccess)
-            std::cerr << "uniform_generation kernel launch failed: " << cudaGetErrorString(launchErr) << "\n";
-
+        CHECK_CUDA_LAUNCH();
         cudaStreamSynchronize(stream);
-        cudaError_t syncErr = cudaGetLastError();
-        if (syncErr != cudaSuccess)
-            std::cerr << "uniform_generation kernel execution failed: " << cudaGetErrorString(syncErr) << "\n";
-
         cudaFree(state);
     }
 
@@ -197,76 +193,79 @@ namespace Deep
 
     void CUDABackend::PrepareForBatchSize(size_t batchSize) noexcept
     {
-        if (onesVector)
-            cudaFree(onesVector);
-        cudaMalloc(&onesVector, batchSize * sizeof(float));
+        if (onesCapacity < batchSize)
+        {
+            if (onesVector)
+                cudaFree(onesVector);
+            cudaMalloc(&onesVector, batchSize * sizeof(float));
+            onesCapacity = batchSize;
+        }
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((batchSize + BLOCK_SIZE - 1) / BLOCK_SIZE);
         FillOnesKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(onesVector, batchSize);
-        cudaStreamSynchronize(stream);
+        CHECK_CUDA_LAUNCH();
     }
 
     void CUDABackend::MatMul(bool transA, bool transB, int M, int N, int K,
-                             float alpha, const float *A, int lda,
-                             const float *B, int ldb,
-                             float beta, float *C, int ldc) noexcept
-    {
-        cublasOperation_t cuTransA = transA ? CUBLAS_OP_T : CUBLAS_OP_N;
-        cublasOperation_t cuTransB = transB ? CUBLAS_OP_T : CUBLAS_OP_N;
-        cublasStatus_t status = cublasSgemm(handle, cuTransB, cuTransA,
-                                            N, M, K, &alpha, B, ldb, A, lda, &beta, C, ldc);
-        if (status != CUBLAS_STATUS_SUCCESS)
-        {
-            cudaError_t cudaErr = cudaGetLastError();
-            std::cerr << "cublasSgemm status: " << status
-                      << ", underlying cudaError: " << cudaErr
-                      << " (" << cudaGetErrorString(cudaErr) << ")\n";
-        }
+                            float alpha, const float *A, int lda,
+                            const float *B, int ldb,
+                            float beta, float *C, int ldc) noexcept
+    {
+        if (!A || !B || !C) return;
+
+        cublasOperation_t opA = transA ? CUBLAS_OP_T : CUBLAS_OP_N;
+        cublasOperation_t opB = transB ? CUBLAS_OP_T : CUBLAS_OP_N;
+
+        cublasSgemm(this->handle, opB, opA, N, M, K, &alpha, B, ldb, A, lda, &beta, C, ldc);
     }
 
     void CUDABackend::SumRows(float *dst, const float *src, size_t batchSize, size_t width) noexcept
     {
+        if (!dst || !src) return;
+
+        if (!onesVector || onesCapacity < batchSize)
+            PrepareForBatchSize(batchSize);
+
         float alpha = 1.0f, beta = 0.0f;
-        cublasStatus_t status = cublasSgemv(handle, CUBLAS_OP_N, width, batchSize,
-                                            &alpha, src, width, onesVector, 1, &beta, dst, 1);
-        if (status != CUBLAS_STATUS_SUCCESS)
-            std::cerr << "cublasSgemv (SumRows) status: " << status << "\n";
+        cublasSgemv(handle, CUBLAS_OP_N, width, batchSize, &alpha, src, width, onesVector, 1, &beta, dst, 1);
     }
 
     void CUDABackend::Scale(float *buf, size_t n, float alpha) noexcept
     {
-        if (cublasSscal(handle, n, &alpha, buf, 1) != CUBLAS_STATUS_SUCCESS)
-            std::cerr << "Failed to perform CUDA Scale.\n";
+        if (!buf || n == 0) return;
+        cublasSscal(handle, n, &alpha, buf, 1);
     }
 
     void CUDABackend::AxpyInto(float *y, const float *x, size_t n, float alpha) noexcept
     {
-        if (cublasSaxpy(handle, n, &alpha, x, 1, y, 1))
-            std::cerr << "Failed to perform CUDA Axpy.\n";
+        if (!x || !y || n == 0) return;
+        cublasSaxpy(handle, n, &alpha, x, 1, y, 1);
     }
-
+    
     __global__ void AddBiasBroadcastKernel(float *buf, const float *bias, size_t batchSize, size_t width)
     {
-        size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < batchSize * width)
-            buf[i] += bias[i % width];
+        size_t col = (size_t)blockIdx.x * blockDim.x + threadIdx.x; // maps to width
+        size_t row = (size_t)blockIdx.y * blockDim.y + threadIdx.y; // maps to batchSize
+
+        if (col < width && row < batchSize)
+            buf[row * width + col] += bias[col];
     }
 
     void CUDABackend::AddBiasBroadcast(float *buf, const float *bias, size_t batchSize, size_t width) noexcept
     {
-        constexpr int BLOCK_SIZE = 256;
-        size_t total = batchSize * width;
-        const int blocks = static_cast<int>((total + BLOCK_SIZE - 1) / BLOCK_SIZE);
-        AddBiasBroadcastKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(buf, bias, batchSize, width);
-    }
+        dim3 threads(32, 8);
+        dim3 blocks(
+            (unsigned int)((width + threads.x - 1) / threads.x),
+            (unsigned int)((batchSize + threads.y - 1) / threads.y));
 
-#pragma region ACTIVATIONS_AND_KERNELS
+        AddBiasBroadcastKernel<<<blocks, threads, 0, stream>>>(buf, bias, batchSize, width);
+        CHECK_CUDA_LAUNCH();
+    }
 
     __global__ void ReluKernelInto(float *dst, const float *src, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            dst[i] = fmaxf(0.0f, src[i]);
+        if (i < n) dst[i] = fmaxf(0.0f, src[i]);
     }
 
     __global__ void GeluKernelInto(float *dst, const float *src, size_t n)
@@ -276,14 +275,12 @@ namespace Deep
         {
             float xi = src[i];
             float inner = MAGIC_GELU_1 * xi * (1.0f + MAGIC_GELU_2 * xi * xi);
-
             float t;
 #if __CUDA_ARCH__ >= 800
             asm("tanh.approx.f32 %0, %1;" : "=f"(t) : "f"(inner));
 #else
             t = tanhf(inner);
 #endif
-
             dst[i] = 0.5f * xi * (1.0f + t);
         }
     }
@@ -323,18 +320,16 @@ namespace Deep
     __global__ void linearKernelInto(float *dst, const float *src, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            dst[i] = src[i];
+        if (i < n) dst[i] = src[i];
     }
 
     __global__ void dReluKernelInto(float *dst, const float *src, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            dst[i] = (float)(src[i] > 0.0f);
+        if (i < n) dst[i] = (float)(src[i] > 0.0f);
     }
 
-    constexpr int MAGIC_GELU_2_3 = 3.0f * MAGIC_GELU_2;
+    constexpr float MAGIC_GELU_2_3 = 3.0f * MAGIC_GELU_2;
 
     __global__ void dGeluKernelInto(float *dst, const float *src, size_t n)
     {
@@ -344,18 +339,15 @@ namespace Deep
             float x = src[i];
             float xsq = x * x;
             float inner = MAGIC_GELU_1 * x * (1.0f + MAGIC_GELU_2 * xsq);
-
             float t;
 #if __CUDA_ARCH__ >= 800
             asm("tanh.approx.f32 %0, %1;" : "=f"(t) : "f"(inner));
 #else
             t = tanhf(inner);
 #endif
-
             float gprime = MAGIC_GELU_1 * (1.0f + MAGIC_GELU_2_3 * xsq);
             float term1 = 0.5f * (1.0f + t);
             float term2 = 0.5f * x * gprime * (1.0f - t * t);
-
             dst[i] = term1 + term2;
         }
     }
@@ -382,16 +374,6 @@ namespace Deep
         }
     }
 
-    __global__ void dSigmoidActivatedKernelInto(float *dst, const float *src, size_t n)
-    {
-        size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-        {
-            float s = src[i];
-            dst[i] = fmaf(-s, s, s);
-        }
-    }
-
     __global__ void d_eSigmoidKernelInto(float *dst, const float *src, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
@@ -402,21 +384,10 @@ namespace Deep
         }
     }
 
-    __global__ void d_eSigmoidActivatedKernelInto(float *dst, const float *src, size_t n)
-    {
-        size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-        {
-            float s = src[i];
-            dst[i] = 2.0f * fmaf(-s, s, s);
-        }
-    }
-
     __global__ void dLinearKernelInto(float *dst, const float *src, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            dst[i] = 1.0f;
+        if (i < n) dst[i] = 1.0f;
     }
 
     __global__ void FusedStateUpdateKernel(float *z, const float *feedback, const float *deriv,
@@ -425,21 +396,19 @@ namespace Deep
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
         if (i < n)
         {
-            z[i] += ir * ((feedback[i] * deriv[i]) - e[i]);
+            float fb = feedback ? feedback[i] : 0.0f;
+            float dev = deriv ? deriv[i] : 1.0f;
+            float err = e ? e[i] : 0.0f;
+            z[i] += ir * ((fb * dev) - err);
         }
     }
 
     __global__ void ComputeErrorKernel(float *e, const float *z, const float *mu, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-        {
-            e[i] = z[i] - mu[i];
-        }
+        if (i < n) e[i] = z[i] - mu[i];
     }
 
-#pragma endregion
-
     void CUDABackend::Activation(ActivationType type, float *buf, size_t n) noexcept
     {
         ActivationInto(type, buf, buf, n);
@@ -447,33 +416,22 @@ namespace Deep
 
     void CUDABackend::ActivationInto(ActivationType type, float *dst, const float *src, size_t n) noexcept
     {
+        if (!dst || !src || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
 
         switch (type)
         {
-        case ActivationType::RELU:
-            ReluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::GELU:
-            GeluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::SIGMOID:
-            sigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::eSIGMOID:
-            eSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::TANH:
-            tanhKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::LINEAR:
-            linearKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
+        case ActivationType::RELU: ReluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::GELU: GeluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::SIGMOID: sigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::eSIGMOID: eSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::TANH: tanhKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::LINEAR: linearKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
         case ActivationType::NONE:
-        default:
-            break;
+        default: break;
         }
+        CHECK_CUDA_LAUNCH();
     }
 
     void CUDABackend::ActivationDerivative(ActivationType type, float *buf, size_t n, bool activated) noexcept
@@ -483,66 +441,67 @@ namespace Deep
 
     void CUDABackend::ActivationDerivativeInto(ActivationType type, float *dst, const float *src, size_t n) noexcept
     {
+        if (!dst || !src || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
 
         switch (type)
         {
-        case ActivationType::dRELU:
-            dReluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::dGELU:
-            dGeluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::dSIGMOID:
-            dSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::d_eSIGMOID:
-            d_eSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::dTANH:
-            dTanhKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
-        case ActivationType::dLINEAR:
-            dLinearKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n);
-            break;
+        case ActivationType::dRELU: dReluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::dGELU: dGeluKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::dSIGMOID: dSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::d_eSIGMOID: d_eSigmoidKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::dTANH: dTanhKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
+        case ActivationType::dLINEAR: dLinearKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
         case ActivationType::NONE:
-        default:
-            break;
+        default: dLinearKernelInto<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, n); break;
         }
+        CHECK_CUDA_LAUNCH();
     }
 
     void CUDABackend::FusedStateUpdate(float *z, const float *feedback, const float *deriv,
                                        const float *e, size_t n, float ir) noexcept
     {
+        if (!z || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
-
         FusedStateUpdateKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(z, feedback, deriv, e, n, ir);
+        CHECK_CUDA_LAUNCH();
     }
 
     float CUDABackend::ComputeErrorAndEnergy(float *e, const float *z, const float *mu, size_t n) noexcept
     {
+        if (!e || !z || !mu || n == 0) return 0.0f;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
 
         ComputeErrorKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(e, z, mu, n);
+        CHECK_CUDA_LAUNCH();
 
         float sum_of_squares = 0.0f;
         cublasSdot(handle, n, e, 1, e, 1, &sum_of_squares);
+        cudaStreamSynchronize(stream);
 
         return 0.5f * sum_of_squares;
     }
 
     void CUDABackend::ComputeError(float *e, const float *z, const float *mu, size_t n) noexcept
     {
+        if (!e || !z || !mu || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
         ComputeErrorKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(e, z, mu, n);
+        CHECK_CUDA_LAUNCH();
     }
 
-    __global__ void IncrementCounterKernel(int *counter) { *counter += 1; }
-
+    __global__ void IncrementCounterKernel(int *counter) { if (counter) *counter += 1; }
+    void CUDABackend::IncrementCounter(int *counter) noexcept
+    {
+        if (!counter) return;
+        IncrementCounterKernel<<<1, 1, 0, stream>>>(counter);
+        CHECK_CUDA_LAUNCH();
+    }
+    
     __global__ void AdamStepKernel(float *param, const float *grad, float *m, float *v,
                                    size_t n, const int *t_ptr, const float *lr_ptr,
                                    float beta1, float beta2, float eps)
@@ -550,14 +509,21 @@ namespace Deep
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
         if (i < n)
         {
-            float beta1_t = 1.0f - powf(beta1, (float)(*t_ptr));
-            float beta2_t = 1.0f - powf(beta2, (float)(*t_ptr));
-            float step_size = *lr_ptr * sqrtf(beta2_t) / beta1_t;
+            int current_t = *t_ptr;
+            if (current_t < 1) current_t = 1;
+            float current_lr = *lr_ptr;
+
+            float beta1_t = 1.0f - powf(beta1, static_cast<float>(current_t));
+            float beta2_t = 1.0f - powf(beta2, static_cast<float>(current_t));
+            float step_size = current_lr * sqrtf(beta2_t) / beta1_t;
 
             float g = grad[i];
-            m[i] = beta1 * m[i] + (1.0f - beta1) * g;
-            v[i] = beta2 * v[i] + (1.0f - beta2) * (g * g);
-            param[i] -= step_size * m[i] / (sqrtf(v[i]) + eps);
+            float m_val = beta1 * m[i] + (1.0f - beta1) * g;
+            float v_val = beta2 * v[i] + (1.0f - beta2) * (g * g);
+
+            m[i] = m_val;
+            v[i] = v_val;
+            param[i] -= step_size * m_val / (sqrtf(v_val) + eps);
         }
     }
 
@@ -568,72 +534,85 @@ namespace Deep
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
         if (i < n)
         {
-            float beta1_t = 1.0f - powf(beta1, (float)(*t_ptr));
-            float beta2_t = 1.0f - powf(beta2, (float)(*t_ptr));
-            float step_size = *lr_ptr * sqrtf(beta2_t) / beta1_t;
+            int current_t = *t_ptr;
+            if (current_t < 1) current_t = 1;
+            float current_lr = *lr_ptr;
+
+            float beta1_t = 1.0f - powf(beta1, static_cast<float>(current_t));
+            float beta2_t = 1.0f - powf(beta2, static_cast<float>(current_t));
+            float step_size = current_lr * sqrtf(beta2_t) / beta1_t;
 
             float g = grad[i];
-            m[i] = beta1 * m[i] + (1.0f - beta1) * g;
-            v[i] = beta2 * v[i] + (1.0f - beta2) * (g * g);
-            param[i] -= *lr_ptr * weightDecay * param[i];
-            param[i] -= step_size * m[i] / (sqrtf(v[i]) + eps);
+            float p = param[i];
+            float m_val = beta1 * m[i] + (1.0f - beta1) * g;
+            float v_val = beta2 * v[i] + (1.0f - beta2) * (g * g);
+
+            m[i] = m_val;
+            v[i] = v_val;
+            p -= current_lr * weightDecay * p;
+            p -= step_size * m_val / (sqrtf(v_val) + eps);
+            param[i] = p;
         }
     }
 
-    void CUDABackend::IncrementCounter(int *counter) noexcept
-    {
-        IncrementCounterKernel<<<1, 1, 0, stream>>>(counter);
-    }
-
     void CUDABackend::AdamStep(float *param, const float *grad, float *m, float *v,
                                size_t n, const int *t, const float *lr,
                                float beta1, float beta2, float eps) noexcept
     {
-        constexpr int BLOCK_SIZE = 256;
-        const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
-        AdamStepKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(param, grad, m, v, n, t, lr, beta1, beta2, eps);
+        if (!param || !grad || !m || !v || !t || !lr || n == 0) return;
+
+        // NO host copy -- t/lr stay as device pointers, dereferenced
+        // inside the kernel itself, so graph capture/replay re-reads the
+        // real, current value every time instead of baking in a
+        // one-time snapshot from whenever capture happened to run.
+        constexpr int blockSize = 256;
+        int numBlocks = static_cast<int>((n + blockSize - 1) / blockSize);
+        AdamStepKernel<<<numBlocks, blockSize, 0, stream>>>(param, grad, m, v, n, t, lr, beta1, beta2, eps);
+        CHECK_CUDA_LAUNCH();
     }
 
     void CUDABackend::AdamWStep(float *param, const float *grad, float *m, float *v,
                                 size_t n, const int *t, const float *lr, float weightDecay,
                                 float beta1, float beta2, float eps) noexcept
     {
-        constexpr int BLOCK_SIZE = 256;
-        const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
-        AdamWStepKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(param, grad, m, v, n, t, lr, weightDecay, beta1, beta2, eps);
+        if (!param || !grad || !m || !v || !t || !lr || n == 0) return;
+
+        constexpr int blockSize = 256;
+        int numBlocks = static_cast<int>((n + blockSize - 1) / blockSize);
+        AdamWStepKernel<<<numBlocks, blockSize, 0, stream>>>(param, grad, m, v, n, t, lr, weightDecay, beta1, beta2, eps);
+        CHECK_CUDA_LAUNCH();
     }
 
     __global__ void MultiplyIntoKernel(float *dst, const float *a, const float *b, size_t n)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            dst[i] = a[i] * b[i];
+        if (i < n) dst[i] = (a && b) ? (a[i] * b[i]) : 0.0f;
     }
 
     void CUDABackend::MultiplyInto(float *dst, const float *a, const float *b, size_t n) noexcept
     {
+        if (!dst || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
         MultiplyIntoKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, a, b, n);
+        CHECK_CUDA_LAUNCH();
     }
 
     __global__ void FillKernel(float *buf, size_t n, float value)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
-        if (i < n)
-            buf[i] = value;
+        if (i < n) buf[i] = value;
     }
 
     void CUDABackend::Fill(float *buf, size_t n, float value) noexcept
     {
+        if (!buf || n == 0) return;
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((n + BLOCK_SIZE - 1) / BLOCK_SIZE);
         FillKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(buf, n, value);
+        CHECK_CUDA_LAUNCH();
     }
 
-    // buf[c*spatialSize + s] += bias[c] for all c,s -- convolutional bias-
-    // add. NOT batch-aware, matching Im2Col/Col2Im's own per-item contract
-    // (caller loops over batch, offsetting buf each time).
     __global__ void AddBiasPerChannelKernel(float *buf, const float *bias, size_t channels, size_t spatialSize)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
@@ -641,33 +620,26 @@ namespace Deep
         if (i < total)
         {
             size_t c = i / spatialSize;
-            buf[i] += bias[c];
+            buf[i] += bias ? bias[c] : 0.0f;
         }
     }
 
     void CUDABackend::AddBiasPerChannel(float *buf, const float *bias, size_t channels, size_t spatialSize) noexcept
     {
+        if (!buf || !bias) return;
         constexpr int BLOCK_SIZE = 256;
         size_t total = channels * spatialSize;
         const int blocks = static_cast<int>((total + BLOCK_SIZE - 1) / BLOCK_SIZE);
         AddBiasPerChannelKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(buf, bias, channels, spatialSize);
+        CHECK_CUDA_LAUNCH();
     }
 
-    // dst[row][batch][:] = src[batch][row][:] -- one thread per output
-    // (dst-indexed) element, decomposing the flat index into
-    // (row, batch, col) to compute the corresponding source offset.
-    // Verified against CPUBackend's identical-purpose loop, same test case
-    // (batchSize=2, rows=3, cols=2), exact match. NOTE: 4 integer div/mod
-    // ops per thread -- correct but not optimal; a future optimization
-    // could use a 3D launch grid to let hardware indexing do this instead,
-    // same "port first, optimize" position as everything else tonight.
     __global__ void RepackForBatchedGemmKernel(float *dst, const float *src,
                                                size_t batchSize, size_t rows, size_t cols)
     {
         size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
         size_t total = rows * batchSize * cols;
-        if (i >= total)
-            return;
+        if (i >= total) return;
 
         size_t row = i / (batchSize * cols);
         size_t rem = i % (batchSize * cols);
@@ -681,17 +653,14 @@ namespace Deep
     void CUDABackend::RepackForBatchedGemm(float *dst, const float *src,
                                            size_t batchSize, size_t rows, size_t cols) noexcept
     {
+        if (!dst || !src) return;
         constexpr int BLOCK_SIZE = 256;
         size_t total = rows * batchSize * cols;
         const int blocks = static_cast<int>((total + BLOCK_SIZE - 1) / BLOCK_SIZE);
         RepackForBatchedGemmKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(dst, src, batchSize, rows, cols);
+        CHECK_CUDA_LAUNCH();
     }
 
-    // Im2Col/Col2Im: one thread per (row, col) element of the
-    // [channels*kH*kW, outH*outW] column matrix, matching Deep::Im2Col/
-    // Deep::Col2Im's exact CPU semantics (NCHW, zero for out-of-bounds
-    // positions). Both NOT batch-aware -- caller loops over batch, offsetting
-    // input/columns each time, same contract as the CPU free functions.
     __global__ void Im2ColKernel(const float *input, int channels, int height, int width,
                                  int kH, int kW, int strideH, int strideW, int padH, int padW,
                                  int outH, int outW, float *columns)
@@ -700,8 +669,7 @@ namespace Deep
         size_t colCols = (size_t)outH * outW;
         size_t colRows = (size_t)channels * kH * kW;
         size_t total = colRows * colCols;
-        if (i >= total)
-            return;
+        if (i >= total) return;
 
         size_t row = i / colCols;
         size_t col = i % colCols;
@@ -726,25 +694,17 @@ namespace Deep
                              int kernelH, int kernelW, int strideH, int strideW, int padH, int padW,
                              float *columns) noexcept
     {
+        if (!input || !columns) return;
         int outH = ConvOutDim(height, kernelH, strideH, padH);
         int outW = ConvOutDim(width, kernelW, strideW, padW);
         size_t total = (size_t)channels * kernelH * kernelW * outH * outW;
 
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((total + BLOCK_SIZE - 1) / BLOCK_SIZE);
-        Im2ColKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(
-            input, channels, height, width,
-            kernelH, kernelW, strideH, strideW, padH, padW,
-            outH, outW, columns);
-    }
-
-    // Col2Im is the adjoint: SCATTERS, accumulating via atomicAdd since
-    // overlapping receptive fields (stride < kernel size) mean multiple
-    // (row,col) source elements can map to the same destination pixel --
-    // a plain, non-atomic write would race. Matches Deep::Col2Im's own
-    // contract exactly: ACCUMULATES into outputImage (does not zero it
-    // first); caller must Zero() the destination first if a fresh result
-    // is wanted.
+        Im2ColKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(input, channels, height, width, kernelH, kernelW, strideH, strideW, padH, padW, outH, outW, columns);
+        CHECK_CUDA_LAUNCH();
+    }
+
     __global__ void Col2ImKernel(const float *columns, int channels, int height, int width,
                                  int kH, int kW, int strideH, int strideW, int padH, int padW,
                                  int outH, int outW, float *outputImage)
@@ -753,8 +713,7 @@ namespace Deep
         size_t colCols = (size_t)outH * outW;
         size_t colRows = (size_t)channels * kH * kW;
         size_t total = colRows * colCols;
-        if (i >= total)
-            return;
+        if (i >= total) return;
 
         size_t row = i / colCols;
         size_t col = i % colCols;
@@ -780,16 +739,15 @@ namespace Deep
                              int kernelH, int kernelW, int strideH, int strideW, int padH, int padW,
                              float *outputImage) noexcept
     {
+        if (!columns || !outputImage) return;
         int outH = ConvOutDim(height, kernelH, strideH, padH);
         int outW = ConvOutDim(width, kernelW, strideW, padW);
         size_t total = (size_t)channels * kernelH * kernelW * outH * outW;
 
         constexpr int BLOCK_SIZE = 256;
         const int blocks = static_cast<int>((total + BLOCK_SIZE - 1) / BLOCK_SIZE);
-        Col2ImKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(
-            columns, channels, height, width,
-            kernelH, kernelW, strideH, strideW, padH, padW,
-            outH, outW, outputImage);
+        Col2ImKernel<<<blocks, BLOCK_SIZE, 0, stream>>>(columns, channels, height, width, kernelH, kernelW, strideH, strideW, padH, padW, outH, outW, outputImage);
+        CHECK_CUDA_LAUNCH();
     }
 }
 
diff --git a/tests/tCUDAFunctionsVerify.cpp b/tests/tCUDAFunctionsVerify.cpp
deleted file mode 100644
index d10b0fa..0000000
--- a/tests/tCUDAFunctionsVerify.cpp
+++ /dev/null
@@ -1,237 +0,0 @@
-/**
- * @file tCUDAFunctionsVerify.cpp
- * @brief Verifies every IComputeBackend method not already covered by
- * tMatMulVerify.cpp -- activations, derivatives, Scale, AxpyInto,
- * FusedStateUpdate, ComputeErrorAndEnergy, AdamStep, AdamWStep -- on
- * both CPUBackend and CUDABackend, against expected values computed
- * independently (via a separate Python script, not derived from this
- * codebase's own formulas).
- *
- * IMPORTANT PRECISION NOTE: CUDABackend's tanh/sigmoid kernels use the
- * hardware-approximate `tanh.approx.f32` PTX instruction, trading
- * precision for speed -- CPUBackend uses SLEEF's high-precision
- * (u10 = <=1.0 ULP error) implementation. These will NOT match to tight
- * precision even when both are correct. Functions using this instruction
- * (SIGMOID, TANH, and their derivatives) use a loose tolerance (1e-3);
- * everything else (exact arithmetic: RELU, LINEAR, eSIGMOID, Scale,
- * AxpyInto, FusedStateUpdate, ComputeErrorAndEnergy, Adam/AdamW) uses a
- * tight one (1e-4), since a loose match there would hide a real bug.
- */
-#include <deepity/backend/Backend.h>
-#include <cstdio>
-#include <cmath>
-#include <vector>
-#include <string>
-
-using namespace Deep;
-
-namespace
-{
-    int g_failures = 0;
-
-    void Check(const std::vector<float> &actual, const std::vector<float> &expected,
-               const char *testName, const char *backendName, float tolerance)
-    {
-        bool ok = true;
-        for (size_t i = 0; i < expected.size(); ++i)
-        {
-            if (std::fabs(actual[i] - expected[i]) > tolerance)
-            {
-                printf("  [%s / %s] MISMATCH at index %zu: got %.6f, expected %.6f (tol %.1e)\n",
-                       backendName, testName, i, actual[i], expected[i], tolerance);
-                ok = false;
-            }
-        }
-        if (ok)
-            printf("  [%s / %s] PASSED\n", backendName, testName);
-        else
-            g_failures++;
-    }
-
-    void CheckScalar(float actual, float expected, const char *testName,
-                     const char *backendName, float tolerance)
-    {
-        Check({actual}, {expected}, testName, backendName, tolerance);
-    }
-
-    void RunAllTests(DeviceType device, const char *backendName)
-    {
-        auto backend = CreateBackend(device);
-        const std::vector<float> xs = {-2.0f, -1.0f, -0.5f, 0.0f, 0.5f, 1.0f, 2.0f};
-        const size_t n = xs.size();
-
-        // --- Activations (independently computed via Python) ---------
-        {
-            Tensor t(backend.get(), device, xs);
-            backend->Activation(ActivationType::RELU, t.Data(), n);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {0, 0, 0, 0, 0.5f, 1.0f, 2.0f}, "relu", backendName, 1e-4f);
-        }
-        {
-            Tensor t(backend.get(), device, xs);
-            backend->Activation(ActivationType::SIGMOID, t.Data(), n);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {0.11920292f, 0.26894142f, 0.37754067f, 0.5f, 0.62245933f, 0.73105858f, 0.88079708f},
-                  "sigmoid", backendName, 1e-3f); // approx tanh-based on GPU
-        }
-        {
-            Tensor t(backend.get(), device, xs);
-            backend->Activation(ActivationType::eSIGMOID, t.Data(), n);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {0.16666667f, 0.25f, 0.33333333f, 0.5f, 0.66666667f, 0.75f, 0.83333333f},
-                  "e_sigmoid", backendName, 1e-4f); // exact arithmetic, no approx instruction
-        }
-        {
-            Tensor t(backend.get(), device, xs);
-            backend->Activation(ActivationType::TANH, t.Data(), n);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {-0.96402758f, -0.76159416f, -0.46211716f, 0.0f, 0.46211716f, 0.76159416f, 0.96402758f},
-                  "tanh", backendName, 1e-3f); // approx instruction on GPU
-        }
-        {
-            Tensor t(backend.get(), device, xs);
-            backend->Activation(ActivationType::LINEAR, t.Data(), n);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {-2, -1, -0.5f, 0, 0.5f, 1, 2}, "linear (identity)", backendName, 1e-4f);
-        }
-
-        // --- Derivatives, raw-input (...Into) variants ----------------
-        {
-            Tensor src(backend.get(), device, xs);
-            Tensor dst(backend.get(), device, n);
-            backend->ActivationDerivativeInto(ActivationType::dRELU, dst.Data(), src.Data(), n);
-            std::vector<float> out;
-            dst.CopyToHost(out);
-            Check(out, {0, 0, 0, 0, 1.0f, 1.0f, 1.0f}, "dRelu", backendName, 1e-4f);
-        }
-        {
-            Tensor src(backend.get(), device, xs);
-            Tensor dst(backend.get(), device, n);
-            backend->ActivationDerivativeInto(ActivationType::dSIGMOID, dst.Data(), src.Data(), n);
-            std::vector<float> out;
-            dst.CopyToHost(out);
-            Check(out, {0.10499359f, 0.19661193f, 0.23500371f, 0.25f, 0.23500371f, 0.19661193f, 0.10499359f},
-                  "dSigmoid", backendName, 1e-3f);
-        }
-        {
-            Tensor src(backend.get(), device, xs);
-            Tensor dst(backend.get(), device, n);
-            backend->ActivationDerivativeInto(ActivationType::d_eSIGMOID, dst.Data(), src.Data(), n);
-            std::vector<float> out;
-            dst.CopyToHost(out);
-            Check(out, {0.05555556f, 0.125f, 0.22222222f, 0.5f, 0.22222222f, 0.125f, 0.05555556f},
-                  "d_eSigmoid", backendName, 1e-4f);
-        }
-        {
-            Tensor src(backend.get(), device, xs);
-            Tensor dst(backend.get(), device, n);
-            backend->ActivationDerivativeInto(ActivationType::dTANH, dst.Data(), src.Data(), n);
-            std::vector<float> out;
-            dst.CopyToHost(out);
-            Check(out, {0.07065082f, 0.41997434f, 0.78644773f, 1.0f, 0.78644773f, 0.41997434f, 0.07065082f},
-                  "dTanh", backendName, 1e-3f);
-        }
-        {
-            Tensor src(backend.get(), device, xs);
-            Tensor dst(backend.get(), device, n);
-            backend->ActivationDerivativeInto(ActivationType::dLINEAR, dst.Data(), src.Data(), n);
-            std::vector<float> out;
-            dst.CopyToHost(out);
-            Check(out, {1, 1, 1, 1, 1, 1, 1}, "dLinear", backendName, 1e-4f);
-        }
-
-        // --- Scale: [1,2,3,4] *= 2.0 -----------------------------------
-        {
-            Tensor t(backend.get(), device, std::vector<float>{1, 2, 3, 4});
-            backend->Scale(t.Data(), 4, 2.0f);
-            std::vector<float> out;
-            t.CopyToHost(out);
-            Check(out, {2, 4, 6, 8}, "Scale", backendName, 1e-4f);
-        }
-
-        // --- AxpyInto: y=[1,1,1,1] += 2.0 * x=[1,2,3,4] -----------------
-        {
-            Tensor y(backend.get(), device, std::vector<float>{1, 1, 1, 1});
-            Tensor x(backend.get(), device, std::vector<float>{1, 2, 3, 4});
-            backend->AxpyInto(y.Data(), x.Data(), 4, 2.0f);
-            std::vector<float> out;
-            y.CopyToHost(out);
-            Check(out, {3, 5, 7, 9}, "AxpyInto", backendName, 1e-4f);
-        }
-
-        // --- FusedStateUpdate: z += ir*(feedback*deriv - e) -------------
-        {
-            Tensor z(backend.get(), device, std::vector<float>{1.0f, 2.0f});
-            Tensor feedback(backend.get(), device, std::vector<float>{0.5f, 1.0f});
-            Tensor deriv(backend.get(), device, std::vector<float>{2.0f, 0.5f});
-            Tensor e(backend.get(), device, std::vector<float>{0.1f, 0.2f});
-            backend->FusedStateUpdate(z.Data(), feedback.Data(), deriv.Data(), e.Data(), 2, 0.1f);
-            std::vector<float> out;
-            z.CopyToHost(out);
-            Check(out, {1.09f, 2.03f}, "FusedStateUpdate", backendName, 1e-4f);
-        }
-
-        // --- ComputeErrorAndEnergy: e=z-mu, returns 0.5*sum(e^2) --------
-        {
-            Tensor z(backend.get(), device, std::vector<float>{3.0f, 5.0f});
-            Tensor mu(backend.get(), device, std::vector<float>{1.0f, 2.0f});
-            Tensor e(backend.get(), device, 2);
-            float energy = backend->ComputeErrorAndEnergy(e.Data(), z.Data(), mu.Data(), 2);
-            std::vector<float> eOut;
-            e.CopyToHost(eOut);
-            Check(eOut, {2.0f, 3.0f}, "ComputeErrorAndEnergy (e)", backendName, 1e-4f);
-            CheckScalar(energy, 6.5f, "ComputeErrorAndEnergy (energy)", backendName, 1e-4f);
-        }
-
-        // --- AdamStep: param=1.0, grad=0.5, m=v=0, t=1 ------------------
-        {
-            Tensor param(backend.get(), device, std::vector<float>{1.0f});
-            Tensor grad(backend.get(), device, std::vector<float>{0.5f});
-            Tensor m(backend.get(), device, 1);
-            Tensor v(backend.get(), device, 1);
-            backend->AdamStep(param.Data(), grad.Data(), m.Data(), v.Data(), 1, 1, 0.1f);
-            std::vector<float> out;
-            param.CopyToHost(out);
-            CheckScalar(out[0], 0.9000000632455132f, "AdamStep", backendName, 1e-4f);
-        }
-
-        // --- AdamWStep: same, plus weightDecay=0.01 ---------------------
-        {
-            Tensor param(backend.get(), device, std::vector<float>{1.0f});
-            Tensor grad(backend.get(), device, std::vector<float>{0.5f});
-            Tensor m(backend.get(), device, 1);
-            Tensor v(backend.get(), device, 1);
-            backend->AdamWStep(param.Data(), grad.Data(), m.Data(), v.Data(), 1, 1, 0.1f, 0.01f);
-            std::vector<float> out;
-            param.CopyToHost(out);
-            CheckScalar(out[0], 0.8990000632455132f, "AdamWStep", backendName, 1e-4f);
-        }
-    }
-}
-
-int main()
-{
-    printf("--- CPUBackend ---\n");
-    RunAllTests(DeviceType::DEVICE_CPU, "CPU");
-
-#ifdef DEEPITY_USE_CUDA
-    printf("\n--- CUDABackend ---\n");
-    RunAllTests(DeviceType::DEVICE_GPU, "CUDA");
-#else
-    printf("\n--- CUDABackend skipped (DEEPITY_USE_CUDA not defined) ---\n");
-#endif
-
-    if (g_failures > 0)
-    {
-        printf("\nFAILED: %d check(s) mismatched.\n", g_failures);
-        return 1;
-    }
-
-    printf("\nPASSED: all backend functions verified correct.\n");
-    return 0;
-}
\ No newline at end of file
diff --git a/tests/tCUDAIm2ColVerify.cpp b/tests/tCUDAIm2ColVerify.cpp
index e882de4..70ccbf4 100644
--- a/tests/tCUDAIm2ColVerify.cpp
+++ b/tests/tCUDAIm2ColVerify.cpp
@@ -1,110 +1,43 @@
 /**
- * @file tCUDAIm2ColVerify.cpp
- * @brief CPU-vs-GPU differential check for CUDABackend::Im2Col/Col2Im --
- * the two genuinely novel, hand-written CUDA kernels added tonight for
- * SimpleConvPCLayer's GPU port. Unlike the other four new backend
- * methods (MultiplyInto, Fill, RepackForBatchedGemm, AddBiasPerChannel,
- * all straightforward translations of existing, simple loops), these
- * two involve real per-thread index decomposition (row -> channel/kh/kw,
- * col -> oh/ow) with no prior verification at all -- "compiles" says
- * nothing about whether that arithmetic is actually correct on real
- * hardware. Same methodology as earlier tonight's CUDA differential
- * tests: identical input run through both backends, compared directly.
- *
- * Config: 3 channels, 6x6 input, 3x3 kernel, stride 1, pad 1 (output
- * stays 6x6) -- large enough to exercise real padding and multi-channel
- * behavior, small enough to run instantly.
+ * @file tCUDAIncrementCounterSyncVerify.cpp
+ * @brief Same as tCUDAIncrementCounterVerify.cpp, but with an explicit
+ * cudaDeviceSynchronize() inserted right after IncrementCounter(),
+ * before the readback. IncrementCounterKernel is launched as
+ * <<<1,1,0,stream>>> -- exactly one block, one thread -- structurally
+ * different from every other kernel tested tonight (all used a real
+ * blocks/BLOCK_SIZE calculation with many threads). Testing whether
+ * this specific, minimal launch configuration has a genuine
+ * stream-timing issue the synchronous CopyToHost isn't actually
+ * catching, despite that same pattern working for every other kernel.
  */
 #include <deepity/backend/Backend.h>
-#include <deepity/utils/Im2Col.h>
+#include <cuda_runtime.h>
 #include <cstdio>
-#include <cmath>
-#include <random>
-#include <vector>
 
 using namespace Deep;
 
 int main()
 {
-    const int channels = 3, height = 6, width = 6;
-    const int kH = 3, kW = 3, strideH = 1, strideW = 1, padH = 1, padW = 1;
-    const int outH = ConvOutDim(height, kH, strideH, padH);
-    const int outW = ConvOutDim(width, kW, strideW, padW);
-
-    printf("outH=%d outW=%d (expected 6, 6)\n", outH, outW);
-
-    const size_t inputSize = (size_t)channels * height * width;
-    const size_t colRows = (size_t)channels * kH * kW;
-    const size_t colCols = (size_t)outH * outW;
-    const size_t colSize = colRows * colCols;
-
-    std::mt19937 rng(42);
-    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
-
-    std::vector<float> hInput(inputSize);
-    for (auto &v : hInput)
-        v = dist(rng);
-
-    // --- Im2Col: CPU ---
-    auto cpuBackend = CreateBackend(DeviceType::DEVICE_CPU);
-    std::vector<float> colsCpu(colSize);
-    cpuBackend->Im2Col(hInput.data(), channels, height, width,
-                       kH, kW, strideH, strideW, padH, padW, colsCpu.data());
-
-    // --- Im2Col: GPU ---
     auto gpuBackend = CreateBackend(DeviceType::DEVICE_GPU);
-    float *dInput = gpuBackend->Allocate(inputSize);
-    float *dCols = gpuBackend->Allocate(colSize);
-    gpuBackend->CopyFromHost(dInput, hInput.data(), inputSize);
-    gpuBackend->Im2Col(dInput, channels, height, width,
-                       kH, kW, strideH, strideW, padH, padW, dCols);
-
-    std::vector<float> colsGpu(colSize);
-    gpuBackend->CopyToHost(colsGpu.data(), dCols, colSize);
-
-    float maxDiffIm2Col = 0.0f;
-    for (size_t i = 0; i < colSize; ++i)
-        maxDiffIm2Col = std::max(maxDiffIm2Col, std::fabs(colsCpu[i] - colsGpu[i]));
-
-    printf("Im2Col max abs diff: %g\n", maxDiffIm2Col);
-    bool im2colPass = maxDiffIm2Col < 1e-5f;
-    printf("Im2Col: %s\n", im2colPass ? "PASS" : "FAIL");
-
-    // --- Col2Im: CPU ---
-    std::vector<float> hColumns(colSize);
-    for (auto &v : hColumns)
-        v = dist(rng);
-
-    std::vector<float> outCpu(inputSize, 0.0f);
-    cpuBackend->Col2Im(hColumns.data(), channels, height, width,
-                       kH, kW, strideH, strideW, padH, padW, outCpu.data());
-
-    // --- Col2Im: GPU ---
-    float *dColumns = gpuBackend->Allocate(colSize);
-    float *dOutput = gpuBackend->Allocate(inputSize);
-    gpuBackend->CopyFromHost(dColumns, hColumns.data(), colSize);
-    gpuBackend->Zero(dOutput, inputSize); // Col2Im ACCUMULATES -- must zero first, matching CPU's own contract
 
-    gpuBackend->Col2Im(dColumns, channels, height, width,
-                       kH, kW, strideH, strideW, padH, padW, dOutput);
+    int *dT = reinterpret_cast<int *>(gpuBackend->Allocate(1));
+    int zero = 0;
+    gpuBackend->CopyFromHost(reinterpret_cast<float *>(dT), reinterpret_cast<float *>(&zero), 1);
 
-    std::vector<float> outGpu(inputSize);
-    gpuBackend->CopyToHost(outGpu.data(), dOutput, inputSize);
+    gpuBackend->IncrementCounter(dT);
 
-    float maxDiffCol2Im = 0.0f;
-    for (size_t i = 0; i < inputSize; ++i)
-        maxDiffCol2Im = std::max(maxDiffCol2Im, std::fabs(outCpu[i] - outGpu[i]));
+    // THE ONLY CHANGE from tCUDAIncrementCounterVerify.cpp:
+    cudaError_t syncErr = cudaDeviceSynchronize();
+    printf("cudaDeviceSynchronize() after IncrementCounter: %s\n", cudaGetErrorString(syncErr));
 
-    printf("Col2Im max abs diff: %g\n", maxDiffCol2Im);
-    bool col2imPass = maxDiffCol2Im < 1e-4f; // slightly looser -- atomicAdd accumulation order can differ from CPU's SIMD order
-    printf("Col2Im: %s\n", col2imPass ? "PASS" : "FAIL");
+    int readback = -999;
+    gpuBackend->CopyToHost(reinterpret_cast<float *>(&readback), reinterpret_cast<float *>(dT), 1);
+    printf("After ONE IncrementCounter call + explicit sync (expect 1): %d\n", readback);
 
-    gpuBackend->Free(dInput);
-    gpuBackend->Free(dCols);
-    gpuBackend->Free(dColumns);
-    gpuBackend->Free(dOutput);
+    bool pass = (readback == 1);
+    printf("\n%s\n", pass ? "PASS (explicit sync fixed it -- real stream-timing bug found)"
+                          : "FAIL (still wrong even with explicit sync -- bug is elsewhere)");
 
-    bool allPassed = im2colPass && col2imPass;
-    printf("\n%s\n", allPassed ? "PASS" : "FAIL");
-    return allPassed ? 0 : 1;
+    gpuBackend->Free(reinterpret_cast<float *>(dT));
+    return pass ? 0 : 1;
 }
\ No newline at end of file