From eb0bbf491a3fe670374f915ea0588aa701741997 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Thu, 23 Apr 2026 10:34:47 +0200
Subject: [PATCH 01/10] start fresh

---
 docs/Project.toml                             |   5 +
 .../tutorials/synthetic_respiration_gpu_2.jl  | 132 ++++++++++++++++++
 docs/setup_gpu_backend.jl                     |  50 +++++++
 3 files changed, 187 insertions(+)
 create mode 100644 docs/literate/tutorials/synthetic_respiration_gpu_2.jl
 create mode 100644 docs/setup_gpu_backend.jl

diff --git a/docs/Project.toml b/docs/Project.toml
index 9b98d244..bf0beca4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -10,9 +11,13 @@ EasyHybrid = "61bb816a-e6af-4913-ab9e-91bff2e122e3"
 Hyperopt = "93e5fe13-2215-51db-baaf-2e9a34fb2712"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
+LuxCUDA = "d0bbae9a-e099-4d5b-a835-1c6931763bda"
 MLDataDevices = "7e8f7934-dd98-4c1a-8fe8-92b47a384d40"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
 
 [sources]
 EasyHybrid = {path = ".."}
diff --git a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
new file mode 100644
index 00000000..35fa1d5d
--- /dev/null
+++ b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
@@ -0,0 +1,132 @@
+# [![CC BY-SA 4.0](https://img.shields.io/badge/License-CC%20BY--SA%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-sa/4.0/)
+#
+# # EasyHybrid Example: Synthetic Data Analysis
+#
+# This example demonstrates how to use EasyHybrid to train a hybrid model
+# on synthetic data for respiration modeling with Q10 temperature sensitivity.
+#
+
+include("../../setup_local_docsrun.jl")
+
+using EasyHybrid
+using MLDataDevices
+
+# Load whichever GPU backend is available on this machine (CUDA / AMDGPU /
+# Metal / oneAPI), or fall back to CPU. Defines `GPU_BACKEND_PKG` and
+# `GPU_DEVICE_TYPE` in the caller's scope.
+include("../../setup_gpu_backend.jl")
+
+cpu_device() isa CPUDevice
+gpu_device() isa GPU_DEVICE_TYPE
+
+# ## Data Loading and Preprocessing
+#
+# Load synthetic dataset from GitHub into DataFrame
+
+df = load_timeseries_netcdf("https://github.com/bask0/q10hybrid/raw/master/data/Synthetic4BookChap.nc");
+
+# Select a subset of data for faster execution
+#df = df[1:20000, :];
+#nothing #hide
+#first(df, 5)
+
+# ## Define the Physical Model
+#
+# **RbQ10 model**: Respiration model with Q10 temperature sensitivity
+#
+# Parameters:
+#   - ta: air temperature [°C]
+#   - Q10: temperature sensitivity factor [-]
+#   - rb: basal respiration rate [μmol/m²/s]
+#   - tref: reference temperature [°C] (default: 15.0)
+
+function RbQ10(; ta, Q10, rb, tref = 15.0f0)
+    reco = rb .* Q10 .^ (0.1f0 .* (ta .- tref))
+    return (; reco, Q10, rb)
+end
+
+# ### Define Model Parameters
+#
+# Parameter specification: (default, lower_bound, upper_bound)
+
+parameters = (
+    ## Parameter name | Default | Lower | Upper      | Description
+    rb = (3.0f0, 0.0f0, 13.0f0),  # Basal respiration [μmol/m²/s]
+    Q10 = (2.0f0, 1.0f0, 4.0f0),   # Temperature sensitivity factor [-]
+)
+
+# ## Configure Hybrid Model Components
+#
+# Define input variables
+forcing = [:ta]                    # Forcing variables (temperature)
+
+# Target variable
+target = [:reco]                   # Target variable (respiration)
+
+# Parameter classification
+global_param_names = [:Q10]        # Global parameters (same for all samples)
+neural_param_names = [:rb]         # Neural network predicted parameters
+
+# ## Single NN Hybrid Model Training
+predictors_single_nn = [:sw_pot, :dsw_pot]   # Predictor variables (solar radiation, and its derivative)
+
+small_nn_hybrid_model = constructHybridModel(
+    predictors_single_nn,              # Input features
+    forcing,                 # Forcing variables
+    target,                  # Target variables
+    RbQ10,                  # Process-based model function
+    parameters,              # Parameter definitions
+    neural_param_names,      # NN-predicted parameters
+    global_param_names,      # Global parameters
+    hidden_layers = [16, 16], # Neural network architecture
+    activation = sigmoid,      # Activation function
+    scale_nn_outputs = true, # Scale neural network outputs
+    input_batchnorm = true   # Apply batch normalization to inputs
+)
+
+large_nn_hybrid_model = constructHybridModel(
+    predictors_single_nn,              # Input features
+    forcing,                 # Forcing variables
+    target,                  # Target variables
+    RbQ10,                  # Process-based model function
+    parameters,              # Parameter definitions
+    neural_param_names,      # NN-predicted parameters
+    global_param_names,      # Global parameters
+    hidden_layers = [512, 512, 512], # Neural network architecture
+    activation = sigmoid,      # Activation function
+    scale_nn_outputs = true, # Scale neural network outputs
+    input_batchnorm = true   # Apply batch normalization to inputs
+)
+
+# ### train on DataFrame
+# Train the hybrid model
+
+cfg = EasyHybrid.TrainConfig(
+    nepochs = 10,
+    batchsize = 512,
+    opt = Adam(0.001),
+    loss_types = [:mse, :nse],
+    show_progress = false,
+    keep_history = false, # set to true to keep per-epoch history, losses, predictions, etc.
+    save_training = false, # Set to true to enable saving training history and checkpoints
+)
+
+using BenchmarkTools
+using Suppressor
+using Printf
+
+# Benchmark `tune` on CPU vs GPU and print the GPU speedup (or slowdown)
+# relative to CPU. Uses a single sample because each run already does a full
+# training loop (`cfg.nepochs` epochs), which is expensive.
+function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, tune_kwargs...)
+    t_cpu = @suppress @belapsed tune($model, $df, $cfg; gdev = cpu_device(), $tune_kwargs...) samples=1 evals=1
+    t_gpu = @suppress @belapsed tune($model, $df, $cfg; gdev = gpu_device(), $tune_kwargs...) samples=1 evals=1
+
+    ratio = t_cpu / t_gpu
+    @printf("%-10s | CPU: %8.3f s | GPU: %8.3f s | with GPU we get %.2fx\n",
+            label, t_cpu, t_gpu, ratio)
+    return (; t_cpu, t_gpu, ratio)
+end
+
+bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN") # on our gpu1-hpc22 0.57x
+bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN") # on our gpu1-hpc22 1.35x
diff --git a/docs/setup_gpu_backend.jl b/docs/setup_gpu_backend.jl
new file mode 100644
index 00000000..e67d659b
--- /dev/null
+++ b/docs/setup_gpu_backend.jl
@@ -0,0 +1,50 @@
+# Load the appropriate GPU backend for `MLDataDevices`, independent of machine.
+#
+# Tries the common GPU trigger packages in order:
+#   - LuxCUDA  -> CUDADevice   (NVIDIA)
+#   - AMDGPU   -> AMDGPUDevice (AMD ROCm)
+#   - Metal    -> MetalDevice  (Apple Silicon)
+#   - oneAPI   -> oneAPIDevice (Intel)
+#
+# The first backend that is both installed AND functional on this machine
+# is kept. If none are available, we silently fall back to CPU and
+# `gpu_device()` will return a `CPUDevice`.
+#
+# After `include`ing this file, the following are defined in the caller's
+# namespace:
+#   - `GPU_BACKEND_PKG`  : `Symbol` of the loaded package (or `nothing`)
+#   - `GPU_DEVICE_TYPE`  : the `MLDataDevices` device type (e.g. `CUDADevice`,
+#                          or `CPUDevice` as fallback)
+
+using MLDataDevices
+
+const _GPU_TRIGGERS = [
+    (:LuxCUDA, :CUDADevice),
+    (:AMDGPU,  :AMDGPUDevice),
+    (:Metal,   :MetalDevice),
+    (:oneAPI,  :oneAPIDevice),
+]
+
+function _load_gpu_backend()
+    for (pkg, dev) in _GPU_TRIGGERS
+        try
+            @eval Main using $pkg
+        catch err
+            @debug "Skipping $pkg (not installed / failed to load)" exception = err
+            continue
+        end
+
+        DeviceT = getfield(MLDataDevices, dev)
+        if MLDataDevices.functional(DeviceT)
+            @info "GPU backend loaded: $pkg → $dev"
+            return pkg, DeviceT
+        else
+            @warn "$pkg loaded but $dev is not functional on this machine; trying next backend."
+        end
+    end
+
+    @info "No functional GPU backend found; falling back to CPU."
+    return nothing, CPUDevice
+end
+
+const GPU_BACKEND_PKG, GPU_DEVICE_TYPE = _load_gpu_backend()

From 5ada65cf644b81a63411bec38d8afaf925d0b6e7 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Thu, 23 Apr 2026 12:12:44 +0200
Subject: [PATCH 02/10] yay speedup when we go wide

---
 .../tutorials/synthetic_respiration_gpu_2.jl  | 124 +++++++++++++++---
 docs/setup_gpu_backend.jl                     |  11 +-
 2 files changed, 118 insertions(+), 17 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
