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