Skip to main content

cortiq_engine/
dit.rs

1//! Lumina2 Next-DiT forward (the Lumina-Image 2.0 denoiser):
2//! noisy latent + caption features + timestep → velocity prediction.
3//!
4//! Third increment of the image-generation runtime
5//! (docs/GENERATIVE.ru.md). Standalone f32 forward loaded from a
6//! diffusers `transformer/` directory, mirroring
7//! Lumina2Transformer2DModel exactly: plain-w RMSNorm (not Gemma's
8//! 1+w), per-head qk-norm, 3-axis complex-interleaved RoPE θ=10000
9//! (caption tokens advance axis 0, image tokens sit at axis0=cap_len
10//! with row/col on axes 1/2), AdaLN modulation with tanh gates from
11//! the 1024-d timestep embedding, sandwich RMSNorms, SwiGLU FFN, and
12//! a final LayerNorm(eps 1e-6, no affine) scaled by (1+scale) before
13//! the patch projection. Attention is full/bidirectional.
14//!
15//! Parity: `python/nextdit_ref.py` + `tests/dit_parity.rs` on the
16//! real Lumina transformer weights.
17
18use crate::pool::Pool;
19use crate::qtensor::QTensor;
20use crate::vae::{StTensor, read_safetensors};
21use cortiq_core::CmfModel;
22use std::collections::HashMap;
23use std::path::Path;
24use std::sync::Arc;
25
26/// A projection weight: exact f32 (diffusers load — Accelerate GEMM)
27/// or a CMF-quantized tensor on the engine's batched dot kernels.
28pub(crate) enum Proj {
29    F32 {
30        w: Vec<f32>,
31        rows: usize,
32        cols: usize,
33    },
34    Q(QTensor),
35}
36
37impl Proj {
38    /// f32 weight `[?, cols]`; rows derived from the data length.
39    pub(crate) fn f32(w: Vec<f32>, cols: usize) -> Self {
40        let rows = w.len() / cols;
41        debug_assert_eq!(w.len(), rows * cols);
42        Proj::F32 { w, rows, cols }
43    }
44
45    /// Load from a CMF directory entry (mmap-resident when quantized;
46    /// an F32 entry dequantizes into the exact-GEMM arm).
47    pub(crate) fn from_model(model: &Arc<CmfModel>, name: &str) -> Result<Self, String> {
48        Ok(match QTensor::from_model(model, name)? {
49            QTensor::F32 { data, rows, cols } => Proj::F32 {
50                w: data,
51                rows,
52                cols,
53            },
54            q => Proj::Q(q),
55        })
56    }
57
58    pub(crate) fn rows(&self) -> usize {
59        match self {
60            Proj::F32 { rows, .. } => *rows,
61            Proj::Q(q) => q.rows(),
62        }
63    }
64
65    /// y[b, rows] = x[b, cols] · Wᵀ.
66    pub(crate) fn matmat(&self, xs: &[f32], b: usize, out: &mut [f32], pool: Option<&Pool>) {
67        match self {
68            Proj::F32 { w, rows, cols } => {
69                crate::fcd_ops::gemm_nt(xs, w, out, b, *cols, *rows, pool)
70            }
71            Proj::Q(q) => q.matmat(xs, b, out, pool),
72        }
73    }
74}
75
76struct Block {
77    /// AdaLN: (linear [4·hidden, 1024], bias [4·hidden]); None = plain norm1.
78    modulation: Option<(Proj, Vec<f32>)>,
79    norm1: Vec<f32>,
80    q: Proj, // [nh·hd, hidden]
81    k: Proj, // [nkv·hd, hidden]
82    v: Proj,
83    o: Proj,          // [hidden, nh·hd]
84    norm_q: Vec<f32>, // [hd]
85    norm_k: Vec<f32>,
86    norm2: Vec<f32>,
87    ffn_norm1: Vec<f32>,
88    w1: Proj, // gate [inter, hidden]
89    w3: Proj, // up
90    w2: Proj, // down [hidden, inter]
91    ffn_norm2: Vec<f32>,
92}
93
94/// Next-DiT: exact f32 from a diffusers directory, or CMF-quantized
95/// (mmap-resident) from a packaged file.
96pub struct NextDit {
97    x_emb: Proj, // [hidden, p·p·c]
98    x_emb_b: Vec<f32>,
99    t_lin1_w: Vec<f32>, // [temb, 256]
100    t_lin1_b: Vec<f32>,
101    t_lin2_w: Vec<f32>, // [temb, temb]
102    t_lin2_b: Vec<f32>,
103    cap_norm: Vec<f32>, // [cap_feat]
104    cap_w: Proj,        // [hidden, cap_feat]
105    cap_b: Vec<f32>,
106    context_refiner: Vec<Block>,
107    noise_refiner: Vec<Block>,
108    layers: Vec<Block>,
109    out_lin1_w: Vec<f32>, // [hidden, temb]
110    out_lin1_b: Vec<f32>,
111    out_lin2: Proj, // [p·p·c, hidden]
112    out_lin2_b: Vec<f32>,
113    pool: Option<Arc<Pool>>,
114    pub hidden: usize,
115    pub in_channels: usize,
116    pub patch: usize,
117    nh: usize,
118    nkv: usize,
119    hd: usize,
120    axes_dim: Vec<usize>,
121    eps: f64,
122}
123
124/// `CMF_DIT_PROF=1`: wall-time totals per forward stage, accumulated
125/// across every block of every step and dumped when the model drops.
126/// Same spirit as `CMF_VAE_PROF` — the knife for "where do the
127/// seconds go" before touching any kernel.
128mod prof {
129    use std::sync::OnceLock;
130    use std::sync::atomic::{AtomicU64, Ordering};
131
132    pub const MODNORM: usize = 0;
133    pub const QKV: usize = 1;
134    pub const ROPE: usize = 2;
135    pub const APACK: usize = 3;
136    pub const AQK: usize = 4;
137    pub const SOFTMAX: usize = 5;
138    pub const APV: usize = 6;
139    pub const OPROJ: usize = 7;
140    pub const FFN: usize = 8;
141    pub const FFNEL: usize = 9;
142    pub const HEADTAIL: usize = 10;
143    pub const GPUBLK: usize = 11;
144    const NAMES: [&str; 12] = [
145        "mod+norms",
146        "qkv-proj",
147        "qknorm+rope",
148        "attn-pack",
149        "attn-qk",
150        "softmax",
151        "attn-pv",
152        "o-proj",
153        "ffn-mm",
154        "ffn-silu",
155        "head+tail",
156        "gpu-block",
157    ];
158    static NS: [AtomicU64; 12] = [const { AtomicU64::new(0) }; 12];
159
160    pub fn on() -> bool {
161        static ON: OnceLock<bool> = OnceLock::new();
162        *ON.get_or_init(|| std::env::var("CMF_DIT_PROF").is_ok_and(|v| v != "0"))
163    }
164
165    /// RAII span: charges its category on drop. Free when prof is off.
166    pub struct Span(Option<(std::time::Instant, usize)>);
167    pub fn span(cat: usize) -> Span {
168        Span(on().then(|| (std::time::Instant::now(), cat)))
169    }
170    impl Drop for Span {
171        fn drop(&mut self) {
172            if let Some((t0, c)) = self.0 {
173                NS[c].fetch_add(t0.elapsed().as_nanos() as u64, Ordering::Relaxed);
174            }
175        }
176    }
177
178    pub fn dump() {
179        if !on() {
180            return;
181        }
182        let total: u64 = NS.iter().map(|a| a.load(Ordering::Relaxed)).sum();
183        if total == 0 {
184            return;
185        }
186        eprintln!("dit prof ({:.1} s total in blocks):", total as f64 / 1e9);
187        for (name, a) in NAMES.iter().zip(&NS) {
188            let ns = a.load(Ordering::Relaxed);
189            eprintln!(
190                "  {name:<12} {:>7.2} s  {:>4.1}%",
191                ns as f64 / 1e9,
192                ns as f64 * 100.0 / total as f64
193            );
194        }
195    }
196}
197
198/// diffusers RMSNorm: x·rsqrt(mean x²+eps) · w (plain w).
199fn rms_norm(x: &[f32], w: &[f32], eps: f64) -> Vec<f32> {
200    let ss = x.iter().map(|&v| (v as f64) * (v as f64)).sum::<f64>() / x.len() as f64;
201    let inv = 1.0 / (ss + eps).sqrt();
202    x.iter()
203        .zip(w)
204        .map(|(&v, &g)| (v as f64 * inv) as f32 * g)
205        .collect()
206}
207
208/// `rms_norm` into a caller buffer — same math, no per-row alloc.
209fn rms_norm_into(x: &[f32], w: &[f32], eps: f64, dst: &mut [f32]) {
210    let ss = x.iter().map(|&v| (v as f64) * (v as f64)).sum::<f64>() / x.len() as f64;
211    let inv = 1.0 / (ss + eps).sqrt();
212    for ((d, &v), &g) in dst.iter_mut().zip(x).zip(w) {
213        *d = (v as f64 * inv) as f32 * g;
214    }
215}
216
217/// `rms_norm` in place — same math, no per-row alloc.
218fn rms_norm_inplace(v: &mut [f32], w: &[f32], eps: f64) {
219    let ss = v.iter().map(|&x| (x as f64) * (x as f64)).sum::<f64>() / v.len() as f64;
220    let inv = 1.0 / (ss + eps).sqrt();
221    for (x, &g) in v.iter_mut().zip(w) {
222        *x = (*x as f64 * inv) as f32 * g;
223    }
224}
225
226/// Rows of `n` items split across pool workers (serial without a pool).
227fn pool_rows(pool: Option<&Pool>, n: usize, f: &(dyn Fn(usize, usize) + Sync)) {
228    match pool {
229        Some(p) => p.run_rows(n, f),
230        None => f(0, n),
231    }
232}
233
234fn silu(v: f32) -> f32 {
235    v / (1.0 + (-v).exp())
236}
237
238/// y = x·Wᵀ + b for a single row.
239fn linear(x: &[f32], w: &[f32], b: &[f32]) -> Vec<f32> {
240    let k = x.len();
241    b.iter()
242        .enumerate()
243        .map(|(o, &bias)| {
244            let row = &w[o * k..(o + 1) * k];
245            bias + row.iter().zip(x).map(|(&a, &c)| a * c).sum::<f32>()
246        })
247        .collect()
248}
249
250/// Row handout for pool workers over one flat buffer (Rust-2021
251/// closures capture the raw pointer field, not the wrapper — hence
252/// the accessor method).
253struct SendRows(*mut f32);
254unsafe impl Send for SendRows {}
255unsafe impl Sync for SendRows {}
256impl SendRows {
257    /// SAFETY: caller guarantees disjoint `[off, off+len)` per worker.
258    #[allow(clippy::mut_from_ref)] // the disjoint-rows contract IS the point
259    unsafe fn row(&self, off: usize, len: usize) -> &mut [f32] {
260        unsafe { std::slice::from_raw_parts_mut(self.0.add(off), len) }
261    }
262
263    /// Single scattered element (transposed stores). SAFETY: as `row`.
264    unsafe fn set(&self, off: usize, v: f32) {
265        unsafe { *self.0.add(off) = v }
266    }
267}
268
269/// Numerically stable in-place softmax of one full row (NEON exp on
270/// aarch64, scalar elsewhere).
271fn softmax_inplace(row: &mut [f32]) {
272    #[cfg(target_arch = "aarch64")]
273    {
274        crate::attention::softmax_row(row);
275    }
276    #[cfg(not(target_arch = "aarch64"))]
277    {
278        let mx = row.iter().cloned().fold(f32::MIN, f32::max);
279        let mut den = 0f32;
280        for r in row.iter_mut() {
281            *r = (*r - mx).exp();
282            den += *r;
283        }
284        if den > 0.0 {
285            let inv = 1.0 / den;
286            for r in row.iter_mut() {
287                *r *= inv;
288            }
289        }
290    }
291}
292
293/// Any CMF directory entry → f32 (norms, biases, f16 conv kernels).
294pub(crate) fn cmf_f32(model: &CmfModel, name: &str) -> Result<Vec<f32>, String> {
295    let entry = model
296        .tensor(name)
297        .ok_or_else(|| format!("missing tensor {name}"))?;
298    let bytes = model.entry_bytes(entry);
299    let mut out = vec![0f32; entry.shape.iter().product()];
300    cortiq_core::quant::dequant_tensor(entry, bytes, &mut out)?;
301    Ok(out)
302}
303
304/// Per-token RoPE table: interleaved-pair rotation angles, f64.
305/// `ids` is [n, 3]; each axis contributes dim/2 frequencies.
306fn rope_table(ids: &[[u32; 3]], axes_dim: &[usize]) -> (Vec<f64>, Vec<f64>) {
307    let pairs: usize = axes_dim.iter().sum::<usize>() / 2;
308    let mut cos = Vec::with_capacity(ids.len() * pairs);
309    let mut sin = Vec::with_capacity(ids.len() * pairs);
310    for id in ids {
311        for (a, &d) in axes_dim.iter().enumerate() {
312            for j in 0..d / 2 {
313                let freq = 1.0 / 10000f64.powf(2.0 * j as f64 / d as f64);
314                let ang = id[a] as f64 * freq;
315                cos.push(ang.cos());
316                sin.push(ang.sin());
317            }
318        }
319    }
320    (cos, sin)
321}
322
323impl Drop for NextDit {
324    fn drop(&mut self) {
325        prof::dump();
326    }
327}
328
329impl NextDit {
330    pub fn load_dir(dir: &Path) -> Result<Self, String> {
331        let cfg: serde_json::Value = serde_json::from_slice(
332            &std::fs::read(dir.join("config.json")).map_err(|e| format!("config.json: {e}"))?,
333        )
334        .map_err(|e| format!("config.json: {e}"))?;
335        let idx: serde_json::Value = serde_json::from_slice(
336            &std::fs::read(dir.join("diffusion_pytorch_model.safetensors.index.json"))
337                .map_err(|e| format!("index: {e}"))?,
338        )
339        .map_err(|e| format!("index: {e}"))?;
340        let mut shards: Vec<String> = idx["weight_map"]
341            .as_object()
342            .ok_or("weight_map")?
343            .values()
344            .filter_map(|v| v.as_str().map(String::from))
345            .collect();
346        shards.sort();
347        shards.dedup();
348        let mut t: HashMap<String, StTensor> = HashMap::new();
349        for sh in &shards {
350            t.extend(read_safetensors(&dir.join(sh))?);
351        }
352        let mut take = |n: String| -> Result<Vec<f32>, String> {
353            t.remove(&n)
354                .map(|v| v.data)
355                .ok_or_else(|| format!("missing tensor {n}"))
356        };
357        let hidden = cfg["hidden_size"].as_u64().ok_or("hidden")? as usize;
358        let mut blocks = |pfx: &str, count: usize, modulated: bool| -> Result<Vec<Block>, String> {
359            (0..count)
360                .map(|l| {
361                    let p = format!("{pfx}.{l}");
362                    let w1 = take(format!("{p}.feed_forward.linear_1.weight"))?;
363                    let inter = w1.len() / hidden;
364                    Ok(Block {
365                        modulation: if modulated {
366                            let mw = take(format!("{p}.norm1.linear.weight"))?;
367                            let cols = mw.len() / (4 * hidden);
368                            Some((Proj::f32(mw, cols), take(format!("{p}.norm1.linear.bias"))?))
369                        } else {
370                            None
371                        },
372                        norm1: if modulated {
373                            take(format!("{p}.norm1.norm.weight"))?
374                        } else {
375                            take(format!("{p}.norm1.weight"))?
376                        },
377                        q: Proj::f32(take(format!("{p}.attn.to_q.weight"))?, hidden),
378                        k: Proj::f32(take(format!("{p}.attn.to_k.weight"))?, hidden),
379                        v: Proj::f32(take(format!("{p}.attn.to_v.weight"))?, hidden),
380                        o: {
381                            let o = take(format!("{p}.attn.to_out.0.weight"))?;
382                            let cols = o.len() / hidden;
383                            Proj::f32(o, cols)
384                        },
385                        norm_q: take(format!("{p}.attn.norm_q.weight"))?,
386                        norm_k: take(format!("{p}.attn.norm_k.weight"))?,
387                        norm2: take(format!("{p}.norm2.weight"))?,
388                        ffn_norm1: take(format!("{p}.ffn_norm1.weight"))?,
389                        w1: Proj::f32(w1, hidden),
390                        w3: Proj::f32(take(format!("{p}.feed_forward.linear_3.weight"))?, hidden),
391                        w2: Proj::f32(take(format!("{p}.feed_forward.linear_2.weight"))?, inter),
392                        ffn_norm2: take(format!("{p}.ffn_norm2.weight"))?,
393                    })
394                })
395                .collect()
396        };
397        let nl = cfg["num_layers"].as_u64().ok_or("num_layers")? as usize;
398        let nr = cfg["num_refiner_layers"].as_u64().unwrap_or(2) as usize;
399        let context_refiner = blocks("context_refiner", nr, false)?;
400        let noise_refiner = blocks("noise_refiner", nr, true)?;
401        let layers = blocks("layers", nl, true)?;
402        let nh = cfg["num_attention_heads"].as_u64().ok_or("nh")? as usize;
403        let in_channels = cfg["in_channels"].as_u64().ok_or("in_channels")? as usize;
404        let patch = cfg["patch_size"].as_u64().unwrap_or(2) as usize;
405        let axes_dim: Vec<usize> = cfg["axes_dim_rope"]
406            .as_array()
407            .ok_or("axes_dim_rope")?
408            .iter()
409            .map(|v| v.as_u64().unwrap_or(0) as usize)
410            .collect();
411        let cap_norm = take("time_caption_embed.caption_embedder.0.weight".into())?;
412        let cap_feat = cap_norm.len();
413        Ok(Self {
414            x_emb: Proj::f32(
415                take("x_embedder.weight".into())?,
416                patch * patch * in_channels,
417            ),
418            x_emb_b: take("x_embedder.bias".into())?,
419            t_lin1_w: take("time_caption_embed.timestep_embedder.linear_1.weight".into())?,
420            t_lin1_b: take("time_caption_embed.timestep_embedder.linear_1.bias".into())?,
421            t_lin2_w: take("time_caption_embed.timestep_embedder.linear_2.weight".into())?,
422            t_lin2_b: take("time_caption_embed.timestep_embedder.linear_2.bias".into())?,
423            cap_norm,
424            cap_w: Proj::f32(
425                take("time_caption_embed.caption_embedder.1.weight".into())?,
426                cap_feat,
427            ),
428            cap_b: take("time_caption_embed.caption_embedder.1.bias".into())?,
429            context_refiner,
430            noise_refiner,
431            layers,
432            out_lin1_w: take("norm_out.linear_1.weight".into())?,
433            out_lin1_b: take("norm_out.linear_1.bias".into())?,
434            out_lin2: Proj::f32(take("norm_out.linear_2.weight".into())?, hidden),
435            out_lin2_b: take("norm_out.linear_2.bias".into())?,
436            pool: Pool::from_env(),
437            hidden,
438            in_channels,
439            patch,
440            nh,
441            nkv: cfg["num_kv_heads"].as_u64().unwrap_or(nh as u64) as usize,
442            hd: hidden / nh,
443            axes_dim,
444            eps: cfg["norm_eps"].as_f64().unwrap_or(1e-5),
445        })
446    }
447
448    /// Load from a packaged imagegen .cmf (`dit.*` tensors +
449    /// `dit.config_json`). Quantized projections stay mmap-resident.
450    pub fn from_cmf(model: &Arc<CmfModel>) -> Result<Self, String> {
451        let cfg: serde_json::Value = serde_json::from_slice(
452            model
453                .tensor_bytes("dit.config_json")
454                .map_err(|e| e.to_string())?,
455        )
456        .map_err(|e| format!("dit.config_json: {e}"))?;
457        let f32v = |n: &str| -> Result<Vec<f32>, String> { cmf_f32(model, n) };
458        let hidden = cfg["hidden_size"].as_u64().ok_or("hidden")? as usize;
459        let blocks = |pfx: &str, count: usize, modulated: bool| -> Result<Vec<Block>, String> {
460            (0..count)
461                .map(|l| {
462                    let p = format!("dit.{pfx}.{l}");
463                    Ok(Block {
464                        modulation: if modulated {
465                            Some((
466                                Proj::from_model(model, &format!("{p}.norm1.linear.weight"))?,
467                                f32v(&format!("{p}.norm1.linear.bias"))?,
468                            ))
469                        } else {
470                            None
471                        },
472                        norm1: if modulated {
473                            f32v(&format!("{p}.norm1.norm.weight"))?
474                        } else {
475                            f32v(&format!("{p}.norm1.weight"))?
476                        },
477                        q: Proj::from_model(model, &format!("{p}.attn.to_q.weight"))?,
478                        k: Proj::from_model(model, &format!("{p}.attn.to_k.weight"))?,
479                        v: Proj::from_model(model, &format!("{p}.attn.to_v.weight"))?,
480                        o: Proj::from_model(model, &format!("{p}.attn.to_out.0.weight"))?,
481                        norm_q: f32v(&format!("{p}.attn.norm_q.weight"))?,
482                        norm_k: f32v(&format!("{p}.attn.norm_k.weight"))?,
483                        norm2: f32v(&format!("{p}.norm2.weight"))?,
484                        ffn_norm1: f32v(&format!("{p}.ffn_norm1.weight"))?,
485                        w1: Proj::from_model(model, &format!("{p}.feed_forward.linear_1.weight"))?,
486                        w3: Proj::from_model(model, &format!("{p}.feed_forward.linear_3.weight"))?,
487                        w2: Proj::from_model(model, &format!("{p}.feed_forward.linear_2.weight"))?,
488                        ffn_norm2: f32v(&format!("{p}.ffn_norm2.weight"))?,
489                    })
490                })
491                .collect()
492        };
493        let nl = cfg["num_layers"].as_u64().ok_or("num_layers")? as usize;
494        let nr = cfg["num_refiner_layers"].as_u64().unwrap_or(2) as usize;
495        let nh = cfg["num_attention_heads"].as_u64().ok_or("nh")? as usize;
496        let axes_dim: Vec<usize> = cfg["axes_dim_rope"]
497            .as_array()
498            .ok_or("axes_dim_rope")?
499            .iter()
500            .map(|v| v.as_u64().unwrap_or(0) as usize)
501            .collect();
502        Ok(Self {
503            x_emb: Proj::from_model(model, "dit.x_embedder.weight")?,
504            x_emb_b: f32v("dit.x_embedder.bias")?,
505            t_lin1_w: f32v("dit.time_caption_embed.timestep_embedder.linear_1.weight")?,
506            t_lin1_b: f32v("dit.time_caption_embed.timestep_embedder.linear_1.bias")?,
507            t_lin2_w: f32v("dit.time_caption_embed.timestep_embedder.linear_2.weight")?,
508            t_lin2_b: f32v("dit.time_caption_embed.timestep_embedder.linear_2.bias")?,
509            cap_norm: f32v("dit.time_caption_embed.caption_embedder.0.weight")?,
510            cap_w: Proj::from_model(model, "dit.time_caption_embed.caption_embedder.1.weight")?,
511            cap_b: f32v("dit.time_caption_embed.caption_embedder.1.bias")?,
512            context_refiner: blocks("context_refiner", nr, false)?,
513            noise_refiner: blocks("noise_refiner", nr, true)?,
514            layers: blocks("layers", nl, true)?,
515            out_lin1_w: f32v("dit.norm_out.linear_1.weight")?,
516            out_lin1_b: f32v("dit.norm_out.linear_1.bias")?,
517            out_lin2: Proj::from_model(model, "dit.norm_out.linear_2.weight")?,
518            out_lin2_b: f32v("dit.norm_out.linear_2.bias")?,
519            pool: Pool::from_env(),
520            hidden,
521            in_channels: cfg["in_channels"].as_u64().ok_or("in_channels")? as usize,
522            patch: cfg["patch_size"].as_u64().unwrap_or(2) as usize,
523            nh,
524            nkv: cfg["num_kv_heads"].as_u64().unwrap_or(nh as u64) as usize,
525            hd: hidden / nh,
526            axes_dim,
527            eps: cfg["norm_eps"].as_f64().unwrap_or(1e-5),
528        })
529    }
530
531    /// Sinusoidal(256, cos-first) → 2-layer MLP → temb [1024].
532    fn time_embed(&self, t: f32) -> Vec<f32> {
533        const HALF: usize = 128;
534        let mut freq = [0f32; 2 * HALF];
535        for i in 0..HALF {
536            let ang = t as f64 * (-(10000f64.ln()) * i as f64 / HALF as f64).exp();
537            freq[i] = ang.cos() as f32;
538            freq[HALF + i] = ang.sin() as f32;
539        }
540        let mut h = linear(&freq, &self.t_lin1_w, &self.t_lin1_b);
541        for v in h.iter_mut() {
542            *v = silu(*v);
543        }
544        linear(&h, &self.t_lin2_w, &self.t_lin2_b)
545    }
546
547    /// Fused on-device SwiGLU FFN (Metal). Taken only once the wide
548    /// GEMM probe has settled on the GPU arm: during probing the
549    /// per-op path feeds the samples, after a CPU verdict (or a
550    /// contention kill) the CPU path is the right one anyway. Carries
551    /// the same work-proportional contention tripwire as the per-op
552    /// route (cold ops exempt).
553    fn gpu_ffn(&self, blk: &Block, xn: &[f32], n: usize, out: &mut [f32]) -> bool {
554        use crate::gpu;
555        if n < 128 || !gpu::enabled_here() || gpu::mm_killed() {
556            return false;
557        }
558        // A fused block is not a wide matmat: see `fused_block_trusted`.
559        if !gpu::fused_block_trusted()
560            && (gpu::probe_deciding(gpu::OpClass::MatmatWide)
561                || !matches!(gpu::probe_arm(gpu::OpClass::MatmatWide), gpu::ProbeArm::Gpu))
562        {
563            return false;
564        }
565        let (Proj::Q(q1), Proj::Q(q3), Proj::Q(q2)) = (&blk.w1, &blk.w3, &blk.w2) else {
566            return false;
567        };
568        // q4t or q4tp — the fused chain exists for both, and picking by dtype
569        // here is what keeps a q4tp image model off the unfused path (which
570        // ships the [b, inter] intermediates across the CPU boundary twice per
571        // layer: 28 s against 14 s on Lumina at 256px).
572        let tp = q1.mapped_q4tp().is_some();
573        let (Some((m, i1)), Some((_, i3)), Some((_, i2))) = (if tp {
574            (q1.mapped_q4tp(), q3.mapped_q4tp(), q2.mapped_q4tp())
575        } else {
576            (q1.mapped_q4t(), q3.mapped_q4t(), q2.mapped_q4t())
577        }) else {
578            return false;
579        };
580        let inter = q1.rows();
581        let t0 = std::time::Instant::now();
582        let ok = if tp {
583            gpu::q4tp_ffn(m, i1, i3, i2, xn, n, self.hidden, inter, out)
584        } else {
585            gpu::q4t_ffn(m, i1, i3, i2, xn, n, self.hidden, inter, out)
586        };
587        if !ok {
588            return false;
589        }
590        let flops = 6.0 * n as f64 * self.hidden as f64 * inter as f64;
591        let budget = std::time::Duration::from_secs_f64(flops / 1.5e12 * 8.0 + 0.020);
592        let el = t0.elapsed();
593        if el > budget && !gpu::probe_was_cold() {
594            tracing::warn!(
595                "gpu ffn took {el:?} (budget {budget:?}) — device contended, \
596                 CPU for the rest of the process"
597            );
598            gpu::mm_kill();
599        }
600        true
601    }
602
603    /// All-heads attention on the device (same probe/kill gating as
604    /// the fused FFN). Packs q/k/v head-major (pool-parallel), runs
605    /// scores→softmax→P·V→unstack in one command buffer, and writes
606    /// straight into the [n][nh·hd] attn layout the O-projection
607    /// consumes. Returns false → caller runs the CPU per-head loop.
608    fn gpu_attention(
609        &self,
610        q_all: &[f32],
611        k_all: &[f32],
612        v_all: &[f32],
613        n: usize,
614        scale: f32,
615        attn: &mut [f32],
616    ) -> bool {
617        use crate::gpu;
618        let (nh, nkv, hd) = (self.nh, self.nkv, self.hd);
619        if n < 128 || !gpu::enabled_here() || gpu::mm_killed() {
620            return false;
621        }
622        // A fused block is not a wide matmat: see `fused_block_trusted`.
623        if !gpu::fused_block_trusted()
624            && (gpu::probe_deciding(gpu::OpClass::MatmatWide)
625                || !matches!(gpu::probe_arm(gpu::OpClass::MatmatWide), gpu::ProbeArm::Gpu))
626        {
627            return false;
628        }
629        let pool = self.pool.as_deref();
630        let mut qh = vec![0f32; nh * n * hd];
631        let mut kh = vec![0f32; nkv * n * hd];
632        let mut vh = vec![0f32; nkv * n * hd];
633        {
634            let _s = prof::span(prof::APACK);
635            let (sq, sk, sv) = (
636                SendRows(qh.as_mut_ptr()),
637                SendRows(kh.as_mut_ptr()),
638                SendRows(vh.as_mut_ptr()),
639            );
640            pool_rows(pool, n, &|start, end| {
641                for p in start..end {
642                    for h in 0..nh {
643                        // SAFETY: workers cover disjoint token ranges.
644                        unsafe { sq.row((h * n + p) * hd, hd) }
645                            .copy_from_slice(&q_all[(p * nh + h) * hd..(p * nh + h + 1) * hd]);
646                    }
647                    for h in 0..nkv {
648                        unsafe { sk.row((h * n + p) * hd, hd) }
649                            .copy_from_slice(&k_all[(p * nkv + h) * hd..(p * nkv + h + 1) * hd]);
650                        unsafe { sv.row((h * n + p) * hd, hd) }
651                            .copy_from_slice(&v_all[(p * nkv + h) * hd..(p * nkv + h + 1) * hd]);
652                    }
653                }
654            });
655        }
656        let _s = prof::span(prof::AQK);
657        let t0 = std::time::Instant::now();
658        if !gpu::dit_attention(&qh, &kh, &vh, nh, nkv, n, hd, scale, attn) {
659            return false;
660        }
661        let flops = 4.0 * nh as f64 * (n as f64) * (n as f64) * hd as f64;
662        let budget = std::time::Duration::from_secs_f64(flops / 1.5e12 * 8.0 + 0.020);
663        let el = t0.elapsed();
664        if el > budget && !gpu::probe_was_cold() {
665            tracing::warn!(
666                "gpu attention took {el:?} (budget {budget:?}) — device contended, \
667                 CPU for the rest of the process"
668            );
669            gpu::mm_kill();
670        }
671        true
672    }
673
674    /// One whole block on the device (norms → qkv → RoPE → attention
675    /// → O → residual → FFN → residual, single command buffer). Same
676    /// gating and contention tripwire as the per-stage GPU arms.
677    fn gpu_block(
678        &self,
679        blk: &Block,
680        x: &mut [f32],
681        n: usize,
682        rope32: &(Vec<f32>, Vec<f32>),
683        m: &[f32],
684    ) -> bool {
685        use crate::gpu;
686        let (hs, nh, nkv, hd) = (self.hidden, self.nh, self.nkv, self.hd);
687        if n < 128 || !gpu::enabled_here() || gpu::mm_killed() {
688            return false;
689        }
690        // A fused block is not a wide matmat: see `fused_block_trusted`.
691        if !gpu::fused_block_trusted()
692            && (gpu::probe_deciding(gpu::OpClass::MatmatWide)
693                || !matches!(gpu::probe_arm(gpu::OpClass::MatmatWide), gpu::ProbeArm::Gpu))
694        {
695            return false;
696        }
697        // The pack kernel assumes the rope table covers the full head
698        // dim (axes_dim sums to hd — true for Lumina; bail otherwise).
699        if rope32.0.len() != n * hd / 2 {
700            return false;
701        }
702        fn q(p: &Proj) -> Option<(&Arc<CmfModel>, usize)> {
703            match p {
704                Proj::Q(q) => q.mapped_q4t(),
705                Proj::F32 { .. } => None,
706            }
707        }
708        let (
709            Some((model, wq)),
710            Some((_, wk)),
711            Some((_, wv)),
712            Some((_, wo)),
713            Some((_, w1)),
714            Some((_, w3)),
715            Some((_, w2)),
716        ) = (
717            q(&blk.q),
718            q(&blk.k),
719            q(&blk.v),
720            q(&blk.o),
721            q(&blk.w1),
722            q(&blk.w3),
723            q(&blk.w2),
724        )
725        else {
726            return false;
727        };
728        let inter = blk.w1.rows();
729        let gate_msa: Vec<f32> = m[hs..2 * hs].iter().map(|&v| v.tanh()).collect();
730        let gate_mlp: Vec<f32> = m[3 * hs..].iter().map(|&v| v.tanh()).collect();
731        let args = gpu::DitBlockArgs {
732            n,
733            hidden: hs,
734            inter,
735            nh,
736            nkv,
737            hd,
738            eps: self.eps as f32,
739            rope_cos: &rope32.0,
740            rope_sin: &rope32.1,
741            norm1: &blk.norm1,
742            norm2: &blk.norm2,
743            ffn_norm1: &blk.ffn_norm1,
744            ffn_norm2: &blk.ffn_norm2,
745            norm_q: &blk.norm_q,
746            norm_k: &blk.norm_k,
747            s_msa: &m[..hs],
748            gate_msa: &gate_msa,
749            s_mlp: &m[2 * hs..3 * hs],
750            gate_mlp: &gate_mlp,
751            wq,
752            wk,
753            wv,
754            wo,
755            w1,
756            w3,
757            w2,
758        };
759        let t0 = std::time::Instant::now();
760        if !gpu::dit_block(model, &args, x) {
761            return false;
762        }
763        let flops = 2.0 * n as f64 * hs as f64 * ((nh + 2 * nkv) * hd) as f64
764            + 4.0 * nh as f64 * (n as f64) * (n as f64) * hd as f64
765            + 2.0 * n as f64 * hs as f64 * (nh * hd) as f64
766            + 6.0 * n as f64 * hs as f64 * inter as f64;
767        let budget = std::time::Duration::from_secs_f64(flops / 1.5e12 * 8.0 + 0.030);
768        let el = t0.elapsed();
769        if el > budget && !gpu::probe_was_cold() {
770            tracing::warn!(
771                "gpu dit block took {el:?} (budget {budget:?}) — device contended, \
772                 CPU for the rest of the process"
773            );
774            gpu::mm_kill();
775        }
776        true
777    }
778
779    fn block_forward(
780        &self,
781        blk: &Block,
782        x: &mut [f32],
783        rope: &(Vec<f64>, Vec<f64>),
784        rope32: Option<&(Vec<f32>, Vec<f32>)>,
785        temb: Option<&[f32]>,
786    ) {
787        let (hs, nh, nkv, hd) = (self.hidden, self.nh, self.nkv, self.hd);
788        let pool = self.pool.as_deref();
789        let n = x.len() / hs;
790        let modv = {
791            let _s = prof::span(prof::MODNORM);
792            blk.modulation.as_ref().zip(temb).map(|((w, b), t)| {
793                let s: Vec<f32> = t.iter().map(|&v| silu(v)).collect();
794                let mut m = vec![0f32; w.rows()];
795                w.matmat(&s, 1, &mut m, pool);
796                for (v, &bias) in m.iter_mut().zip(b) {
797                    *v += bias;
798                }
799                m
800            })
801        };
802        if let (Some(m), Some(r32)) = (&modv, rope32) {
803            let _s = prof::span(prof::GPUBLK);
804            if self.gpu_block(blk, x, n, r32, m) {
805                return;
806            }
807        }
808        let modnorm = prof::span(prof::MODNORM);
809        let (s_msa, g_msa, s_mlp, g_mlp) = match &modv {
810            Some(m) => (
811                Some(&m[..hs]),
812                Some(&m[hs..2 * hs]),
813                Some(&m[2 * hs..3 * hs]),
814                Some(&m[3 * hs..]),
815            ),
816            None => (None, None, None, None),
817        };
818        // Gates: tanh once per block — every row shares the same gate
819        // vector, and the naive per-element tanh in the residual loop
820        // was billions of repeated evaluations per render.
821        let gate_msa: Option<Vec<f32>> = g_msa.map(|g| g.iter().map(|&v| v.tanh()).collect());
822        let gate_mlp: Option<Vec<f32>> = g_mlp.map(|g| g.iter().map(|&v| v.tanh()).collect());
823        // Pool-parallel row helpers: dst = rms(src)·w · (1+s)  and
824        // x += gate ⊙ rms(src)·w. Same math and per-row summation
825        // order as the serial loops — rows are independent, so the
826        // parallel split is bit-exact.
827        let norm_scaled = |src: &[f32], w: &[f32], s: Option<&[f32]>, dst: &mut [f32]| {
828            let sr = SendRows(dst.as_mut_ptr());
829            pool_rows(pool, n, &|start, end| {
830                for p in start..end {
831                    // SAFETY: workers cover disjoint row ranges.
832                    let row = unsafe { sr.row(p * hs, hs) };
833                    rms_norm_into(&src[p * hs..(p + 1) * hs], w, self.eps, row);
834                    if let Some(s) = s {
835                        for (r, &sc) in row.iter_mut().zip(s) {
836                            *r *= 1.0 + sc;
837                        }
838                    }
839                }
840            });
841        };
842        let residual = |src: &[f32], w: &[f32], gate: Option<&[f32]>, x: &mut [f32]| {
843            let sr = SendRows(x.as_mut_ptr());
844            pool_rows(pool, n, &|start, end| {
845                let mut tmp = vec![0f32; hs];
846                for p in start..end {
847                    rms_norm_into(&src[p * hs..(p + 1) * hs], w, self.eps, &mut tmp);
848                    // SAFETY: workers cover disjoint row ranges.
849                    let dst = unsafe { sr.row(p * hs, hs) };
850                    match gate {
851                        Some(g) => {
852                            for ((d, &v), &gt) in dst.iter_mut().zip(&tmp).zip(g) {
853                                *d += gt * v;
854                            }
855                        }
856                        None => {
857                            for (d, &v) in dst.iter_mut().zip(&tmp) {
858                                *d += v;
859                            }
860                        }
861                    }
862                }
863            });
864        };
865        // ── attention ──
866        let mut xn = vec![0f32; n * hs];
867        norm_scaled(x, &blk.norm1, s_msa, &mut xn);
868        drop(modnorm);
869        let mut q_all = vec![0f32; n * nh * hd];
870        let mut k_all = vec![0f32; n * nkv * hd];
871        let mut v_all = vec![0f32; n * nkv * hd];
872        {
873            let _s = prof::span(prof::QKV);
874            blk.q.matmat(&xn, n, &mut q_all, pool);
875            blk.k.matmat(&xn, n, &mut k_all, pool);
876            blk.v.matmat(&xn, n, &mut v_all, pool);
877        }
878        let rope_span = prof::span(prof::ROPE);
879        // per-head qk-norm, then interleaved-pair RoPE
880        let (cos, sin) = rope;
881        let pairs = hd / 2;
882        for (all, heads, w) in [
883            (&mut q_all, nh, &blk.norm_q),
884            (&mut k_all, nkv, &blk.norm_k),
885        ] {
886            let sr = SendRows(all.as_mut_ptr());
887            pool_rows(pool, n, &|start, end| {
888                for p in start..end {
889                    for hh in 0..heads {
890                        // SAFETY: workers cover disjoint token ranges.
891                        let v = unsafe { sr.row((p * heads + hh) * hd, hd) };
892                        rms_norm_inplace(v, w, 1e-5);
893                        for j in 0..pairs {
894                            let (c, s) = (cos[p * pairs + j], sin[p * pairs + j]);
895                            let (a, b) = (v[2 * j] as f64, v[2 * j + 1] as f64);
896                            v[2 * j] = (a * c - b * s) as f32;
897                            v[2 * j + 1] = (a * s + b * c) as f32;
898                        }
899                    }
900                }
901            });
902        }
903        drop(rope_span);
904        // full (bidirectional) softmax attention, GQA — per head:
905        // scores = (Q·s)·Kᵀ and P·V as GEMMs (Accelerate/blocked),
906        // pool-parallel row softmax between them. The naive
907        // per-position loop was the depth wall: at 512px (1064
908        // tokens) attention alone cost hundreds of serial GFLOP.
909        let scale = 1.0 / (hd as f32).sqrt();
910        let hpk = nh / nkv;
911        let mut attn = vec![0f32; n * nh * hd];
912        if !self.gpu_attention(&q_all, &k_all, &v_all, n, scale, &mut attn) {
913            let mut qh = vec![0f32; n * hd];
914            let mut kh = vec![0f32; n * hd];
915            let mut vt = vec![0f32; hd * n]; // V transposed: gemm_nt's W layout
916            let mut scores = vec![0f32; n * n];
917            let mut oh = vec![0f32; n * hd];
918            for hh in 0..nh {
919                let kv = hh / hpk;
920                {
921                    let _s = prof::span(prof::APACK);
922                    let (sq, sk, sv) = (
923                        SendRows(qh.as_mut_ptr()),
924                        SendRows(kh.as_mut_ptr()),
925                        SendRows(vt.as_mut_ptr()),
926                    );
927                    pool_rows(pool, n, &|start, end| {
928                        for p in start..end {
929                            let qsrc = &q_all[(p * nh + hh) * hd..(p * nh + hh + 1) * hd];
930                            // SAFETY: workers cover disjoint token ranges
931                            // (`vt` columns are indexed by token too).
932                            let qd = unsafe { sq.row(p * hd, hd) };
933                            for (d, &v) in qsrc.iter().enumerate() {
934                                qd[d] = v * scale;
935                            }
936                            unsafe { sk.row(p * hd, hd) }.copy_from_slice(
937                                &k_all[(p * nkv + kv) * hd..(p * nkv + kv + 1) * hd],
938                            );
939                            let vv = &v_all[(p * nkv + kv) * hd..(p * nkv + kv + 1) * hd];
940                            for (d, &val) in vv.iter().enumerate() {
941                                unsafe { sv.set(d * n + p, val) };
942                            }
943                        }
944                    });
945                }
946                {
947                    let _s = prof::span(prof::AQK);
948                    crate::fcd_ops::gemm_nt(&qh, &kh, &mut scores, n, hd, n, pool);
949                }
950                {
951                    let _s = prof::span(prof::SOFTMAX);
952                    let sp = SendRows(scores.as_mut_ptr());
953                    let soft = |start: usize, end: usize| {
954                        for r in start..end {
955                            // SAFETY: workers cover disjoint row ranges.
956                            softmax_inplace(unsafe { sp.row(r * n, n) });
957                        }
958                    };
959                    match pool {
960                        Some(p) => p.run_rows(n, &soft),
961                        None => soft(0, n),
962                    }
963                }
964                {
965                    let _s = prof::span(prof::APV);
966                    crate::fcd_ops::gemm_nt(&scores, &vt, &mut oh, n, n, hd, pool);
967                }
968                let _s = prof::span(prof::APACK);
969                let sa = SendRows(attn.as_mut_ptr());
970                pool_rows(pool, n, &|start, end| {
971                    for p in start..end {
972                        // SAFETY: workers cover disjoint token ranges.
973                        unsafe { sa.row((p * nh + hh) * hd, hd) }
974                            .copy_from_slice(&oh[p * hd..(p + 1) * hd]);
975                    }
976                });
977            }
978        }
979        let mut proj = vec![0f32; n * hs];
980        {
981            let _s = prof::span(prof::OPROJ);
982            blk.o.matmat(&attn, n, &mut proj, pool);
983        }
984        let modnorm = prof::span(prof::MODNORM);
985        residual(&proj, &blk.norm2, gate_msa.as_deref(), x);
986        // ── SwiGLU FFN ──
987        norm_scaled(x, &blk.ffn_norm1, s_mlp, &mut xn);
988        drop(modnorm);
989        let mut d_all = vec![0f32; n * hs];
990        let fused = {
991            let _s = prof::span(prof::FFN);
992            self.gpu_ffn(blk, &xn, n, &mut d_all)
993        };
994        if !fused {
995            let inter = blk.w1.rows();
996            let mut g_all = vec![0f32; n * inter];
997            let mut u_all = vec![0f32; n * inter];
998            {
999                let _s = prof::span(prof::FFN);
1000                blk.w1.matmat(&xn, n, &mut g_all, pool);
1001                blk.w3.matmat(&xn, n, &mut u_all, pool);
1002            }
1003            {
1004                let _s = prof::span(prof::FFNEL);
1005                let sg = SendRows(g_all.as_mut_ptr());
1006                pool_rows(pool, n, &|start, end| {
1007                    for p in start..end {
1008                        // SAFETY: workers cover disjoint token ranges.
1009                        let g = unsafe { sg.row(p * inter, inter) };
1010                        for (gv, &uv) in g.iter_mut().zip(&u_all[p * inter..(p + 1) * inter]) {
1011                            *gv = silu(*gv) * uv;
1012                        }
1013                    }
1014                });
1015            }
1016            {
1017                let _s = prof::span(prof::FFN);
1018                blk.w2.matmat(&g_all, n, &mut d_all, pool);
1019            }
1020        }
1021        let _modnorm = prof::span(prof::MODNORM);
1022        residual(&d_all, &blk.ffn_norm2, gate_mlp.as_deref(), x);
1023    }
1024
1025    /// One denoising forward: latent `[c, h, w]` (NCHW), caption
1026    /// features `[cap_n, cap_feat]`, timestep `t` ∈ [0,1] (the
1027    /// pipeline's `1 − σ`). Returns the velocity prediction `[c, h, w]`.
1028    pub fn forward(
1029        &self,
1030        latent: &[f32],
1031        h: usize,
1032        w: usize,
1033        cap: &[f32],
1034        cap_n: usize,
1035        t: f32,
1036    ) -> Vec<f32> {
1037        self.forward_with_cap(latent, h, w, &self.refine_caption(cap, cap_n), cap_n, t)
1038    }
1039
1040    /// Caption features → hidden, through the context refiner. Depends on
1041    /// NOTHING that moves during denoising — not the timestep, not the
1042    /// latents — so the whole thing is a constant of the prompt. The
1043    /// denoise loop hoists it out and hands the result to
1044    /// `forward_with_cap`; it used to be recomputed on every model call,
1045    /// which for 30 steps under CFG meant 60 evaluations of a value with
1046    /// two distinct instances.
1047    pub fn refine_caption(&self, cap: &[f32], cap_n: usize) -> Vec<f32> {
1048        let hs = self.hidden;
1049        let cap_feat = self.cap_norm.len();
1050        let mut cap_n_all = vec![0f32; cap_n * cap_feat];
1051        for i in 0..cap_n {
1052            cap_n_all[i * cap_feat..(i + 1) * cap_feat].copy_from_slice(&rms_norm(
1053                &cap[i * cap_feat..(i + 1) * cap_feat],
1054                &self.cap_norm,
1055                self.eps,
1056            ));
1057        }
1058        let mut cap_e = vec![0f32; cap_n * hs];
1059        self.cap_w
1060            .matmat(&cap_n_all, cap_n, &mut cap_e, self.pool.as_deref());
1061        for i in 0..cap_n {
1062            for (v, &b) in cap_e[i * hs..(i + 1) * hs].iter_mut().zip(&self.cap_b) {
1063                *v += b;
1064            }
1065        }
1066        let cap_ids: Vec<[u32; 3]> = (0..cap_n).map(|i| [i as u32, 0, 0]).collect();
1067        let cap_rope = rope_table(&cap_ids, &self.axes_dim);
1068        for blk in &self.context_refiner {
1069            self.block_forward(blk, &mut cap_e, &cap_rope, None, None);
1070        }
1071        cap_e
1072    }
1073
1074    /// The rest of the forward, from an already-refined caption.
1075    pub fn forward_with_cap(
1076        &self,
1077        latent: &[f32],
1078        h: usize,
1079        w: usize,
1080        cap_e_in: &[f32],
1081        cap_n: usize,
1082        t: f32,
1083    ) -> Vec<f32> {
1084        let (c, p, hs) = (self.in_channels, self.patch, self.hidden);
1085        assert_eq!(latent.len(), c * h * w);
1086        let (hp, wp) = (h / p, w / p);
1087        let n_img = hp * wp;
1088        let head = prof::span(prof::HEADTAIL);
1089        let temb = self.time_embed(t);
1090        let mut cap_e = cap_e_in.to_vec();
1091
1092        // patchify (dy, dx, ch inner order) + x_embedder
1093        let pv = p * p * c;
1094        let mut tok = vec![0f32; n_img * pv];
1095        for ph in 0..hp {
1096            for pw in 0..wp {
1097                let dst = &mut tok[(ph * wp + pw) * pv..(ph * wp + pw + 1) * pv];
1098                for dy in 0..p {
1099                    for dx in 0..p {
1100                        for ch in 0..c {
1101                            dst[(dy * p + dx) * c + ch] =
1102                                latent[ch * h * w + (ph * p + dy) * w + pw * p + dx];
1103                        }
1104                    }
1105                }
1106            }
1107        }
1108        let mut img = vec![0f32; n_img * hs];
1109        self.x_emb
1110            .matmat(&tok, n_img, &mut img, self.pool.as_deref());
1111        for i in 0..n_img {
1112            for (v, &b) in img[i * hs..(i + 1) * hs].iter_mut().zip(&self.x_emb_b) {
1113                *v += b;
1114            }
1115        }
1116
1117        // 3-axis position ids: caption (i,0,0), image (cap_n, row, col)
1118        let cap_ids: Vec<[u32; 3]> = (0..cap_n).map(|i| [i as u32, 0, 0]).collect();
1119        let img_ids: Vec<[u32; 3]> = (0..n_img)
1120            .map(|i| [cap_n as u32, (i / wp) as u32, (i % wp) as u32])
1121            .collect();
1122        let cap_rope = rope_table(&cap_ids, &self.axes_dim);
1123        let img_rope = rope_table(&img_ids, &self.axes_dim);
1124        // f32 twins for the on-device block (values stay f64-derived).
1125        let to32 = |r: &(Vec<f64>, Vec<f64>)| {
1126            (
1127                r.0.iter().map(|&v| v as f32).collect::<Vec<f32>>(),
1128                r.1.iter().map(|&v| v as f32).collect::<Vec<f32>>(),
1129            )
1130        };
1131        let img_rope32 = to32(&img_rope);
1132        drop(head);
1133
1134        // The context refiner already ran in `refine_caption` — it is a
1135        // constant of the prompt, not of this call.
1136        for blk in &self.noise_refiner {
1137            self.block_forward(blk, &mut img, &img_rope, Some(&img_rope32), Some(&temb));
1138        }
1139
1140        // joint sequence: caption first
1141        let n = cap_n + n_img;
1142        let mut x = cap_e;
1143        x.extend_from_slice(&img);
1144        let joint_rope = (
1145            [cap_rope.0, img_rope.0].concat(),
1146            [cap_rope.1, img_rope.1].concat(),
1147        );
1148        let joint_rope32 = to32(&joint_rope);
1149        for blk in &self.layers {
1150            self.block_forward(blk, &mut x, &joint_rope, Some(&joint_rope32), Some(&temb));
1151        }
1152
1153        // norm_out: LayerNorm(eps 1e-6, no affine) · (1+scale), project
1154        let _tail = prof::span(prof::HEADTAIL);
1155        let silu_t: Vec<f32> = temb.iter().map(|&v| silu(v)).collect();
1156        let scale = linear(&silu_t, &self.out_lin1_w, &self.out_lin1_b);
1157        for row in x.chunks_exact_mut(hs) {
1158            let mean = row.iter().map(|&v| v as f64).sum::<f64>() / hs as f64;
1159            let var = row
1160                .iter()
1161                .map(|&v| (v as f64 - mean) * (v as f64 - mean))
1162                .sum::<f64>()
1163                / hs as f64;
1164            let inv = 1.0 / (var + 1e-6).sqrt();
1165            for (v, &s) in row.iter_mut().zip(&scale) {
1166                *v = ((*v as f64 - mean) * inv) as f32 * (1.0 + s);
1167            }
1168        }
1169        let mut out = vec![0f32; n * pv];
1170        self.out_lin2.matmat(&x, n, &mut out, self.pool.as_deref());
1171        for i in 0..n {
1172            for (v, &b) in out[i * pv..(i + 1) * pv].iter_mut().zip(&self.out_lin2_b) {
1173                *v += b;
1174            }
1175        }
1176
1177        // unpatchify image tokens → [c, h, w]
1178        let mut pred = vec![0f32; c * h * w];
1179        for ph in 0..hp {
1180            for pw in 0..wp {
1181                let src = &out[(cap_n + ph * wp + pw) * pv..(cap_n + ph * wp + pw + 1) * pv];
1182                for dy in 0..p {
1183                    for dx in 0..p {
1184                        for ch in 0..c {
1185                            pred[ch * h * w + (ph * p + dy) * w + pw * p + dx] =
1186                                src[(dy * p + dx) * c + ch];
1187                        }
1188                    }
1189                }
1190            }
1191        }
1192        pred
1193    }
1194}