index 35fa1d5d..e0195545 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
@@ -6,7 +6,8 @@
 # on synthetic data for respiration modeling with Q10 temperature sensitivity.
 #
 
-include("../../setup_local_docsrun.jl")
+## for development, use local setup
+# include("../../setup_local_docsrun.jl")
 
 using EasyHybrid
 using MLDataDevices
@@ -109,24 +110,117 @@ cfg = EasyHybrid.TrainConfig(
     show_progress = false,
     keep_history = false, # set to true to keep per-epoch history, losses, predictions, etc.
     save_training = false, # Set to true to enable saving training history and checkpoints
+    plotting = false,
 )
 
-using BenchmarkTools
 using Suppressor
 using Printf
 
-# Benchmark `tune` on CPU vs GPU and print the GPU speedup (or slowdown)
-# relative to CPU. Uses a single sample because each run already does a full
-# training loop (`cfg.nepochs` epochs), which is expensive.
-function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, tune_kwargs...)
-    t_cpu = @suppress @belapsed tune($model, $df, $cfg; gdev = cpu_device(), $tune_kwargs...) samples=1 evals=1
-    t_gpu = @suppress @belapsed tune($model, $df, $cfg; gdev = gpu_device(), $tune_kwargs...) samples=1 evals=1
-
-    ratio = t_cpu / t_gpu
-    @printf("%-10s | CPU: %8.3f s | GPU: %8.3f s | with GPU we get %.2fx\n",
-            label, t_cpu, t_gpu, ratio)
-    return (; t_cpu, t_gpu, ratio)
+function _pure(s)
+    return s.time - s.gctime - s.compile_time - s.recompile_time
 end
 
-bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN") # on our gpu1-hpc22 0.57x
-bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN") # on our gpu1-hpc22 1.35x
+function _warn_overhead(tag, label, s, warn_frac)
+    gc_frac        = s.gctime         / s.time
+    compile_frac   = s.compile_time   / s.time
+    recompile_frac = s.recompile_time / s.time
+    overhead_frac  = gc_frac + compile_frac + recompile_frac
+    if overhead_frac > warn_frac
+        pct(x) = round(100 * x; digits = 1)
+        @warn """[$label/$tag] $(pct(overhead_frac))% of elapsed time was overhead:
+                 gc=$(pct(gc_frac))% ($(round(s.gctime; digits=3))s), \
+                 compile=$(pct(compile_frac))% ($(round(s.compile_time; digits=3))s), \
+                 recompile=$(pct(recompile_frac))% ($(round(s.recompile_time; digits=3))s).
+                 Elapsed ratio will be misleading; prefer the pure ratio."""
+    end
+end
+
+# Benchmark `tune` on CPU vs GPU. Uses `@timed` so we can separate out
+# compile/recompile/GC overhead from steady-state work, and reports BOTH:
+#   - elapsed ratio (wall clock CPU / wall clock GPU)
+#   - pure ratio    (elapsed minus GC/compile/recompile, CPU/GPU)
+#
+# `warn_frac` is the fraction of total time spent in GC / (re)compile above
+# which a warning is emitted for that run; defaults to 5%.
+function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool = true, warn_frac::Real = 0.05, tune_kwargs...)
+    # A short 1-epoch run on each device triggers Julia/Zygote specialization
+    # and GPU kernel compilation, so the timed runs below measure steady-state
+    # training rather than first-call compile latency.
+    if warmup
+        warm_cfg = EasyHybrid.TrainConfig(;
+            nepochs = 1,
+            batchsize = cfg.batchsize,
+            opt = cfg.opt,
+            show_progress = false,
+            save_training = false,
+            keep_history = false,
+            plotting = false,
+        )
+        @suppress tune(model, df, warm_cfg; gdev = cpu_device(), tune_kwargs...);
+        @suppress tune(model, df, warm_cfg; gdev = gpu_device(), tune_kwargs...);
+    end
+
+    cpu_stats = @suppress @timed tune(model, df, cfg; gdev = cpu_device(), tune_kwargs...);
+    gpu_stats = @suppress @timed tune(model, df, cfg; gdev = gpu_device(), tune_kwargs...);
+
+    _warn_overhead("CPU", label, cpu_stats, warn_frac)
+    _warn_overhead("GPU", label, gpu_stats, warn_frac)
+
+    cpu_pure = _pure(cpu_stats)
+    gpu_pure = _pure(gpu_stats)
+    elapsed_ratio = cpu_stats.time / gpu_stats.time
+    pure_ratio    = cpu_pure      / gpu_pure
+
+    @printf("[%s] CPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+            label, cpu_stats.time, cpu_stats.compile_time, cpu_stats.recompile_time, cpu_stats.gctime, cpu_pure)
+    @printf("[%s] GPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+            label, gpu_stats.time, gpu_stats.compile_time, gpu_stats.recompile_time, gpu_stats.gctime, gpu_pure)
+    @printf("[%s] with GPU we get: elapsed=%.2fx  pure=%.2fx\n",
+            label, elapsed_ratio, pure_ratio)
+
+    # Return a slim summary. We deliberately drop `.value` (the full TrainResults)
+    # and `.gcstats` so that REPL/Literate printing of the return value stays compact.
+    slim(s) = (; s.time, s.bytes, s.gctime, s.compile_time, s.recompile_time)
+    return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
+end
+
+bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN"); # on our gpu1-hpc22 0.57x
+bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our gpu1-hpc22 2.73x
+
+# ## Sweep: GPU speedup vs hidden-layer width
+#
+# Build a hybrid model for a range of hidden-layer widths (fixed depth of 3),
+# benchmark each on CPU and GPU, and plot the elapsed and pure ratios.
+# A ratio > 1 means GPU is faster than CPU; a ratio < 1 means GPU is slower.
+
+widths = [16, 64, 256, 512, 1024]
+
+results = map(widths) do w
+    r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
+        label = "w=$w", hidden_layers = [w, w, w])
+    (; width = w, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+       cpu_time = r.cpu.time, gpu_time = r.gpu.time)
+end
+
+using CairoMakie
+
+fig = Figure(size = (720, 480))
+ax = Axis(fig[1, 1];
+    xlabel = "Hidden-layer width (depth = 3)",
+    ylabel = "GPU speedup (CPU / GPU)",
+    title  = "GPU speedup vs hidden-layer width",
+    xscale = log2,
+    xticks = (widths, string.(widths)),
+)
+
+ws        = [r.width for r in results]
+elapsed_x = [r.elapsed_ratio for r in results]
+pure_x    = [r.pure_ratio    for r in results]
+
+hlines!(ax, [1.0]; color = :gray, linestyle = :dash, label = "CPU = GPU")
+scatterlines!(ax, ws, elapsed_x; label = "elapsed", marker = :circle)
+scatterlines!(ax, ws, pure_x;    label = "pure",    marker = :diamond)
+
+axislegend(ax; position = :lt)
+fig
+
diff --git a/docs/setup_gpu_backend.jl b/docs/setup_gpu_backend.jl
index e67d659b..679dfcf0 100644
--- a/docs/setup_gpu_backend.jl
+++ b/docs/setup_gpu_backend.jl
@@ -34,8 +34,15 @@ function _load_gpu_backend()
             continue
         end
 
-        DeviceT = getfield(MLDataDevices, dev)
-        if MLDataDevices.functional(DeviceT)
+        # `using $pkg` may load an `MLDataDevices` package extension that
+        # adds new methods (e.g. a real `functional(::Type{CUDADevice})`).
+        # Those methods live in a newer world age than this function, so
+        # we must dispatch via `invokelatest` or we'll get the stale
+        # (pre-extension) method that reports `false`.
+        DeviceT = Base.invokelatest(getfield, MLDataDevices, dev)
+        is_functional = Base.invokelatest(MLDataDevices.functional, DeviceT)
+
+        if is_functional
             @info "GPU backend loaded: $pkg → $dev"
             return pkg, DeviceT
         else

From 1e722f328062e8057203fbc860f0dee2aa3bcd06 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Thu, 23 Apr 2026 13:07:10 +0200
Subject: [PATCH 03/10] two factors depth and width

---
 .../tutorials/synthetic_respiration_gpu_2.jl  | 85 +++++++++++++++----
 1 file changed, 68 insertions(+), 17 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
index e0195545..acebaa15 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
@@ -194,33 +194,84 @@ bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our g
 # A ratio > 1 means GPU is faster than CPU; a ratio < 1 means GPU is slower.
 
 widths = [16, 64, 256, 512, 1024]
-
-results = map(widths) do w
-    r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
-        label = "w=$w", hidden_layers = [w, w, w])
-    (; width = w, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-       cpu_time = r.cpu.time, gpu_time = r.gpu.time)
-end
+depths = 2:5
+
+results = [
+    begin
+        r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
+            label = "d=$d w=$w", hidden_layers = fill(w, d))
+        (; depth = d, width = w,
+           elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+           cpu_time = r.cpu.time, gpu_time = r.gpu.time)
+    end
+    for d in depths for w in widths
+]
 
 using CairoMakie
 
-fig = Figure(size = (720, 480))
+fig = Figure(size = (820, 520))
 ax = Axis(fig[1, 1];
-    xlabel = "Hidden-layer width (depth = 3)",
+    xlabel = "Hidden-layer width",
     ylabel = "GPU speedup (CPU / GPU)",
-    title  = "GPU speedup vs hidden-layer width",
+    title  = "GPU speedup vs hidden-layer width, per depth",
     xscale = log2,
     xticks = (widths, string.(widths)),
 )
 
-ws        = [r.width for r in results]
-elapsed_x = [r.elapsed_ratio for r in results]
-pure_x    = [r.pure_ratio    for r in results]
+hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
+
+# One colour per depth; `pure` is solid, `elapsed` is dashed.
+palette = Makie.wong_colors()
+for (i, d) in enumerate(depths)
+    rs = filter(r -> r.depth == d, results)
+    ws = [r.width for r in rs]
+    el = [r.elapsed_ratio for r in rs]
+    pr = [r.pure_ratio    for r in rs]
+    c  = palette[mod1(i, length(palette))]
+    # Only the solid (pure) line gets a label → legend entries are per-depth.
+    scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
+    scatterlines!(ax, ws, el; color = c, linestyle = :dash,  marker = :xcross)
+end
+
+# Legend 1: depth → colour (auto-picked up from the labeled solid lines).
+axislegend(ax, "depth"; position = :lt)
 
-hlines!(ax, [1.0]; color = :gray, linestyle = :dash, label = "CPU = GPU")
-scatterlines!(ax, ws, elapsed_x; label = "elapsed", marker = :circle)
-scatterlines!(ax, ws, pure_x;    label = "pure",    marker = :diamond)
+# Legend 2: line style → metric (built manually).
+style_elements = [
+    LineElement(color = :black, linestyle = :solid),
+    LineElement(color = :black, linestyle = :dash),
+]
+Legend(fig[1, 2], style_elements, ["pure", "elapsed"], "metric")
 
-axislegend(ax; position = :lt)
 fig
 
+
+batchsizes = [64, 128, 256, 512, 1024, 2048, 4096, 8192]
+bs_hidden  = [64, 64, 64]
+
+bs_results = map(batchsizes) do b
+    r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
+        label = "b=$b", batchsize = b, hidden_layers = bs_hidden)
+    (; batchsize = b, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+       cpu_time = r.cpu.time, gpu_time = r.gpu.time)
+end
+
+fig_bs = Figure(size = (720, 480))
+ax_bs  = Axis(fig_bs[1, 1];
+    xlabel = "Batch size (hidden = $bs_hidden)",
+    ylabel = "GPU speedup (CPU / GPU)",
+    title  = "GPU speedup vs batch size",
+    xscale = log2,
+    xticks = (batchsizes, string.(batchsizes)),
+)
+
+bs         = [r.batchsize     for r in bs_results]
+bs_elapsed = [r.elapsed_ratio for r in bs_results]
+bs_pure    = [r.pure_ratio    for r in bs_results]
+
+hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
+scatterlines!(ax_bs, bs, bs_pure;    color = :steelblue, linestyle = :solid, marker = :circle,  label = "pure")
+scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash,  marker = :xcross, label = "elapsed")
+
+axislegend(ax_bs; position = :lt)
+fig_bs
\ No newline at end of file

From 8f5192fe5c44231b6621aceb5a47d6a0b10c8404 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Thu, 23 Apr 2026 13:13:06 +0200
Subject: [PATCH 04/10] local dev setup

---
 .gitignore                  | 2 --
 docs/setup_local_docsrun.jl | 4 ++++
 2 files changed, 4 insertions(+), 2 deletions(-)
 create mode 100644 docs/setup_local_docsrun.jl

diff --git a/.gitignore b/.gitignore
index 0aa43f4e..35c143e8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,5 +32,3 @@ coverage/lcov.info
 !MPI_BGC_wide_E_neg_cmyk.png
 !logo.png
 !logo_sq.png
-test/setup_local_testrun.jl
-docs/setup_local_docsrun.jl
diff --git a/docs/setup_local_docsrun.jl b/docs/setup_local_docsrun.jl
new file mode 100644
index 00000000..75bfdfa8
--- /dev/null
+++ b/docs/setup_local_docsrun.jl
@@ -0,0 +1,4 @@
+using Pkg
+Pkg.activate("docs")
+Pkg.develop(path = pwd())
+Pkg.instantiate()

From f55824d57ee0c83e2c9d7cd436fbef37cc567e32 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Thu, 23 Apr 2026 15:15:13 +0200
Subject: [PATCH 05/10] works with Metal?

---
 .../tutorials/synthetic_respiration_gpu.jl    | 223 +++++++++++---
 .../tutorials/synthetic_respiration_gpu_2.jl  | 277 ------------------
 2 files changed, 177 insertions(+), 323 deletions(-)
 delete mode 100644 docs/literate/tutorials/synthetic_respiration_gpu_2.jl

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index 88479d78..44a17d1c 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -6,10 +6,19 @@
 # on synthetic data for respiration modeling with Q10 temperature sensitivity.
 #
 
+## for development, use local setup
+# include("../../setup_local_docsrun.jl")
+
 using EasyHybrid
-using Metal
-## using CUDA
-using Lux
+using MLDataDevices
+
+# Load whichever GPU backend is available on this machine (CUDA / AMDGPU /
+# Metal / oneAPI), or fall back to CPU. Defines `GPU_BACKEND_PKG` and
+# `GPU_DEVICE_TYPE` in the caller's scope.
+include("../../setup_gpu_backend.jl")
+
+cpu_device() isa CPUDevice
+gpu_device() isa GPU_DEVICE_TYPE
 
 # ## Data Loading and Preprocessing
 #
@@ -18,9 +27,9 @@ using Lux
 df = load_timeseries_netcdf("https://github.com/bask0/q10hybrid/raw/master/data/Synthetic4BookChap.nc");
 
 # Select a subset of data for faster execution
-df = df[1:5000, :];
-nothing #hide
-first(df, 5)
+#df = df[1:20000, :];
+#nothing #hide
+#first(df, 5)
 
 # ## Define the Physical Model
 #
@@ -84,7 +93,7 @@ large_nn_hybrid_model = constructHybridModel(
     parameters,              # Parameter definitions
     neural_param_names,      # NN-predicted parameters
     global_param_names,      # Global parameters
-    hidden_layers = [1024, 512, 256, 128, 64], # Neural network architecture
+    hidden_layers = [512, 512, 512], # Neural network architecture
     activation = sigmoid,      # Activation function
     scale_nn_outputs = true, # Scale neural network outputs
     input_batchnorm = true   # Apply batch normalization to inputs
@@ -94,53 +103,175 @@ large_nn_hybrid_model = constructHybridModel(
 # Train the hybrid model
 
 cfg = EasyHybrid.TrainConfig(
-    nepochs = 20,
-    batchsize = 64,
-    opt = RMSProp(0.01),
+    nepochs = 10,
+    batchsize = 512,
+    opt = Adam(0.001),
     loss_types = [:mse, :nse],
     show_progress = false,
     keep_history = false, # set to true to keep per-epoch history, losses, predictions, etc.
     save_training = false, # Set to true to enable saving training history and checkpoints
+    plotting = false,
 )
 
-using BenchmarkTools
+using Suppressor
+using Printf
 
-gpu_small_nn() = tune(
-    small_nn_hybrid_model, df, cfg;
-    gdev = gpu_device(), model_name = "small_nn_gpu"
-)
+function _pure(s)
+    return s.time - s.gctime - s.compile_time - s.recompile_time
+end
 
-cpu_small_nn() = tune(
-    small_nn_hybrid_model, df, cfg;
-    gdev = cpu_device(), model_name = "small_nn_cpu"
-)
+function _warn_overhead(tag, label, s, warn_frac)
+    gc_frac        = s.gctime         / s.time
+    compile_frac   = s.compile_time   / s.time
+    recompile_frac = s.recompile_time / s.time
+    overhead_frac  = gc_frac + compile_frac + recompile_frac
+    if overhead_frac > warn_frac
+        pct(x) = round(100 * x; digits = 1)
+        @warn """[$label/$tag] $(pct(overhead_frac))% of elapsed time was overhead:
+                 gc=$(pct(gc_frac))% ($(round(s.gctime; digits=3))s), \
+                 compile=$(pct(compile_frac))% ($(round(s.compile_time; digits=3))s), \
+                 recompile=$(pct(recompile_frac))% ($(round(s.recompile_time; digits=3))s).
+                 Elapsed ratio will be misleading; prefer the pure ratio."""
+    end
+end
 
-gpu_large_nn() = tune(
-    large_nn_hybrid_model, df, cfg;
-    gdev = gpu_device(), model_name = "large_nn_gpu"
-)
+# Benchmark `tune` on CPU vs GPU. Uses `@timed` so we can separate out
+# compile/recompile/GC overhead from steady-state work, and reports BOTH:
+#   - elapsed ratio (wall clock CPU / wall clock GPU)
+#   - pure ratio    (elapsed minus GC/compile/recompile, CPU/GPU)
+#
+# `warn_frac` is the fraction of total time spent in GC / (re)compile above
+# which a warning is emitted for that run; defaults to 5%.
+function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool = true, warn_frac::Real = 0.05, tune_kwargs...)
+    # A short 1-epoch run on each device triggers Julia/Zygote specialization
+    # and GPU kernel compilation, so the timed runs below measure steady-state
+    # training rather than first-call compile latency.
+    if warmup
+        warm_cfg = EasyHybrid.TrainConfig(;
+            nepochs = 1,
+            batchsize = cfg.batchsize,
+            opt = cfg.opt,
+            show_progress = false,
+            save_training = false,
+            keep_history = false,
+            plotting = false,
+        )
+        @suppress tune(model, df, warm_cfg; gdev = cpu_device(), tune_kwargs...);
+        @suppress tune(model, df, warm_cfg; gdev = gpu_device(), tune_kwargs...);
+    end
 
-cpu_large_nn() = tune(
-    large_nn_hybrid_model, df, cfg;
-    gdev = cpu_device(), model_name = "large_nn_cpu"
-)
+    cpu_stats = @suppress @timed tune(model, df, cfg; gdev = cpu_device(), tune_kwargs...);
+    gpu_stats = @suppress @timed tune(model, df, cfg; gdev = gpu_device(), tune_kwargs...);
+
+    _warn_overhead("CPU", label, cpu_stats, warn_frac)
+    _warn_overhead("GPU", label, gpu_stats, warn_frac)
+
+    cpu_pure = _pure(cpu_stats)
+    gpu_pure = _pure(gpu_stats)
+    elapsed_ratio = cpu_stats.time / gpu_stats.time
+    pure_ratio    = cpu_pure      / gpu_pure
+
+    @printf("[%s] CPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+            label, cpu_stats.time, cpu_stats.compile_time, cpu_stats.recompile_time, cpu_stats.gctime, cpu_pure)
+    @printf("[%s] GPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+            label, gpu_stats.time, gpu_stats.compile_time, gpu_stats.recompile_time, gpu_stats.gctime, gpu_pure)
+    @printf("[%s] with GPU we get: elapsed=%.2fx  pure=%.2fx\n",
+            label, elapsed_ratio, pure_ratio)
+
+    # Return a slim summary. We deliberately drop `.value` (the full TrainResults)
+    # and `.gcstats` so that REPL/Literate printing of the return value stays compact.
+    slim(s) = (; s.time, s.bytes, s.gctime, s.compile_time, s.recompile_time)
+    return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
+end
+
+bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN"); # on our gpu1-hpc22 0.57x
+bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our gpu1-hpc22 2.73x
+
+# ## Sweep: GPU speedup vs hidden-layer width
+#
+# Build a hybrid model for a range of hidden-layer widths (fixed depth of 3),
+# benchmark each on CPU and GPU, and plot the elapsed and pure ratios.
+# A ratio > 1 means GPU is faster than CPU; a ratio < 1 means GPU is slower.
+
+# widths = [16, 64, 256, 512, 1024]
+# depths = 2:5
+
+# results = [
+#     begin
+#         r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
+#             label = "d=$d w=$w", hidden_layers = fill(w, d))
+#         (; depth = d, width = w,
+#            elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+#            cpu_time = r.cpu.time, gpu_time = r.gpu.time)
+#     end
+#     for d in depths for w in widths
+# ]
+
+# using CairoMakie
+
+# fig = Figure(size = (820, 520))
+# ax = Axis(fig[1, 1];
+#     xlabel = "Hidden-layer width",
+#     ylabel = "GPU speedup (CPU / GPU)",
+#     title  = "GPU speedup vs hidden-layer width, per depth",
+#     xscale = log2,
+#     xticks = (widths, string.(widths)),
+# )
+
+# hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
+
+# # One colour per depth; `pure` is solid, `elapsed` is dashed.
+# palette = Makie.wong_colors()
+# for (i, d) in enumerate(depths)
+#     rs = filter(r -> r.depth == d, results)
+#     ws = [r.width for r in rs]
+#     el = [r.elapsed_ratio for r in rs]
+#     pr = [r.pure_ratio    for r in rs]
+#     c  = palette[mod1(i, length(palette))]
+#     # Only the solid (pure) line gets a label → legend entries are per-depth.
+#     scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
+#     scatterlines!(ax, ws, el; color = c, linestyle = :dash,  marker = :xcross)
+# end
+
+# # Legend 1: depth → colour (auto-picked up from the labeled solid lines).
+# axislegend(ax, "depth"; position = :lt)
+
+# # Legend 2: line style → metric (built manually).
+# style_elements = [
+#     LineElement(color = :black, linestyle = :solid),
+#     LineElement(color = :black, linestyle = :dash),
+# ]
+# Legend(fig[1, 2], style_elements, ["pure", "elapsed"], "metric")
+
+# fig
+
+
+# batchsizes = [64, 128, 256, 512, 1024, 2048, 4096, 8192]
+# bs_hidden  = [64, 64, 64]
+
+# bs_results = map(batchsizes) do b
+#     r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
+#         label = "b=$b", batchsize = b, hidden_layers = bs_hidden)
+#     (; batchsize = b, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+#        cpu_time = r.cpu.time, gpu_time = r.gpu.time)
+# end
+
+# fig_bs = Figure(size = (720, 480))
+# ax_bs  = Axis(fig_bs[1, 1];
+#     xlabel = "Batch size (hidden = $bs_hidden)",
+#     ylabel = "GPU speedup (CPU / GPU)",
+#     title  = "GPU speedup vs batch size",
+#     xscale = log2,
+#     xticks = (batchsizes, string.(batchsizes)),
+# )
+
+# bs         = [r.batchsize     for r in bs_results]
+# bs_elapsed = [r.elapsed_ratio for r in bs_results]
+# bs_pure    = [r.pure_ratio    for r in bs_results]
+
+# hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
+# scatterlines!(ax_bs, bs, bs_pure;    color = :steelblue, linestyle = :solid, marker = :circle,  label = "pure")
+# scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash,  marker = :xcross, label = "elapsed")
 
-# warm-up to pay compilation once
-gpu_small_nn();
-cpu_small_nn();
-gpu_large_nn();
-cpu_large_nn();
-nothing # hide
-
-# ## With Large NN CPU is slower than GPU
-# Large NN on GPU
-@benchmark gpu_large_nn() samples = 4 evals = 1
-# Large NN on CPU
-@benchmark cpu_large_nn() samples = 4 evals = 1
-
-# ## With Small NN CPU and GPU are on par
-using BenchmarkTools
-# Small NN on GPU
-@benchmark gpu_small_nn() samples = 4 evals = 1
-# Small NN on CPU
-@benchmark cpu_small_nn() samples = 4 evals = 1
+# axislegend(ax_bs; position = :lt)
+# fig_bs
\ No newline at end of file
diff --git a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl b/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
deleted file mode 100644
index acebaa15..00000000
--- a/docs/literate/tutorials/synthetic_respiration_gpu_2.jl
+++ /dev/null
@@ -1,277 +0,0 @@
-# [![CC BY-SA 4.0](https://img.shields.io/badge/License-CC%20BY--SA%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-sa/4.0/)
-#
-# # EasyHybrid Example: Synthetic Data Analysis
-#
-# This example demonstrates how to use EasyHybrid to train a hybrid model
-# on synthetic data for respiration modeling with Q10 temperature sensitivity.
-#
-
-## for development, use local setup
-# include("../../setup_local_docsrun.jl")
-
-using EasyHybrid
-using MLDataDevices
-
-# Load whichever GPU backend is available on this machine (CUDA / AMDGPU /
-# Metal / oneAPI), or fall back to CPU. Defines `GPU_BACKEND_PKG` and
-# `GPU_DEVICE_TYPE` in the caller's scope.
-include("../../setup_gpu_backend.jl")
-
-cpu_device() isa CPUDevice
-gpu_device() isa GPU_DEVICE_TYPE
-
-# ## Data Loading and Preprocessing
-#
-# Load synthetic dataset from GitHub into DataFrame
-
-df = load_timeseries_netcdf("https://github.com/bask0/q10hybrid/raw/master/data/Synthetic4BookChap.nc");
-
-# Select a subset of data for faster execution
-#df = df[1:20000, :];
-#nothing #hide
-#first(df, 5)
-
-# ## Define the Physical Model
-#
-# **RbQ10 model**: Respiration model with Q10 temperature sensitivity
-#
-# Parameters:
-#   - ta: air temperature [°C]
-#   - Q10: temperature sensitivity factor [-]
-#   - rb: basal respiration rate [μmol/m²/s]
-#   - tref: reference temperature [°C] (default: 15.0)
-
-function RbQ10(; ta, Q10, rb, tref = 15.0f0)
-    reco = rb .* Q10 .^ (0.1f0 .* (ta .- tref))
-    return (; reco, Q10, rb)
-end
-
-# ### Define Model Parameters
-#
-# Parameter specification: (default, lower_bound, upper_bound)
-
-parameters = (
-    ## Parameter name | Default | Lower | Upper      | Description
-    rb = (3.0f0, 0.0f0, 13.0f0),  # Basal respiration [μmol/m²/s]
-    Q10 = (2.0f0, 1.0f0, 4.0f0),   # Temperature sensitivity factor [-]
-)
-
-# ## Configure Hybrid Model Components
-#
-# Define input variables
-forcing = [:ta]                    # Forcing variables (temperature)
-
-# Target variable
-target = [:reco]                   # Target variable (respiration)
-
-# Parameter classification
-global_param_names = [:Q10]        # Global parameters (same for all samples)
-neural_param_names = [:rb]         # Neural network predicted parameters
-
-# ## Single NN Hybrid Model Training
-predictors_single_nn = [:sw_pot, :dsw_pot]   # Predictor variables (solar radiation, and its derivative)
-
-small_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn,              # Input features
-    forcing,                 # Forcing variables
-    target,                  # Target variables
-    RbQ10,                  # Process-based model function
-    parameters,              # Parameter definitions
-    neural_param_names,      # NN-predicted parameters
-    global_param_names,      # Global parameters
-    hidden_layers = [16, 16], # Neural network architecture
-    activation = sigmoid,      # Activation function
-    scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true   # Apply batch normalization to inputs
-)
-
-large_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn,              # Input features
-    forcing,                 # Forcing variables
-    target,                  # Target variables
-    RbQ10,                  # Process-based model function
-    parameters,              # Parameter definitions
-    neural_param_names,      # NN-predicted parameters
-    global_param_names,      # Global parameters
-    hidden_layers = [512, 512, 512], # Neural network architecture
-    activation = sigmoid,      # Activation function
-    scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true   # Apply batch normalization to inputs
-)
-
-# ### train on DataFrame
-# Train the hybrid model
-
-cfg = EasyHybrid.TrainConfig(
-    nepochs = 10,
-    batchsize = 512,
-    opt = Adam(0.001),
-    loss_types = [:mse, :nse],
-    show_progress = false,
-    keep_history = false, # set to true to keep per-epoch history, losses, predictions, etc.
-    save_training = false, # Set to true to enable saving training history and checkpoints
-    plotting = false,
-)
-
-using Suppressor
-using Printf
-
-function _pure(s)
-    return s.time - s.gctime - s.compile_time - s.recompile_time
-end
-
-function _warn_overhead(tag, label, s, warn_frac)
-    gc_frac        = s.gctime         / s.time
-    compile_frac   = s.compile_time   / s.time
-    recompile_frac = s.recompile_time / s.time
-    overhead_frac  = gc_frac + compile_frac + recompile_frac
-    if overhead_frac > warn_frac
-        pct(x) = round(100 * x; digits = 1)
-        @warn """[$label/$tag] $(pct(overhead_frac))% of elapsed time was overhead:
-                 gc=$(pct(gc_frac))% ($(round(s.gctime; digits=3))s), \
-                 compile=$(pct(compile_frac))% ($(round(s.compile_time; digits=3))s), \
-                 recompile=$(pct(recompile_frac))% ($(round(s.recompile_time; digits=3))s).
-                 Elapsed ratio will be misleading; prefer the pure ratio."""
-    end
-end
-
-# Benchmark `tune` on CPU vs GPU. Uses `@timed` so we can separate out
-# compile/recompile/GC overhead from steady-state work, and reports BOTH:
-#   - elapsed ratio (wall clock CPU / wall clock GPU)
-#   - pure ratio    (elapsed minus GC/compile/recompile, CPU/GPU)
-#
-# `warn_frac` is the fraction of total time spent in GC / (re)compile above
-# which a warning is emitted for that run; defaults to 5%.
-function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool = true, warn_frac::Real = 0.05, tune_kwargs...)
-    # A short 1-epoch run on each device triggers Julia/Zygote specialization
-    # and GPU kernel compilation, so the timed runs below measure steady-state
-    # training rather than first-call compile latency.
-    if warmup
-        warm_cfg = EasyHybrid.TrainConfig(;
-            nepochs = 1,
-            batchsize = cfg.batchsize,
-            opt = cfg.opt,
-            show_progress = false,
-            save_training = false,
-            keep_history = false,
-            plotting = false,
-        )
-        @suppress tune(model, df, warm_cfg; gdev = cpu_device(), tune_kwargs...);
-        @suppress tune(model, df, warm_cfg; gdev = gpu_device(), tune_kwargs...);
-    end
-
-    cpu_stats = @suppress @timed tune(model, df, cfg; gdev = cpu_device(), tune_kwargs...);
-    gpu_stats = @suppress @timed tune(model, df, cfg; gdev = gpu_device(), tune_kwargs...);
-
-    _warn_overhead("CPU", label, cpu_stats, warn_frac)
-    _warn_overhead("GPU", label, gpu_stats, warn_frac)
-
-    cpu_pure = _pure(cpu_stats)
-    gpu_pure = _pure(gpu_stats)
-    elapsed_ratio = cpu_stats.time / gpu_stats.time
-    pure_ratio    = cpu_pure      / gpu_pure
-
-    @printf("[%s] CPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
-            label, cpu_stats.time, cpu_stats.compile_time, cpu_stats.recompile_time, cpu_stats.gctime, cpu_pure)
-    @printf("[%s] GPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
-            label, gpu_stats.time, gpu_stats.compile_time, gpu_stats.recompile_time, gpu_stats.gctime, gpu_pure)
-    @printf("[%s] with GPU we get: elapsed=%.2fx  pure=%.2fx\n",
-            label, elapsed_ratio, pure_ratio)
-
-    # Return a slim summary. We deliberately drop `.value` (the full TrainResults)
-    # and `.gcstats` so that REPL/Literate printing of the return value stays compact.
-    slim(s) = (; s.time, s.bytes, s.gctime, s.compile_time, s.recompile_time)
-    return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
-end
-
-bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN"); # on our gpu1-hpc22 0.57x
-bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our gpu1-hpc22 2.73x
-
-# ## Sweep: GPU speedup vs hidden-layer width
-#
-# Build a hybrid model for a range of hidden-layer widths (fixed depth of 3),
-# benchmark each on CPU and GPU, and plot the elapsed and pure ratios.
-# A ratio > 1 means GPU is faster than CPU; a ratio < 1 means GPU is slower.
-
-widths = [16, 64, 256, 512, 1024]
-depths = 2:5
-
-results = [
-    begin
-        r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
-            label = "d=$d w=$w", hidden_layers = fill(w, d))
-        (; depth = d, width = w,
-           elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-           cpu_time = r.cpu.time, gpu_time = r.gpu.time)
-    end
-    for d in depths for w in widths
-]
-
-using CairoMakie
-
-fig = Figure(size = (820, 520))
-ax = Axis(fig[1, 1];
-    xlabel = "Hidden-layer width",
-    ylabel = "GPU speedup (CPU / GPU)",
-    title  = "GPU speedup vs hidden-layer width, per depth",
-    xscale = log2,
-    xticks = (widths, string.(widths)),
-)
-
-hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
-
-# One colour per depth; `pure` is solid, `elapsed` is dashed.
-palette = Makie.wong_colors()
-for (i, d) in enumerate(depths)
-    rs = filter(r -> r.depth == d, results)
-    ws = [r.width for r in rs]
-    el = [r.elapsed_ratio for r in rs]
-    pr = [r.pure_ratio    for r in rs]
-    c  = palette[mod1(i, length(palette))]
-    # Only the solid (pure) line gets a label → legend entries are per-depth.
-    scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
-    scatterlines!(ax, ws, el; color = c, linestyle = :dash,  marker = :xcross)
-end
-
-# Legend 1: depth → colour (auto-picked up from the labeled solid lines).
-axislegend(ax, "depth"; position = :lt)
-
-# Legend 2: line style → metric (built manually).
-style_elements = [
-    LineElement(color = :black, linestyle = :solid),
-    LineElement(color = :black, linestyle = :dash),
-]
-Legend(fig[1, 2], style_elements, ["pure", "elapsed"], "metric")
-
-fig
-
-
-batchsizes = [64, 128, 256, 512, 1024, 2048, 4096, 8192]
-bs_hidden  = [64, 64, 64]
-
-bs_results = map(batchsizes) do b
-    r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
-        label = "b=$b", batchsize = b, hidden_layers = bs_hidden)
-    (; batchsize = b, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-       cpu_time = r.cpu.time, gpu_time = r.gpu.time)
-end
-
-fig_bs = Figure(size = (720, 480))
-ax_bs  = Axis(fig_bs[1, 1];
-    xlabel = "Batch size (hidden = $bs_hidden)",
-    ylabel = "GPU speedup (CPU / GPU)",
-    title  = "GPU speedup vs batch size",
-    xscale = log2,
-    xticks = (batchsizes, string.(batchsizes)),
-)
-
-bs         = [r.batchsize     for r in bs_results]
-bs_elapsed = [r.elapsed_ratio for r in bs_results]
-bs_pure    = [r.pure_ratio    for r in bs_results]
-
-hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
-scatterlines!(ax_bs, bs, bs_pure;    color = :steelblue, linestyle = :solid, marker = :circle,  label = "pure")
-scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash,  marker = :xcross, label = "elapsed")
-
-axislegend(ax_bs; position = :lt)
-fig_bs
\ No newline at end of file

From f4cbef85c28112a8441c907d89c5de5eb6ea1bba Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Wed, 29 Apr 2026 14:28:50 +0200
Subject: [PATCH 06/10] runic

---
 .../tutorials/synthetic_respiration_gpu.jl    | 220 ++++++++++--------
 docs/setup_gpu_backend.jl                     |   6 +-
 2 files changed, 121 insertions(+), 105 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index 44a17d1c..3c23ff1e 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -121,17 +121,17 @@ function _pure(s)
 end
 
 function _warn_overhead(tag, label, s, warn_frac)
-    gc_frac        = s.gctime         / s.time
-    compile_frac   = s.compile_time   / s.time
+    gc_frac = s.gctime / s.time
+    compile_frac = s.compile_time / s.time
     recompile_frac = s.recompile_time / s.time
-    overhead_frac  = gc_frac + compile_frac + recompile_frac
-    if overhead_frac > warn_frac
+    overhead_frac = gc_frac + compile_frac + recompile_frac
+    return if overhead_frac > warn_frac
         pct(x) = round(100 * x; digits = 1)
         @warn """[$label/$tag] $(pct(overhead_frac))% of elapsed time was overhead:
-                 gc=$(pct(gc_frac))% ($(round(s.gctime; digits=3))s), \
-                 compile=$(pct(compile_frac))% ($(round(s.compile_time; digits=3))s), \
-                 recompile=$(pct(recompile_frac))% ($(round(s.recompile_time; digits=3))s).
-                 Elapsed ratio will be misleading; prefer the pure ratio."""
+        gc=$(pct(gc_frac))% ($(round(s.gctime; digits = 3))s), \
+        compile=$(pct(compile_frac))% ($(round(s.compile_time; digits = 3))s), \
+        recompile=$(pct(recompile_frac))% ($(round(s.recompile_time; digits = 3))s).
+        Elapsed ratio will be misleading; prefer the pure ratio."""
     end
 end
 
@@ -156,12 +156,12 @@ function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool =
             keep_history = false,
             plotting = false,
         )
-        @suppress tune(model, df, warm_cfg; gdev = cpu_device(), tune_kwargs...);
-        @suppress tune(model, df, warm_cfg; gdev = gpu_device(), tune_kwargs...);
+        @suppress tune(model, df, warm_cfg; gdev = cpu_device(), tune_kwargs...)
+        @suppress tune(model, df, warm_cfg; gdev = gpu_device(), tune_kwargs...)
     end
 
-    cpu_stats = @suppress @timed tune(model, df, cfg; gdev = cpu_device(), tune_kwargs...);
-    gpu_stats = @suppress @timed tune(model, df, cfg; gdev = gpu_device(), tune_kwargs...);
+    cpu_stats = @suppress @timed tune(model, df, cfg; gdev = cpu_device(), tune_kwargs...)
+    gpu_stats = @suppress @timed tune(model, df, cfg; gdev = gpu_device(), tune_kwargs...)
 
     _warn_overhead("CPU", label, cpu_stats, warn_frac)
     _warn_overhead("GPU", label, gpu_stats, warn_frac)
@@ -169,14 +169,20 @@ function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool =
     cpu_pure = _pure(cpu_stats)
     gpu_pure = _pure(gpu_stats)
     elapsed_ratio = cpu_stats.time / gpu_stats.time
-    pure_ratio    = cpu_pure      / gpu_pure
-
-    @printf("[%s] CPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
-            label, cpu_stats.time, cpu_stats.compile_time, cpu_stats.recompile_time, cpu_stats.gctime, cpu_pure)
-    @printf("[%s] GPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
-            label, gpu_stats.time, gpu_stats.compile_time, gpu_stats.recompile_time, gpu_stats.gctime, gpu_pure)
-    @printf("[%s] with GPU we get: elapsed=%.2fx  pure=%.2fx\n",
-            label, elapsed_ratio, pure_ratio)
+    pure_ratio = cpu_pure / gpu_pure
+
+    @printf(
+        "[%s] CPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+        label, cpu_stats.time, cpu_stats.compile_time, cpu_stats.recompile_time, cpu_stats.gctime, cpu_pure
+    )
+    @printf(
+        "[%s] GPU: total=%.3fs compile=%.3fs recompile=%.3fs gc=%.3fs pure=%.3fs\n",
+        label, gpu_stats.time, gpu_stats.compile_time, gpu_stats.recompile_time, gpu_stats.gctime, gpu_pure
+    )
+    @printf(
+        "[%s] with GPU we get: elapsed=%.2fx  pure=%.2fx\n",
+        label, elapsed_ratio, pure_ratio
+    )
 
     # Return a slim summary. We deliberately drop `.value` (the full TrainResults)
     # and `.gcstats` so that REPL/Literate printing of the return value stays compact.
@@ -193,85 +199,95 @@ bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our g
 # benchmark each on CPU and GPU, and plot the elapsed and pure ratios.
 # A ratio > 1 means GPU is faster than CPU; a ratio < 1 means GPU is slower.
 
-# widths = [16, 64, 256, 512, 1024]
-# depths = 2:5
-
-# results = [
-#     begin
-#         r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
-#             label = "d=$d w=$w", hidden_layers = fill(w, d))
-#         (; depth = d, width = w,
-#            elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-#            cpu_time = r.cpu.time, gpu_time = r.gpu.time)
-#     end
-#     for d in depths for w in widths
-# ]
-
-# using CairoMakie
-
-# fig = Figure(size = (820, 520))
-# ax = Axis(fig[1, 1];
-#     xlabel = "Hidden-layer width",
-#     ylabel = "GPU speedup (CPU / GPU)",
-#     title  = "GPU speedup vs hidden-layer width, per depth",
-#     xscale = log2,
-#     xticks = (widths, string.(widths)),
-# )
-
-# hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
-
-# # One colour per depth; `pure` is solid, `elapsed` is dashed.
-# palette = Makie.wong_colors()
-# for (i, d) in enumerate(depths)
-#     rs = filter(r -> r.depth == d, results)
-#     ws = [r.width for r in rs]
-#     el = [r.elapsed_ratio for r in rs]
-#     pr = [r.pure_ratio    for r in rs]
-#     c  = palette[mod1(i, length(palette))]
-#     # Only the solid (pure) line gets a label → legend entries are per-depth.
-#     scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
-#     scatterlines!(ax, ws, el; color = c, linestyle = :dash,  marker = :xcross)
-# end
-
-# # Legend 1: depth → colour (auto-picked up from the labeled solid lines).
-# axislegend(ax, "depth"; position = :lt)
-
-# # Legend 2: line style → metric (built manually).
-# style_elements = [
-#     LineElement(color = :black, linestyle = :solid),
-#     LineElement(color = :black, linestyle = :dash),
-# ]
-# Legend(fig[1, 2], style_elements, ["pure", "elapsed"], "metric")
-
-# fig
-
-
-# batchsizes = [64, 128, 256, 512, 1024, 2048, 4096, 8192]
-# bs_hidden  = [64, 64, 64]
-
-# bs_results = map(batchsizes) do b
-#     r = bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg;
-#         label = "b=$b", batchsize = b, hidden_layers = bs_hidden)
-#     (; batchsize = b, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-#        cpu_time = r.cpu.time, gpu_time = r.gpu.time)
-# end
-
-# fig_bs = Figure(size = (720, 480))
-# ax_bs  = Axis(fig_bs[1, 1];
-#     xlabel = "Batch size (hidden = $bs_hidden)",
-#     ylabel = "GPU speedup (CPU / GPU)",
-#     title  = "GPU speedup vs batch size",
-#     xscale = log2,
-#     xticks = (batchsizes, string.(batchsizes)),
-# )
-
-# bs         = [r.batchsize     for r in bs_results]
-# bs_elapsed = [r.elapsed_ratio for r in bs_results]
-# bs_pure    = [r.pure_ratio    for r in bs_results]
-
-# hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
-# scatterlines!(ax_bs, bs, bs_pure;    color = :steelblue, linestyle = :solid, marker = :circle,  label = "pure")
-# scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash,  marker = :xcross, label = "elapsed")
-
-# axislegend(ax_bs; position = :lt)
-# fig_bs
\ No newline at end of file
+widths = [16, 64, 256, 512, 1024]
+depths = 2:5
+
+results = [
+    begin
+            r = bench_cpu_vs_gpu(
+                small_nn_hybrid_model, df, cfg;
+                label = "d=$d w=$w", hidden_layers = fill(w, d)
+            )
+            (;
+                depth = d, width = w,
+                elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+                cpu_time = r.cpu.time, gpu_time = r.gpu.time,
+            )
+        end
+        for d in depths for w in widths
+]
+
+using CairoMakie
+
+fig = Figure(size = (820, 520))
+ax = Axis(
+    fig[1, 1];
+    xlabel = "Hidden-layer width",
+    ylabel = "GPU speedup (CPU / GPU)",
+    title = "GPU speedup vs hidden-layer width, per depth",
+    xscale = log2,
+    xticks = (widths, string.(widths)),
+)
+
+hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
+
+# One colour per depth; `pure` is solid, `elapsed` is dashed.
+palette = Makie.wong_colors()
+for (i, d) in enumerate(depths)
+    rs = filter(r -> r.depth == d, results)
+    ws = [r.width for r in rs]
+    el = [r.elapsed_ratio for r in rs]
+    pr = [r.pure_ratio    for r in rs]
+    c = palette[mod1(i, length(palette))]
+    # Only the solid (pure) line gets a label → legend entries are per-depth.
+    scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
+    scatterlines!(ax, ws, el; color = c, linestyle = :dash, marker = :xcross)
+end
+
+# Legend 1: depth → colour (auto-picked up from the labeled solid lines).
+axislegend(ax, "depth"; position = :lt)
+
+# Legend 2: line style → metric (built manually).
+style_elements = [
+    LineElement(color = :black, linestyle = :solid),
+    LineElement(color = :black, linestyle = :dash),
+]
+Legend(fig[1, 2], style_elements, ["pure", "elapsed"], "metric")
+
+fig
+
+
+batchsizes = [64, 128, 256, 512, 1024, 2048, 4096, 8192]
+bs_hidden = [64, 64, 64]
+
+bs_results = map(batchsizes) do b
+    r = bench_cpu_vs_gpu(
+        small_nn_hybrid_model, df, cfg;
+        label = "b=$b", batchsize = b, hidden_layers = bs_hidden
+    )
+    (;
+        batchsize = b, elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
+        cpu_time = r.cpu.time, gpu_time = r.gpu.time,
+    )
+end
+
+fig_bs = Figure(size = (720, 480))
+ax_bs = Axis(
+    fig_bs[1, 1];
+    xlabel = "Batch size (hidden = $bs_hidden)",
+    ylabel = "GPU speedup (CPU / GPU)",
+    title = "GPU speedup vs batch size",
+    xscale = log2,
+    xticks = (batchsizes, string.(batchsizes)),
+)
+
+bs = [r.batchsize     for r in bs_results]
+bs_elapsed = [r.elapsed_ratio for r in bs_results]
+bs_pure = [r.pure_ratio    for r in bs_results]
+
+hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
+scatterlines!(ax_bs, bs, bs_pure; color = :steelblue, linestyle = :solid, marker = :circle, label = "pure")
+scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash, marker = :xcross, label = "elapsed")
+
+axislegend(ax_bs; position = :lt)
+fig_bs
diff --git a/docs/setup_gpu_backend.jl b/docs/setup_gpu_backend.jl
index 679dfcf0..52252e90 100644
--- a/docs/setup_gpu_backend.jl
+++ b/docs/setup_gpu_backend.jl
@@ -20,9 +20,9 @@ using MLDataDevices
 
 const _GPU_TRIGGERS = [
     (:LuxCUDA, :CUDADevice),
-    (:AMDGPU,  :AMDGPUDevice),
-    (:Metal,   :MetalDevice),
-    (:oneAPI,  :oneAPIDevice),
+    (:AMDGPU, :AMDGPUDevice),
+    (:Metal, :MetalDevice),
+    (:oneAPI, :oneAPIDevice),
 ]
 
 function _load_gpu_backend()

From d735f2153654631860f38ec6417a6b308ba35550 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Wed, 29 Apr 2026 17:47:36 +0200
Subject: [PATCH 07/10] literate

---
 .../tutorials/synthetic_respiration_gpu.jl    | 84 ++++++++++---------
 1 file changed, 43 insertions(+), 41 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index 3c23ff1e..ac712f2e 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -6,7 +6,7 @@
 # on synthetic data for respiration modeling with Q10 temperature sensitivity.
 #
 
-## for development, use local setup
+# ## For development (local docs run)
 # include("../../setup_local_docsrun.jl")
 
 using EasyHybrid
@@ -19,17 +19,16 @@ include("../../setup_gpu_backend.jl")
 
 cpu_device() isa CPUDevice
 gpu_device() isa GPU_DEVICE_TYPE
+nothing #hide
 
 # ## Data Loading and Preprocessing
 #
-# Load synthetic dataset from GitHub into DataFrame
+# Load synthetic dataset from GitHub into a `DataFrame`.
 
-df = load_timeseries_netcdf("https://github.com/bask0/q10hybrid/raw/master/data/Synthetic4BookChap.nc");
+df = load_timeseries_netcdf("https://github.com/bask0/q10hybrid/raw/master/data/Synthetic4BookChap.nc")
 
-# Select a subset of data for faster execution
-#df = df[1:20000, :];
-#nothing #hide
-#first(df, 5)
+# Select a subset of data for faster execution (especially when benchmarking).
+df = df[1:20000, :]
 
 # ## Define the Physical Model
 #
@@ -50,10 +49,10 @@ end
 #
 # Parameter specification: (default, lower_bound, upper_bound)
 
+# Parameter name | Default | Lower | Upper
 parameters = (
-    ## Parameter name | Default | Lower | Upper      | Description
-    rb = (3.0f0, 0.0f0, 13.0f0),  # Basal respiration [μmol/m²/s]
-    Q10 = (2.0f0, 1.0f0, 4.0f0),   # Temperature sensitivity factor [-]
+    rb = (3.0f0, 0.0f0, 13.0f0), # Basal respiration [μmol/m²/s]
+    Q10 = (2.0f0, 1.0f0, 4.0f0), # Temperature sensitivity factor [-]
 )
 
 # ## Configure Hybrid Model Components
@@ -72,31 +71,31 @@ neural_param_names = [:rb]         # Neural network predicted parameters
 predictors_single_nn = [:sw_pot, :dsw_pot]   # Predictor variables (solar radiation, and its derivative)
 
 small_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn,              # Input features
-    forcing,                 # Forcing variables
-    target,                  # Target variables
-    RbQ10,                  # Process-based model function
-    parameters,              # Parameter definitions
-    neural_param_names,      # NN-predicted parameters
-    global_param_names,      # Global parameters
+    predictors_single_nn, # Input features
+    forcing, # Forcing variables
+    target, # Target variables
+    RbQ10, # Process-based model function
+    parameters, # Parameter definitions
+    neural_param_names, # NN-predicted parameters
+    global_param_names, # Global parameters
     hidden_layers = [16, 16], # Neural network architecture
-    activation = sigmoid,      # Activation function
+    activation = sigmoid, # Activation function
     scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true   # Apply batch normalization to inputs
+    input_batchnorm = true, # Apply batch normalization to inputs
 )
 
 large_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn,              # Input features
-    forcing,                 # Forcing variables
-    target,                  # Target variables
-    RbQ10,                  # Process-based model function
-    parameters,              # Parameter definitions
-    neural_param_names,      # NN-predicted parameters
-    global_param_names,      # Global parameters
+    predictors_single_nn, # Input features
+    forcing, # Forcing variables
+    target, # Target variables
+    RbQ10, # Process-based model function
+    parameters, # Parameter definitions
+    neural_param_names, # NN-predicted parameters
+    global_param_names, # Global parameters
     hidden_layers = [512, 512, 512], # Neural network architecture
-    activation = sigmoid,      # Activation function
+    activation = sigmoid, # Activation function
     scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true   # Apply batch normalization to inputs
+    input_batchnorm = true, # Apply batch normalization to inputs
 )
 
 # ### train on DataFrame
@@ -190,8 +189,8 @@ function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool =
     return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
 end
 
-bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN"); # on our gpu1-hpc22 0.57x
-bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN"); # on our gpu1-hpc22 2.73x
+bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN") # on our gpu1-hpc22 0.57x
+bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN") # on our gpu1-hpc22 2.73x
 
 # ## Sweep: GPU speedup vs hidden-layer width
 #
@@ -204,17 +203,20 @@ depths = 2:5
 
 results = [
     begin
-            r = bench_cpu_vs_gpu(
-                small_nn_hybrid_model, df, cfg;
-                label = "d=$d w=$w", hidden_layers = fill(w, d)
-            )
-            (;
-                depth = d, width = w,
-                elapsed_ratio = r.elapsed_ratio, pure_ratio = r.pure_ratio,
-                cpu_time = r.cpu.time, gpu_time = r.gpu.time,
-            )
-        end
-        for d in depths for w in widths
+        r = bench_cpu_vs_gpu(
+            small_nn_hybrid_model, df, cfg;
+            label = "d=$d w=$w", hidden_layers = fill(w, d),
+        )
+        (;
+            depth = d,
+            width = w,
+            elapsed_ratio = r.elapsed_ratio,
+            pure_ratio = r.pure_ratio,
+            cpu_time = r.cpu.time,
+            gpu_time = r.gpu.time,
+        )
+    end
+    for d in depths for w in widths
 ]
 
 using CairoMakie

From 11484daf57b3be41a8099fbd2235b1332f762e3e Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Tue, 5 May 2026 13:19:30 +0200
Subject: [PATCH 08/10] comments outside functions for Literate.jl

---
 .../tutorials/synthetic_respiration_gpu.jl    | 40 ++++++++++---------
 docs/setup_gpu_backend.jl                     | 12 +++---
 2 files changed, 27 insertions(+), 25 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index ac712f2e..7e0e814e 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -141,10 +141,14 @@ end
 #
 # `warn_frac` is the fraction of total time spent in GC / (re)compile above
 # which a warning is emitted for that run; defaults to 5%.
+# A short 1-epoch run on each device triggers Julia/Zygote specialization
+# and GPU kernel compilation, so the timed runs below measure steady-state
+# training rather than first-call compile latency.
+#
+# Return a slim summary. We deliberately drop `.value` (the full TrainResults)
+# and `.gcstats` so that REPL/Literate printing of the return value stays compact.
 function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool = true, warn_frac::Real = 0.05, tune_kwargs...)
-    # A short 1-epoch run on each device triggers Julia/Zygote specialization
-    # and GPU kernel compilation, so the timed runs below measure steady-state
-    # training rather than first-call compile latency.
+
     if warmup
         warm_cfg = EasyHybrid.TrainConfig(;
             nepochs = 1,
@@ -183,8 +187,6 @@ function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool =
         label, elapsed_ratio, pure_ratio
     )
 
-    # Return a slim summary. We deliberately drop `.value` (the full TrainResults)
-    # and `.gcstats` so that REPL/Literate printing of the return value stays compact.
     slim(s) = (; s.time, s.bytes, s.gctime, s.compile_time, s.recompile_time)
     return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
 end
@@ -203,20 +205,20 @@ depths = 2:5
 
 results = [
     begin
-        r = bench_cpu_vs_gpu(
-            small_nn_hybrid_model, df, cfg;
-            label = "d=$d w=$w", hidden_layers = fill(w, d),
-        )
-        (;
-            depth = d,
-            width = w,
-            elapsed_ratio = r.elapsed_ratio,
-            pure_ratio = r.pure_ratio,
-            cpu_time = r.cpu.time,
-            gpu_time = r.gpu.time,
-        )
-    end
-    for d in depths for w in widths
+            r = bench_cpu_vs_gpu(
+                small_nn_hybrid_model, df, cfg;
+                label = "d=$d w=$w", hidden_layers = fill(w, d),
+            )
+            (;
+                depth = d,
+                width = w,
+                elapsed_ratio = r.elapsed_ratio,
+                pure_ratio = r.pure_ratio,
+                cpu_time = r.cpu.time,
+                gpu_time = r.gpu.time,
+            )
+        end
+        for d in depths for w in widths
 ]
 
 using CairoMakie
diff --git a/docs/setup_gpu_backend.jl b/docs/setup_gpu_backend.jl
index 52252e90..0e8d09fc 100644
--- a/docs/setup_gpu_backend.jl
+++ b/docs/setup_gpu_backend.jl
@@ -25,6 +25,12 @@ const _GPU_TRIGGERS = [
     (:oneAPI, :oneAPIDevice),
 ]
 
+# Note for Literate.jl:
+# `using $pkg` may load an `MLDataDevices` package extension that adds new
+# methods (e.g. a real `functional(::Type{CUDADevice})`). Those methods live in
+# a newer world age than `_load_gpu_backend()`, so we must dispatch via
+# `invokelatest` or we'll get the stale (pre-extension) method that reports
+# `false`.
 function _load_gpu_backend()
     for (pkg, dev) in _GPU_TRIGGERS
         try
@@ -33,12 +39,6 @@ function _load_gpu_backend()
             @debug "Skipping $pkg (not installed / failed to load)" exception = err
             continue
         end
-
-        # `using $pkg` may load an `MLDataDevices` package extension that
-        # adds new methods (e.g. a real `functional(::Type{CUDADevice})`).
-        # Those methods live in a newer world age than this function, so
-        # we must dispatch via `invokelatest` or we'll get the stale
-        # (pre-extension) method that reports `false`.
         DeviceT = Base.invokelatest(getfield, MLDataDevices, dev)
         is_functional = Base.invokelatest(MLDataDevices.functional, DeviceT)
 

From 47c861edc564d8837955c8cdf21f2181a33d83ac Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Wed, 6 May 2026 10:10:35 +0200
Subject: [PATCH 09/10] more inline loop comments - map problem?

---
 .../tutorials/synthetic_respiration_gpu.jl    | 95 +++++++++++--------
 1 file changed, 55 insertions(+), 40 deletions(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index 7e0e814e..ed47cb2a 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -51,55 +51,65 @@ end
 
 # Parameter name | Default | Lower | Upper
 parameters = (
-    rb = (3.0f0, 0.0f0, 13.0f0), # Basal respiration [μmol/m²/s]
-    Q10 = (2.0f0, 1.0f0, 4.0f0), # Temperature sensitivity factor [-]
+    rb = (3.0f0, 0.0f0, 13.0f0),
+    Q10 = (2.0f0, 1.0f0, 4.0f0),
 )
 
 # ## Configure Hybrid Model Components
 #
-# Define input variables
-forcing = [:ta]                    # Forcing variables (temperature)
+# Define input variables (forcing = temperature).
+forcing = [:ta]
 
-# Target variable
-target = [:reco]                   # Target variable (respiration)
+# Target variable (respiration).
+target = [:reco]
 
-# Parameter classification
-global_param_names = [:Q10]        # Global parameters (same for all samples)
-neural_param_names = [:rb]         # Neural network predicted parameters
+# Parameter classification: global parameters are shared across all samples,
+# neural-network-predicted parameters are produced by the NN per sample.
+global_param_names = [:Q10]
+neural_param_names = [:rb]
 
 # ## Single NN Hybrid Model Training
-predictors_single_nn = [:sw_pot, :dsw_pot]   # Predictor variables (solar radiation, and its derivative)
+#
+# Predictor variables for the NN: solar radiation potential and its derivative.
+predictors_single_nn = [:sw_pot, :dsw_pot]
 
+# `constructHybridModel` arguments: predictors, forcing, targets, the
+# process-based model function, parameter definitions, NN-predicted parameters
+# and global parameters. Keyword arguments configure the neural network
+# architecture, activation, output scaling and input batch normalization.
 small_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn, # Input features
-    forcing, # Forcing variables
-    target, # Target variables
-    RbQ10, # Process-based model function
-    parameters, # Parameter definitions
-    neural_param_names, # NN-predicted parameters
-    global_param_names, # Global parameters
-    hidden_layers = [16, 16], # Neural network architecture
-    activation = sigmoid, # Activation function
-    scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true, # Apply batch normalization to inputs
+    predictors_single_nn,
+    forcing,
+    target,
+    RbQ10,
+    parameters,
+    neural_param_names,
+    global_param_names,
+    hidden_layers = [16, 16],
+    activation = sigmoid,
+    scale_nn_outputs = true,
+    input_batchnorm = true,
 )
 
 large_nn_hybrid_model = constructHybridModel(
-    predictors_single_nn, # Input features
-    forcing, # Forcing variables
-    target, # Target variables
-    RbQ10, # Process-based model function
-    parameters, # Parameter definitions
-    neural_param_names, # NN-predicted parameters
-    global_param_names, # Global parameters
-    hidden_layers = [512, 512, 512], # Neural network architecture
-    activation = sigmoid, # Activation function
-    scale_nn_outputs = true, # Scale neural network outputs
-    input_batchnorm = true, # Apply batch normalization to inputs
+    predictors_single_nn,
+    forcing,
+    target,
+    RbQ10,
+    parameters,
+    neural_param_names,
+    global_param_names,
+    hidden_layers = [512, 512, 512],
+    activation = sigmoid,
+    scale_nn_outputs = true,
+    input_batchnorm = true,
 )
 
-# ### train on DataFrame
-# Train the hybrid model
+# ### Train on DataFrame
+#
+# Configure training. Set `keep_history = true` to keep per-epoch history
+# (losses, predictions, etc.), and `save_training = true` to persist training
+# history and checkpoints to disk.
 
 cfg = EasyHybrid.TrainConfig(
     nepochs = 10,
@@ -107,8 +117,8 @@ cfg = EasyHybrid.TrainConfig(
     opt = Adam(0.001),
     loss_types = [:mse, :nse],
     show_progress = false,
-    keep_history = false, # set to true to keep per-epoch history, losses, predictions, etc.
-    save_training = false, # Set to true to enable saving training history and checkpoints
+    keep_history = false,
+    save_training = false,
     plotting = false,
 )
 
@@ -191,8 +201,11 @@ function bench_cpu_vs_gpu(model, df, cfg; label::AbstractString, warmup::Bool =
     return (; cpu = slim(cpu_stats), gpu = slim(gpu_stats), elapsed_ratio, pure_ratio)
 end
 
-bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN") # on our gpu1-hpc22 0.57x
-bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN") # on our gpu1-hpc22 2.73x
+# On our `gpu1-hpc22` machine these typically come out around 0.57x for the
+# small NN and 2.73x for the large NN.
+
+bench_cpu_vs_gpu(small_nn_hybrid_model, df, cfg; label = "small NN")
+bench_cpu_vs_gpu(large_nn_hybrid_model, df, cfg; label = "large NN")
 
 # ## Sweep: GPU speedup vs hidden-layer width
 #
@@ -233,7 +246,8 @@ ax = Axis(
     xticks = (widths, string.(widths)),
 )
 
-hlines!(ax, [1.0]; color = :gray, linestyle = :dash) # break-even
+# Break-even line (ratio = 1).
+hlines!(ax, [1.0]; color = :gray, linestyle = :dash)
 
 # One colour per depth; `pure` is solid, `elapsed` is dashed.
 palette = Makie.wong_colors()
@@ -289,7 +303,8 @@ bs = [r.batchsize     for r in bs_results]
 bs_elapsed = [r.elapsed_ratio for r in bs_results]
 bs_pure = [r.pure_ratio    for r in bs_results]
 
-hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash) # break-even
+# Break-even line (ratio = 1).
+hlines!(ax_bs, [1.0]; color = :gray, linestyle = :dash)
 scatterlines!(ax_bs, bs, bs_pure; color = :steelblue, linestyle = :solid, marker = :circle, label = "pure")
 scatterlines!(ax_bs, bs, bs_elapsed; color = :steelblue, linestyle = :dash, marker = :xcross, label = "elapsed")
 

From ff460885780a9c971c0c58a9a8923ec829ce1b91 Mon Sep 17 00:00:00 2001
From: Bernhard Ahrens <bahrens@bgc-jena.mpg.de>
Date: Wed, 6 May 2026 11:17:49 +0200
Subject: [PATCH 10/10] leftover

---
 docs/literate/tutorials/synthetic_respiration_gpu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/literate/tutorials/synthetic_respiration_gpu.jl b/docs/literate/tutorials/synthetic_respiration_gpu.jl
index ed47cb2a..18fb7ab0 100644
--- a/docs/literate/tutorials/synthetic_respiration_gpu.jl
+++ b/docs/literate/tutorials/synthetic_respiration_gpu.jl
@@ -250,6 +250,7 @@ ax = Axis(
 hlines!(ax, [1.0]; color = :gray, linestyle = :dash)
 
 # One colour per depth; `pure` is solid, `elapsed` is dashed.
+# Only the solid (pure) line gets a label → legend entries are per-depth.
 palette = Makie.wong_colors()
 for (i, d) in enumerate(depths)
     rs = filter(r -> r.depth == d, results)
@@ -257,7 +258,6 @@ for (i, d) in enumerate(depths)
     el = [r.elapsed_ratio for r in rs]
     pr = [r.pure_ratio    for r in rs]
     c = palette[mod1(i, length(palette))]
-    # Only the solid (pure) line gets a label → legend entries are per-depth.
     scatterlines!(ax, ws, pr; color = c, linestyle = :solid, marker = :circle, label = "depth=$d")
     scatterlines!(ax, ws, el; color = c, linestyle = :dash, marker = :xcross)
 end