diff --git "a/hexstate_quantize.c" "b/hexstate_quantize.c" --- "a/hexstate_quantize.c" +++ "b/hexstate_quantize.c" @@ -253,7 +253,7 @@ static ConfigJson parse_config_json(const char *path) } static void detect_architecture(const STMultiFile *mf, ModelArchitecture *arch, - const char *config_json_path) + const char *config_json_path) { memset(arch, 0, sizeof(*arch)); @@ -275,41 +275,41 @@ static void detect_architecture(const STMultiFile *mf, ModelArchitecture *arch, if (strcmp(cfg.model_type, "llama") == 0 || strcmp(cfg.model_type, "mistral") == 0) { strcpy(arch->architecture, "llama"); - } else if (strcmp(cfg.model_type, "qwen2") == 0) { - strcpy(arch->architecture, "qwen2"); - } else if (strcmp(cfg.model_type, "qwen2_moe") == 0) { - strcpy(arch->architecture, "qwen2moe"); - } else if (strcmp(cfg.model_type, "qwen3_5") == 0 || - strcmp(cfg.model_type, "qwen3_5_text") == 0 || - strcmp(cfg.model_type, "qwen3_5_moe") == 0) { - strcpy(arch->architecture, "qwen2"); /* GGUF arch: qwen2 compat */ - } else if (strcmp(cfg.model_type, "phi3") == 0 || - strcmp(cfg.model_type, "phi") == 0) { - strcpy(arch->architecture, "phi3"); - } else if (strcmp(cfg.model_type, "gemma") == 0 || - strcmp(cfg.model_type, "gemma2") == 0) { - strcpy(arch->architecture, "gemma"); - } else if (strcmp(cfg.model_type, "deepseek_v2") == 0) { - strcpy(arch->architecture, "llama"); - } else if (strcmp(cfg.model_type, "gpt_neox") == 0) { - strcpy(arch->architecture, "gpt_neox"); - } else if (strcmp(cfg.model_type, "falcon") == 0) { - strcpy(arch->architecture, "falcon"); - } else if (cfg.model_type[0]) { - /* Unknown — try llama as fallback */ - strcpy(arch->architecture, "llama"); - } + } else if (strcmp(cfg.model_type, "qwen2") == 0) { + strcpy(arch->architecture, "qwen2"); + } else if (strcmp(cfg.model_type, "qwen2_moe") == 0) { + strcpy(arch->architecture, "qwen2moe"); + } else if (strcmp(cfg.model_type, "qwen3_5") == 0 || + strcmp(cfg.model_type, "qwen3_5_text") == 0 || + strcmp(cfg.model_type, "qwen3_5_moe") == 0) { + strcpy(arch->architecture, "qwen2"); /* GGUF arch: qwen2 compat */ + } else if (strcmp(cfg.model_type, "phi3") == 0 || + strcmp(cfg.model_type, "phi") == 0) { + strcpy(arch->architecture, "phi3"); + } else if (strcmp(cfg.model_type, "gemma") == 0 || + strcmp(cfg.model_type, "gemma2") == 0) { + strcpy(arch->architecture, "gemma"); + } else if (strcmp(cfg.model_type, "deepseek_v2") == 0) { + strcpy(arch->architecture, "llama"); + } else if (strcmp(cfg.model_type, "gpt_neox") == 0) { + strcpy(arch->architecture, "gpt_neox"); + } else if (strcmp(cfg.model_type, "falcon") == 0) { + strcpy(arch->architecture, "falcon"); + } else if (cfg.model_type[0]) { + /* Unknown — try llama as fallback */ + strcpy(arch->architecture, "llama"); + } - if (cfg.hidden_size) arch->embedding_length = cfg.hidden_size; - if (cfg.intermediate_size) arch->feed_forward_length = cfg.intermediate_size; - if (cfg.num_attention_heads) arch->head_count = cfg.num_attention_heads; - if (cfg.num_key_value_heads) arch->head_count_kv = cfg.num_key_value_heads; - if (cfg.num_hidden_layers) arch->block_count = cfg.num_hidden_layers; - if (cfg.vocab_size) arch->vocab_size = cfg.vocab_size; - if (cfg.max_position_embeddings) arch->context_length = cfg.max_position_embeddings; - if (cfg.rope_theta > 0) arch->rope_freq_base = cfg.rope_theta; - if (cfg.rms_norm_eps > 0) arch->rms_norm_eps = cfg.rms_norm_eps; - arch->tie_word_embeddings = cfg.tie_word_embeddings; + if (cfg.hidden_size) arch->embedding_length = cfg.hidden_size; + if (cfg.intermediate_size) arch->feed_forward_length = cfg.intermediate_size; + if (cfg.num_attention_heads) arch->head_count = cfg.num_attention_heads; + if (cfg.num_key_value_heads) arch->head_count_kv = cfg.num_key_value_heads; + if (cfg.num_hidden_layers) arch->block_count = cfg.num_hidden_layers; + if (cfg.vocab_size) arch->vocab_size = cfg.vocab_size; + if (cfg.max_position_embeddings) arch->context_length = cfg.max_position_embeddings; + if (cfg.rope_theta > 0) arch->rope_freq_base = cfg.rope_theta; + if (cfg.rms_norm_eps > 0) arch->rms_norm_eps = cfg.rms_norm_eps; + arch->tie_word_embeddings = cfg.tie_word_embeddings; printf(" Architecture determined from config.json: %s\n", cfg.model_type); } @@ -586,11 +586,11 @@ static int should_quantize(const STTensorInfo *ti, const char *gguf_name) /* Never quantize DeltaNet state-space parameters (1D or small) */ if (strstr(gguf_name, "ssm_a") != NULL) return 0; /* A_log */ - if (strstr(gguf_name, "ssm_dt") != NULL) return 0; /* dt_bias */ - if (strstr(gguf_name, "ssm_conv1d") != NULL) return 0; /* conv kernel */ + if (strstr(gguf_name, "ssm_dt") != NULL) return 0; /* dt_bias */ + if (strstr(gguf_name, "ssm_conv1d") != NULL) return 0; /* conv kernel */ - /* Quantize everything else (attention projections, FFN weights, SSM projections) */ - return 1; + /* Quantize everything else (attention projections, FFN weights, SSM projections) */ + return 1; } /* Detect attention Q/K/V/O projection tensors. @@ -727,8 +727,8 @@ static inline float fast_pow_2_4(float x) /* Compute the Q2_K sub-block reconstruction error for a block at a given * scale multiplier, optionally weighted by importance vector */ static float compute_block_error_q2k(const float *weights, int block_size, - float scale_mult, - const float *importance, int imp_offset) + float scale_mult, + const float *importance, int imp_offset) { float min_val = weights[0]; float max_val = weights[0]; @@ -767,10 +767,10 @@ static float compute_block_error_q2k(const float *weights, int block_size, * Inter-block: CZ(2i, 2(i+1)) — coarse↔coarse neighbor * CZ(2i+1, 2(i+1)+1) — fine↔fine neighbor */ static HPCGraph *build_sensitivity_graph(const float *weights, - int64_t n_elements, - int block_size, - float temperature, - const float *importance) + int64_t n_elements, + int block_size, + float temperature, + const float *importance) { int64_t n_blocks = n_elements / block_size; if (n_blocks < 2) return NULL; @@ -798,9 +798,9 @@ static HPCGraph *build_sensitivity_graph(const float *weights, float min_err = 1e30f; for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) { errors[c] = compute_block_error_q2k(block_weights, block_size, - SCALE_TABLE[c], - importance, - (int)(block_idx * block_size)); + SCALE_TABLE[c], + importance, + (int)(block_idx * block_size)); if (errors[c] < min_err) min_err = errors[c]; } @@ -814,7 +814,7 @@ static HPCGraph *build_sensitivity_graph(const float *weights, for (int v1 = 0; v1 < 6; v1++) { int idx = v0 * 6 + v1; coarse_re[v0] += exp(-(double)(errors[idx] - min_err) / - (2.0 * (double)temperature)); + (2.0 * (double)temperature)); } coarse_norm += coarse_re[v0] * coarse_re[v0]; } @@ -833,7 +833,7 @@ static HPCGraph *build_sensitivity_graph(const float *weights, for (int v0 = 0; v0 < 6; v0++) { int idx = v0 * 6 + v1; fine_re[v1] += exp(-(double)(errors[idx] - min_err) / - (2.0 * (double)temperature)); + (2.0 * (double)temperature)); } fine_norm += fine_re[v1] * fine_re[v1]; } @@ -915,9 +915,9 @@ static const MSEGridConfig MSE_DEFAULT_CONFIG = { * Returns optimized scale; stores absolute min in *out_min. * importance: per-element weights (can be NULL for uniform). */ static float mse_grid_search_q2k_subblock(const float *x, int n, int nmax, - uint8_t *L, float *out_min, - const float *importance, - const MSEGridConfig *cfg) + uint8_t *L, float *out_min, + const float *importance, + const MSEGridConfig *cfg) { float min_val = x[0], max_val = x[0]; for (int i = 1; i < n; i++) { @@ -1029,13 +1029,30 @@ static float mse_grid_search_q2k_subblock(const float *x, int n, int nmax, * constructive interference of per-sub-block phase coherence. * ═════════════════���═════════════════════════════════════════════════════════ */ +/* ── Seed solver search grid ── + * ggml's make_qkx2_quants scans 16 iscale candidates on nmax + [-0.5, +1.0] + * in steps of 0.1. The seed (scales[], mins[]) it produces is what every + * Phase-2b candidate derives its 4-bit Ls/Lm codes from, so seed quality + * propagates through the whole pipeline. Halving the step and extending the + * lower tail costs ~2.4x in a phase that is not the bottleneck. */ +#define QKX2_RMIN (-0.9f) +#define QKX2_RDELTA (0.05f) +#define QKX2_NSTEP 38 /* inclusive: 39 candidates */ + /* Weighted least-squares quantization for a sub-block (ggml make_qkx2_quants). - * Finds optimal (scale, min) by searching 16 candidate iscale values - * and solving weighted least-squares for each. - * Returns scale; *the_min is set to the negative of the optimal min. */ + * Finds optimal (scale, min) by searching QKX2_NSTEP+1 candidate iscale + * values and solving weighted least-squares for each. + * Returns scale; *the_min is set to the negative of the optimal min. + * + * OBJECTIVE: candidate selection uses the weighted SUM OF SQUARES, not the + * weighted mean absolute deviation ggml uses. ggml picks MAD for robustness + * against outliers, but the metric this quantizer is tuned against is RMSE, + * and selecting a seed on L1 while scoring the result on L2 leaves the seed + * systematically off the squared-error optimum. Scoring on L2 here aligns + * the seed with every downstream acceptance test. */ static float hpc_make_qkx2_quants(int n, int nmax, const float *x, - const float *w, uint8_t *L, - float *the_min, uint8_t *Laux) + const float *w, uint8_t *L, + float *the_min, uint8_t *Laux) { float xmin = x[0], xmax = x[0]; float sum_w = w[0], sum_x = w[0] * x[0]; @@ -1054,19 +1071,20 @@ static float hpc_make_qkx2_quants(int n, int nmax, const float *x, float iscale = (float)nmax / (xmax - xmin); float scale = 1.0f / iscale; - float best_mad = 0; + float best_sse = 0; for (int i = 0; i < n; i++) { int l = gguf_nearest_int(iscale * (x[i] - xmin)); if (l < 0) l = 0; if (l > nmax) l = nmax; L[i] = (uint8_t)l; float diff = scale * (float)l + xmin - x[i]; - best_mad += w[i] * fabsf(diff); + best_sse += w[i] * diff * diff; } - /* 16 candidate iscale values: search [-0.5, -0.5 + 0.1*15] + nmax */ - for (int is = 0; is <= 15; is++) { - float try_iscale = (-0.5f + 0.1f * (float)is + (float)nmax) / (xmax - xmin); + /* QKX2_NSTEP+1 candidate iscale values on nmax + [RMIN, RMIN+RDELTA*NSTEP] */ + for (int is = 0; is <= QKX2_NSTEP; is++) { + float try_iscale = (QKX2_RMIN + QKX2_RDELTA * (float)is + (float)nmax) + / (xmax - xmin); float sl = 0, sl2 = 0, sxl = 0; for (int i = 0; i < n; i++) { int l = gguf_nearest_int(try_iscale * (x[i] - xmin)); @@ -1085,14 +1103,14 @@ static float hpc_make_qkx2_quants(int n, int nmax, const float *x, this_min = 0; this_scale = sxl / sl2; } - float mad = 0; + float sse = 0; for (int i = 0; i < n; i++) { float diff = this_scale * (float)Laux[i] + this_min - x[i]; - mad += w[i] * fabsf(diff); + sse += w[i] * diff * diff; } - if (mad < best_mad) { + if (sse < best_sse) { for (int i = 0; i < n; i++) L[i] = Laux[i]; - best_mad = mad; + best_sse = sse; scale = this_scale; xmin = this_min; } @@ -1105,16 +1123,16 @@ static float hpc_make_qkx2_quants(int n, int nmax, const float *x, /* Quantize the scale/min arrays into 4-bit values: make_qp_quants equivalent. * Returns the optimal d such that scales[j] ≈ d × Ls[j]. */ static float hpc_make_qp_quants(int n, int nmax, const float *x, - uint8_t *L, const float *sw) + uint8_t *L, const float *sw) { float xmax = 0; for (int i = 0; i < n; i++) if (x[i] > xmax) xmax = x[i]; - if (xmax < 1e-15f) { - for (int i = 0; i < n; i++) L[i] = 0; - return 0.0f; - } - float iscale = (float)nmax / xmax; + if (xmax < 1e-15f) { + for (int i = 0; i < n; i++) L[i] = 0; + return 0.0f; + } + float iscale = (float)nmax / xmax; for (int i = 0; i < n; i++) { int l = gguf_nearest_int(iscale * x[i]); if (l < 0) l = 0; @@ -1235,7 +1253,7 @@ static void shor_collapse_site(HPCGraph *graph, int target_site, int outcome) uint64_t eid = adj->edge_ids[ei]; HPCEdge *edge = &graph->edges[eid]; uint64_t partner = (edge->site_a == (uint64_t)target_site) ? - edge->site_b : edge->site_a; + edge->site_b : edge->site_a; TrialityQuhit *pq = &graph->locals[partner]; for (int d = 0; d < 6; d++) { @@ -1268,7 +1286,7 @@ static void shor_collapse_site(HPCGraph *graph, int target_site, int outcome) uint64_t eid = adj->edge_ids[ei]; HPCEdge *edge = &graph->edges[eid]; uint64_t partner = (edge->site_a == (uint64_t)target_site) ? - edge->site_b : edge->site_a; + edge->site_b : edge->site_a; /* Remove this edge from partner's adj list */ HPCAdjList *padj = &graph->adj[partner]; @@ -1288,8 +1306,8 @@ static void shor_collapse_site(HPCGraph *graph, int target_site, int outcome) * Kept for reference only; all call sites now use sieve_measure_graph. * ═══════════════════════════════════════════════════════════════════════════ */ static void shor_measure_graph(HPCGraph *graph, int64_t n_sites, - double (*marg_out)[6], int *measured_out, - int deterministic) + double (*marg_out)[6], int *measured_out, + int deterministic) { /* Measure sites from last to first (MSB→LSB, same as Griffiths-Niu) */ for (int64_t k = n_sites - 1; k >= 0; k--) { @@ -1320,8 +1338,8 @@ static void shor_measure_graph(HPCGraph *graph, int64_t n_sites, uint64_t eid = adj->edge_ids[ei]; const HPCEdge *edge = &graph->edges[eid]; if (edge->fidelity < 0.0) continue; /* Skip dead edges */ - uint64_t partner = (edge->site_a == (uint64_t)site_k) ? - edge->site_b : edge->site_a; + uint64_t partner = (edge->site_a == (uint64_t)site_k) ? + edge->site_b : edge->site_a; const TrialityQuhit *pq = &graph->locals[partner]; for (int d = 0; d < 6; d++) { @@ -1393,7 +1411,7 @@ static void shor_measure_graph(HPCGraph *graph, int64_t n_sites, double total = 0.0; for (int v = 0; v < 6; v++) { probs[v] = graph->locals[site_k].edge_re[v] * graph->locals[site_k].edge_re[v] + - graph->locals[site_k].edge_im[v] * graph->locals[site_k].edge_im[v]; + graph->locals[site_k].edge_im[v] * graph->locals[site_k].edge_im[v]; total += probs[v]; } if (total > 1e-30) { @@ -1485,8 +1503,8 @@ static void sieve_collapse_site(HPCGraph *graph, int target_site, int outcome) uint64_t eid = adj->edge_ids[ei]; HPCEdge *edge = &graph->edges[eid]; if (edge->fidelity < 0.0) continue; /* skip dead edges */ - uint64_t partner = (edge->site_a == (uint64_t)target_site) ? - edge->site_b : edge->site_a; + uint64_t partner = (edge->site_a == (uint64_t)target_site) ? + edge->site_b : edge->site_a; TrialityQuhit *pq = &graph->locals[partner]; double norm2 = 0.0; @@ -1495,7 +1513,7 @@ static void sieve_collapse_site(HPCGraph *graph, int target_site, int outcome) pq->edge_re[d] *= f; pq->edge_im[d] *= f; norm2 += pq->edge_re[d] * pq->edge_re[d] - + pq->edge_im[d] * pq->edge_im[d]; + + pq->edge_im[d] * pq->edge_im[d]; } if (norm2 > 1e-30) { double inv = 1.0 / sqrt(norm2); @@ -1513,7 +1531,7 @@ static void sieve_collapse_site(HPCGraph *graph, int target_site, int outcome) uint64_t eid = adj->edge_ids[ei]; HPCEdge *edge = &graph->edges[eid]; uint64_t partner = (edge->site_a == (uint64_t)target_site) ? - edge->site_b : edge->site_a; + edge->site_b : edge->site_a; HPCAdjList *padj = &graph->adj[partner]; for (uint64_t pi = 0; pi < padj->count; pi++) { @@ -1540,7 +1558,7 @@ static void sieve_measure_graph(HPCGraph *graph, int64_t n_sites, double p0[6], total0 = 0.0; for (int d = 0; d < 6; d++) { p0[d] = graph->locals[site_k].edge_re[d] * graph->locals[site_k].edge_re[d] - + graph->locals[site_k].edge_im[d] * graph->locals[site_k].edge_im[d]; + + graph->locals[site_k].edge_im[d] * graph->locals[site_k].edge_im[d]; total0 += p0[d]; } if (total0 > 1e-30) { @@ -1563,13 +1581,13 @@ static void sieve_measure_graph(HPCGraph *graph, int64_t n_sites, const HPCEdge *edge = &graph->edges[eid]; if (edge->fidelity < 0.0) continue; uint64_t partner = (edge->site_a == (uint64_t)site_k) ? - edge->site_b : edge->site_a; + edge->site_b : edge->site_a; const TrialityQuhit *pq = &graph->locals[partner]; double q[6], qt = 0.0; for (int w = 0; w < 6; w++) { q[w] = pq->edge_re[w] * pq->edge_re[w] - + pq->edge_im[w] * pq->edge_im[w]; + + pq->edge_im[w] * pq->edge_im[w]; qt += q[w]; } if (qt > 1e-30) { @@ -1591,7 +1609,7 @@ static void sieve_measure_graph(HPCGraph *graph, int64_t n_sites, for (int d = 1; d < 6; d++) if (log_score[d] > best) best = log_score[d]; - double probs[6], total = 0.0; + double probs[6], total = 0.0; for (int d = 0; d < 6; d++) { if (best - log_score[d] > SIEVE_LOG_SLACK) probs[d] = 0.0; @@ -1690,7 +1708,37 @@ static const int Q4_CAND_TO_QUHIT[Q4_N_CAND] = { * polish (move 3). Larger radii let coordinate descent escape shallower * local minima at O(radius²) extra cost per polish iteration. */ #ifndef HEX_POLISH_ULP -#define HEX_POLISH_ULP 4 +#define HEX_POLISH_ULP 8 + +/* ── Phase-4 multi-start: how many of the best Phase-2 candidates to run + * the full assembly + polish from, keeping the lowest-error result. + * + * The assembly chain (WLS/sieve → ULP search → shaping → diffusion → + * refit → coordinate-descent polish) is monotone but LOCAL: it descends + * into whatever basin its starting (d, dmin) candidate sits in. Restarting + * from the K best candidates and keeping the winner is the single largest + * remaining RMSE lever, since candidate rank by raw scored error is only a + * weak predictor of which basin polishes best. + * + * Cost is roughly linear in K (see HEX_TOPK_PRUNE below, which cuts most of + * it in practice). K = 1 reproduces the previous single-start behaviour. */ +#ifndef HEX_TOPK_RESTARTS +#define HEX_TOPK_RESTARTS 4 +#endif + +/* Skip a restart whose raw candidate error already exceeds the best + * FINISHED block error found so far by this ratio — such a candidate is + * very unlikely to polish into the winner. 1e30f disables pruning. */ +#ifndef HEX_TOPK_PRUNE +#define HEX_TOPK_PRUNE 1.25f +#endif + +/* Maximum coordinate-descent sweeps in the Phase-4.6 polish. The loop + * breaks early on convergence (typically 2-3 sweeps), so raising this + * only costs time on the blocks that are still improving. */ +#ifndef HEX_POLISH_ITERS +#define HEX_POLISH_ITERS 12 +#endif #endif /* ── DC + vesica/wave extended objective (dot-product error cancellation) ── @@ -1729,19 +1777,29 @@ static const int Q4_CAND_TO_QUHIT[Q4_N_CAND] = { * terms are a proxy for row-level structure (the API sees a flat stream); * the Phase-3.9 rolling-DC pass handles cross-block linkage. */ #ifndef HEX_DC_LAMBDA -#define HEX_DC_LAMBDA 1.0f +#define HEX_DC_LAMBDA 0.0f #endif #ifndef HEX_VW_LAMBDA -#define HEX_VW_LAMBDA 1.0f +#define HEX_VW_LAMBDA 0.0f #endif -/* Default (1, 1): unit-strength spectral prior. Empirically (synthetic +/* DEFAULT CHANGED (1,1) → (0,0): PURE-RMSE MODE. + * + * (1, 1) is the unit-strength spectral prior. Empirically (synthetic * benchmark, identical inputs): lowers dot-product output error ~0.8-1.4% * on both mean-only and fold-correlated activation models for ~+0.05% * weight RMSE. The theoretically optimal λ grows with the deployment * model's activation mean energy and row length (the per-block term * under-counts cross-block row coupling); the synthetic sweep kept - * improving monotonically through λ = 4 at ~+0.1% RMSE. Set both to - * 0.0f to recover the exact pure weighted-SSE / minimum-RMSE pipeline. */ + * improving monotonically through λ = 4 at ~+0.1% RMSE. + * + * λ = 0 on both knobs reduces exactly to the pure weighted-SSE objective, + * which is what reported RMSE measures. Because every selection and + * acceptance stage scores with hex_spectral_penalty(), zeroing these two + * constants makes the ENTIRE pipeline minimise exactly the reported metric + * — no stage gives back RMSE for output-error cancellation any more. + * + * Rebuild with -DHEX_DC_LAMBDA=1.0f -DHEX_VW_LAMBDA=1.0f to restore the + * output-error-optimised behaviour. */ /* Spectral penalty of the extended objective for one block: residuals e[n], * fold at n/2. Negative values are possible (anti-phase credit) — the total @@ -1756,15 +1814,16 @@ static inline float hex_spectral_penalty(const float *e, int n) cross += e[i] * e[i + half]; } return (HEX_DC_LAMBDA / (float)n) * dc * dc - + (HEX_VW_LAMBDA / (float)n) * 4.0f * cross; + + (HEX_VW_LAMBDA / (float)n) * 4.0f * cross; } static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, - BlockQ4_0 *output, float *out_total_error, - const float *imat_importance, int verbose) + BlockQ4_0 *output, float *out_total_error, + const float *imat_importance, int verbose) { int64_t n_blocks = n_elements / QK4_0; - float total_err = 0.0f; + /* double: same float-saturation hazard as the Q2_K accumulator. */ + double total_err = 0.0; (void)verbose; /* kept for API symmetry with the Q2_K path */ /* ── Phase 1: Greedy seed — compute scale per block ── */ float *greedy_d = (float *)calloc(n_blocks, sizeof(float)); @@ -1784,9 +1843,9 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, * First find the true optimal d* via 3-iteration WLS, * then generate candidates centered on d* with tight spacing. */ float (*cand_errors)[Q4_N_CAND] = (float (*)[Q4_N_CAND]) - calloc(n_blocks, sizeof(float[Q4_N_CAND])); + calloc(n_blocks, sizeof(float[Q4_N_CAND])); uint16_t (*cand_d16)[Q4_N_CAND] = (uint16_t (*)[Q4_N_CAND]) - calloc(n_blocks, sizeof(uint16_t[Q4_N_CAND])); + calloc(n_blocks, sizeof(uint16_t[Q4_N_CAND])); #pragma omp parallel for schedule(dynamic, 64) for (int64_t blk = 0; blk < n_blocks; blk++) { @@ -1805,7 +1864,7 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, if (q < 0) q = 0; if (q > 15) q = 15; float qc = (float)q - 8.0f; float w = (imat_importance) ? - imat_importance[blk * QK4_0 + j] : 1.0f; + imat_importance[blk * QK4_0 + j] : 1.0f; num += w * bw[j] * qc; den += w * qc * qc; dcS += bw[j]; @@ -1824,7 +1883,7 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, } uint16_t cur_wls_d16 = gguf_fp32_to_fp16(wls_d); if (cur_wls_d16 == prev_wls_d16) break; /* converged in FP16 */ - prev_wls_d16 = cur_wls_d16; + prev_wls_d16 = cur_wls_d16; } /* ── Step 2b: Generate candidates centered on WLS optimum ── */ @@ -1861,572 +1920,572 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, for (int64_t i = 0; i < n_blocks; i++) best_candidate[i] = 11; /* Q4_NEIGHBOR_MULTS[11] = 1.00 */ - if (n_blocks >= 2) { - float temperature = 0.5f; - int64_t graph_blocks = (n_blocks > 200) ? 200 : n_blocks; - int64_t stride = n_blocks / graph_blocks; - int64_t n_sites = graph_blocks; /* 1 quhit per block */ + if (n_blocks >= 2) { + float temperature = 0.5f; + int64_t graph_blocks = (n_blocks > 200) ? 200 : n_blocks; + int64_t stride = n_blocks / graph_blocks; + int64_t n_sites = graph_blocks; /* 1 quhit per block */ - HPCGraph *graph = hpc_create(n_sites); - if (graph) { - hpc_ran_q4 = 1; - for (int64_t i = 0; i < n_sites; i++) - triality_dft(&graph->locals[i]); + HPCGraph *graph = hpc_create(n_sites); + if (graph) { + hpc_ran_q4 = 1; + for (int64_t i = 0; i < n_sites; i++) + triality_dft(&graph->locals[i]); - /* Adaptive temperature from error landscape */ - { - double err_accum = 0.0; - int err_count = 0; - for (int64_t gi = 0; gi < graph_blocks && gi < 100; gi++) { - int64_t blk = gi * stride; - float max_e = 0.0f; - for (int c = 0; c < Q4_N_CAND; c++) - if (cand_errors[blk][c] > max_e) - max_e = cand_errors[blk][c]; - err_accum += (double)max_e; - err_count++; - } - if (err_count > 0) { - temperature = (float)(err_accum / err_count) * 0.1f; - if (temperature < 1e-10f) temperature = 1e-10f; + /* Adaptive temperature from error landscape */ + { + double err_accum = 0.0; + int err_count = 0; + for (int64_t gi = 0; gi < graph_blocks && gi < 100; gi++) { + int64_t blk = gi * stride; + float max_e = 0.0f; + for (int c = 0; c < Q4_N_CAND; c++) + if (cand_errors[blk][c] > max_e) + max_e = cand_errors[blk][c]; + err_accum += (double)max_e; + err_count++; + } + if (err_count > 0) { + temperature = (float)(err_accum / err_count) * 0.1f; + if (temperature < 1e-10f) temperature = 1e-10f; + } } - } - /* Encode stride-group AGGREGATED candidate errors as Boltzmann amplitudes */ - for (int64_t i = 0; i < graph_blocks; i++) { - /* Aggregate errors across stride group */ - float agg_errors[Q4_N_CAND]; - for (int c = 0; c < Q4_N_CAND; c++) - agg_errors[c] = 0.0f; + /* Encode stride-group AGGREGATED candidate errors as Boltzmann amplitudes */ + for (int64_t i = 0; i < graph_blocks; i++) { + /* Aggregate errors across stride group */ + float agg_errors[Q4_N_CAND]; + for (int c = 0; c < Q4_N_CAND; c++) + agg_errors[c] = 0.0f; - int64_t blk_start = i * stride; - int64_t blk_end = blk_start + stride; - if (blk_end > n_blocks) blk_end = n_blocks; - int64_t group_size = blk_end - blk_start; + int64_t blk_start = i * stride; + int64_t blk_end = blk_start + stride; + if (blk_end > n_blocks) blk_end = n_blocks; + int64_t group_size = blk_end - blk_start; - for (int64_t b = blk_start; b < blk_end; b++) { - for (int c = 0; c < Q4_N_CAND; c++) - agg_errors[c] += cand_errors[b][c]; - } - if (group_size > 1) { - float inv_gs = 1.0f / (float)group_size; - for (int c = 0; c < Q4_N_CAND; c++) - agg_errors[c] *= inv_gs; - } + for (int64_t b = blk_start; b < blk_end; b++) { + for (int c = 0; c < Q4_N_CAND; c++) + agg_errors[c] += cand_errors[b][c]; + } + if (group_size > 1) { + float inv_gs = 1.0f / (float)group_size; + for (int c = 0; c < Q4_N_CAND; c++) + agg_errors[c] *= inv_gs; + } - float min_err = 1e30f; - for (int c = 0; c < Q4_N_CAND; c++) - if (agg_errors[c] < min_err) - min_err = agg_errors[c]; + float min_err = 1e30f; + for (int c = 0; c < Q4_N_CAND; c++) + if (agg_errors[c] < min_err) + min_err = agg_errors[c]; - double amp_re[6]; - double amp_norm = 0.0; - for (int qi = 0; qi < 6; qi++) amp_re[qi] = 0.0; - for (int ci = 0; ci < Q4_N_CAND; ci++) { - int qi = Q4_CAND_TO_QUHIT[ci]; - amp_re[qi] += exp(-(double)(agg_errors[ci] - min_err) / - (2.0 * (double)temperature)); - } - for (int qi = 0; qi < 6; qi++) - amp_norm += amp_re[qi] * amp_re[qi]; - if (amp_norm > 1e-30) { - double inv = 1.0 / sqrt(amp_norm); - for (int v = 0; v < 6; v++) amp_re[v] *= inv; - } + double amp_re[6]; + double amp_norm = 0.0; + for (int qi = 0; qi < 6; qi++) amp_re[qi] = 0.0; + for (int ci = 0; ci < Q4_N_CAND; ci++) { + int qi = Q4_CAND_TO_QUHIT[ci]; + amp_re[qi] += exp(-(double)(agg_errors[ci] - min_err) / + (2.0 * (double)temperature)); + } + for (int qi = 0; qi < 6; qi++) + amp_norm += amp_re[qi] * amp_re[qi]; + if (amp_norm > 1e-30) { + double inv = 1.0 / sqrt(amp_norm); + for (int v = 0; v < 6; v++) amp_re[v] *= inv; + } - for (int v = 0; v < 6; v++) { - graph->locals[i].edge_re[v] = amp_re[v]; - graph->locals[i].edge_im[v] = 0.0; + for (int v = 0; v < 6; v++) { + graph->locals[i].edge_re[v] = amp_re[v]; + graph->locals[i].edge_im[v] = 0.0; + } + graph->locals[i].primary = VIEW_EDGE; + graph->locals[i].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; + graph->locals[i].delta_valid = 0; + triality_update_mask(&graph->locals[i]); } - graph->locals[i].primary = VIEW_EDGE; - graph->locals[i].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; - graph->locals[i].delta_valid = 0; - triality_update_mask(&graph->locals[i]); - } - - /* Neighbor edges */ - for (int64_t i = 0; i < graph_blocks - 1; i++) - hpc_cz(graph, i, i + 1); - /* ── Sieve sequential selection ── - * Replaces Shor/BP with real-domain log-sieve + parity - * back-action (see sieve_measure_graph). - * Single pass, no iteration, no message damping. */ - double (*marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); - int *shor_measured = (int *)calloc(graph_blocks, sizeof(int)); + /* Neighbor edges */ + for (int64_t i = 0; i < graph_blocks - 1; i++) + hpc_cz(graph, i, i + 1); - sieve_measure_graph(graph, graph_blocks, marg, shor_measured, 1); + /* ── Sieve sequential selection ── + * Replaces Shor/BP with real-domain log-sieve + parity + * back-action (see sieve_measure_graph). + * Single pass, no iteration, no message damping. */ + double (*marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); + int *shor_measured = (int *)calloc(graph_blocks, sizeof(int)); - free(shor_measured); + sieve_measure_graph(graph, graph_blocks, marg, shor_measured, 1); - /* Beam search over candidates */ - typedef struct { double acc_error; int history_idx; } Q4Beam; - typedef struct { int cand_idx; int parent_idx; } Q4BeamHistory; + free(shor_measured); - Q4Beam beams[Q4_N_BEAMS]; - int active_beams = 1; - Q4BeamHistory *history = (Q4BeamHistory *)malloc(n_blocks * Q4_N_BEAMS * sizeof(Q4BeamHistory)); + /* Beam search over candidates */ + typedef struct { double acc_error; int history_idx; } Q4Beam; + typedef struct { int cand_idx; int parent_idx; } Q4BeamHistory; - for (int b = 0; b < Q4_N_BEAMS; b++) { - beams[b].acc_error = 0.0; - beams[b].history_idx = -1; - } + Q4Beam beams[Q4_N_BEAMS]; + int active_beams = 1; + Q4BeamHistory *history = (Q4BeamHistory *)malloc(n_blocks * Q4_N_BEAMS * sizeof(Q4BeamHistory)); - for (int64_t i = 0; i < graph_blocks; i++) { - double m_total = 0.0; - for (int v = 0; v < 6; v++) m_total += marg[i][v]; - - double cand_score[Q4_N_CAND]; - int64_t blk = i * stride; - /* Count candidates per quhit bin for normalization */ - int q4_bin_count[6] = {0}; - for (int ci = 0; ci < Q4_N_CAND; ci++) - q4_bin_count[Q4_CAND_TO_QUHIT[ci]]++; - /* Per-block error normalization: divide by block mean error - * so small-weight blocks don't dominate beam selection */ - float blk_mean_err = 0.0f; - for (int ci = 0; ci < Q4_N_CAND; ci++) - blk_mean_err += cand_errors[blk][ci]; - blk_mean_err /= (float)Q4_N_CAND; - if (blk_mean_err < 1e-30f) blk_mean_err = 1e-30f; - for (int ci = 0; ci < Q4_N_CAND; ci++) { - int qi = Q4_CAND_TO_QUHIT[ci]; - double p = (m_total > 1e-30) ? marg[i][qi] / m_total : 1.0/6.0; - p /= (double)q4_bin_count[qi]; /* normalize by bin occupancy */ - cand_score[ci] = p / (cand_errors[blk][ci] / blk_mean_err + 1e-15); + for (int b = 0; b < Q4_N_BEAMS; b++) { + beams[b].acc_error = 0.0; + beams[b].history_idx = -1; } - typedef struct { double score; int beam_idx; int cand_idx; } Q4Ext; - Q4Ext extensions[Q4_N_BEAMS * Q4_N_CAND]; - int n_ext = 0; - for (int b = 0; b < active_beams; b++) { - for (int c = 0; c < Q4_N_CAND; c++) { - double ext_err = beams[b].acc_error + cand_errors[blk][c]; - extensions[n_ext].score = cand_score[c] / (ext_err + 1e-15); - extensions[n_ext].beam_idx = b; - extensions[n_ext].cand_idx = c; - n_ext++; + for (int64_t i = 0; i < graph_blocks; i++) { + double m_total = 0.0; + for (int v = 0; v < 6; v++) m_total += marg[i][v]; + + double cand_score[Q4_N_CAND]; + int64_t blk = i * stride; + /* Count candidates per quhit bin for normalization */ + int q4_bin_count[6] = {0}; + for (int ci = 0; ci < Q4_N_CAND; ci++) + q4_bin_count[Q4_CAND_TO_QUHIT[ci]]++; + /* Per-block error normalization: divide by block mean error + * so small-weight blocks don't dominate beam selection */ + float blk_mean_err = 0.0f; + for (int ci = 0; ci < Q4_N_CAND; ci++) + blk_mean_err += cand_errors[blk][ci]; + blk_mean_err /= (float)Q4_N_CAND; + if (blk_mean_err < 1e-30f) blk_mean_err = 1e-30f; + for (int ci = 0; ci < Q4_N_CAND; ci++) { + int qi = Q4_CAND_TO_QUHIT[ci]; + double p = (m_total > 1e-30) ? marg[i][qi] / m_total : 1.0/6.0; + p /= (double)q4_bin_count[qi]; /* normalize by bin occupancy */ + cand_score[ci] = p / (cand_errors[blk][ci] / blk_mean_err + 1e-15); } - } - int top_k = (n_ext < Q4_N_BEAMS) ? n_ext : Q4_N_BEAMS; - int top_indices[Q4_N_BEAMS]; - for (int k = 0; k < top_k; k++) { - int best = -1; double best_s = -1e30; - for (int e = 0; e < n_ext; e++) { - if (extensions[e].score > best_s) { - best_s = extensions[e].score; best = e; + typedef struct { double score; int beam_idx; int cand_idx; } Q4Ext; + Q4Ext extensions[Q4_N_BEAMS * Q4_N_CAND]; + int n_ext = 0; + for (int b = 0; b < active_beams; b++) { + for (int c = 0; c < Q4_N_CAND; c++) { + double ext_err = beams[b].acc_error + cand_errors[blk][c]; + extensions[n_ext].score = cand_score[c] / (ext_err + 1e-15); + extensions[n_ext].beam_idx = b; + extensions[n_ext].cand_idx = c; + n_ext++; } } - top_indices[k] = best; - extensions[best].score = -2e30; - } - Q4Beam new_beams[Q4_N_BEAMS]; - for (int k = 0; k < top_k; k++) { - int ei = top_indices[k]; - int sb = extensions[ei].beam_idx; - int cand = extensions[ei].cand_idx; + int top_k = (n_ext < Q4_N_BEAMS) ? n_ext : Q4_N_BEAMS; + int top_indices[Q4_N_BEAMS]; + for (int k = 0; k < top_k; k++) { + int best = -1; double best_s = -1e30; + for (int e = 0; e < n_ext; e++) { + if (extensions[e].score > best_s) { + best_s = extensions[e].score; best = e; + } + } + top_indices[k] = best; + extensions[best].score = -2e30; + } - int hist_idx = i * Q4_N_BEAMS + k; - history[hist_idx].cand_idx = cand; - history[hist_idx].parent_idx = beams[sb].history_idx; + Q4Beam new_beams[Q4_N_BEAMS]; + for (int k = 0; k < top_k; k++) { + int ei = top_indices[k]; + int sb = extensions[ei].beam_idx; + int cand = extensions[ei].cand_idx; - new_beams[k].history_idx = hist_idx; - new_beams[k].acc_error = beams[sb].acc_error + cand_errors[blk][cand]; - } - for (int k = 0; k < top_k; k++) beams[k] = new_beams[k]; - active_beams = top_k; - } + int hist_idx = i * Q4_N_BEAMS + k; + history[hist_idx].cand_idx = cand; + history[hist_idx].parent_idx = beams[sb].history_idx; - int curr_hist = beams[0].history_idx; - for (int64_t i = graph_blocks - 1; i >= 0; i--) { - int group_cidx; - if (curr_hist >= 0) { - group_cidx = history[curr_hist].cand_idx; - curr_hist = history[curr_hist].parent_idx; - } else { - group_cidx = 11; + new_beams[k].history_idx = hist_idx; + new_beams[k].acc_error = beams[sb].acc_error + cand_errors[blk][cand]; + } + for (int k = 0; k < top_k; k++) beams[k] = new_beams[k]; + active_beams = top_k; } - if (stride <= 1) { - best_candidate[i] = group_cidx; - } else { - /* Per-block local optimization within stride group. - * Beam picks the quhit bin; each block picks its best - * candidate in that bin from its own error landscape. */ - int target_bin = Q4_CAND_TO_QUHIT[group_cidx]; - - for (int64_t b = i * stride; b < (i+1) * stride && b < n_blocks; b++) { - float best_err = 1e30f; - int best_c = group_cidx; - for (int c = 0; c < Q4_N_CAND; c++) { - if (Q4_CAND_TO_QUHIT[c] != target_bin) continue; - if (cand_errors[b][c] < best_err) { - best_err = cand_errors[b][c]; - best_c = c; + int curr_hist = beams[0].history_idx; + for (int64_t i = graph_blocks - 1; i >= 0; i--) { + int group_cidx; + if (curr_hist >= 0) { + group_cidx = history[curr_hist].cand_idx; + curr_hist = history[curr_hist].parent_idx; + } else { + group_cidx = 11; + } + + if (stride <= 1) { + best_candidate[i] = group_cidx; + } else { + /* Per-block local optimization within stride group. + * Beam picks the quhit bin; each block picks its best + * candidate in that bin from its own error landscape. */ + int target_bin = Q4_CAND_TO_QUHIT[group_cidx]; + + for (int64_t b = i * stride; b < (i+1) * stride && b < n_blocks; b++) { + float best_err = 1e30f; + int best_c = group_cidx; + for (int c = 0; c < Q4_N_CAND; c++) { + if (Q4_CAND_TO_QUHIT[c] != target_bin) continue; + if (cand_errors[b][c] < best_err) { + best_err = cand_errors[b][c]; + best_c = c; + } } - } - /* Greedy override if global best is >5% better */ - float global_best = 1e30f; - int global_best_c = group_cidx; - for (int c = 0; c < Q4_N_CAND; c++) { - if (cand_errors[b][c] < global_best) { - global_best = cand_errors[b][c]; - global_best_c = c; + /* Greedy override if global best is >5% better */ + float global_best = 1e30f; + int global_best_c = group_cidx; + for (int c = 0; c < Q4_N_CAND; c++) { + if (cand_errors[b][c] < global_best) { + global_best = cand_errors[b][c]; + global_best_c = c; + } } + if (global_best < best_err * HEX_GREEDY_OVERRIDE_RATIO) + best_candidate[b] = global_best_c; + else + best_candidate[b] = best_c; } - if (global_best < best_err * HEX_GREEDY_OVERRIDE_RATIO) - best_candidate[b] = global_best_c; - else - best_candidate[b] = best_c; } } - } - free(history); - - /* ══════════════════════════════════════════════════════════════ - * Phase 3.5: Born-Rule Multi-Shot Scale Refinement - * - * The beam search found the MAP candidate sequence. But the - * triality marginals encode quantum phase-coherent structure - * that a greedy beam can miss. - * ══════════════════════════════════════════════════════════════ */ - { - #define Q4_BORN_SHOTS 128 - - /* Build per-block CDFs from triality marginals */ - unsigned int born_rng = 314159; - - /* Compute tail error once (blocks beyond graph coverage) */ - float tail_err_q4 = 0.0f; - for (int64_t bi = graph_blocks * stride; bi < n_blocks; bi++) - tail_err_q4 += cand_errors[bi][best_candidate[bi]]; - - /* Beam-search baseline over the SAME set of blocks a Born - * shot covers: stride representatives + tail. The previous - * code summed the baseline over ALL blocks (including - * mid-stride blocks the shots never touch), making shot_err - * systematically smaller than the baseline and letting - * strictly worse configurations be adopted whenever - * stride > 1. */ - float beam_total_err = tail_err_q4; - for (int64_t gi = 0; gi < graph_blocks; gi++) { - int64_t rep = gi * stride; - beam_total_err += cand_errors[rep][best_candidate[rep]]; - } - - /* Sparse shot buffer: only track stride-sampled blocks */ - int *shot_sparse_q4 = (int *)malloc(graph_blocks * sizeof(int)); - - for (int shot = 0; shot < Q4_BORN_SHOTS; shot++) { - float shot_err = tail_err_q4; - + free(history); + + /* ══════════════════════════════════════════════════════════════ + * Phase 3.5: Born-Rule Multi-Shot Scale Refinement + * + * The beam search found the MAP candidate sequence. But the + * triality marginals encode quantum phase-coherent structure + * that a greedy beam can miss. + * ══════════════════════════════════════════════════════════════ */ + { + #define Q4_BORN_SHOTS 128 + + /* Build per-block CDFs from triality marginals */ + unsigned int born_rng = 314159; + + /* Compute tail error once (blocks beyond graph coverage) */ + float tail_err_q4 = 0.0f; + for (int64_t bi = graph_blocks * stride; bi < n_blocks; bi++) + tail_err_q4 += cand_errors[bi][best_candidate[bi]]; + + /* Beam-search baseline over the SAME set of blocks a Born + * shot covers: stride representatives + tail. The previous + * code summed the baseline over ALL blocks (including + * mid-stride blocks the shots never touch), making shot_err + * systematically smaller than the baseline and letting + * strictly worse configurations be adopted whenever + * stride > 1. */ + float beam_total_err = tail_err_q4; for (int64_t gi = 0; gi < graph_blocks; gi++) { - /* Normalize marginals to CDF */ - double m_total = 0.0; - for (int v = 0; v < 6; v++) m_total += marg[gi][v]; - - /* Born sample: CDF inversion (same as born_sample) */ - born_rng = born_rng * 1664525u + 1013904223u; - double rnd = (double)(born_rng >> 8) / 16777216.0; - double target = rnd * m_total; - double cum = 0.0; - int sampled_qi = 5; - for (int v = 0; v < 6; v++) { - cum += marg[gi][v]; - if (cum > target) { sampled_qi = v; break; } - } + int64_t rep = gi * stride; + beam_total_err += cand_errors[rep][best_candidate[rep]]; + } - /* Find the best candidate WITHIN this quhit bin */ - int64_t blk = gi * stride; - float best_bin_err = 1e30f; - int best_bin_cand = 11; /* default */ - for (int ci = 0; ci < Q4_N_CAND; ci++) { - if (Q4_CAND_TO_QUHIT[ci] == sampled_qi) { - if (cand_errors[blk][ci] < best_bin_err) { - best_bin_err = cand_errors[blk][ci]; - best_bin_cand = ci; + /* Sparse shot buffer: only track stride-sampled blocks */ + int *shot_sparse_q4 = (int *)malloc(graph_blocks * sizeof(int)); + + for (int shot = 0; shot < Q4_BORN_SHOTS; shot++) { + float shot_err = tail_err_q4; + + for (int64_t gi = 0; gi < graph_blocks; gi++) { + /* Normalize marginals to CDF */ + double m_total = 0.0; + for (int v = 0; v < 6; v++) m_total += marg[gi][v]; + + /* Born sample: CDF inversion (same as born_sample) */ + born_rng = born_rng * 1664525u + 1013904223u; + double rnd = (double)(born_rng >> 8) / 16777216.0; + double target = rnd * m_total; + double cum = 0.0; + int sampled_qi = 5; + for (int v = 0; v < 6; v++) { + cum += marg[gi][v]; + if (cum > target) { sampled_qi = v; break; } + } + + /* Find the best candidate WITHIN this quhit bin */ + int64_t blk = gi * stride; + float best_bin_err = 1e30f; + int best_bin_cand = 11; /* default */ + for (int ci = 0; ci < Q4_N_CAND; ci++) { + if (Q4_CAND_TO_QUHIT[ci] == sampled_qi) { + if (cand_errors[blk][ci] < best_bin_err) { + best_bin_err = cand_errors[blk][ci]; + best_bin_cand = ci; + } } } + + shot_sparse_q4[gi] = best_bin_cand; + shot_err += cand_errors[blk][best_bin_cand]; } - shot_sparse_q4[gi] = best_bin_cand; - shot_err += cand_errors[blk][best_bin_cand]; + /* Metropolis acceptance: adopt if better than current best */ + if (shot_err < beam_total_err) { + for (int64_t gi = 0; gi < graph_blocks; gi++) + best_candidate[gi * stride] = shot_sparse_q4[gi]; + beam_total_err = shot_err; + } } - /* Metropolis acceptance: adopt if better than current best */ - if (shot_err < beam_total_err) { - for (int64_t gi = 0; gi < graph_blocks; gi++) - best_candidate[gi * stride] = shot_sparse_q4[gi]; - beam_total_err = shot_err; - } + free(shot_sparse_q4); } - free(shot_sparse_q4); - } - - /* Born refinement pass: non-stride blocks were set during beam - * traceback and never revisited by Born shots. For each such block - * pick the lowest-error candidate within the same quhit bin that - * the winning Born shot chose for its stride-representative. */ - if (stride > 1) { - for (int64_t b = 0; b < n_blocks; b++) { - if (b % stride == 0) continue; - int64_t rep = (b / stride) * stride; - int target_bin = Q4_CAND_TO_QUHIT[best_candidate[rep]]; - float best_b_err = 1e30f; - int best_b_cand = best_candidate[rep]; - for (int ci = 0; ci < Q4_N_CAND; ci++) { - if (Q4_CAND_TO_QUHIT[ci] != target_bin) continue; - if (cand_errors[b][ci] < best_b_err) { - best_b_err = cand_errors[b][ci]; - best_b_cand = ci; + /* Born refinement pass: non-stride blocks were set during beam + * traceback and never revisited by Born shots. For each such block + * pick the lowest-error candidate within the same quhit bin that + * the winning Born shot chose for its stride-representative. */ + if (stride > 1) { + for (int64_t b = 0; b < n_blocks; b++) { + if (b % stride == 0) continue; + int64_t rep = (b / stride) * stride; + int target_bin = Q4_CAND_TO_QUHIT[best_candidate[rep]]; + float best_b_err = 1e30f; + int best_b_cand = best_candidate[rep]; + for (int ci = 0; ci < Q4_N_CAND; ci++) { + if (Q4_CAND_TO_QUHIT[ci] != target_bin) continue; + if (cand_errors[b][ci] < best_b_err) { + best_b_err = cand_errors[b][ci]; + best_b_cand = ci; + } } + best_candidate[b] = best_b_cand; } - best_candidate[b] = best_b_cand; } - } - free(marg); - hpc_destroy(graph); + free(marg); + hpc_destroy(graph); + } } - } - /* Fallback when the HPC graph never ran (single block, or hpc_create - * failure): pick the per-block argmin over the candidate grid instead - * of silently leaving every block on the neutral ×1.00 candidate. */ - if (!hpc_ran_q4) { - #pragma omp parallel for schedule(static) - for (int64_t blk = 0; blk < n_blocks; blk++) { - float best_e = cand_errors[blk][0]; - int best_c = 0; - for (int c = 1; c < Q4_N_CAND; c++) { - if (cand_errors[blk][c] < best_e) { - best_e = cand_errors[blk][c]; - best_c = c; + /* Fallback when the HPC graph never ran (single block, or hpc_create + * failure): pick the per-block argmin over the candidate grid instead + * of silently leaving every block on the neutral ×1.00 candidate. */ + if (!hpc_ran_q4) { + #pragma omp parallel for schedule(static) + for (int64_t blk = 0; blk < n_blocks; blk++) { + float best_e = cand_errors[blk][0]; + int best_c = 0; + for (int c = 1; c < Q4_N_CAND; c++) { + if (cand_errors[blk][c] < best_e) { + best_e = cand_errors[blk][c]; + best_c = c; + } } + best_candidate[blk] = best_c; } - best_candidate[blk] = best_c; } - } - /* ══════════════════════════════════════════════════════════════════ - * PHASE 4: Assemble blocks via least-squares scale extraction - * ══════════════════════════════════════════════════════════════════ */ - - #pragma omp parallel for schedule(dynamic, 64) reduction(+:total_err) - for (int64_t blk = 0; blk < n_blocks; blk++) { - const float *bw = weights + blk * QK4_0; - int cidx = best_candidate[blk]; + /* ══════════════════════════════════════════════════════════════════ + * PHASE 4: Assemble blocks via least-squares scale extraction + * ══════════════════════════════════════════════════════════════════ */ - /* Start from the grid-selected scale (the "assembled frequency") */ - float d_current = gguf_fp16_to_fp32(cand_d16[blk][cidx]); + #pragma omp parallel for schedule(dynamic, 64) reduction(+:total_err) + for (int64_t blk = 0; blk < n_blocks; blk++) { + const float *bw = weights + blk * QK4_0; + int cidx = best_candidate[blk]; - /* Analog assembly: iterate to full convergence. */ - for (int ls_iter = 0; ls_iter < 5; ls_iter++) { - if (d_current < 1e-15f) break; - float id = 1.0f / d_current; + /* Start from the grid-selected scale (the "assembled frequency") */ + float d_current = gguf_fp16_to_fp32(cand_d16[blk][cidx]); - int qs_tmp[QK4_0]; - for (int j = 0; j < QK4_0; j++) { - int q = (int)(bw[j] * id + 8.5f); - if (q < 0) q = 0; if (q > 15) q = 15; - qs_tmp[j] = q; - } + /* Analog assembly: iterate to full convergence. */ + for (int ls_iter = 0; ls_iter < 5; ls_iter++) { + if (d_current < 1e-15f) break; + float id = 1.0f / d_current; - float num = 0.0f, den = 0.0f; - float dc4S = 0.0f, dc4Q = 0.0f; - for (int j = 0; j < QK4_0; j++) { - float q_centered = (float)qs_tmp[j] - 8.0f; - float w = (imat_importance) ? - imat_importance[blk * QK4_0 + j] : 1.0f; - num += w * bw[j] * q_centered; - den += w * q_centered * q_centered; - dc4S += bw[j]; - dc4Q += q_centered; - } - num += (HEX_DC_LAMBDA / (float)QK4_0) * dc4S * dc4Q; - den += (HEX_DC_LAMBDA / (float)QK4_0) * dc4Q * dc4Q; + int qs_tmp[QK4_0]; + for (int j = 0; j < QK4_0; j++) { + int q = (int)(bw[j] * id + 8.5f); + if (q < 0) q = 0; if (q > 15) q = 15; + qs_tmp[j] = q; + } - if (den > 1e-15f) { - float d_new = num / den; - float d_seed = gguf_fp16_to_fp32(cand_d16[blk][cidx]); - if (fabsf(d_new) < 4.0f * (fabsf(d_seed) + 1e-10f)) { - uint16_t d16 = gguf_fp32_to_fp16(d_new); - d_current = gguf_fp16_to_fp32(d16); + float num = 0.0f, den = 0.0f; + float dc4S = 0.0f, dc4Q = 0.0f; + for (int j = 0; j < QK4_0; j++) { + float q_centered = (float)qs_tmp[j] - 8.0f; + float w = (imat_importance) ? + imat_importance[blk * QK4_0 + j] : 1.0f; + num += w * bw[j] * q_centered; + den += w * q_centered * q_centered; + dc4S += bw[j]; + dc4Q += q_centered; + } + num += (HEX_DC_LAMBDA / (float)QK4_0) * dc4S * dc4Q; + den += (HEX_DC_LAMBDA / (float)QK4_0) * dc4Q * dc4Q; + + if (den > 1e-15f) { + float d_new = num / den; + float d_seed = gguf_fp16_to_fp32(cand_d16[blk][cidx]); + if (fabsf(d_new) < 4.0f * (fabsf(d_seed) + 1e-10f)) { + uint16_t d16 = gguf_fp32_to_fp16(d_new); + d_current = gguf_fp16_to_fp32(d16); + } } } - } - /* ── FP16 ULP neighborhood search + sign-flip exploration ── */ - { - uint16_t base_d16 = gguf_fp32_to_fp16(d_current); - uint16_t best_d16 = base_d16; - float best_ulp_err = 1e30f; - - /* Try ±8 ULP neighborhood + sign flip = up to 34 candidates */ - uint16_t ulp_candidates[35]; - int n_ulp = 0; - for (int delta = -8; delta <= 8; delta++) { - int cand16 = (int)base_d16 + delta; - if (cand16 >= 0 && cand16 <= 0x7BFF) - ulp_candidates[n_ulp++] = (uint16_t)cand16; - } + /* ── FP16 ULP neighborhood search + sign-flip exploration ── */ { - float neg_d = -d_current; - uint16_t neg_d16 = gguf_fp32_to_fp16(neg_d); + uint16_t base_d16 = gguf_fp32_to_fp16(d_current); + uint16_t best_d16 = base_d16; + float best_ulp_err = 1e30f; + + /* Try ±8 ULP neighborhood + sign flip = up to 34 candidates */ + uint16_t ulp_candidates[35]; + int n_ulp = 0; for (int delta = -8; delta <= 8; delta++) { - int cand16 = (int)neg_d16 + delta; + int cand16 = (int)base_d16 + delta; if (cand16 >= 0 && cand16 <= 0x7BFF) ulp_candidates[n_ulp++] = (uint16_t)cand16; } - } - - for (int ui = 0; ui < n_ulp; ui++) { - float trial_d = gguf_fp16_to_fp32(ulp_candidates[ui]); - float trial_id = (fabsf(trial_d) > 1e-15f) ? 1.0f / trial_d : 0.0f; - float err = 0.0f; - float e_ulp[QK4_0]; - for (int j = 0; j < QK4_0; j++) { - int q = (int)(bw[j] * trial_id + 8.5f); - if (q < 0) q = 0; if (q > 15) q = 15; - float deq = ((float)q - 8.0f) * trial_d; - float w = (imat_importance) ? imat_importance[blk * QK4_0 + j] : 1.0f; - e_ulp[j] = bw[j] - deq; - err += e_ulp[j] * e_ulp[j] * w; + { + float neg_d = -d_current; + uint16_t neg_d16 = gguf_fp32_to_fp16(neg_d); + for (int delta = -8; delta <= 8; delta++) { + int cand16 = (int)neg_d16 + delta; + if (cand16 >= 0 && cand16 <= 0x7BFF) + ulp_candidates[n_ulp++] = (uint16_t)cand16; + } } - err += hex_spectral_penalty(e_ulp, QK4_0); - if (err < best_ulp_err) { - best_ulp_err = err; - best_d16 = ulp_candidates[ui]; + + for (int ui = 0; ui < n_ulp; ui++) { + float trial_d = gguf_fp16_to_fp32(ulp_candidates[ui]); + float trial_id = (fabsf(trial_d) > 1e-15f) ? 1.0f / trial_d : 0.0f; + float err = 0.0f; + float e_ulp[QK4_0]; + for (int j = 0; j < QK4_0; j++) { + int q = (int)(bw[j] * trial_id + 8.5f); + if (q < 0) q = 0; if (q > 15) q = 15; + float deq = ((float)q - 8.0f) * trial_d; + float w = (imat_importance) ? imat_importance[blk * QK4_0 + j] : 1.0f; + e_ulp[j] = bw[j] - deq; + err += e_ulp[j] * e_ulp[j] * w; + } + err += hex_spectral_penalty(e_ulp, QK4_0); + if (err < best_ulp_err) { + best_ulp_err = err; + best_d16 = ulp_candidates[ui]; + } } + d_current = gguf_fp16_to_fp32(best_d16); } - d_current = gguf_fp16_to_fp32(best_d16); - } - output[blk].d = gguf_fp32_to_fp16(d_current); - float actual_d = d_current; - float id = (fabsf(actual_d) > 1e-15f) ? 1.0f / actual_d : 0.0f; - - /* ── D₆ Hadamard Error Shaping with Simulated Annealing ── */ - int q_base[QK4_0], q_shaped[QK4_0]; - float q_cont[QK4_0]; - for (int j = 0; j < QK4_0; j++) { - q_cont[j] = bw[j] * id + 8.0f; - q_base[j] = (int)(q_cont[j] + 0.5f); - if (q_base[j] < 0) q_base[j] = 0; - if (q_base[j] > 15) q_base[j] = 15; - } - memcpy(q_shaped, q_base, QK4_0 * sizeof(int)); + output[blk].d = gguf_fp32_to_fp16(d_current); + float actual_d = d_current; + float id = (fabsf(actual_d) > 1e-15f) ? 1.0f / actual_d : 0.0f; - { - float e_live[QK4_0]; + /* ── D₆ Hadamard Error Shaping with Simulated Annealing ── */ + int q_base[QK4_0], q_shaped[QK4_0]; + float q_cont[QK4_0]; for (int j = 0; j < QK4_0; j++) { - float deq = ((float)q_shaped[j] - 8.0f) * actual_d; - e_live[j] = bw[j] - deq; + q_cont[j] = bw[j] * id + 8.0f; + q_base[j] = (int)(q_cont[j] + 0.5f); + if (q_base[j] < 0) q_base[j] = 0; + if (q_base[j] > 15) q_base[j] = 15; } + memcpy(q_shaped, q_base, QK4_0 * sizeof(int)); - float v_live[QK4_0 / 2]; - float vesica_cur = 0.0f, dc_cur = 0.0f; - for (int j = 0; j < QK4_0 / 2; j++) { - v_live[j] = e_live[j] + e_live[j + QK4_0 / 2]; - vesica_cur += v_live[j] * v_live[j]; - } - for (int j = 0; j < QK4_0; j++) dc_cur += e_live[j]; - float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; + { + float e_live[QK4_0]; + for (int j = 0; j < QK4_0; j++) { + float deq = ((float)q_shaped[j] - 8.0f) * actual_d; + e_live[j] = bw[j] - deq; + } - /* Deterministic greedy descent: only strict improvements. - * The previous SA acceptance called rand() inside an OpenMP - * parallel region (data race in the shared PRNG state, and - * non-reproducible output). Uphill moves were pointless anyway: - * the base-vs-shaped MSE guard below discards any shaped result - * that ends up worse, so accepted uphill excursions could only - * waste the pass budget or strand the descent. */ - for (int pass = 0; pass < QK4_0; pass++) { - int best_k = -1; - int best_q_alt = 0; - float best_delta = 0.0f; /* strictly positive threshold */ - - for (int k = 0; k < QK4_0; k++) { - int q_cur = q_shaped[k]; - int q_try = (q_cont[k] - (float)q_cur >= 0.0f) - ? q_cur + 1 : q_cur - 1; - if (q_try < 0 || q_try > 15) continue; - - float deq_try = ((float)q_try - 8.0f) * actual_d; - float e_new = bw[k] - deq_try; - float de = e_new - e_live[k]; - - int pi = (k < QK4_0 / 2) ? k : k - QK4_0 / 2; - float v_old = v_live[pi]; - float v_new = v_old + de; - - float vesica_alt = vesica_cur - v_old * v_old + v_new * v_new; - float dc_alt = dc_cur + de; - float metric_alt = 4.0f * vesica_alt + dc_alt * dc_alt; - - float delta = metric_cur - metric_alt; - if (delta > best_delta) { - best_delta = delta; - best_k = k; - best_q_alt = q_try; - } + float v_live[QK4_0 / 2]; + float vesica_cur = 0.0f, dc_cur = 0.0f; + for (int j = 0; j < QK4_0 / 2; j++) { + v_live[j] = e_live[j] + e_live[j + QK4_0 / 2]; + vesica_cur += v_live[j] * v_live[j]; } + for (int j = 0; j < QK4_0; j++) dc_cur += e_live[j]; + float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; + + /* Deterministic greedy descent: only strict improvements. + * The previous SA acceptance called rand() inside an OpenMP + * parallel region (data race in the shared PRNG state, and + * non-reproducible output). Uphill moves were pointless anyway: + * the base-vs-shaped MSE guard below discards any shaped result + * that ends up worse, so accepted uphill excursions could only + * waste the pass budget or strand the descent. */ + for (int pass = 0; pass < QK4_0; pass++) { + int best_k = -1; + int best_q_alt = 0; + float best_delta = 0.0f; /* strictly positive threshold */ + + for (int k = 0; k < QK4_0; k++) { + int q_cur = q_shaped[k]; + int q_try = (q_cont[k] - (float)q_cur >= 0.0f) + ? q_cur + 1 : q_cur - 1; + if (q_try < 0 || q_try > 15) continue; + + float deq_try = ((float)q_try - 8.0f) * actual_d; + float e_new = bw[k] - deq_try; + float de = e_new - e_live[k]; + + int pi = (k < QK4_0 / 2) ? k : k - QK4_0 / 2; + float v_old = v_live[pi]; + float v_new = v_old + de; + + float vesica_alt = vesica_cur - v_old * v_old + v_new * v_new; + float dc_alt = dc_cur + de; + float metric_alt = 4.0f * vesica_alt + dc_alt * dc_alt; + + float delta = metric_cur - metric_alt; + if (delta > best_delta) { + best_delta = delta; + best_k = k; + best_q_alt = q_try; + } + } - if (best_k < 0) break; /* converged — no improving flip */ + if (best_k < 0) break; /* converged — no improving flip */ - q_shaped[best_k] = best_q_alt; - { - float deq_commit = ((float)best_q_alt - 8.0f) * actual_d; - float e_new_commit = bw[best_k] - deq_commit; - float de_commit = e_new_commit - e_live[best_k]; + q_shaped[best_k] = best_q_alt; + { + float deq_commit = ((float)best_q_alt - 8.0f) * actual_d; + float e_new_commit = bw[best_k] - deq_commit; + float de_commit = e_new_commit - e_live[best_k]; - int pi_commit = (best_k < QK4_0 / 2) ? best_k : best_k - QK4_0 / 2; - float v_old_commit = v_live[pi_commit]; - float v_new_commit = v_old_commit + de_commit; + int pi_commit = (best_k < QK4_0 / 2) ? best_k : best_k - QK4_0 / 2; + float v_old_commit = v_live[pi_commit]; + float v_new_commit = v_old_commit + de_commit; - vesica_cur += v_new_commit * v_new_commit - v_old_commit * v_old_commit; - dc_cur += de_commit; - metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; + vesica_cur += v_new_commit * v_new_commit - v_old_commit * v_old_commit; + dc_cur += de_commit; + metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; - v_live[pi_commit] = v_new_commit; - e_live[best_k] = e_new_commit; + v_live[pi_commit] = v_new_commit; + e_live[best_k] = e_new_commit; + } } } - } - float err_base = 0.0f, err_shaped = 0.0f; - float e_gb[QK4_0], e_gs[QK4_0]; - for (int j = 0; j < QK4_0; j++) { - float w = (imat_importance) ? imat_importance[blk * QK4_0 + j] : 1.0f; - float deq_b = ((float)q_base[j] - 8.0f) * actual_d; - float deq_s = ((float)q_shaped[j] - 8.0f) * actual_d; - e_gb[j] = bw[j] - deq_b; - e_gs[j] = bw[j] - deq_s; - err_base += e_gb[j] * e_gb[j] * w; - err_shaped += e_gs[j] * e_gs[j] * w; - } - err_base += hex_spectral_penalty(e_gb, QK4_0); - err_shaped += hex_spectral_penalty(e_gs, QK4_0); - int *q_final = (err_shaped <= err_base) ? q_shaped : q_base; - - for (int j = 0; j < QK4_0 / 2; j++) { - int q0 = q_final[j]; - int q1 = q_final[j + QK4_0/2]; - output[blk].qs[j] = (uint8_t)(q0 | (q1 << 4)); - - float deq0 = ((float)q0 - 8.0f) * actual_d; - float deq1 = ((float)q1 - 8.0f) * actual_d; - total_err += (bw[j] - deq0) * (bw[j] - deq0) + (bw[j + QK4_0/2] - deq1) * (bw[j + QK4_0/2] - deq1); + float err_base = 0.0f, err_shaped = 0.0f; + float e_gb[QK4_0], e_gs[QK4_0]; + for (int j = 0; j < QK4_0; j++) { + float w = (imat_importance) ? imat_importance[blk * QK4_0 + j] : 1.0f; + float deq_b = ((float)q_base[j] - 8.0f) * actual_d; + float deq_s = ((float)q_shaped[j] - 8.0f) * actual_d; + e_gb[j] = bw[j] - deq_b; + e_gs[j] = bw[j] - deq_s; + err_base += e_gb[j] * e_gb[j] * w; + err_shaped += e_gs[j] * e_gs[j] * w; + } + err_base += hex_spectral_penalty(e_gb, QK4_0); + err_shaped += hex_spectral_penalty(e_gs, QK4_0); + int *q_final = (err_shaped <= err_base) ? q_shaped : q_base; + + for (int j = 0; j < QK4_0 / 2; j++) { + int q0 = q_final[j]; + int q1 = q_final[j + QK4_0/2]; + output[blk].qs[j] = (uint8_t)(q0 | (q1 << 4)); + + float deq0 = ((float)q0 - 8.0f) * actual_d; + float deq1 = ((float)q1 - 8.0f) * actual_d; + total_err += (bw[j] - deq0) * (bw[j] - deq0) + (bw[j + QK4_0/2] - deq1) * (bw[j + QK4_0/2] - deq1); + } } - } - *out_total_error = total_err; - free(greedy_d); - free(cand_errors); - free(cand_d16); - free(best_candidate); + *out_total_error = (float)total_err; + free(greedy_d); + free(cand_errors); + free(cand_d16); + free(best_candidate); } /* ════════════════════════════════════════════════════════════════════════ @@ -2486,13 +2545,14 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, const float *imat_importance, int verbose) { int64_t n_blocks = n_elements / QK8_0; - float total_err = 0.0f; + /* double: same float-saturation hazard as the Q2_K accumulator. */ + double total_err = 0.0; (void)verbose; float (*cand_errors)[Q8_N_CAND] = (float (*)[Q8_N_CAND]) - calloc(n_blocks, sizeof(float[Q8_N_CAND])); + calloc(n_blocks, sizeof(float[Q8_N_CAND])); uint16_t (*cand_d16)[Q8_N_CAND] = (uint16_t (*)[Q8_N_CAND]) - calloc(n_blocks, sizeof(uint16_t[Q8_N_CAND])); + calloc(n_blocks, sizeof(uint16_t[Q8_N_CAND])); int *best_candidate = (int *)malloc(n_blocks * sizeof(int)); if (!cand_errors || !cand_d16 || !best_candidate) { free(cand_errors); free(cand_d16); free(best_candidate); @@ -2590,12 +2650,12 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, for (int c = 0; c < Q8_N_CAND; c++) if (agg_errors[c] < min_err) min_err = agg_errors[c]; - double amp_re[6] = {0,0,0,0,0,0}; + double amp_re[6] = {0,0,0,0,0,0}; double amp_norm = 0.0; for (int ci = 0; ci < Q8_N_CAND; ci++) amp_re[Q8_CAND_TO_QUHIT[ci]] += - exp(-(double)(agg_errors[ci] - min_err) / - (2.0 * (double)temperature)); + exp(-(double)(agg_errors[ci] - min_err) / + (2.0 * (double)temperature)); for (int v = 0; v < 6; v++) amp_norm += amp_re[v] * amp_re[v]; if (amp_norm > 1e-30) { double inv = 1.0 / sqrt(amp_norm); @@ -2660,7 +2720,7 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, if (cand_errors[blk][c] < g_best) { g_best = cand_errors[blk][c]; g_cand = c; } - best_candidate[blk] = g_cand; + best_candidate[blk] = g_cand; } } @@ -2709,14 +2769,14 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, float best_delta = 0.0f; for (int k = 0; k < QK8_0; k++) { int q_try = (e_live[k] >= 0.0f) ? qs_shaped[k] + 1 - : qs_shaped[k] - 1; + : qs_shaped[k] - 1; if (q_try < -127 || q_try > 127) continue; float e_new = bw[k] - (float)q_try * d; float de = e_new - e_live[k]; int pi = (k < QK8_0 / 2) ? k : k - QK8_0 / 2; float v_new = v_live[pi] + de; float ves_a = vesica_cur - v_live[pi] * v_live[pi] - + v_new * v_new; + + v_new * v_new; float dc_a = dc_cur + de; float delta = metric_cur - (4.0f * ves_a + dc_a * dc_a); if (delta > best_delta) { @@ -2728,7 +2788,7 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, float e_new = bw[best_k] - (float)best_q_alt * d; float de = e_new - e_live[best_k]; int pi = (best_k < QK8_0 / 2) ? best_k - : best_k - QK8_0 / 2; + : best_k - QK8_0 / 2; float v_new = v_live[pi] + de; vesica_cur += v_new * v_new - v_live[pi] * v_live[pi]; dc_cur += de; @@ -2764,7 +2824,7 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, free(cand_errors); free(cand_d16); free(best_candidate); - if (out_total_error) *out_total_error = total_err; + if (out_total_error) *out_total_error = (float)total_err; } @@ -2789,42 +2849,393 @@ static inline void hex_derive_subscales(const float *scales, const float *mins, } } -static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, - BlockQ2K *output, float *out_total_error, - OptimizerMode opt_mode, - const float *imat_importance, - int verbose) +/* Extended-objective score of a fully-assembled Q2_K block state, measured + * against the ORIGINAL weights. This is the single scoring function used to + * compare finished candidates in the Phase-4 multi-start loop and to test + * the Phase-4.7 floor, so no two stages can disagree about which block is + * better. Writes the per-weight residuals into e_out (caller-owned, QK_K). + * At HEX_DC_LAMBDA = HEX_VW_LAMBDA = 0 this is exactly the importance- + * weighted SSE, i.e. the quantity the reported RMSE is built from. */ +static inline float hex_block_ext_err(const float *block_x, const float *bw, + float dm, float mm, + const uint8_t *scales, const uint8_t *L, + float *e_out) { - int64_t n_blocks = n_elements / QK_K; - float total_err = 0.0f; const int N_SUB = QK_K / 16; - - /* ── Outlier Clamping for WLS Seeds ── - * Protects the Phase 1 greedy seed from being violently warped by extreme - * >4.0 sigma outliers, which creates better centering for the grid search. */ - double t_sum_sq = 0.0, t_sum_4 = 0.0; - for (int64_t i = 0; i < n_elements; i++) { - double w2 = (double)weights[i] * (double)weights[i]; - t_sum_sq += w2; - t_sum_4 += w2 * w2; + float err = 0.0f; + for (int j = 0; j < N_SUB; j++) { + float d_sub = dm * (float)(scales[j] & 0xF); + float m_sub = mm * (float)(scales[j] >> 4); + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float w = bw ? bw[idx] : 1.0f; + /* deq = d*ls*q - dmin*lm; equals -m_sub when ls == 0 */ + float e = block_x[idx] - (d_sub * (float)L[idx] - m_sub); + e_out[idx] = e; + err += e * e * w; + } } - float w_sigma = sqrtf((float)(t_sum_sq / (double)n_elements)); - - /* ── Adaptive outlier clamp (kurtosis-driven) ── - * The fixed 3.5σ clamp suppressed the heavy-tail mass that dominates - * reconstruction error, inflating RMSE on near-Gaussian tensors that did - * not need clamping at all. Instead, gate the clamp on the tensor's raw - * kurtosis (Gaussian = 3): leave near-Gaussian tensors untouched and only - * apply a stabilising clamp to genuinely heavy-tailed tensors, where the - * final (d, dmin) refit later recovers fidelity against the UNCLIPPED - * weights anyway. */ - double t_var = t_sum_sq / (double)n_elements; - double t_kurt = (t_var > 1e-30) ? (t_sum_4 / (double)n_elements) / (t_var * t_var) : 3.0; - float clamp_sigma; - if (t_kurt <= 6.0) clamp_sigma = 1.0e9f; /* ~Gaussian: effectively no clamp */ - else if (t_kurt <= 20.0) clamp_sigma = 6.0f; /* moderately heavy tails */ - else clamp_sigma = 4.0f; /* very heavy tails: stabilise seed */ - float clamp_val = w_sigma * clamp_sigma; + return err + hex_spectral_penalty(e_out, QK_K); +} + +/* ══ PHASE 4.6 CORE: MONOTONE COORDINATE-DESCENT POLISH (RMSE-guaranteed) ══ + * + * Extracted from the Phase-4 assembly loop so it can be invoked from more + * than one place: once per multi-start candidate, and again after the + * Phase-4.7 candidate floor replaces the block wholesale (previously the + * floor adopted a RAW candidate that had never been polished at all). + * + * Runs alternating coordinate descent on the TRUE objective (importance- + * weighted SSE against the ORIGINAL weights, plus the spectral terms when + * HEX_DC_LAMBDA / HEX_VW_LAMBDA are non-zero): + * + * (1) For each 16-weight sub-block, an exact joint re-search of (ls, lm) + * over the full 16x16 grid with per-weight optimal q in {0..3}, + * committed only on strict improvement of the extended objective E. + * With lambda_dc = lambda_vw = 0 sub-blocks are independent given + * (d, dmin); with spectral terms active the coupling (DC: all subs; + * fold: sub j <-> sub j^8) is handled exactly via live residual + * bookkeeping. + * (2) Closed-form weighted LS refit of the two fp16 scalars (d, dmin) + * with all codes held fixed, committed only on strict improvement. + * (3) On stall, a +/-HEX_POLISH_ULP joint fp16 micro-search of (d, dmin) + * on the true objective, then one more sweep so (1) can re-optimise + * the codes for the new scalars. + * + * All moves are accept-only-if-better on E, so E is monotonically + * non-increasing and the routine can never make a block worse than the + * state it was handed. The state space is finite (4-bit codes, fp16 + * scalars), so the loop terminates; it typically converges in 2-3 sweeps + * and breaks out well before HEX_POLISH_ITERS. + * + * block_x : QK_K original (unclipped, un-DC-shifted) weights + * bw : QK_K importance weights, or NULL for uniform + * dm_io : in/out d scalar (fp16-representable value) + * mm_io : in/out dmin scalar (fp16-representable value) + * pl_Ls : in/out 16 four-bit sub-block scale codes + * pl_Lm : in/out 16 four-bit sub-block min codes + * L : in/out QK_K two-bit weight codes + */ +static void hex_polish_block(const float *block_x, const float *bw, + float *dm_io, float *mm_io, + uint8_t *pl_Ls, uint8_t *pl_Lm, uint8_t *L) +{ + const int N_SUB = QK_K / 16; + float dm = *dm_io; + float mm = *mm_io; + + for (int pol_iter = 0; pol_iter < HEX_POLISH_ITERS; pol_iter++) { + int pol_improved = 0; + + /* ── (1) Exact per-sub-block (ls, lm, q) re-search on the + * EXTENDED objective. Under the spectral terms sub-blocks + * are no longer independent: every sub couples to all others + * through the DC term and to its fold partner (sub j ⊕ 8, + * i.e. weights i ↔ i+128) through vesica² − wave². The + * search therefore keeps live residuals pe[] and scores each + * candidate against the whole-block penalty with the partner + * residuals held fixed — exact coordinate descent on E. */ + float pe[QK_K]; + float sub_sse[16], sub_dc[16], pair_cross[8]; + float dc_tot = 0.0f, cross_tot = 0.0f; + for (int j = 0; j < N_SUB; j++) { + float d_sub = dm * (float)pl_Ls[j]; + float m_sub = mm * (float)pl_Lm[j]; + sub_sse[j] = 0.0f; + sub_dc[j] = 0.0f; + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float w = (bw) ? + bw[idx] : 1.0f; + /* deq = d·ls·q − dmin·lm; equals −m_sub at ls==0 */ + float e = block_x[idx] - (d_sub * (float)L[idx] - m_sub); + pe[idx] = e; + sub_sse[j] += e * e * w; + sub_dc[j] += e; + } + dc_tot += sub_dc[j]; + } + for (int p = 0; p < 8; p++) { + pair_cross[p] = 0.0f; + for (int k = 0; k < 16; k++) + pair_cross[p] += pe[16*p + k] * pe[16*(p+8) + k]; + cross_tot += pair_cross[p]; + } + + for (int j = 0; j < N_SUB; j++) { + const float *sx = block_x + 16 * j; + int pi = j & 7; /* fold-pair index */ + int pj = j ^ 8; /* partner sub-block */ + const float *ppe = pe + 16 * pj; /* partner residuals */ + float dc_rest = dc_tot - sub_dc[j]; + float cross_rest = cross_tot - pair_cross[pi]; + + /* Extended score of the CURRENT committed state */ + float best_sub = sub_sse[j] + + (HEX_DC_LAMBDA / (float)QK_K) * dc_tot * dc_tot + + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f * cross_tot; + int best_ls = -1, best_lm = 0; + uint8_t best_q[16]; + float best_e[16]; + float best_sse = 0.0f, best_dcc = 0.0f, best_cxc = 0.0f; + + for (int try_ls = 0; try_ls <= 15; try_ls++) { + float d_sub = dm * (float)try_ls; + for (int try_lm = 0; try_lm <= 15; try_lm++) { + float m_sub = mm * (float)try_lm; + float sub_err = 0.0f, dcc = 0.0f, cxc = 0.0f; + uint8_t q_loc[16]; + float e_loc[16]; + int aborted = 0; + for (int k = 0; k < 16; k++) { + float x = sx[k]; + float w = (bw) ? + bw[16*j + k] : 1.0f; + int q = 0; + if (d_sub >= 1e-15f) { + q = gguf_nearest_int((x + m_sub) / d_sub); + if (q < 0) q = 0; if (q > 3) q = 3; + } + q_loc[k] = (uint8_t)q; + /* deq = d·ls·q − dmin·lm; −m_sub at ls==0 */ + float e = x - (d_sub * (float)q - m_sub); + e_loc[k] = e; + sub_err += e * e * w; + dcc += e; + cxc += e * ppe[k]; + /* SSE-partial prune is a valid lower bound + * only while the spectral terms are ≥ 0, + * i.e. when the (signable) vw credit is off */ + if (HEX_VW_LAMBDA == 0.0f && + sub_err >= best_sub) { aborted = 1; break; } + } + if (aborted) continue; + float score = sub_err + + (HEX_DC_LAMBDA / (float)QK_K) + * (dc_rest + dcc) * (dc_rest + dcc) + + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f + * (cross_rest + cxc); + if (score < best_sub) { + best_sub = score; + best_ls = try_ls; + best_lm = try_lm; + memcpy(best_q, q_loc, 16); + memcpy(best_e, e_loc, sizeof(e_loc)); + best_sse = sub_err; + best_dcc = dcc; + best_cxc = cxc; + } + } + } + + if (best_ls >= 0) { /* strict improvement in E found */ + pl_Ls[j] = (uint8_t)best_ls; + pl_Lm[j] = (uint8_t)best_lm; + memcpy(L + 16 * j, best_q, 16); + memcpy(pe + 16 * j, best_e, sizeof(best_e)); + sub_sse[j] = best_sse; + sub_dc[j] = best_dcc; + pair_cross[pi] = best_cxc; + dc_tot = dc_rest + best_dcc; + cross_tot = cross_rest + best_cxc; + pol_improved = 1; + } + } + + /* ── (2) Closed-form (d, dmin) refit vs ORIGINAL, codes fixed ── */ + { + double pSaa = 0, pSab = 0, pSbb = 0, pSxa = 0, pSxb = 0; + double pA = 0, pB = 0, pS = 0; /* DC rank-1 augmentation */ + for (int j = 0; j < N_SUB; j++) { + float ls_f = (float)pl_Ls[j]; + float lm_f = (float)pl_Lm[j]; + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float x = block_x[idx]; + float w = (bw) ? + bw[idx] : 1.0f; + float a = ls_f * (float)L[idx]; + float b = lm_f; + pSaa += (double)w * a * a; + pSab += (double)w * a * b; + pSbb += (double)w * b * b; + pSxa += (double)w * x * a; + pSxb += (double)w * x * b; + pA += a; pB += b; pS += x; + } + } + { + double pw = (double)HEX_DC_LAMBDA / (double)QK_K; + pSaa += pw * pA * pA; pSab += pw * pA * pB; + pSbb += pw * pB * pB; pSxa += pw * pS * pA; + pSxb += pw * pS * pB; + } + double pdet = pSaa * pSbb - pSab * pSab; + if (fabs(pdet) > 1e-30) { + double d_ref = (pSbb * pSxa - pSab * pSxb) / pdet; + double m_ref = (pSab * pSxa - pSaa * pSxb) / pdet; + if (d_ref > 0.0) { + float dm_try = gguf_fp16_to_fp32( + gguf_fp32_to_fp16((float)d_ref)); + float mm_try = (m_ref > 0.0) + ? gguf_fp16_to_fp32( + gguf_fp32_to_fp16((float)m_ref)) + : mm; + float err_cur = 0.0f, err_try = 0.0f; + float e_pc[QK_K], e_pt[QK_K]; + for (int j = 0; j < N_SUB; j++) { + float ls_f = (float)pl_Ls[j]; + float lm_f = (float)pl_Lm[j]; + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float x = block_x[idx]; + float w = (bw) ? + bw[idx] : 1.0f; + float qf = (float)L[idx]; + float dc = dm * ls_f * qf - mm * lm_f; + float dt = dm_try * ls_f * qf - mm_try * lm_f; + e_pc[idx] = x - dc; + e_pt[idx] = x - dt; + err_cur += e_pc[idx] * e_pc[idx] * w; + err_try += e_pt[idx] * e_pt[idx] * w; + } + } + err_cur += hex_spectral_penalty(e_pc, QK_K); + err_try += hex_spectral_penalty(e_pt, QK_K); + if (err_try < err_cur) { + dm = dm_try; + mm = mm_try; + pol_improved = 1; + } + } + } + } + + if (!pol_improved) { + /* ── (3) ±2 ULP joint (d, dmin) micro-search vs ORIGINAL ── + * The closed-form refit rounds its real-valued optimum to + * fp16, which can land 1–2 ULP away from the best + * representable pair (and the earlier ±8 ULP search ran + * against the DC-shifted objective). With codes fixed, + * scan the (2·HEX_POLISH_ULP+1)² fp16 neighborhood on the + * true objective; + * accept only strict improvement, then loop once more so + * move (1) can re-optimise codes for the new scalars. + * Monotone ⇒ final RMSE can only decrease. */ + uint16_t base_d16 = gguf_fp32_to_fp16(dm); + uint16_t base_m16 = gguf_fp32_to_fp16(mm); + + float cur_err = 0.0f; + float e_u[QK_K]; + for (int j = 0; j < N_SUB; j++) { + float d_sub = dm * (float)pl_Ls[j]; + float m_sub = mm * (float)pl_Lm[j]; + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float w = (bw) ? + bw[idx] : 1.0f; + e_u[idx] = block_x[idx] - + (d_sub * (float)L[idx] - m_sub); + cur_err += e_u[idx] * e_u[idx] * w; + } + } + cur_err += hex_spectral_penalty(e_u, QK_K); + + float best_err = cur_err; + uint16_t best_d16 = base_d16, best_m16 = base_m16; + for (int dd = -HEX_POLISH_ULP; dd <= HEX_POLISH_ULP; dd++) { + int cd16 = (int)base_d16 + dd; + if (cd16 < 0 || cd16 > 0x7BFF) continue; + float t_dm = gguf_fp16_to_fp32((uint16_t)cd16); + for (int dmm = -HEX_POLISH_ULP; dmm <= HEX_POLISH_ULP; dmm++) { + if (dd == 0 && dmm == 0) continue; + int cm16 = (int)base_m16 + dmm; + if (cm16 < 0 || cm16 > 0x7BFF) continue; + float t_mm = gguf_fp16_to_fp32((uint16_t)cm16); + + float err = 0.0f; + /* SSE-partial prune valid only without the + * signable vesica/wave credit */ + for (int j = 0; + j < N_SUB && (HEX_VW_LAMBDA != 0.0f || + err < best_err); j++) { + float d_sub = t_dm * (float)pl_Ls[j]; + float m_sub = t_mm * (float)pl_Lm[j]; + for (int k = 0; k < 16; k++) { + int idx = 16 * j + k; + float w = (bw) ? + bw[idx] : 1.0f; + e_u[idx] = block_x[idx] - + (d_sub * (float)L[idx] - m_sub); + err += e_u[idx] * e_u[idx] * w; + } + } + if (HEX_DC_LAMBDA != 0.0f || HEX_VW_LAMBDA != 0.0f) + err = (err < best_err || HEX_VW_LAMBDA != 0.0f) + ? err + hex_spectral_penalty(e_u, QK_K) + : err; + if (err < best_err) { + best_err = err; + best_d16 = (uint16_t)cd16; + best_m16 = (uint16_t)cm16; + } + } + } + if (best_d16 != base_d16 || best_m16 != base_m16) { + dm = gguf_fp16_to_fp32(best_d16); + mm = gguf_fp16_to_fp32(best_m16); + pol_improved = 1; + } + } + + if (!pol_improved) break; /* converged on true objective */ + } + + *dm_io = dm; + *mm_io = mm; +} + +static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, + BlockQ2K *output, float *out_total_error, + OptimizerMode opt_mode, + const float *imat_importance, + int verbose) +{ + int64_t n_blocks = n_elements / QK_K; + /* double, not float: this accumulates one SSE term per 256-weight block + * over tensors with 10^5-10^6 blocks. In float, once the running total is + * ~2^24 times a typical per-block term the additions round to no-ops and + * the reported error silently saturates (under-reporting RMSE, and doing + * so worst on exactly the largest tensors). */ + double total_err = 0.0; + const int N_SUB = QK_K / 16; + + /* ── Outlier Clamping for WLS Seeds ── + * Protects the Phase 1 greedy seed from being violently warped by extreme + * >4.0 sigma outliers, which creates better centering for the grid search. */ + double t_sum_sq = 0.0, t_sum_4 = 0.0; + for (int64_t i = 0; i < n_elements; i++) { + double w2 = (double)weights[i] * (double)weights[i]; + t_sum_sq += w2; + t_sum_4 += w2 * w2; + } + float w_sigma = sqrtf((float)(t_sum_sq / (double)n_elements)); + + /* ── Adaptive outlier clamp (kurtosis-driven) ── + * The fixed 3.5σ clamp suppressed the heavy-tail mass that dominates + * reconstruction error, inflating RMSE on near-Gaussian tensors that did + * not need clamping at all. Instead, gate the clamp on the tensor's raw + * kurtosis (Gaussian = 3): leave near-Gaussian tensors untouched and only + * apply a stabilising clamp to genuinely heavy-tailed tensors, where the + * final (d, dmin) refit later recovers fidelity against the UNCLIPPED + * weights anyway. */ + double t_var = t_sum_sq / (double)n_elements; + double t_kurt = (t_var > 1e-30) ? (t_sum_4 / (double)n_elements) / (t_var * t_var) : 3.0; + float clamp_sigma; + if (t_kurt <= 6.0) clamp_sigma = 1.0e9f; /* ~Gaussian: effectively no clamp */ + else if (t_kurt <= 20.0) clamp_sigma = 6.0f; /* moderately heavy tails */ + else clamp_sigma = 4.0f; /* very heavy tails: stabilise seed */ + float clamp_val = w_sigma * clamp_sigma; /* ══════════════════════════════════════════════════════════════════ * PHASE 1: Greedy quantization — produce seed (d, dmin) per block @@ -2867,18 +3278,18 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, * double-counts magnitude. Without an imatrix, fall back to the * magnitude-relative heuristic. */ wt[l] = (imat_importance) - ? imp - : sqrtf(sigma2 + sx_clipped[l] * sx_clipped[l]); + ? imp + : sqrtf(sigma2 + sx_clipped[l] * sx_clipped[l]); seeds[blk].sw[j] += wt[l]; } seeds[blk].scales[j] = hpc_make_qkx2_quants(16, 3, sx_clipped, wt, - L + 16 * j, &seeds[blk].mins[j], Laux); + L + 16 * j, &seeds[blk].mins[j], Laux); } seeds[blk].dm = hpc_make_qp_quants(N_SUB, 15, seeds[blk].scales, - seeds[blk].Ls, seeds[blk].sw); + seeds[blk].Ls, seeds[blk].sw); seeds[blk].mm = hpc_make_qp_quants(N_SUB, 15, seeds[blk].mins, - seeds[blk].Lm, seeds[blk].sw); + seeds[blk].Lm, seeds[blk].sw); seeds[blk].d_fp16 = gguf_fp32_to_fp16(seeds[blk].dm); seeds[blk].dmin_fp16 = gguf_fp32_to_fp16(seeds[blk].mm); } @@ -2914,11 +3325,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, uint16_t (*candidate_dmin)[TOTAL_SCALE_CANDIDATES] = NULL; candidate_errors = (float (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(float[TOTAL_SCALE_CANDIDATES])); + sizeof(float[TOTAL_SCALE_CANDIDATES])); candidate_d = (uint16_t (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); + sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); candidate_dmin = (uint16_t (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); + sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); /* NOTE: the per-candidate sub-scale codes (Ls/Lm) are NOT stored. * They are a pure function of (seeds[blk].scales/mins, candidate fp16 * d/dmin) and are re-derived where needed. Storing them cost @@ -2968,7 +3379,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, for (int k = 0; k < 16; k++) { float x = clipped_block_x[16*j+k]; float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j+k] : 1.0f; + imat_importance[blk * QK_K + 16*j+k] : 1.0f; float a = ls_f * (float)L_wls[16*j+k]; float b = lm_f; Saa += w * a * a; @@ -3056,7 +3467,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, err += e * e * w; } candidate_errors[blk][cidx] = - err + hex_spectral_penalty(e_arr, QK_K); + err + hex_spectral_penalty(e_arr, QK_K); } } } @@ -3069,682 +3480,821 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, for (int64_t i = 0; i < n_blocks; i++) best_candidate[i] = 11 * N_CAND_M + 11; /* index 11 = 1.0 multiplier */ - if (opt_mode != OPT_MSE && n_blocks >= 2) { - int64_t graph_blocks = (n_blocks > 2000) ? 2000 : n_blocks; - int64_t stride = n_blocks / graph_blocks; - float temperature = 0.5f; - int64_t n_sites = graph_blocks * QUHITS_PER_BLOCK; + if (opt_mode != OPT_MSE && n_blocks >= 2) { + int64_t graph_blocks = (n_blocks > 2000) ? 2000 : n_blocks; + int64_t stride = n_blocks / graph_blocks; + float temperature = 0.5f; + int64_t n_sites = graph_blocks * QUHITS_PER_BLOCK; - HPCGraph *graph = hpc_create(n_sites); - if (graph) { - for (int64_t i = 0; i < n_sites; i++) - triality_dft(&graph->locals[i]); + HPCGraph *graph = hpc_create(n_sites); + if (graph) { + for (int64_t i = 0; i < n_sites; i++) + triality_dft(&graph->locals[i]); - { - double err_accum = 0.0; - int err_count = 0; - for (int64_t gi = 0; gi < graph_blocks && gi < 100; gi++) { - int64_t blk = gi * stride; - float max_e = 0.0f; - for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) - if (candidate_errors[blk][c] > max_e) - max_e = candidate_errors[blk][c]; - err_accum += (double)max_e; - err_count++; - } - if (err_count > 0) { - float median_err = (float)(err_accum / err_count); - temperature = median_err * 0.1f; - if (temperature < 1e-10f) temperature = 1e-10f; + { + double err_accum = 0.0; + int err_count = 0; + for (int64_t gi = 0; gi < graph_blocks && gi < 100; gi++) { + int64_t blk = gi * stride; + float max_e = 0.0f; + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) + if (candidate_errors[blk][c] > max_e) + max_e = candidate_errors[blk][c]; + err_accum += (double)max_e; + err_count++; + } + if (err_count > 0) { + float median_err = (float)(err_accum / err_count); + temperature = median_err * 0.1f; + if (temperature < 1e-10f) temperature = 1e-10f; + } } - } - for (int64_t i = 0; i < graph_blocks; i++) { - float agg_errors[TOTAL_SCALE_CANDIDATES]; - for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) agg_errors[c] = 0.0f; + for (int64_t i = 0; i < graph_blocks; i++) { + float agg_errors[TOTAL_SCALE_CANDIDATES]; + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) agg_errors[c] = 0.0f; - int64_t blk_start = i * stride; - int64_t blk_end = blk_start + stride; - if (blk_end > n_blocks) blk_end = n_blocks; - int64_t group_size = blk_end - blk_start; + int64_t blk_start = i * stride; + int64_t blk_end = blk_start + stride; + if (blk_end > n_blocks) blk_end = n_blocks; + int64_t group_size = blk_end - blk_start; - for (int64_t b = blk_start; b < blk_end; b++) { - for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) - agg_errors[c] += candidate_errors[b][c]; - } - if (group_size > 1) { - float inv_gs = 1.0f / (float)group_size; + for (int64_t b = blk_start; b < blk_end; b++) { + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) + agg_errors[c] += candidate_errors[b][c]; + } + if (group_size > 1) { + float inv_gs = 1.0f / (float)group_size; + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) + agg_errors[c] *= inv_gs; + } + + float min_err = 1e30f; for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) - agg_errors[c] *= inv_gs; - } + if (agg_errors[c] < min_err) + min_err = agg_errors[c]; - float min_err = 1e30f; - for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) - if (agg_errors[c] < min_err) - min_err = agg_errors[c]; - - double coarse_re[6]; - double coarse_norm = 0.0; - for (int qi = 0; qi < 6; qi++) coarse_re[qi] = 0.0; - for (int di = 0; di < N_CAND_D; di++) { - int qi = CAND_TO_QUHIT[di]; - for (int mi = 0; mi < N_CAND_M; mi++) { - int cidx = di * N_CAND_M + mi; - coarse_re[qi] += exp(-(double)(agg_errors[cidx] - min_err) / - (2.0 * (double)temperature)); + double coarse_re[6]; + double coarse_norm = 0.0; + for (int qi = 0; qi < 6; qi++) coarse_re[qi] = 0.0; + for (int di = 0; di < N_CAND_D; di++) { + int qi = CAND_TO_QUHIT[di]; + for (int mi = 0; mi < N_CAND_M; mi++) { + int cidx = di * N_CAND_M + mi; + coarse_re[qi] += exp(-(double)(agg_errors[cidx] - min_err) / + (2.0 * (double)temperature)); + } + } + for (int qi = 0; qi < 6; qi++) coarse_norm += coarse_re[qi] * coarse_re[qi]; + if (coarse_norm > 1e-30) { + double inv = 1.0 / sqrt(coarse_norm); + for (int v = 0; v < 6; v++) coarse_re[v] *= inv; } - } - for (int qi = 0; qi < 6; qi++) coarse_norm += coarse_re[qi] * coarse_re[qi]; - if (coarse_norm > 1e-30) { - double inv = 1.0 / sqrt(coarse_norm); - for (int v = 0; v < 6; v++) coarse_re[v] *= inv; - } - double fine_re[6]; - double fine_norm = 0.0; - for (int qi = 0; qi < 6; qi++) fine_re[qi] = 0.0; - for (int mi = 0; mi < N_CAND_M; mi++) { - int qi = CAND_TO_QUHIT[mi]; - for (int di = 0; di < N_CAND_D; di++) { - int cidx = di * N_CAND_M + mi; - fine_re[qi] += exp(-(double)(agg_errors[cidx] - min_err) / - (2.0 * (double)temperature)); + double fine_re[6]; + double fine_norm = 0.0; + for (int qi = 0; qi < 6; qi++) fine_re[qi] = 0.0; + for (int mi = 0; mi < N_CAND_M; mi++) { + int qi = CAND_TO_QUHIT[mi]; + for (int di = 0; di < N_CAND_D; di++) { + int cidx = di * N_CAND_M + mi; + fine_re[qi] += exp(-(double)(agg_errors[cidx] - min_err) / + (2.0 * (double)temperature)); + } + } + for (int qi = 0; qi < 6; qi++) fine_norm += fine_re[qi] * fine_re[qi]; + if (fine_norm > 1e-30) { + double inv = 1.0 / sqrt(fine_norm); + for (int v = 0; v < 6; v++) fine_re[v] *= inv; } - } - for (int qi = 0; qi < 6; qi++) fine_norm += fine_re[qi] * fine_re[qi]; - if (fine_norm > 1e-30) { - double inv = 1.0 / sqrt(fine_norm); - for (int v = 0; v < 6; v++) fine_re[v] *= inv; - } - int64_t s0 = 2 * i, s1 = 2 * i + 1; - for (int v = 0; v < 6; v++) { - graph->locals[s0].edge_re[v] = coarse_re[v]; - graph->locals[s0].edge_im[v] = 0.0; - graph->locals[s1].edge_re[v] = fine_re[v]; - graph->locals[s1].edge_im[v] = 0.0; + int64_t s0 = 2 * i, s1 = 2 * i + 1; + for (int v = 0; v < 6; v++) { + graph->locals[s0].edge_re[v] = coarse_re[v]; + graph->locals[s0].edge_im[v] = 0.0; + graph->locals[s1].edge_re[v] = fine_re[v]; + graph->locals[s1].edge_im[v] = 0.0; + } + graph->locals[s0].primary = VIEW_EDGE; + graph->locals[s0].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; + graph->locals[s0].delta_valid = 0; + triality_update_mask(&graph->locals[s0]); + graph->locals[s1].primary = VIEW_EDGE; + graph->locals[s1].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; + graph->locals[s1].delta_valid = 0; + triality_update_mask(&graph->locals[s1]); } - graph->locals[s0].primary = VIEW_EDGE; - graph->locals[s0].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; - graph->locals[s0].delta_valid = 0; - triality_update_mask(&graph->locals[s0]); - graph->locals[s1].primary = VIEW_EDGE; - graph->locals[s1].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; - graph->locals[s1].delta_valid = 0; - triality_update_mask(&graph->locals[s1]); - } - for (int64_t i = 0; i < graph_blocks; i++) { - hpc_cz(graph, 2 * i, 2 * i + 1); - if (i + 1 < graph_blocks) { - hpc_cz(graph, 2 * i, 2 * (i + 1)); - hpc_cz(graph, 2 * i + 1, 2 * (i + 1) + 1); + for (int64_t i = 0; i < graph_blocks; i++) { + hpc_cz(graph, 2 * i, 2 * i + 1); + if (i + 1 < graph_blocks) { + hpc_cz(graph, 2 * i, 2 * (i + 1)); + hpc_cz(graph, 2 * i + 1, 2 * (i + 1) + 1); + } } - } - double (*shor_marg)[6] = (double (*)[6])calloc(n_sites, sizeof(double[6])); - int *shor_measured = (int *)calloc(n_sites, sizeof(int)); + double (*shor_marg)[6] = (double (*)[6])calloc(n_sites, sizeof(double[6])); + int *shor_measured = (int *)calloc(n_sites, sizeof(int)); - sieve_measure_graph(graph, n_sites, shor_marg, shor_measured, 1); + sieve_measure_graph(graph, n_sites, shor_marg, shor_measured, 1); - double (*coarse_marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); - double (*fine_marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); + double (*coarse_marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); + double (*fine_marg)[6] = (double (*)[6])calloc(graph_blocks, sizeof(double[6])); - for (int64_t i = 0; i < graph_blocks; i++) { - for (int v = 0; v < 6; v++) { - coarse_marg[i][v] = shor_marg[2 * i][v]; - fine_marg[i][v] = shor_marg[2 * i + 1][v]; + for (int64_t i = 0; i < graph_blocks; i++) { + for (int v = 0; v < 6; v++) { + coarse_marg[i][v] = shor_marg[2 * i][v]; + fine_marg[i][v] = shor_marg[2 * i + 1][v]; + } } - } - - free(shor_marg); - free(shor_measured); - - /* ══════════════════════════════════════════════════════════════ - * PHASE 3 — DETERMINISTIC VITERBI DP - * - * Replaces the probabilistic beam-search + Born-rule Monte-Carlo - * shots with an exact, fully-deterministic DP over the 36-state - * sieve quhit space (6 coarse bins × 6 fine bins). - * - * For each graph block i and combined state s = qi_d*6 + qi_m: - * - * bin_best_err[i][s] = min candidate error in that (d,m)-bin - * aggregated over the stride group - * bin_log_prior[i][s] = log P_coarse(qi_d) + log P_fine(qi_m) - * from sieve marginals → HPC prior bonus - * - * Local Viterbi cost (lower = better): - * vcost[i][s] = bin_best_err[i][s] - * − VITERBI_BETA × scale_err × bin_log_prior[i][s] - * - * Transition cost (cross-block smoothness prior): - * trans(s′→s) = VITERBI_ALPHA × scale_err - * × (|qi_d − qi_d′| + |qi_m − qi_m′|) - * - * DP recurrence: - * dp[0][s] = vcost[0][s] - * dp[i][s] = vcost[i][s] + min_{s′}(dp[i-1][s′] + trans(s′→s)) - * - * Traceback yields the globally optimal sequence of bin choices, - * which is then mapped to per-block best_candidate[] indices. - * A 5%-threshold greedy override rescues blocks where the local - * MSE-optimal candidate is meaningfully better than the bin winner. - * ══════════════════════════════════════════════════════════════ */ - #define VIT_N_STATES 36 /* 6 coarse × 6 fine quhit bins */ - #define VITERBI_BETA 0.25f /* log-prior bonus weight */ - #define VITERBI_ALPHA 0.08f /* cross-block smoothness penalty weight */ + free(shor_marg); + free(shor_measured); + + /* ══════════════════════════════════════════════════════════════ + * PHASE 3 — DETERMINISTIC VITERBI DP + * + * Replaces the probabilistic beam-search + Born-rule Monte-Carlo + * shots with an exact, fully-deterministic DP over the 36-state + * sieve quhit space (6 coarse bins × 6 fine bins). + * + * For each graph block i and combined state s = qi_d*6 + qi_m: + * + * bin_best_err[i][s] = min candidate error in that (d,m)-bin + * aggregated over the stride group + * bin_log_prior[i][s] = log P_coarse(qi_d) + log P_fine(qi_m) + * from sieve marginals → HPC prior bonus + * + * Local Viterbi cost (lower = better): + * vcost[i][s] = bin_best_err[i][s] + * − VITERBI_BETA × scale_err × bin_log_prior[i][s] + * + * Transition cost (cross-block smoothness prior): + * trans(s′→s) = VITERBI_ALPHA × scale_err + * × (|qi_d − qi_d′| + |qi_m − qi_m′|) + * + * DP recurrence: + * dp[0][s] = vcost[0][s] + * dp[i][s] = vcost[i][s] + min_{s′}(dp[i-1][s′] + trans(s′→s)) + * + * Traceback yields the globally optimal sequence of bin choices, + * which is then mapped to per-block best_candidate[] indices. + * A 5%-threshold greedy override rescues blocks where the local + * MSE-optimal candidate is meaningfully better than the bin winner. + * ══════════════════════════════════════════════════════════════ */ + + #define VIT_N_STATES 36 /* 6 coarse × 6 fine quhit bins */ + #define VITERBI_BETA 0.25f /* log-prior bonus weight */ + #define VITERBI_ALPHA 0.08f /* cross-block smoothness penalty weight */ - { - int64_t vit_gi, vit_b; - int vit_s, vit_sp; + { + int64_t vit_gi, vit_b; + int vit_s, vit_sp; - /* Per-graph-block per-state workspace */ - float (*vit_bin_err )[VIT_N_STATES] = + /* Per-graph-block per-state workspace */ + float (*vit_bin_err )[VIT_N_STATES] = (float (*)[VIT_N_STATES])malloc(graph_blocks * sizeof(float[VIT_N_STATES])); - int (*vit_bin_cand)[VIT_N_STATES] = + int (*vit_bin_cand)[VIT_N_STATES] = (int (*)[VIT_N_STATES])malloc(graph_blocks * sizeof(int [VIT_N_STATES])); - float (*vit_log_pri )[VIT_N_STATES] = + float (*vit_log_pri )[VIT_N_STATES] = (float (*)[VIT_N_STATES])malloc(graph_blocks * sizeof(float[VIT_N_STATES])); - float (*vit_dp )[VIT_N_STATES] = + float (*vit_dp )[VIT_N_STATES] = (float (*)[VIT_N_STATES])malloc(graph_blocks * sizeof(float[VIT_N_STATES])); - int (*vit_back )[VIT_N_STATES] = + int (*vit_back )[VIT_N_STATES] = (int (*)[VIT_N_STATES])malloc(graph_blocks * sizeof(int [VIT_N_STATES])); - /* ── Step A: build per-block per-bin statistics ── */ - for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { - double c_tot = 0.0, f_tot = 0.0; + /* ── Step A: build per-block per-bin statistics ── */ + for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { + double c_tot = 0.0, f_tot = 0.0; - for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { - vit_bin_err [vit_gi][vit_s] = 1e30f; - vit_bin_cand[vit_gi][vit_s] = -1; - } + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + vit_bin_err [vit_gi][vit_s] = 1e30f; + vit_bin_cand[vit_gi][vit_s] = -1; + } - /* Best candidate per (qi_d, qi_m) bin over stride group */ - for (vit_b = vit_gi * stride; - vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; - vit_b++) { - int vit_c; - for (vit_c = 0; vit_c < TOTAL_SCALE_CANDIDATES; vit_c++) { - int qi_d = CAND_TO_QUHIT[vit_c / N_CAND_M]; - int qi_m = CAND_TO_QUHIT[vit_c % N_CAND_M]; - vit_s = qi_d * 6 + qi_m; - float e = candidate_errors[vit_b][vit_c]; - if (e < vit_bin_err[vit_gi][vit_s]) { - vit_bin_err[vit_gi][vit_s] = e; - /* Canonical candidate = stride-rep block's best */ - if (vit_b == vit_gi * stride) - vit_bin_cand[vit_gi][vit_s] = vit_c; + /* Best candidate per (qi_d, qi_m) bin over stride group */ + for (vit_b = vit_gi * stride; + vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; + vit_b++) { + int vit_c; + for (vit_c = 0; vit_c < TOTAL_SCALE_CANDIDATES; vit_c++) { + int qi_d = CAND_TO_QUHIT[vit_c / N_CAND_M]; + int qi_m = CAND_TO_QUHIT[vit_c % N_CAND_M]; + vit_s = qi_d * 6 + qi_m; + float e = candidate_errors[vit_b][vit_c]; + if (e < vit_bin_err[vit_gi][vit_s]) { + vit_bin_err[vit_gi][vit_s] = e; + /* Canonical candidate = stride-rep block's best */ + if (vit_b == vit_gi * stride) + vit_bin_cand[vit_gi][vit_s] = vit_c; + } } } - } - /* HPC log-prior from sieve marginals */ - for (int v = 0; v < 6; v++) { - c_tot += coarse_marg[vit_gi][v]; - f_tot += fine_marg [vit_gi][v]; - } - for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { - int qi_d = vit_s / 6, qi_m = vit_s % 6; - double pc = (c_tot > 1e-30) - ? coarse_marg[vit_gi][qi_d] / c_tot : 1.0/6.0; - double pf = (f_tot > 1e-30) - ? fine_marg [vit_gi][qi_m] / f_tot : 1.0/6.0; - vit_log_pri[vit_gi][vit_s] = + /* HPC log-prior from sieve marginals */ + for (int v = 0; v < 6; v++) { + c_tot += coarse_marg[vit_gi][v]; + f_tot += fine_marg [vit_gi][v]; + } + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + int qi_d = vit_s / 6, qi_m = vit_s % 6; + double pc = (c_tot > 1e-30) + ? coarse_marg[vit_gi][qi_d] / c_tot : 1.0/6.0; + double pf = (f_tot > 1e-30) + ? fine_marg [vit_gi][qi_m] / f_tot : 1.0/6.0; + vit_log_pri[vit_gi][vit_s] = (float)(log(pc + 1e-30) + log(pf + 1e-30)); + } } - } - /* ── Step B: scale_err normaliser for transition cost ── */ - float vit_scale_err = 0.0f; - int vit_scale_cnt = 0; - for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { - for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { - if (vit_bin_err[vit_gi][vit_s] < 1e29f) { - vit_scale_err += vit_bin_err[vit_gi][vit_s]; - vit_scale_cnt++; + /* ── Step B: scale_err normaliser for transition cost ── */ + float vit_scale_err = 0.0f; + int vit_scale_cnt = 0; + for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + if (vit_bin_err[vit_gi][vit_s] < 1e29f) { + vit_scale_err += vit_bin_err[vit_gi][vit_s]; + vit_scale_cnt++; + } } } - } - vit_scale_err = (vit_scale_cnt > 0) - ? vit_scale_err / (float)vit_scale_cnt : 1e-10f; - if (vit_scale_err < 1e-20f) vit_scale_err = 1e-20f; + vit_scale_err = (vit_scale_cnt > 0) + ? vit_scale_err / (float)vit_scale_cnt : 1e-10f; + if (vit_scale_err < 1e-20f) vit_scale_err = 1e-20f; - /* ── Step C: Forward Viterbi pass ── */ + /* ── Step C: Forward Viterbi pass ── */ - /* Block 0 — no predecessor */ - for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { - float local = (vit_bin_err[0][vit_s] < 1e29f) + /* Block 0 — no predecessor */ + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + float local = (vit_bin_err[0][vit_s] < 1e29f) ? vit_bin_err[0][vit_s] - - VITERBI_BETA * vit_scale_err * vit_log_pri[0][vit_s] + - VITERBI_BETA * vit_scale_err * vit_log_pri[0][vit_s] : 1e30f; - vit_dp [0][vit_s] = local; - vit_back[0][vit_s] = -1; - } - - /* Blocks 1..graph_blocks-1 */ - for (vit_gi = 1; vit_gi < graph_blocks; vit_gi++) { - for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { - float local; - float best_pred = 1e30f; - int best_sp = 0; - int qi_d = vit_s / 6; - int qi_m = vit_s % 6; + vit_dp [0][vit_s] = local; + vit_back[0][vit_s] = -1; + } - if (vit_bin_err[vit_gi][vit_s] > 1e29f) { - vit_dp [vit_gi][vit_s] = 1e30f; - vit_back[vit_gi][vit_s] = 0; - continue; + /* Blocks 1..graph_blocks-1 */ + for (vit_gi = 1; vit_gi < graph_blocks; vit_gi++) { + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + float local; + float best_pred = 1e30f; + int best_sp = 0; + int qi_d = vit_s / 6; + int qi_m = vit_s % 6; + + if (vit_bin_err[vit_gi][vit_s] > 1e29f) { + vit_dp [vit_gi][vit_s] = 1e30f; + vit_back[vit_gi][vit_s] = 0; + continue; + } + local = vit_bin_err[vit_gi][vit_s] + - VITERBI_BETA * vit_scale_err * vit_log_pri[vit_gi][vit_s]; + + /* Min-cost predecessor with Manhattan transition penalty */ + for (vit_sp = 0; vit_sp < VIT_N_STATES; vit_sp++) { + float prev = vit_dp[vit_gi - 1][vit_sp]; + if (prev > 1e29f) continue; + int td = abs(qi_d - (vit_sp / 6)); + int tm = abs(qi_m - (vit_sp % 6)); + float trans = VITERBI_ALPHA * vit_scale_err * (float)(td + tm); + float total = prev + trans; + if (total < best_pred) { + best_pred = total; + best_sp = vit_sp; + } + } + vit_dp [vit_gi][vit_s] = (best_pred < 1e29f) + ? best_pred + local : 1e30f; + vit_back[vit_gi][vit_s] = best_sp; } - local = vit_bin_err[vit_gi][vit_s] - - VITERBI_BETA * vit_scale_err * vit_log_pri[vit_gi][vit_s]; - - /* Min-cost predecessor with Manhattan transition penalty */ - for (vit_sp = 0; vit_sp < VIT_N_STATES; vit_sp++) { - float prev = vit_dp[vit_gi - 1][vit_sp]; - if (prev > 1e29f) continue; - int td = abs(qi_d - (vit_sp / 6)); - int tm = abs(qi_m - (vit_sp % 6)); - float trans = VITERBI_ALPHA * vit_scale_err * (float)(td + tm); - float total = prev + trans; - if (total < best_pred) { - best_pred = total; - best_sp = vit_sp; + } + + /* ── Step D: Traceback ── */ + int *vit_path = (int *)malloc(graph_blocks * sizeof(int)); + { + int best_s = 0; + float best_f = vit_dp[graph_blocks - 1][0]; + for (vit_s = 1; vit_s < VIT_N_STATES; vit_s++) { + if (vit_dp[graph_blocks - 1][vit_s] < best_f) { + best_f = vit_dp[graph_blocks - 1][vit_s]; + best_s = vit_s; } } - vit_dp [vit_gi][vit_s] = (best_pred < 1e29f) - ? best_pred + local : 1e30f; - vit_back[vit_gi][vit_s] = best_sp; + vit_path[graph_blocks - 1] = best_s; + for (vit_gi = graph_blocks - 2; vit_gi >= 0; vit_gi--) + vit_path[vit_gi] = vit_back[vit_gi + 1][vit_path[vit_gi + 1]]; } - } - /* ── Step D: Traceback ── */ - int *vit_path = (int *)malloc(graph_blocks * sizeof(int)); - { - int best_s = 0; - float best_f = vit_dp[graph_blocks - 1][0]; - for (vit_s = 1; vit_s < VIT_N_STATES; vit_s++) { - if (vit_dp[graph_blocks - 1][vit_s] < best_f) { - best_f = vit_dp[graph_blocks - 1][vit_s]; - best_s = vit_s; + /* ── Step E: Map Viterbi path → best_candidate[] ── */ + for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { + vit_s = vit_path[vit_gi]; + int qi_d = vit_s / 6; + int qi_m = vit_s % 6; + int64_t blk_rep = vit_gi * stride; + + /* Stride-representative block: use precomputed bin winner */ + if (vit_bin_cand[vit_gi][vit_s] >= 0) + best_candidate[blk_rep] = vit_bin_cand[vit_gi][vit_s]; + + /* Non-representative blocks in the stride group */ + for (vit_b = blk_rep + 1; + vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; + vit_b++) { + int vit_c; + float best_e = 1e30f; + int best_c = best_candidate[blk_rep]; + for (vit_c = 0; vit_c < TOTAL_SCALE_CANDIDATES; vit_c++) { + if (CAND_TO_QUHIT[vit_c / N_CAND_M] != qi_d) continue; + if (CAND_TO_QUHIT[vit_c % N_CAND_M] != qi_m) continue; + if (candidate_errors[vit_b][vit_c] < best_e) { + best_e = candidate_errors[vit_b][vit_c]; + best_c = vit_c; + } + } + best_candidate[vit_b] = best_c; } } - vit_path[graph_blocks - 1] = best_s; - for (vit_gi = graph_blocks - 2; vit_gi >= 0; vit_gi--) - vit_path[vit_gi] = vit_back[vit_gi + 1][vit_path[vit_gi + 1]]; - } - /* ── Step E: Map Viterbi path → best_candidate[] ── */ - for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { - vit_s = vit_path[vit_gi]; - int qi_d = vit_s / 6; - int qi_m = vit_s % 6; - int64_t blk_rep = vit_gi * stride; - - /* Stride-representative block: use precomputed bin winner */ - if (vit_bin_cand[vit_gi][vit_s] >= 0) - best_candidate[blk_rep] = vit_bin_cand[vit_gi][vit_s]; - - /* Non-representative blocks in the stride group */ - for (vit_b = blk_rep + 1; - vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; - vit_b++) { + /* ── Step F: 5 % greedy override (pure MSE safety net) ── */ + for (vit_b = 0; vit_b < n_blocks; vit_b++) { int vit_c; - float best_e = 1e30f; - int best_c = best_candidate[blk_rep]; + float cur_err = candidate_errors[vit_b][best_candidate[vit_b]]; + float g_best = cur_err; + int g_cand = best_candidate[vit_b]; for (vit_c = 0; vit_c < TOTAL_SCALE_CANDIDATES; vit_c++) { - if (CAND_TO_QUHIT[vit_c / N_CAND_M] != qi_d) continue; - if (CAND_TO_QUHIT[vit_c % N_CAND_M] != qi_m) continue; - if (candidate_errors[vit_b][vit_c] < best_e) { - best_e = candidate_errors[vit_b][vit_c]; - best_c = vit_c; + if (candidate_errors[vit_b][vit_c] < g_best) { + g_best = candidate_errors[vit_b][vit_c]; + g_cand = vit_c; } } - best_candidate[vit_b] = best_c; + if (g_best < cur_err * HEX_GREEDY_OVERRIDE_RATIO) + best_candidate[vit_b] = g_cand; } - } - /* ── Step F: 5 % greedy override (pure MSE safety net) ── */ - for (vit_b = 0; vit_b < n_blocks; vit_b++) { - int vit_c; - float cur_err = candidate_errors[vit_b][best_candidate[vit_b]]; - float g_best = cur_err; - int g_cand = best_candidate[vit_b]; - for (vit_c = 0; vit_c < TOTAL_SCALE_CANDIDATES; vit_c++) { - if (candidate_errors[vit_b][vit_c] < g_best) { - g_best = candidate_errors[vit_b][vit_c]; - g_cand = vit_c; - } - } - if (g_best < cur_err * HEX_GREEDY_OVERRIDE_RATIO) - best_candidate[vit_b] = g_cand; + free(vit_path); + free(vit_dp); + free(vit_back); + free(vit_bin_err); + free(vit_bin_cand); + free(vit_log_pri); } - free(vit_path); - free(vit_dp); - free(vit_back); - free(vit_bin_err); - free(vit_bin_cand); - free(vit_log_pri); + free(coarse_marg); + free(fine_marg); + hpc_destroy(graph); } - - free(coarse_marg); - free(fine_marg); - hpc_destroy(graph); - } - } else { - for (int64_t blk = 0; blk < n_blocks; blk++) { - float best_err = candidate_errors[blk][0]; - int best_idx = 0; - for (int c = 1; c < TOTAL_SCALE_CANDIDATES; c++) { - if (candidate_errors[blk][c] < best_err) { - best_err = candidate_errors[blk][c]; - best_idx = c; + } else { + for (int64_t blk = 0; blk < n_blocks; blk++) { + float best_err = candidate_errors[blk][0]; + int best_idx = 0; + for (int c = 1; c < TOTAL_SCALE_CANDIDATES; c++) { + if (candidate_errors[blk][c] < best_err) { + best_err = candidate_errors[blk][c]; + best_idx = c; + } } + best_candidate[blk] = best_idx; } - best_candidate[blk] = best_idx; } - } - - /* ═══���══════════════════════════════════════════════════════════════ - * PHASE 3.9 — ROLLING DC BOUNDARY CONDITION PRE-PASS - * - * Transforms the tensor from a collection of isolated 256-element - * Q2_K superblocks into a single, continuous error-cancelling waveform. - * - * After Phase 3 has selected the optimal (d, dmin) candidate for every - * block, this sequential pass computes the net DC residual left by each - * block using a cheap round-nearest forward quantization, then feeds the - * negated, exponentially-decayed residual as a correction bias into the - * WLS solver of the immediately following block. - * - * Mathematically, for block N with final DC residual R_N = Σ εᵢ: - * - * dc_bias[N+1] = −DC_DECAY × R_N / QK_K (per-element offset) - * - * Block N+1's WLS targets become x′ᵢ = xᵢ − dc_bias[N+1], steering the - * quantizer toward codes whose reconstruction deq ≈ x′, so that - * - * Σ (xᵢ − deqᵢ) ≈ dc_bias[N+1] × QK_K = −DC_DECAY × R_N - * - * The accumulated cross-block DC collapses geometrically: - * - * R₀, DC_DECAY·R₀, DC_DECAY²·R₀, … → 0 - * - * The result is written into block_dc_bias[n_blocks]. Phase 4 reads - * this array (safe: written sequentially before the parallel loop). - * ══════════════════════════════════════════════════════════════════ */ - - #define DC_DECAY 0.85f /* Boundary-condition leak factor (0 = isolated, 1 = full) */ - float *block_dc_bias = (float *)calloc(n_blocks, sizeof(float)); + /* ══════════════════════════════════════════════════════════════════ + * PHASE 3.9 — ROLLING DC BOUNDARY CONDITION PRE-PASS + * + * Transforms the tensor from a collection of isolated 256-element + * Q2_K superblocks into a single, continuous error-cancelling waveform. + * + * After Phase 3 has selected the optimal (d, dmin) candidate for every + * block, this sequential pass computes the net DC residual left by each + * block using a cheap round-nearest forward quantization, then feeds the + * negated, exponentially-decayed residual as a correction bias into the + * WLS solver of the immediately following block. + * + * Mathematically, for block N with final DC residual R_N = Σ εᵢ: + * + * dc_bias[N+1] = −DC_DECAY × R_N / QK_K (per-element offset) + * + * Block N+1's WLS targets become x′ᵢ = xᵢ − dc_bias[N+1], steering the + * quantizer toward codes whose reconstruction deq ≈ x′, so that + * + * Σ (xᵢ − deqᵢ) ≈ dc_bias[N+1] × QK_K = −DC_DECAY × R_N + * + * The accumulated cross-block DC collapses geometrically: + * + * R₀, DC_DECAY·R₀, DC_DECAY²·R₀, … → 0 + * + * The result is written into block_dc_bias[n_blocks]. Phase 4 reads + * this array (safe: written sequentially before the parallel loop). + * ══════════════════════════════════════════════════════════════════ */ - if (block_dc_bias) { - float rolling_dc = 0.0f; + /* Boundary-condition leak factor (0 = isolated, 1 = full). + * + * DEFAULT CHANGED 0.85f → 0.0f: OBJECTIVE-MISMATCH FIX. + * + * With DC_DECAY > 0 every stage that commits 2-bit codes — the + * ls_iter WLS/sieve loop, the ±ULP (d, dmin) search, the 16x16 + * (ls, lm) search, the greedy shaping descent and the Floyd- + * Steinberg diffusion — minimises error against adj_block_x, the + * DC-SHIFTED target. Reported RMSE is measured against block_x. + * The codes therefore settle at the optimum of a shifted objective, + * and only the scalar (d, dmin) refits target the true one (with all + * codes frozen). The Phase-4.6 polish repairs this partially, but it + * is coordinate descent from an already-displaced starting point. + * + * At 0.0f, adj_block_x == block_x identically and every stage in the + * pipeline optimises the same objective the RMSE report measures. + * Rebuild with -DDC_DECAY=0.85f to restore cross-block DC cancellation + * (which buys layer-output accuracy, not RMSE). */ + #ifndef DC_DECAY + #define DC_DECAY 0.0f + #endif + + /* At DC_DECAY == 0 the bias array is identically zero, so the whole + * sequential pre-pass is dead weight — a full extra forward + * quantisation of the tensor. Leave the pointer NULL; Phase 4 already + * treats NULL as "no adjustment". */ + float *block_dc_bias = (DC_DECAY != 0.0f) + ? (float *)calloc(n_blocks, sizeof(float)) + : NULL; + + if (block_dc_bias) { + float rolling_dc = 0.0f; + + for (int64_t blk = 0; blk < n_blocks; blk++) { + const float *bx = weights + blk * QK_K; + int cidx = best_candidate[blk]; + float dm0 = gguf_fp16_to_fp32(candidate_d [blk][cidx]); + float mm0 = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); + + uint8_t dc_Ls[16], dc_Lm[16]; + hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, + dm0, mm0, dc_Ls, dc_Lm); + + /* Bias applied to THIS block's WLS targets */ + float dc_bias = (DC_DECAY * rolling_dc) / (float)QK_K; + block_dc_bias[blk] = dc_bias; + + /* Quick round-nearest quant to estimate DC residual for NEXT block. + * We quantize the adjusted target x′ = x − dc_bias, then measure + * the residual of the ORIGINAL weight against the chosen code. */ + float dc_res = 0.0f; + int j, k; + for (j = 0; j < N_SUB; j++) { + float d_sub = dm0 * (float)dc_Ls[j]; + float m_sub = mm0 * (float)dc_Lm[j]; + for (k = 0; k < 16; k++) { + float x_adj = bx[16*j + k] - dc_bias; + int q = 0; + if (d_sub >= 1e-15f) { + q = gguf_nearest_int((x_adj + m_sub) / d_sub); + if (q < 0) q = 0; + if (q > 3) q = 3; + } + float deq = d_sub * (float)q - m_sub; + /* Residual against ORIGINAL weight (not adjusted) */ + dc_res += bx[16*j + k] - deq; + } + } + rolling_dc = dc_res; + } + } - for (int64_t blk = 0; blk < n_blocks; blk++) { - const float *bx = weights + blk * QK_K; - int cidx = best_candidate[blk]; - float dm0 = gguf_fp16_to_fp32(candidate_d [blk][cidx]); - float mm0 = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); + /* ══════════════════════════════════════════════════════════════════ + * PHASE 4: Assemble blocks via least-squares (d, dmin) extraction + * ══════════════════════════════════════════════════════════════════ */ + + int _n_omp_threads = 1; + #ifdef _OPENMP + _n_omp_threads = omp_get_max_threads(); + #endif + HPCGraph **_tl_graphs = (HPCGraph **)calloc(_n_omp_threads, sizeof(HPCGraph *)); + for (int _ti = 0; _ti < _n_omp_threads; _ti++) + _tl_graphs[_ti] = hpc_create(N_SUB); + + /* ── Sub-scale saturation diagnostics ── + * Ls / Lm are 4-bit codes: scales[j] = round(scale_j / d) clamped to + * [0, 15]. A sub-block that CLAMPS AT 15 is encoded with a scale that is + * too small for its dynamic range and will dominate the block's error; + * one that clamps at 0 is encoded as flat. Both are pathologies the + * (d, dmin) candidate grid is supposed to prevent, and if either is + * common the grid's tails are too narrow. Counted on the FINAL committed + * codes, which is what actually ships. */ + int64_t sat_ls_hi = 0, sat_lm_hi = 0, sat_ls_lo = 0, sat_subs = 0; + + #pragma omp parallel for schedule(dynamic, 64) \ + reduction(+:total_err, sat_ls_hi, sat_lm_hi, sat_ls_lo, sat_subs) + for (int64_t blk = 0; blk < n_blocks; blk++) { + const float *block_x = weights + blk * QK_K; + const float *blk_imat = imat_importance ? imat_importance + blk * QK_K + : NULL; - uint8_t dc_Ls[16], dc_Lm[16]; - hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, - dm0, mm0, dc_Ls, dc_Lm); - - /* Bias applied to THIS block's WLS targets */ - float dc_bias = (DC_DECAY * rolling_dc) / (float)QK_K; - block_dc_bias[blk] = dc_bias; - - /* Quick round-nearest quant to estimate DC residual for NEXT block. - * We quantize the adjusted target x′ = x − dc_bias, then measure - * the residual of the ORIGINAL weight against the chosen code. */ - float dc_res = 0.0f; - int j, k; - for (j = 0; j < N_SUB; j++) { - float d_sub = dm0 * (float)dc_Ls[j]; - float m_sub = mm0 * (float)dc_Lm[j]; - for (k = 0; k < 16; k++) { - float x_adj = bx[16*j + k] - dc_bias; - int q = 0; - if (d_sub >= 1e-15f) { - q = gguf_nearest_int((x_adj + m_sub) / d_sub); - if (q < 0) q = 0; - if (q > 3) q = 3; - } - float deq = d_sub * (float)q - m_sub; - /* Residual against ORIGINAL weight (not adjusted) */ - dc_res += bx[16*j + k] - deq; + /* ══════════════════════════════════════════════════════════════════ + * PHASE 4.0: MULTI-START CANDIDATE SHORTLIST + * + * The assembly chain below is monotone but LOCAL — it descends into + * whatever basin its starting (d, dmin) candidate occupies. Raw + * candidate score is only a weak predictor of which basin polishes + * best, because the candidates are scored with seed-derived Ls/Lm and + * nearest rounding, while the finished block gets an exact 16x16 + * (ls, lm) re-search, an fp16 micro-search and a closed-form refit. + * + * So: shortlist the K best-scoring candidates, run the FULL assembly + * from each, and keep whichever finishes lowest. HEX_TOPK_RESTARTS = 1 + * reproduces the previous single-start behaviour exactly. + * ══════════════════════════════════════════════════════════════════ */ + int topk_c[HEX_TOPK_RESTARTS]; + float topk_e[HEX_TOPK_RESTARTS]; + int n_topk = 0; + { + /* Insertion-ordered partial selection: topk_e stays ascending. */ + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) { + float e = candidate_errors[blk][c]; + int pos; + if (n_topk < HEX_TOPK_RESTARTS) { + pos = n_topk++; + } else if (e < topk_e[n_topk - 1]) { + pos = n_topk - 1; + } else { + continue; + } + while (pos > 0 && topk_e[pos - 1] > e) { + topk_e[pos] = topk_e[pos - 1]; + topk_c[pos] = topk_c[pos - 1]; + pos--; + } + topk_e[pos] = e; + topk_c[pos] = c; + } + /* Always include the Phase-3 selection. With + * HEX_GREEDY_OVERRIDE_RATIO = 1.0 it is the per-block argmin and + * therefore already at topk_c[0], but do not rely on that. */ + { + int bc = best_candidate[blk], found = 0; + for (int t = 0; t < n_topk; t++) + if (topk_c[t] == bc) { found = 1; break; } + if (!found && n_topk > 0) { + topk_c[n_topk - 1] = bc; + topk_e[n_topk - 1] = candidate_errors[blk][bc]; } } - rolling_dc = dc_res; } - } - /* ══════════════════════════════════════════════════════════════════ - * PHASE 4: Assemble blocks via least-squares (d, dmin) extraction - * ══════════════════════════════════════════════════════════════════ */ + /* Best FINISHED state seen across restarts */ + float ms_best_err = 1e30f; + uint16_t ms_best_d16 = 0, ms_best_m16 = 0; + uint8_t ms_best_scales[16]; + uint8_t ms_best_L[QK_K]; + int ms_have = 0; - int _n_omp_threads = 1; - #ifdef _OPENMP - _n_omp_threads = omp_get_max_threads(); - #endif - HPCGraph **_tl_graphs = (HPCGraph **)calloc(_n_omp_threads, sizeof(HPCGraph *)); - for (int _ti = 0; _ti < _n_omp_threads; _ti++) - _tl_graphs[_ti] = hpc_create(N_SUB); + for (int kstart = 0; kstart < n_topk; kstart++) { - #pragma omp parallel for schedule(dynamic, 64) reduction(+:total_err) - for (int64_t blk = 0; blk < n_blocks; blk++) { - const float *block_x = weights + blk * QK_K; - int cidx = best_candidate[blk]; - uint8_t Ls_blk[16], Lm_blk[16]; + /* topk_e is ascending, so once one candidate is too far above the + * best finished result, every later one is too — break, not continue. */ + if (ms_have && topk_e[kstart] > ms_best_err * HEX_TOPK_PRUNE) break; - /* ── Rolling DC boundary condition ────────────────────���───────── - * dc_adj shifts every WLS target in this block so that the net - * quantisation error steers toward cancelling the previous block's - * DC residual (written by the sequential Phase 3.9 pre-pass). */ - float dc_adj = (block_dc_bias) ? block_dc_bias[blk] : 0.0f; + int cidx = topk_c[kstart]; + uint8_t Ls_blk[16], Lm_blk[16]; - /* Adjusted weight view — WLS and sieve work on this array; - * the final error is always reported against the original block_x. */ - float adj_block_x[QK_K]; - { - int _i; - for (_i = 0; _i < QK_K; _i++) - adj_block_x[_i] = block_x[_i] - dc_adj; - } + /* ── Rolling DC boundary condition ────────────────────────────── + * dc_adj shifts every WLS target in this block so that the net + * quantisation error steers toward cancelling the previous block's + * DC residual (written by the sequential Phase 3.9 pre-pass). */ + float dc_adj = (block_dc_bias) ? block_dc_bias[blk] : 0.0f; - float dm = gguf_fp16_to_fp32(candidate_d[blk][cidx]); - float mm = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); + /* Adjusted weight view — WLS and sieve work on this array; + * the final error is always reported against the original block_x. */ + float adj_block_x[QK_K]; + { + int _i; + for (_i = 0; _i < QK_K; _i++) + adj_block_x[_i] = block_x[_i] - dc_adj; + } - hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, - dm, mm, Ls_blk, Lm_blk); + float dm = gguf_fp16_to_fp32(candidate_d[blk][cidx]); + float mm = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); - uint16_t prev_dm16 = 0, prev_mm16 = 0; - for (int ls_iter = 0; ls_iter < 5; ls_iter++) { + hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, + dm, mm, Ls_blk, Lm_blk); - uint8_t state_ls[N_SUB][6]; - uint8_t state_lm[N_SUB][6]; - float state_err[N_SUB][6]; + uint16_t prev_dm16 = 0, prev_mm16 = 0; + for (int ls_iter = 0; ls_iter < 5; ls_iter++) { - for (int j = 0; j < N_SUB; j++) { - const float *sx = adj_block_x + 16 * j; - for (int v = 0; v < 6; v++) state_err[j][v] = 1e30f; + uint8_t state_ls[N_SUB][6]; + uint8_t state_lm[N_SUB][6]; + float state_err[N_SUB][6]; - for (int try_ls = 0; try_ls <= 15; try_ls++) { - float d_sub = dm * (float)try_ls; - for (int try_lm = 0; try_lm <= 15; try_lm++) { - float m_sub = mm * (float)try_lm; - float sub_err = 0.0f; + for (int j = 0; j < N_SUB; j++) { + const float *sx = adj_block_x + 16 * j; + for (int v = 0; v < 6; v++) state_err[j][v] = 1e30f; - for (int k = 0; k < 16; k++) { - float x = sx[k]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j + k] : 1.0f; - int q = 0; - if (d_sub >= 1e-15f) { - q = gguf_nearest_int((x + m_sub) / d_sub); - if (q < 0) q = 0; if (q > 3) q = 3; + for (int try_ls = 0; try_ls <= 15; try_ls++) { + float d_sub = dm * (float)try_ls; + for (int try_lm = 0; try_lm <= 15; try_lm++) { + float m_sub = mm * (float)try_lm; + float sub_err = 0.0f; + + for (int k = 0; k < 16; k++) { + float x = sx[k]; + float w = (imat_importance) ? + imat_importance[blk * QK_K + 16*j + k] : 1.0f; + int q = 0; + if (d_sub >= 1e-15f) { + q = gguf_nearest_int((x + m_sub) / d_sub); + if (q < 0) q = 0; if (q > 3) q = 3; + } + float deq = d_sub * (float)q - m_sub; + float diff = x - deq; + sub_err += diff * diff * w; } - float deq = d_sub * (float)q - m_sub; - float diff = x - deq; - sub_err += diff * diff * w; - } - for (int v = 0; v < 6; v++) { - if (sub_err < state_err[j][v]) { - for (int u = 5; u > v; u--) { - state_err[j][u] = state_err[j][u-1]; - state_ls[j][u] = state_ls[j][u-1]; - state_lm[j][u] = state_lm[j][u-1]; + for (int v = 0; v < 6; v++) { + if (sub_err < state_err[j][v]) { + for (int u = 5; u > v; u--) { + state_err[j][u] = state_err[j][u-1]; + state_ls[j][u] = state_ls[j][u-1]; + state_lm[j][u] = state_lm[j][u-1]; + } + state_err[j][v] = sub_err; + state_ls[j][v] = (uint8_t)try_ls; + state_lm[j][v] = (uint8_t)try_lm; + break; } - state_err[j][v] = sub_err; - state_ls[j][v] = (uint8_t)try_ls; - state_lm[j][v] = (uint8_t)try_lm; - break; } } } } - } - int _tid = 0; - #ifdef _OPENMP - _tid = omp_get_thread_num(); - #endif - HPCGraph *sg = _tl_graphs[_tid]; - hpc_reset_for_subblock(sg, N_SUB); - { - float min_sub_err[N_SUB]; - for (int j = 0; j < N_SUB; j++) min_sub_err[j] = state_err[j][0]; + int _tid = 0; + #ifdef _OPENMP + _tid = omp_get_thread_num(); + #endif + HPCGraph *sg = _tl_graphs[_tid]; + hpc_reset_for_subblock(sg, N_SUB); + { + float min_sub_err[N_SUB]; + for (int j = 0; j < N_SUB; j++) min_sub_err[j] = state_err[j][0]; - for (int j = 0; j < N_SUB; j++) { - triality_dft(&sg->locals[j]); - double amp_re[6]; - double amp_norm = 0.0; - for (int v = 0; v < 6; v++) { - float err_spread = state_err[j][5] - state_err[j][0]; - float sub_temp = (err_spread > 1e-15f) ? err_spread * 0.3f : 0.1f; - if (sub_temp < 1e-12f) sub_temp = 1e-12f; - amp_re[v] = exp(-(double)(state_err[j][v] - min_sub_err[j]) / (double)sub_temp); - amp_norm += amp_re[v] * amp_re[v]; - } - if (amp_norm > 1e-30) { - double inv = 1.0 / sqrt(amp_norm); - for (int v = 0; v < 6; v++) amp_re[v] *= inv; - } - for (int v = 0; v < 6; v++) { - sg->locals[j].edge_re[v] = amp_re[v]; - sg->locals[j].edge_im[v] = 0.0; + for (int j = 0; j < N_SUB; j++) { + triality_dft(&sg->locals[j]); + double amp_re[6]; + double amp_norm = 0.0; + for (int v = 0; v < 6; v++) { + float err_spread = state_err[j][5] - state_err[j][0]; + float sub_temp = (err_spread > 1e-15f) ? err_spread * 0.3f : 0.1f; + if (sub_temp < 1e-12f) sub_temp = 1e-12f; + amp_re[v] = exp(-(double)(state_err[j][v] - min_sub_err[j]) / (double)sub_temp); + amp_norm += amp_re[v] * amp_re[v]; + } + if (amp_norm > 1e-30) { + double inv = 1.0 / sqrt(amp_norm); + for (int v = 0; v < 6; v++) amp_re[v] *= inv; + } + for (int v = 0; v < 6; v++) { + sg->locals[j].edge_re[v] = amp_re[v]; + sg->locals[j].edge_im[v] = 0.0; + } + sg->locals[j].primary = VIEW_EDGE; + sg->locals[j].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; + sg->locals[j].delta_valid = 0; + triality_update_mask(&sg->locals[j]); } - sg->locals[j].primary = VIEW_EDGE; - sg->locals[j].dirty = DIRTY_VERTEX | DIRTY_DIAGONAL | DIRTY_FOLDED; - sg->locals[j].delta_valid = 0; - triality_update_mask(&sg->locals[j]); - } - for (int j = 0; j < N_SUB - 1; j++) - hpc_cz(sg, j, j + 1); + for (int j = 0; j < N_SUB - 1; j++) + hpc_cz(sg, j, j + 1); - double sub_marg[N_SUB][6]; - int sub_measured[N_SUB]; - memset(sub_marg, 0, sizeof(sub_marg)); - memset(sub_measured, 0, sizeof(sub_measured)); + double sub_marg[N_SUB][6]; + int sub_measured[N_SUB]; + memset(sub_marg, 0, sizeof(sub_marg)); + memset(sub_measured, 0, sizeof(sub_measured)); - sieve_measure_graph(sg, N_SUB, sub_marg, sub_measured, 1); + sieve_measure_graph(sg, N_SUB, sub_marg, sub_measured, 1); - for (int j = 0; j < N_SUB; j++) { - double best_prob = -1.0; - int best_v = 0; - for (int v = 0; v < 6; v++) { - if (sub_marg[j][v] > best_prob) { - best_prob = sub_marg[j][v]; - best_v = v; + for (int j = 0; j < N_SUB; j++) { + double best_prob = -1.0; + int best_v = 0; + for (int v = 0; v < 6; v++) { + if (sub_marg[j][v] > best_prob) { + best_prob = sub_marg[j][v]; + best_v = v; + } } + Ls_blk[j] = state_ls[j][best_v]; + Lm_blk[j] = state_lm[j][best_v]; } - Ls_blk[j] = state_ls[j][best_v]; - Lm_blk[j] = state_lm[j][best_v]; } - } - uint8_t L[QK_K]; - for (int j = 0; j < N_SUB; j++) { - float d_sub = dm * (float)Ls_blk[j]; - float m_sub = mm * (float)Lm_blk[j]; - if (d_sub < 1e-15f) { - for (int k = 0; k < 16; k++) L[16*j+k] = 0; - continue; + uint8_t L[QK_K]; + for (int j = 0; j < N_SUB; j++) { + float d_sub = dm * (float)Ls_blk[j]; + float m_sub = mm * (float)Lm_blk[j]; + if (d_sub < 1e-15f) { + for (int k = 0; k < 16; k++) L[16*j+k] = 0; + continue; + } + for (int k = 0; k < 16; k++) { + int q = gguf_nearest_int((adj_block_x[16*j+k] + m_sub) / d_sub); + if (q < 0) q = 0; if (q > 3) q = 3; + L[16*j+k] = (uint8_t)q; + } } - for (int k = 0; k < 16; k++) { - int q = gguf_nearest_int((adj_block_x[16*j+k] + m_sub) / d_sub); - if (q < 0) q = 0; if (q > 3) q = 3; - L[16*j+k] = (uint8_t)q; + + double Saa = 0, Sab = 0, Sbb = 0, Sxa = 0, Sxb = 0; + for (int j = 0; j < N_SUB; j++) { + float ls_f = (float)Ls_blk[j]; + float lm_f = (float)Lm_blk[j]; + for (int k = 0; k < 16; k++) { + float x = adj_block_x[16*j+k]; + float w = (imat_importance) ? + imat_importance[blk * QK_K + 16*j+k] : 1.0f; + float a = ls_f * (float)L[16*j+k]; + float b = lm_f; + Saa += w * a * a; + Sab += w * a * b; + Sbb += w * b * b; + Sxa += w * x * a; + Sxb += w * x * b; + } } - } - double Saa = 0, Sab = 0, Sbb = 0, Sxa = 0, Sxb = 0; - for (int j = 0; j < N_SUB; j++) { - float ls_f = (float)Ls_blk[j]; - float lm_f = (float)Lm_blk[j]; - for (int k = 0; k < 16; k++) { - float x = adj_block_x[16*j+k]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j+k] : 1.0f; - float a = ls_f * (float)L[16*j+k]; - float b = lm_f; - Saa += w * a * a; - Sab += w * a * b; - Sbb += w * b * b; - Sxa += w * x * a; - Sxb += w * x * b; + double det = Saa * Sbb - Sab * Sab; + if (fabs(det) > 1e-30) { + double d_new = (Sbb * Sxa - Sab * Sxb) / det; + double dm_new = (Sab * Sxa - Saa * Sxb) / det; + float d_seed = gguf_fp16_to_fp32(candidate_d[blk][cidx]); + float m_seed = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); + if (d_new > 0.0 && d_new < 4.0 * (d_seed + 1e-10)) + dm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)d_new)); + if (dm_new > 0.0 && dm_new < 4.0 * (m_seed + 1e-10)) + mm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)dm_new)); } - } - double det = Saa * Sbb - Sab * Sab; - if (fabs(det) > 1e-30) { - double d_new = (Sbb * Sxa - Sab * Sxb) / det; - double dm_new = (Sab * Sxa - Saa * Sxb) / det; - float d_seed = gguf_fp16_to_fp32(candidate_d[blk][cidx]); - float m_seed = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); - if (d_new > 0.0 && d_new < 4.0 * (d_seed + 1e-10)) - dm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)d_new)); - if (dm_new > 0.0 && dm_new < 4.0 * (m_seed + 1e-10)) - mm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)dm_new)); + uint16_t cur_dm16 = gguf_fp32_to_fp16(dm); + uint16_t cur_mm16 = gguf_fp32_to_fp16(mm); + if (cur_dm16 == prev_dm16 && cur_mm16 == prev_mm16) break; + prev_dm16 = cur_dm16; + prev_mm16 = cur_mm16; } - uint16_t cur_dm16 = gguf_fp32_to_fp16(dm); - uint16_t cur_mm16 = gguf_fp32_to_fp16(mm); - if (cur_dm16 == prev_dm16 && cur_mm16 == prev_mm16) break; - prev_dm16 = cur_dm16; - prev_mm16 = cur_mm16; - } - - /* ── FP16 ULP neighborhood search for (d, dmin) — Expanded to ±8 ── */ - { - uint16_t base_d16 = gguf_fp32_to_fp16(dm); - uint16_t base_m16 = gguf_fp32_to_fp16(mm); - uint16_t best_d16 = base_d16, best_m16 = base_m16; - float best_ulp_err = 1e30f; - - for (int dd = -8; dd <= 8; dd++) { - int cd16 = (int)base_d16 + dd; - if (cd16 < 0 || cd16 > 0x7BFF) continue; - float trial_dm = gguf_fp16_to_fp32((uint16_t)cd16); - - for (int dm_delta = -8; dm_delta <= 8; dm_delta++) { - int cm16 = (int)base_m16 + dm_delta; - if (cm16 < 0 || cm16 > 0x7BFF) continue; - float trial_mm = gguf_fp16_to_fp32((uint16_t)cm16); + /* ── FP16 ULP neighborhood search for (d, dmin) — Expanded to ±8 ── */ + { + uint16_t base_d16 = gguf_fp32_to_fp16(dm); + uint16_t base_m16 = gguf_fp32_to_fp16(mm); + uint16_t best_d16 = base_d16, best_m16 = base_m16; + float best_ulp_err = 1e30f; + + for (int dd = -8; dd <= 8; dd++) { + int cd16 = (int)base_d16 + dd; + if (cd16 < 0 || cd16 > 0x7BFF) continue; + float trial_dm = gguf_fp16_to_fp32((uint16_t)cd16); + + for (int dm_delta = -8; dm_delta <= 8; dm_delta++) { + int cm16 = (int)base_m16 + dm_delta; + if (cm16 < 0 || cm16 > 0x7BFF) continue; + float trial_mm = gguf_fp16_to_fp32((uint16_t)cm16); + + float err = 0.0f; + for (int j = 0; j < N_SUB; j++) { + float d_sub = trial_dm * (float)Ls_blk[j]; + float m_sub = trial_mm * (float)Lm_blk[j]; + for (int k = 0; k < 16; k++) { + float x = adj_block_x[16*j+k]; + float w = (imat_importance) ? + imat_importance[blk * QK_K + 16*j+k] : 1.0f; + int q; + if (d_sub < 1e-15f) { q = 0; } + else { + q = gguf_nearest_int((x + m_sub) / d_sub); + if (q < 0) q = 0; if (q > 3) q = 3; + } + float deq = d_sub * (float)q - m_sub; + float diff = x - deq; + err += diff * diff * w; + } + } + if (err < best_ulp_err) { + best_ulp_err = err; + best_d16 = (uint16_t)cd16; + best_m16 = (uint16_t)cm16; + } + } + } + dm = gguf_fp16_to_fp32(best_d16); + mm = gguf_fp16_to_fp32(best_m16); + } - float err = 0.0f; - for (int j = 0; j < N_SUB; j++) { - float d_sub = trial_dm * (float)Ls_blk[j]; - float m_sub = trial_mm * (float)Lm_blk[j]; + for (int j = 0; j < N_SUB; j++) { + const float *sx = adj_block_x + 16 * j; + float best_sub_err = 1e30f; + uint8_t best_ls = Ls_blk[j], best_lm = Lm_blk[j]; + for (int try_ls = 0; try_ls <= 15; try_ls++) { + float d_sub = dm * (float)try_ls; + for (int try_lm = 0; try_lm <= 15; try_lm++) { + float m_sub = mm * (float)try_lm; + float sub_err = 0.0f; for (int k = 0; k < 16; k++) { - float x = adj_block_x[16*j+k]; + float x = sx[k]; float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j+k] : 1.0f; + imat_importance[blk * QK_K + 16*j + k] : 1.0f; int q; if (d_sub < 1e-15f) { q = 0; } else { @@ -3753,730 +4303,491 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, } float deq = d_sub * (float)q - m_sub; float diff = x - deq; - err += diff * diff * w; + sub_err += diff * diff * w; } - } - if (err < best_ulp_err) { - best_ulp_err = err; - best_d16 = (uint16_t)cd16; - best_m16 = (uint16_t)cm16; - } - } - } - dm = gguf_fp16_to_fp32(best_d16); - mm = gguf_fp16_to_fp32(best_m16); - } - - for (int j = 0; j < N_SUB; j++) { - const float *sx = adj_block_x + 16 * j; - float best_sub_err = 1e30f; - uint8_t best_ls = Ls_blk[j], best_lm = Lm_blk[j]; - for (int try_ls = 0; try_ls <= 15; try_ls++) { - float d_sub = dm * (float)try_ls; - for (int try_lm = 0; try_lm <= 15; try_lm++) { - float m_sub = mm * (float)try_lm; - float sub_err = 0.0f; - for (int k = 0; k < 16; k++) { - float x = sx[k]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j + k] : 1.0f; - int q; - if (d_sub < 1e-15f) { q = 0; } - else { - q = gguf_nearest_int((x + m_sub) / d_sub); - if (q < 0) q = 0; if (q > 3) q = 3; + if (sub_err < best_sub_err) { + best_sub_err = sub_err; + best_ls = (uint8_t)try_ls; + best_lm = (uint8_t)try_lm; } - float deq = d_sub * (float)q - m_sub; - float diff = x - deq; - sub_err += diff * diff * w; - } - if (sub_err < best_sub_err) { - best_sub_err = sub_err; - best_ls = (uint8_t)try_ls; - best_lm = (uint8_t)try_lm; } } + Ls_blk[j] = best_ls; + Lm_blk[j] = best_lm; } - Ls_blk[j] = best_ls; - Lm_blk[j] = best_lm; - } - - output[blk].d = gguf_fp32_to_fp16(dm); - output[blk].dmin = gguf_fp32_to_fp16(mm); - - for (int j = 0; j < N_SUB; j++) - output[blk].scales[j] = Ls_blk[j] | (Lm_blk[j] << 4); - /* ── Final quantization: D₆ Hadamard Greedy Descent (deterministic) ── - * - * The original Simulated Annealing acceptance rule is replaced by a - * strict greedy descent: only accept a flip if it strictly reduces the - * D₆ Hadamard metric (4·‖vesica‖² + DC²). This makes error shaping - * fully deterministic and thread-safe (no rand() inside omp parallel), - * consistent with the Viterbi philosophy applied in Phase 3. - * - * The metric measures both: - * - Vesica Piscis term: correlated error between weights i and i+QK_K/2 - * (targets the first non-DC harmonic — halfwave symmetry) - * - DC term: total signed error across the 256-weight superblock - * (captured and propagated to the next block by Phase 3.9) - */ - uint8_t L[QK_K]; - { - float q_cont_all[QK_K]; - int q_base_all[QK_K]; - int q_shaped_all[QK_K]; - - for (int i = 0; i < QK_K; i++) { - int jj = i >> 4; - float d_s = dm * (float)(output[blk].scales[jj] & 0xF); - float m_s = mm * (float)(output[blk].scales[jj] >> 4); - if (d_s < 1e-15f) { - q_cont_all[i] = 0.0f; - q_base_all[i] = 0; - } else { - /* Quantize the DC-adjusted target */ - float qc = (adj_block_x[i] + m_s) / d_s; - q_cont_all[i] = qc; - int qr = gguf_nearest_int(qc); - if (qr < 0) qr = 0; if (qr > 3) qr = 3; - q_base_all[i] = qr; - } - } - memcpy(q_shaped_all, q_base_all, QK_K * sizeof(int)); - - float e_live[QK_K]; - for (int i = 0; i < QK_K; i++) { - int jj = i >> 4; - float d_s = dm * (float)(output[blk].scales[jj] & 0xF); - float m_s = mm * (float)(output[blk].scales[jj] >> 4); - /* Decoder semantics: deq = d_s·q − m_s, which is −m_s when - * d_s == 0 (NOT 0 — the −dmin·lm term always applies). */ - float deq = d_s * (float)q_shaped_all[i] - m_s; - /* Residual against the adjusted target (DC-corrected view) */ - e_live[i] = adj_block_x[i] - deq; - } - - float v_live[QK_K / 2]; - float vesica_cur = 0.0f, dc_cur = 0.0f; - for (int i = 0; i < QK_K / 2; i++) { - v_live[i] = e_live[i] + e_live[i + QK_K / 2]; - vesica_cur += v_live[i] * v_live[i]; - } - for (int i = 0; i < QK_K; i++) dc_cur += e_live[i]; - float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; - - /* Deterministic greedy descent: accept only strict improvements */ - for (int pass = 0; pass < QK_K; pass++) { - int best_k = -1; - int best_q_alt = 0; - float best_delta = 0.0f; /* strictly positive threshold */ + output[blk].d = gguf_fp32_to_fp16(dm); + output[blk].dmin = gguf_fp32_to_fp16(mm); - for (int k = 0; k < QK_K; k++) { - int jj = k >> 4; - float d_s = dm * (float)(output[blk].scales[jj] & 0xF); - if (d_s < 1e-15f) continue; + for (int j = 0; j < N_SUB; j++) + output[blk].scales[j] = Ls_blk[j] | (Lm_blk[j] << 4); - int q_cur = q_shaped_all[k]; - int q_try = (q_cont_all[k] - (float)q_cur >= 0.0f) - ? q_cur + 1 : q_cur - 1; - if (q_try < 0 || q_try > 3) continue; + /* ── Final quantization: D₆ Hadamard Greedy Descent (deterministic) ── + * + * The original Simulated Annealing acceptance rule is replaced by a + * strict greedy descent: only accept a flip if it strictly reduces the + * D₆ Hadamard metric (4·‖vesica‖² + DC²). This makes error shaping + * fully deterministic and thread-safe (no rand() inside omp parallel), + * consistent with the Viterbi philosophy applied in Phase 3. + * + * The metric measures both: + * - Vesica Piscis term: correlated error between weights i and i+QK_K/2 + * (targets the first non-DC harmonic — halfwave symmetry) + * - DC term: total signed error across the 256-weight superblock + * (captured and propagated to the next block by Phase 3.9) + */ + uint8_t L[QK_K]; + { + float q_cont_all[QK_K]; + int q_base_all[QK_K]; + int q_shaped_all[QK_K]; - float m_s = mm * (float)(output[blk].scales[jj] >> 4); - float e_new = adj_block_x[k] - (d_s * (float)q_try - m_s); - float de = e_new - e_live[k]; + for (int i = 0; i < QK_K; i++) { + int jj = i >> 4; + float d_s = dm * (float)(output[blk].scales[jj] & 0xF); + float m_s = mm * (float)(output[blk].scales[jj] >> 4); + if (d_s < 1e-15f) { + q_cont_all[i] = 0.0f; + q_base_all[i] = 0; + } else { + /* Quantize the DC-adjusted target */ + float qc = (adj_block_x[i] + m_s) / d_s; + q_cont_all[i] = qc; + int qr = gguf_nearest_int(qc); + if (qr < 0) qr = 0; if (qr > 3) qr = 3; + q_base_all[i] = qr; + } + } + memcpy(q_shaped_all, q_base_all, QK_K * sizeof(int)); + + /* The descent minimises the PURELY SPECTRAL metric + * 4*||vesica||^2 + DC^2, and its result is accepted only if the + * extended objective E improves. At HEX_DC_LAMBDA = + * HEX_VW_LAMBDA = 0, E is plain weighted SSE, for which + * q_base_all (nearest rounding at fixed d, dmin, ls, lm) is + * already per-weight optimal — so no shaped state can ever win + * and the whole pass is provably dead weight. Skip it and spend + * the time on the Phase-4.0 restarts instead. Compiled out + * entirely at lambda = 0; fully restored when either lambda is + * rebuilt non-zero. */ + if (HEX_DC_LAMBDA != 0.0f || HEX_VW_LAMBDA != 0.0f) { + + float e_live[QK_K]; + for (int i = 0; i < QK_K; i++) { + int jj = i >> 4; + float d_s = dm * (float)(output[blk].scales[jj] & 0xF); + float m_s = mm * (float)(output[blk].scales[jj] >> 4); + /* Decoder semantics: deq = d_s·q − m_s, which is −m_s when + * d_s == 0 (NOT 0 — the −dmin·lm term always applies). */ + float deq = d_s * (float)q_shaped_all[i] - m_s; + /* Residual against the adjusted target (DC-corrected view) */ + e_live[i] = adj_block_x[i] - deq; + } - int pi = (k < QK_K / 2) ? k : k - QK_K / 2; - float v_new = v_live[pi] + de; + float v_live[QK_K / 2]; + float vesica_cur = 0.0f, dc_cur = 0.0f; + for (int i = 0; i < QK_K / 2; i++) { + v_live[i] = e_live[i] + e_live[i + QK_K / 2]; + vesica_cur += v_live[i] * v_live[i]; + } + for (int i = 0; i < QK_K; i++) dc_cur += e_live[i]; + float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; + + /* Deterministic greedy descent: accept only strict improvements */ + for (int pass = 0; pass < QK_K; pass++) { + int best_k = -1; + int best_q_alt = 0; + float best_delta = 0.0f; /* strictly positive threshold */ + + for (int k = 0; k < QK_K; k++) { + int jj = k >> 4; + float d_s = dm * (float)(output[blk].scales[jj] & 0xF); + if (d_s < 1e-15f) continue; + + int q_cur = q_shaped_all[k]; + int q_try = (q_cont_all[k] - (float)q_cur >= 0.0f) + ? q_cur + 1 : q_cur - 1; + if (q_try < 0 || q_try > 3) continue; + + float m_s = mm * (float)(output[blk].scales[jj] >> 4); + float e_new = adj_block_x[k] - (d_s * (float)q_try - m_s); + float de = e_new - e_live[k]; + + int pi = (k < QK_K / 2) ? k : k - QK_K / 2; + float v_new = v_live[pi] + de; + + float vesica_alt = vesica_cur - v_live[pi]*v_live[pi] + v_new*v_new; + float dc_alt = dc_cur + de; + float delta = metric_cur - (4.0f * vesica_alt + dc_alt * dc_alt); + + if (delta > best_delta) { + best_delta = delta; + best_k = k; + best_q_alt = q_try; + } + } - float vesica_alt = vesica_cur - v_live[pi]*v_live[pi] + v_new*v_new; - float dc_alt = dc_cur + de; - float delta = metric_cur - (4.0f * vesica_alt + dc_alt * dc_alt); + if (best_k < 0) break; /* converged — no further improvement */ - if (delta > best_delta) { - best_delta = delta; - best_k = k; - best_q_alt = q_try; + q_shaped_all[best_k] = best_q_alt; + { + int jj_c = best_k >> 4; + float d_c = dm * (float)(output[blk].scales[jj_c] & 0xF); + float m_c = mm * (float)(output[blk].scales[jj_c] >> 4); + float e_new_c = adj_block_x[best_k] - (d_c * (float)best_q_alt - m_c); + float de_c = e_new_c - e_live[best_k]; + int pi_c = (best_k < QK_K / 2) ? best_k : best_k - QK_K / 2; + float v_new_c = v_live[pi_c] + de_c; + vesica_cur += v_new_c * v_new_c - v_live[pi_c] * v_live[pi_c]; + dc_cur += de_c; + metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; + v_live[pi_c] = v_new_c; + e_live[best_k]= e_new_c; } } - if (best_k < 0) break; /* converged — no further improvement */ - - q_shaped_all[best_k] = best_q_alt; + /* Choose base vs shaped on the EXTENDED objective vs originals */ + float err_base = 0.0f, err_shaped = 0.0f; + float e_qb[QK_K], e_qs[QK_K]; + for (int i = 0; i < QK_K; i++) { + int jj = i >> 4; + float d_s = dm * (float)(output[blk].scales[jj] & 0xF); + float m_s = mm * (float)(output[blk].scales[jj] >> 4); + float w = (imat_importance) ? imat_importance[blk * QK_K + i] : 1.0f; + float deq_b = d_s * (float)q_base_all[i] - m_s; /* −m_s when d_s==0 */ + float deq_s = d_s * (float)q_shaped_all[i] - m_s; + float xv = block_x[i]; /* original weight for error report */ + e_qb[i] = xv - deq_b; + e_qs[i] = xv - deq_s; + err_base += e_qb[i] * e_qb[i] * w; + err_shaped += e_qs[i] * e_qs[i] * w; + } + err_base += hex_spectral_penalty(e_qb, QK_K); + err_shaped += hex_spectral_penalty(e_qs, QK_K); { - int jj_c = best_k >> 4; - float d_c = dm * (float)(output[blk].scales[jj_c] & 0xF); - float m_c = mm * (float)(output[blk].scales[jj_c] >> 4); - float e_new_c = adj_block_x[best_k] - (d_c * (float)best_q_alt - m_c); - float de_c = e_new_c - e_live[best_k]; - int pi_c = (best_k < QK_K / 2) ? best_k : best_k - QK_K / 2; - float v_new_c = v_live[pi_c] + de_c; - vesica_cur += v_new_c * v_new_c - v_live[pi_c] * v_live[pi_c]; - dc_cur += de_c; - metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur; - v_live[pi_c] = v_new_c; - e_live[best_k]= e_new_c; + int use_shaped = (err_shaped <= err_base); + for (int i = 0; i < QK_K; i++) + L[i] = (uint8_t)(use_shaped ? q_shaped_all[i] : q_base_all[i]); } - } - /* Choose base vs shaped on the EXTENDED objective vs originals */ - float err_base = 0.0f, err_shaped = 0.0f; - float e_qb[QK_K], e_qs[QK_K]; - for (int i = 0; i < QK_K; i++) { - int jj = i >> 4; - float d_s = dm * (float)(output[blk].scales[jj] & 0xF); - float m_s = mm * (float)(output[blk].scales[jj] >> 4); - float w = (imat_importance) ? imat_importance[blk * QK_K + i] : 1.0f; - float deq_b = d_s * (float)q_base_all[i] - m_s; /* −m_s when d_s==0 */ - float deq_s = d_s * (float)q_shaped_all[i] - m_s; - float xv = block_x[i]; /* original weight for error report */ - e_qb[i] = xv - deq_b; - e_qs[i] = xv - deq_s; - err_base += e_qb[i] * e_qb[i] * w; - err_shaped += e_qs[i] * e_qs[i] * w; - } - err_base += hex_spectral_penalty(e_qb, QK_K); - err_shaped += hex_spectral_penalty(e_qs, QK_K); - { - int use_shaped = (err_shaped <= err_base); + } else { + /* lambda = 0 path: commit nearest rounding directly. */ for (int i = 0; i < QK_K; i++) - L[i] = (uint8_t)(use_shaped ? q_shaped_all[i] : q_base_all[i]); + L[i] = (uint8_t)q_base_all[i]; + } } - } - /* ── Cross-weight error diffusion — intra-sub-block Floyd-Steinberg ── - * - * Implements cross-weight error diffusion within each 16-weight sub-block. - * After the greedy descent has committed quantisation codes, the residual - * of each weight is partially propagated forward to the next position in - * the same sub-block (7/16 of the error), re-quantising if the diffused - * target falls in a different bin. - * - * This is the "cross-weight" dimension of the error-diffusion request: - * neighbouring weights share and partially absorb each other's rounding - * error, shaping the within-block spectrum away from the DC component - * that Phase 3.9 already propagates between blocks. - * - * Staying within sub-blocks avoids scale-mismatch artefacts that would - * arise from diffusing across the dm * Ls[j] boundary between sub-blocks. - * - * The diffused codes are accepted only when they reduce the weighted MSE - * against the ORIGINAL weight (not the adjusted target), so the diffusion - * cannot increase the total reconstruction error. - */ - { - int fs_j, fs_k; - for (fs_j = 0; fs_j < N_SUB; fs_j++) { - int base = fs_j * 16; - float d_s = dm * (float)(output[blk].scales[fs_j] & 0xF); - float m_s = mm * (float)(output[blk].scales[fs_j] >> 4); - if (d_s < 1e-15f) continue; - - float carry = 0.0f; /* FS carry from position k-1 */ - - for (fs_k = 0; fs_k < 16; fs_k++) { - int idx = base + fs_k; - float x_orig = block_x[idx]; - float x_adj = adj_block_x[idx] + carry; /* adjusted + diffused */ - - /* Propose new code from diffused target */ - int q_fs = gguf_nearest_int((x_adj + m_s) / d_s); - if (q_fs < 0) q_fs = 0; if (q_fs > 3) q_fs = 3; - - if (q_fs != (int)L[idx]) { - /* Accept only when MSE against original weight improves */ - float w_imp = (imat_importance) - ? imat_importance[blk * QK_K + idx] : 1.0f; - float deq_old = d_s * (float)L[idx] - m_s; - float deq_new = d_s * (float)q_fs - m_s; - float e_old = (x_orig - deq_old) * (x_orig - deq_old) * w_imp; - float e_new = (x_orig - deq_new) * (x_orig - deq_new) * w_imp; - if (e_new < e_old) - L[idx] = (uint8_t)q_fs; - } + /* ── Cross-weight error diffusion — intra-sub-block Floyd-Steinberg ── + * + * Implements cross-weight error diffusion within each 16-weight sub-block. + * After the greedy descent has committed quantisation codes, the residual + * of each weight is partially propagated forward to the next position in + * the same sub-block (7/16 of the error), re-quantising if the diffused + * target falls in a different bin. + * + * This is the "cross-weight" dimension of the error-diffusion request: + * neighbouring weights share and partially absorb each other's rounding + * error, shaping the within-block spectrum away from the DC component + * that Phase 3.9 already propagates between blocks. + * + * Staying within sub-blocks avoids scale-mismatch artefacts that would + * arise from diffusing across the dm * Ls[j] boundary between sub-blocks. + * + * The diffused codes are accepted only when they reduce the weighted MSE + * against the ORIGINAL weight (not the adjusted target), so the diffusion + * cannot increase the total reconstruction error. + * + * At DC_DECAY = 0 AND lambda = 0 this pass is provably a no-op: L + * already holds nearest-of-original rounding, which minimises the + * per-weight squared error the accept guard tests, so no proposal + * can ever be accepted. Skip it under those settings; it stays + * fully live whenever either knob is rebuilt non-zero. + */ + if (DC_DECAY != 0.0f || HEX_DC_LAMBDA != 0.0f || + HEX_VW_LAMBDA != 0.0f) { + int fs_j, fs_k; + for (fs_j = 0; fs_j < N_SUB; fs_j++) { + int base = fs_j * 16; + float d_s = dm * (float)(output[blk].scales[fs_j] & 0xF); + float m_s = mm * (float)(output[blk].scales[fs_j] >> 4); + if (d_s < 1e-15f) continue; - /* Propagate 7/16 of the residual (adj target vs committed code) */ - { - float deq_final = d_s * (float)L[idx] - m_s; - float residual = (adj_block_x[idx] - deq_final); - carry = (fs_k < 15) ? residual * (7.0f / 16.0f) : 0.0f; - } - } - } - } + float carry = 0.0f; /* FS carry from position k-1 */ + + for (fs_k = 0; fs_k < 16; fs_k++) { + int idx = base + fs_k; + float x_orig = block_x[idx]; + float x_adj = adj_block_x[idx] + carry; /* adjusted + diffused */ + + /* Propose new code from diffused target */ + int q_fs = gguf_nearest_int((x_adj + m_s) / d_s); + if (q_fs < 0) q_fs = 0; if (q_fs > 3) q_fs = 3; + + if (q_fs != (int)L[idx]) { + /* Accept only when MSE against original weight improves */ + float w_imp = (imat_importance) + ? imat_importance[blk * QK_K + idx] : 1.0f; + float deq_old = d_s * (float)L[idx] - m_s; + float deq_new = d_s * (float)q_fs - m_s; + float e_old = (x_orig - deq_old) * (x_orig - deq_old) * w_imp; + float e_new = (x_orig - deq_new) * (x_orig - deq_new) * w_imp; + if (e_new < e_old) + L[idx] = (uint8_t)q_fs; + } - /* ── Final closed-form (d, dmin) refit against the UNCLIPPED weights ── - * (issues #2 / #5) - * - * Every earlier (d, dmin) solve fits the DC-adjusted, soft-clipped - * target and runs BEFORE the greedy descent and Floyd-Steinberg passes - * mutate the committed 2-bit codes. Once L[], and the 4-bit sub-block - * scale codes (Ls = scales & 0xF, Lm = scales >> 4), are final, the two - * fp16 scalars (d, dmin) that minimise the importance-weighted SSE - * against the ORIGINAL weights have a closed form. Solve it and adopt it - * only when it lowers the weighted block error — so it can never raise - * RMSE, and because the integer codes are held fixed, the vesica/wave/DC - * error shaping baked into them is preserved intact. */ - { - double rSaa = 0, rSab = 0, rSbb = 0, rSxa = 0, rSxb = 0; - double rA = 0, rB = 0, rS = 0; /* DC rank-1 augmentation */ - for (int j = 0; j < N_SUB; j++) { - float ls_f = (float)(output[blk].scales[j] & 0xF); - float lm_f = (float)(output[blk].scales[j] >> 4); - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float x = block_x[idx]; /* unclipped original */ - float w = (imat_importance) ? imat_importance[blk * QK_K + idx] : 1.0f; - float a = ls_f * (float)L[idx]; - float b = lm_f; - rSaa += (double)w * a * a; - rSab += (double)w * a * b; - rSbb += (double)w * b * b; - rSxa += (double)w * x * a; - rSxb += (double)w * x * b; - rA += a; rB += b; rS += x; - } - } - /* DC term as one augmented observation (S ~ A·d − B·m), weight - * λ_dc/n; vesica/wave handled by the extended-E acceptance. */ - { - double rw = (double)HEX_DC_LAMBDA / (double)QK_K; - rSaa += rw * rA * rA; rSab += rw * rA * rB; - rSbb += rw * rB * rB; rSxa += rw * rS * rA; - rSxb += rw * rS * rB; - } - double rdet = rSaa * rSbb - rSab * rSab; - if (fabs(rdet) > 1e-30) { - double d_ref = (rSbb * rSxa - rSab * rSxb) / rdet; - double m_ref = (rSab * rSxa - rSaa * rSxb) / rdet; - if (d_ref > 0.0) { - float dm_try = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)d_ref)); - float mm_try = (m_ref > 0.0) - ? gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)m_ref)) - : mm; - /* Extended-objective acceptance test vs original weights. */ - float err_cur = 0.0f, err_try = 0.0f; - float e_rc[QK_K], e_rt[QK_K]; - for (int j = 0; j < N_SUB; j++) { - float ls_f = (float)(output[blk].scales[j] & 0xF); - float lm_f = (float)(output[blk].scales[j] >> 4); - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float x = block_x[idx]; - float w = (imat_importance) ? imat_importance[blk * QK_K + idx] : 1.0f; - float qf = (float)L[idx]; - float dc = dm * ls_f * qf - mm * lm_f; - float dt = dm_try * ls_f * qf - mm_try * lm_f; - e_rc[idx] = x - dc; - e_rt[idx] = x - dt; - err_cur += e_rc[idx] * e_rc[idx] * w; - err_try += e_rt[idx] * e_rt[idx] * w; + /* Propagate 7/16 of the residual (adj target vs committed code) */ + { + float deq_final = d_s * (float)L[idx] - m_s; + float residual = (adj_block_x[idx] - deq_final); + carry = (fs_k < 15) ? residual * (7.0f / 16.0f) : 0.0f; } } - err_cur += hex_spectral_penalty(e_rc, QK_K); - err_try += hex_spectral_penalty(e_rt, QK_K); - if (err_try < err_cur) { dm = dm_try; mm = mm_try; } } } - output[blk].d = gguf_fp32_to_fp16(dm); - output[blk].dmin = gguf_fp32_to_fp16(mm); - } - - /* ══ PHASE 4.6: MONOTONE COORDINATE-DESCENT POLISH (RMSE-guaranteed) ══ - * - * Objective-function mismatch fix: the final passes that commit the - * 2-bit codes — the 16×16 (ls, lm) sub-block search, the ±8 ULP - * (d, dmin) neighborhood search, and the greedy-descent error shaping - * — all minimise error against the DC-ADJUSTED target adj_block_x. - * The reported RMSE, however, is measured against the ORIGINAL - * weights. The codes are therefore stranded at the optimum of a - * SHIFTED objective, while only the scalar (d, dmin) refit above - * targets the true one (and it holds all codes frozen). - * - * This polish runs alternating coordinate descent on the TRUE - * objective (importance-weighted SSE vs the original weights): - * - * (1) For each 16-weight sub-block, an exact joint re-search of - * (ls, lm) over the full 16×16 grid with per-weight optimal - * q ∈ {0..3}, committed only on strict improvement of the - * extended objective E. With λ_dc = λ_vw = 0 sub-blocks are - * independent given (d, dmin); with spectral terms active the - * coupling (DC: all subs; fold: sub j ↔ sub j⊕8) is handled - * exactly via live residual bookkeeping. - * (2) Closed-form weighted LS refit of the two fp16 scalars - * (d, dmin) with all codes held fixed, committed only on - * strict improvement (same guard as the refit above). - * - * All moves are accept-only-if-better on E ⇒ the extended block - * objective is monotonically non-increasing; at λ = 0 this reduces - * to RMSE-monotone (final RMSE can only go DOWN relative to the - * unpatched pipeline), at λ > 0 small RMSE giveback is permitted - * exactly where it buys dot-product error cancellation. The state space is finite - * (4-bit codes, fp16 scalars), so the loop terminates; in practice - * it converges in 2–3 sweeps. The vesica/DC spectral shaping baked - * into L survives wherever it is SSE-neutral, and is overridden - * only where it was costing true reconstruction error. */ - { - uint8_t pl_Ls[16], pl_Lm[16]; - for (int j = 0; j < N_SUB; j++) { - pl_Ls[j] = output[blk].scales[j] & 0xF; - pl_Lm[j] = output[blk].scales[j] >> 4; - } - for (int pol_iter = 0; pol_iter < 6; pol_iter++) { - int pol_improved = 0; - - /* ── (1) Exact per-sub-block (ls, lm, q) re-search on the - * EXTENDED objective. Under the spectral terms sub-blocks - * are no longer independent: every sub couples to all others - * through the DC term and to its fold partner (sub j ⊕ 8, - * i.e. weights i ↔ i+128) through vesica² − wave². The - * search therefore keeps live residuals pe[] and scores each - * candidate against the whole-block penalty with the partner - * residuals held fixed — exact coordinate descent on E. */ - float pe[QK_K]; - float sub_sse[16], sub_dc[16], pair_cross[8]; - float dc_tot = 0.0f, cross_tot = 0.0f; + /* ── Final closed-form (d, dmin) refit against the UNCLIPPED weights ── + * (issues #2 / #5) + * + * Every earlier (d, dmin) solve fits the DC-adjusted, soft-clipped + * target and runs BEFORE the greedy descent and Floyd-Steinberg passes + * mutate the committed 2-bit codes. Once L[], and the 4-bit sub-block + * scale codes (Ls = scales & 0xF, Lm = scales >> 4), are final, the two + * fp16 scalars (d, dmin) that minimise the importance-weighted SSE + * against the ORIGINAL weights have a closed form. Solve it and adopt it + * only when it lowers the weighted block error — so it can never raise + * RMSE, and because the integer codes are held fixed, the vesica/wave/DC + * error shaping baked into them is preserved intact. */ + { + double rSaa = 0, rSab = 0, rSbb = 0, rSxa = 0, rSxb = 0; + double rA = 0, rB = 0, rS = 0; /* DC rank-1 augmentation */ for (int j = 0; j < N_SUB; j++) { - float d_sub = dm * (float)pl_Ls[j]; - float m_sub = mm * (float)pl_Lm[j]; - sub_sse[j] = 0.0f; - sub_dc[j] = 0.0f; + float ls_f = (float)(output[blk].scales[j] & 0xF); + float lm_f = (float)(output[blk].scales[j] >> 4); for (int k = 0; k < 16; k++) { int idx = 16 * j + k; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - /* deq = d·ls·q − dmin·lm; equals −m_sub at ls==0 */ - float e = block_x[idx] - (d_sub * (float)L[idx] - m_sub); - pe[idx] = e; - sub_sse[j] += e * e * w; - sub_dc[j] += e; + float x = block_x[idx]; /* unclipped original */ + float w = (imat_importance) ? imat_importance[blk * QK_K + idx] : 1.0f; + float a = ls_f * (float)L[idx]; + float b = lm_f; + rSaa += (double)w * a * a; + rSab += (double)w * a * b; + rSbb += (double)w * b * b; + rSxa += (double)w * x * a; + rSxb += (double)w * x * b; + rA += a; rB += b; rS += x; } - dc_tot += sub_dc[j]; } - for (int p = 0; p < 8; p++) { - pair_cross[p] = 0.0f; - for (int k = 0; k < 16; k++) - pair_cross[p] += pe[16*p + k] * pe[16*(p+8) + k]; - cross_tot += pair_cross[p]; + /* DC term as one augmented observation (S ~ A·d − B·m), weight + * λ_dc/n; vesica/wave handled by the extended-E acceptance. */ + { + double rw = (double)HEX_DC_LAMBDA / (double)QK_K; + rSaa += rw * rA * rA; rSab += rw * rA * rB; + rSbb += rw * rB * rB; rSxa += rw * rS * rA; + rSxb += rw * rS * rB; } - - for (int j = 0; j < N_SUB; j++) { - const float *sx = block_x + 16 * j; - int pi = j & 7; /* fold-pair index */ - int pj = j ^ 8; /* partner sub-block */ - const float *ppe = pe + 16 * pj; /* partner residuals */ - float dc_rest = dc_tot - sub_dc[j]; - float cross_rest = cross_tot - pair_cross[pi]; - - /* Extended score of the CURRENT committed state */ - float best_sub = sub_sse[j] - + (HEX_DC_LAMBDA / (float)QK_K) * dc_tot * dc_tot - + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f * cross_tot; - int best_ls = -1, best_lm = 0; - uint8_t best_q[16]; - float best_e[16]; - float best_sse = 0.0f, best_dcc = 0.0f, best_cxc = 0.0f; - - for (int try_ls = 0; try_ls <= 15; try_ls++) { - float d_sub = dm * (float)try_ls; - for (int try_lm = 0; try_lm <= 15; try_lm++) { - float m_sub = mm * (float)try_lm; - float sub_err = 0.0f, dcc = 0.0f, cxc = 0.0f; - uint8_t q_loc[16]; - float e_loc[16]; - int aborted = 0; + double rdet = rSaa * rSbb - rSab * rSab; + if (fabs(rdet) > 1e-30) { + double d_ref = (rSbb * rSxa - rSab * rSxb) / rdet; + double m_ref = (rSab * rSxa - rSaa * rSxb) / rdet; + if (d_ref > 0.0) { + float dm_try = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)d_ref)); + float mm_try = (m_ref > 0.0) + ? gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)m_ref)) + : mm; + /* Extended-objective acceptance test vs original weights. */ + float err_cur = 0.0f, err_try = 0.0f; + float e_rc[QK_K], e_rt[QK_K]; + for (int j = 0; j < N_SUB; j++) { + float ls_f = (float)(output[blk].scales[j] & 0xF); + float lm_f = (float)(output[blk].scales[j] >> 4); for (int k = 0; k < 16; k++) { - float x = sx[k]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + 16*j + k] : 1.0f; - int q = 0; - if (d_sub >= 1e-15f) { - q = gguf_nearest_int((x + m_sub) / d_sub); - if (q < 0) q = 0; if (q > 3) q = 3; - } - q_loc[k] = (uint8_t)q; - /* deq = d·ls·q − dmin·lm; −m_sub at ls==0 */ - float e = x - (d_sub * (float)q - m_sub); - e_loc[k] = e; - sub_err += e * e * w; - dcc += e; - cxc += e * ppe[k]; - /* SSE-partial prune is a valid lower bound - * only while the spectral terms are ≥ 0, - * i.e. when the (signable) vw credit is off */ - if (HEX_VW_LAMBDA == 0.0f && - sub_err >= best_sub) { aborted = 1; break; } - } - if (aborted) continue; - float score = sub_err - + (HEX_DC_LAMBDA / (float)QK_K) - * (dc_rest + dcc) * (dc_rest + dcc) - + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f - * (cross_rest + cxc); - if (score < best_sub) { - best_sub = score; - best_ls = try_ls; - best_lm = try_lm; - memcpy(best_q, q_loc, 16); - memcpy(best_e, e_loc, sizeof(e_loc)); - best_sse = sub_err; - best_dcc = dcc; - best_cxc = cxc; + int idx = 16 * j + k; + float x = block_x[idx]; + float w = (imat_importance) ? imat_importance[blk * QK_K + idx] : 1.0f; + float qf = (float)L[idx]; + float dc = dm * ls_f * qf - mm * lm_f; + float dt = dm_try * ls_f * qf - mm_try * lm_f; + e_rc[idx] = x - dc; + e_rt[idx] = x - dt; + err_cur += e_rc[idx] * e_rc[idx] * w; + err_try += e_rt[idx] * e_rt[idx] * w; } } + err_cur += hex_spectral_penalty(e_rc, QK_K); + err_try += hex_spectral_penalty(e_rt, QK_K); + if (err_try < err_cur) { dm = dm_try; mm = mm_try; } } + } + output[blk].d = gguf_fp32_to_fp16(dm); + output[blk].dmin = gguf_fp32_to_fp16(mm); + } - if (best_ls >= 0) { /* strict improvement in E found */ - pl_Ls[j] = (uint8_t)best_ls; - pl_Lm[j] = (uint8_t)best_lm; - memcpy(L + 16 * j, best_q, 16); - memcpy(pe + 16 * j, best_e, sizeof(best_e)); - sub_sse[j] = best_sse; - sub_dc[j] = best_dcc; - pair_cross[pi] = best_cxc; - dc_tot = dc_rest + best_dcc; - cross_tot = cross_rest + best_cxc; - pol_improved = 1; - } + /* ══ PHASE 4.6: MONOTONE COORDINATE-DESCENT POLISH (RMSE-guaranteed) ══ + * + * Objective-function mismatch fix: the final passes that commit the + * 2-bit codes — the 16×16 (ls, lm) sub-block search, the ±8 ULP + * (d, dmin) neighborhood search, and the greedy-descent error shaping + * — all minimise error against the DC-ADJUSTED target adj_block_x. + * The reported RMSE, however, is measured against the ORIGINAL + * weights. The codes are therefore stranded at the optimum of a + * SHIFTED objective, while only the scalar (d, dmin) refit above + * targets the true one (and it holds all codes frozen). + * + * This polish runs alternating coordinate descent on the TRUE + * objective (importance-weighted SSE vs the original weights): + * + * (1) For each 16-weight sub-block, an exact joint re-search of + * (ls, lm) over the full 16×16 grid with per-weight optimal + * q ∈ {0..3}, committed only on strict improvement of the + * extended objective E. With λ_dc = λ_vw = 0 sub-blocks are + * independent given (d, dmin); with spectral terms active the + * coupling (DC: all subs; fold: sub j ↔ sub j⊕8) is handled + * exactly via live residual bookkeeping. + * (2) Closed-form weighted LS refit of the two fp16 scalars + * (d, dmin) with all codes held fixed, committed only on + * strict improvement (same guard as the refit above). + * + * All moves are accept-only-if-better on E ⇒ the extended block + * objective is monotonically non-increasing; at λ = 0 this reduces + * to RMSE-monotone (final RMSE can only go DOWN relative to the + * unpatched pipeline), at λ > 0 small RMSE giveback is permitted + * exactly where it buys dot-product error cancellation. The state space is finite + * (4-bit codes, fp16 scalars), so the loop terminates; in practice + * it converges in 2–3 sweeps. The vesica/DC spectral shaping baked + * into L survives wherever it is SSE-neutral, and is overridden + * only where it was costing true reconstruction error. */ + { + uint8_t pl_Ls[16], pl_Lm[16]; + for (int j = 0; j < N_SUB; j++) { + pl_Ls[j] = output[blk].scales[j] & 0xF; + pl_Lm[j] = output[blk].scales[j] >> 4; } - /* ── (2) Closed-form (d, dmin) refit vs ORIGINAL, codes fixed ── */ - { - double pSaa = 0, pSab = 0, pSbb = 0, pSxa = 0, pSxb = 0; - double pA = 0, pB = 0, pS = 0; /* DC rank-1 augmentation */ - for (int j = 0; j < N_SUB; j++) { - float ls_f = (float)pl_Ls[j]; - float lm_f = (float)pl_Lm[j]; - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float x = block_x[idx]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - float a = ls_f * (float)L[idx]; - float b = lm_f; - pSaa += (double)w * a * a; - pSab += (double)w * a * b; - pSbb += (double)w * b * b; - pSxa += (double)w * x * a; - pSxb += (double)w * x * b; - pA += a; pB += b; pS += x; - } - } - { - double pw = (double)HEX_DC_LAMBDA / (double)QK_K; - pSaa += pw * pA * pA; pSab += pw * pA * pB; - pSbb += pw * pB * pB; pSxa += pw * pS * pA; - pSxb += pw * pS * pB; - } - double pdet = pSaa * pSbb - pSab * pSab; - if (fabs(pdet) > 1e-30) { - double d_ref = (pSbb * pSxa - pSab * pSxb) / pdet; - double m_ref = (pSab * pSxa - pSaa * pSxb) / pdet; - if (d_ref > 0.0) { - float dm_try = gguf_fp16_to_fp32( - gguf_fp32_to_fp16((float)d_ref)); - float mm_try = (m_ref > 0.0) - ? gguf_fp16_to_fp32( - gguf_fp32_to_fp16((float)m_ref)) - : mm; - float err_cur = 0.0f, err_try = 0.0f; - float e_pc[QK_K], e_pt[QK_K]; - for (int j = 0; j < N_SUB; j++) { - float ls_f = (float)pl_Ls[j]; - float lm_f = (float)pl_Lm[j]; - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float x = block_x[idx]; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - float qf = (float)L[idx]; - float dc = dm * ls_f * qf - mm * lm_f; - float dt = dm_try * ls_f * qf - mm_try * lm_f; - e_pc[idx] = x - dc; - e_pt[idx] = x - dt; - err_cur += e_pc[idx] * e_pc[idx] * w; - err_try += e_pt[idx] * e_pt[idx] * w; - } - } - err_cur += hex_spectral_penalty(e_pc, QK_K); - err_try += hex_spectral_penalty(e_pt, QK_K); - if (err_try < err_cur) { - dm = dm_try; - mm = mm_try; - pol_improved = 1; - } - } + hex_polish_block(block_x, blk_imat, &dm, &mm, pl_Ls, pl_Lm, L); + + /* Write back polished codes and scalars */ + for (int j = 0; j < N_SUB; j++) + output[blk].scales[j] = pl_Ls[j] | (pl_Lm[j] << 4); + output[blk].d = gguf_fp32_to_fp16(dm); + output[blk].dmin = gguf_fp32_to_fp16(mm); + } + + /* ══ PHASE 4.7: CANDIDATE FLOOR (worst-case bound) ══ + * + * candidate_errors[blk][c] is the EXACT weighted SSE of a directly + * encodable configuration (fp16 d/dmin + derived Ls/Lm + nearest + * rounding vs the original weights). The multi-stage assembly + * usually improves on its seed, but each stage optimises a slightly + * different objective and coordinate descent can land in a worse + * basin. Compare the finished block against the best raw candidate + * and fall back when the pipeline ended up worse — guaranteeing + * final weighted SSE ≤ min_c candidate_errors[blk][c]. + * + * CHANGED: the fallback used to adopt the raw candidate as-is, so a + * block that hit the floor shipped a configuration that had never + * been through the polish at all — strictly worse than it needed to + * be. Now the floor state is handed to the same monotone coordinate + * descent, which can only lower its error further. The bound above + * still holds (polish is accept-only-if-better). */ + { + float e_f[QK_K]; + float fin_err = hex_block_ext_err(block_x, blk_imat, dm, mm, + output[blk].scales, L, e_f); + + float g_best = candidate_errors[blk][0]; + int g_cand = 0; + for (int c = 1; c < TOTAL_SCALE_CANDIDATES; c++) { + if (candidate_errors[blk][c] < g_best) { + g_best = candidate_errors[blk][c]; + g_cand = c; } } - if (!pol_improved) { - /* ── (3) ±2 ULP joint (d, dmin) micro-search vs ORIGINAL ── - * The closed-form refit rounds its real-valued optimum to - * fp16, which can land 1–2 ULP away from the best - * representable pair (and the earlier ±8 ULP search ran - * against the DC-shifted objective). With codes fixed, - * scan the (2·HEX_POLISH_ULP+1)² fp16 neighborhood on the - * true objective; - * accept only strict improvement, then loop once more so - * move (1) can re-optimise codes for the new scalars. - * Monotone ⇒ final RMSE can only decrease. */ - uint16_t base_d16 = gguf_fp32_to_fp16(dm); - uint16_t base_m16 = gguf_fp32_to_fp16(mm); - - float cur_err = 0.0f; - float e_u[QK_K]; + if (g_best < fin_err) { + /* Rebuild the block exactly as the candidate was scored */ + float c_dm = gguf_fp16_to_fp32(candidate_d [blk][g_cand]); + float c_mm = gguf_fp16_to_fp32(candidate_dmin[blk][g_cand]); + uint8_t c_Ls[16], c_Lm[16]; + hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, + c_dm, c_mm, c_Ls, c_Lm); for (int j = 0; j < N_SUB; j++) { - float d_sub = dm * (float)pl_Ls[j]; - float m_sub = mm * (float)pl_Lm[j]; + float d_sub = c_dm * (float)c_Ls[j]; + float m_sub = c_mm * (float)c_Lm[j]; for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - e_u[idx] = block_x[idx] - - (d_sub * (float)L[idx] - m_sub); - cur_err += e_u[idx] * e_u[idx] * w; - } - } - cur_err += hex_spectral_penalty(e_u, QK_K); - - float best_err = cur_err; - uint16_t best_d16 = base_d16, best_m16 = base_m16; - for (int dd = -HEX_POLISH_ULP; dd <= HEX_POLISH_ULP; dd++) { - int cd16 = (int)base_d16 + dd; - if (cd16 < 0 || cd16 > 0x7BFF) continue; - float t_dm = gguf_fp16_to_fp32((uint16_t)cd16); - for (int dmm = -HEX_POLISH_ULP; dmm <= HEX_POLISH_ULP; dmm++) { - if (dd == 0 && dmm == 0) continue; - int cm16 = (int)base_m16 + dmm; - if (cm16 < 0 || cm16 > 0x7BFF) continue; - float t_mm = gguf_fp16_to_fp32((uint16_t)cm16); - - float err = 0.0f; - /* SSE-partial prune valid only without the - * signable vesica/wave credit */ - for (int j = 0; - j < N_SUB && (HEX_VW_LAMBDA != 0.0f || - err < best_err); j++) { - float d_sub = t_dm * (float)pl_Ls[j]; - float m_sub = t_mm * (float)pl_Lm[j]; - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - e_u[idx] = block_x[idx] - - (d_sub * (float)L[idx] - m_sub); - err += e_u[idx] * e_u[idx] * w; - } - } - if (HEX_DC_LAMBDA != 0.0f || HEX_VW_LAMBDA != 0.0f) - err = (err < best_err || HEX_VW_LAMBDA != 0.0f) - ? err + hex_spectral_penalty(e_u, QK_K) - : err; - if (err < best_err) { - best_err = err; - best_d16 = (uint16_t)cd16; - best_m16 = (uint16_t)cm16; + int idx = 16 * j + k; + int q = 0; + if (d_sub >= 1e-15f) { + q = gguf_nearest_int((block_x[idx] + m_sub) / d_sub); + if (q < 0) q = 0; if (q > 3) q = 3; } + L[idx] = (uint8_t)q; } } - if (best_d16 != base_d16 || best_m16 != base_m16) { - dm = gguf_fp16_to_fp32(best_d16); - mm = gguf_fp16_to_fp32(best_m16); - pol_improved = 1; - } + dm = c_dm; mm = c_mm; + + /* Polish the floor state on the true objective */ + hex_polish_block(block_x, blk_imat, &dm, &mm, c_Ls, c_Lm, L); + + for (int j = 0; j < N_SUB; j++) + output[blk].scales[j] = c_Ls[j] | (c_Lm[j] << 4); + output[blk].d = gguf_fp32_to_fp16(dm); + output[blk].dmin = gguf_fp32_to_fp16(mm); } + } - if (!pol_improved) break; /* converged on true objective */ + /* ── Multi-start bookkeeping: score this restart, keep the best ── */ + { + float e_ms[QK_K]; + float err_ms = hex_block_ext_err(block_x, blk_imat, dm, mm, + output[blk].scales, L, e_ms); + if (!ms_have || err_ms < ms_best_err) { + ms_best_err = err_ms; + ms_best_d16 = output[blk].d; + ms_best_m16 = output[blk].dmin; + memcpy(ms_best_scales, output[blk].scales, N_SUB); + memcpy(ms_best_L, L, QK_K); + ms_have = 1; + } } - /* Write back polished codes and scalars */ - for (int j = 0; j < N_SUB; j++) - output[blk].scales[j] = pl_Ls[j] | (pl_Lm[j] << 4); - output[blk].d = gguf_fp32_to_fp16(dm); - output[blk].dmin = gguf_fp32_to_fp16(mm); - } + } /* ── end multi-start restart loop ── */ - /* ══ PHASE 4.7: CANDIDATE FLOOR (worst-case bound) ══ - * - * candidate_errors[blk][c] is the EXACT weighted SSE of a directly - * encodable configuration (fp16 d/dmin + derived Ls/Lm + nearest - * rounding vs the original weights). The multi-stage assembly - * (DC-shifted WLS, shaping, diffusion, polish) usually improves on - * its seed, but each stage optimises a slightly different objective - * and coordinate descent can land in a worse basin. Compare the - * finished block against the best raw candidate and fall back when - * the pipeline ended up worse — guaranteeing - * final weighted SSE ≤ min_c candidate_errors[blk][c]. */ + /* ── Commit the winning restart ── */ { - float fin_err = 0.0f; - float e_f[QK_K]; + output[blk].d = ms_best_d16; + output[blk].dmin = ms_best_m16; + memcpy(output[blk].scales, ms_best_scales, N_SUB); + for (int j = 0; j < N_SUB; j++) { - float d_sub = dm * (float)(output[blk].scales[j] & 0xF); - float m_sub = mm * (float)(output[blk].scales[j] >> 4); - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - float w = (imat_importance) ? - imat_importance[blk * QK_K + idx] : 1.0f; - e_f[idx] = block_x[idx] - - (d_sub * (float)L[idx] - m_sub); - fin_err += e_f[idx] * e_f[idx] * w; - } - } - fin_err += hex_spectral_penalty(e_f, QK_K); - - float g_best = candidate_errors[blk][0]; - int g_cand = 0; - for (int c = 1; c < TOTAL_SCALE_CANDIDATES; c++) { - if (candidate_errors[blk][c] < g_best) { - g_best = candidate_errors[blk][c]; - g_cand = c; - } + int ls = ms_best_scales[j] & 0xF; + int lm = ms_best_scales[j] >> 4; + sat_subs++; + if (ls == 15) sat_ls_hi++; + if (lm == 15) sat_lm_hi++; + if (ls == 0) sat_ls_lo++; } - if (g_best < fin_err) { - /* Rebuild the block exactly as the candidate was scored */ - float c_dm = gguf_fp16_to_fp32(candidate_d [blk][g_cand]); - float c_mm = gguf_fp16_to_fp32(candidate_dmin[blk][g_cand]); - uint8_t c_Ls[16], c_Lm[16]; - hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, - c_dm, c_mm, c_Ls, c_Lm); - for (int j = 0; j < N_SUB; j++) { - float d_sub = c_dm * (float)c_Ls[j]; - float m_sub = c_mm * (float)c_Lm[j]; - for (int k = 0; k < 16; k++) { - int idx = 16 * j + k; - int q = 0; - if (d_sub >= 1e-15f) { - q = gguf_nearest_int((block_x[idx] + m_sub) / d_sub); - if (q < 0) q = 0; if (q > 3) q = 3; - } - L[idx] = (uint8_t)q; - } - output[blk].scales[j] = c_Ls[j] | (c_Lm[j] << 4); + for (int j = 0; j < QK_K; j += 128) { + for (int l = 0; l < 32; l++) { + output[blk].qs[j / 4 + l] = ms_best_L[j + l] + | (ms_best_L[j + l + 32] << 2) + | (ms_best_L[j + l + 64] << 4) + | (ms_best_L[j + l + 96] << 6); } - dm = c_dm; mm = c_mm; - output[blk].d = candidate_d [blk][g_cand]; - output[blk].dmin = candidate_dmin[blk][g_cand]; } - } - for (int j = 0; j < QK_K; j += 128) { - for (int l = 0; l < 32; l++) { - output[blk].qs[j / 4 + l] = L[j + l] - | (L[j + l + 32] << 2) - | (L[j + l + 64] << 4) - | (L[j + l + 96] << 6); + float berr = gguf_q2_k_block_error(block_x, &output[blk]); + if (isnan(berr)) { + printf("NaN block error at blk %ld! d16=%u dmin16=%u\n", + (long)blk, (unsigned)output[blk].d, + (unsigned)output[blk].dmin); + for (int j = 0; j < N_SUB; j++) + printf("Ls[%d]=%d Lm[%d]=%d\n", j, + ms_best_scales[j] & 0xF, j, ms_best_scales[j] >> 4); + exit(1); } + total_err += berr; } - - float berr = gguf_q2_k_block_error(block_x, &output[blk]); - if (isnan(berr)) { - printf("NaN block error at blk %ld! dm=%f mm=%f\n", (long)blk, dm, mm); - for (int j=0; j<16; j++) printf("Ls[%d]=%d Lm[%d]=%d\n", j, Ls_blk[j], j, Lm_blk[j]); - exit(1); - } - total_err += berr; } for (int _ti = 0; _ti < _n_omp_threads; _ti++) @@ -4489,10 +4800,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, free(candidate_d); free(candidate_dmin); free(best_candidate); - if (out_total_error) *out_total_error = total_err; - + if (out_total_error) *out_total_error = (float)total_err; + if (verbose) { - float rmse = sqrtf(total_err / (float)n_elements); + /* Divide in double, narrow once at the end. */ + float rmse = (float)sqrt(total_err / (double)n_elements); double w_sum2 = 0.0; for (int64_t i = 0; i < n_elements; i++) @@ -4521,14 +4833,30 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, (long long)n_elements, (long long)(n_elements / QK_K)); printf(" │ Weight σ: %-12.4e Range: [%.4e, %.4e] │\n", w_sigma, w_sigma * -4.0f, w_sigma * 4.0f); - printf(" │ Total MSE: %-12.6f │\n", total_err); + printf(" │ Total MSE: %-12.6g │\n", total_err); printf(" │ RMSE: %-12.4e RMSE/σ: %-8.4f │\n", rmse, rmse_over_sigma); printf(" │ Fidelity: %s %-14s │\n", fidelity_icon, fidelity_class); - printf(" │ Engine: Sieve sequential (log-sieve + parity) │\n"); + printf(" │ Engine: Sieve sequential + %d-start assembly │\n", + (int)HEX_TOPK_RESTARTS); + printf(" │ Objective: lam_dc=%.2f lam_vw=%.2f dc_decay=%.2f%s │\n", + (double)HEX_DC_LAMBDA, (double)HEX_VW_LAMBDA, (double)DC_DECAY, + (HEX_DC_LAMBDA == 0.0f && HEX_VW_LAMBDA == 0.0f && + DC_DECAY == 0.0f) ? " [pure RMSE]" : " "); + if (sat_subs > 0) { + double p_hi = 100.0 * (double)sat_ls_hi / (double)sat_subs; + double p_mi = 100.0 * (double)sat_lm_hi / (double)sat_subs; + double p_lo = 100.0 * (double)sat_ls_lo / (double)sat_subs; + printf(" │ Sub-scale sat: ls=15 %6.3f%% lm=15 %6.3f%% ls=0 %6.3f%% │\n", + p_hi, p_mi, p_lo); + if (p_hi > 0.10 || p_mi > 0.10) + printf(" │ ^ clamped sub-scales dominate block error: widen the │\n" + " │ upper tail of NEIGHBOR_MULTS_D / NEIGHBOR_MULTS_M. │\n"); + } printf(" └─────────────────────────────────────────────────────────────────┘\n"); } + (void)sat_ls_hi; (void)sat_lm_hi; (void)sat_ls_lo; (void)sat_subs; } @@ -4537,7 +4865,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, * ═══════════════════════════════════════════════════════════════════════════ */ static void print_progress_bar(int current, int total, const char *label, - time_t start_time) + time_t start_time) { if (total <= 0) return; float pct = (float)current / (float)total; @@ -4567,11 +4895,11 @@ static void print_progress_bar(int current, int total, const char *label, * ═══════════════════════════════════════════════════════════════════════════ */ static int write_gguf(const char *output_path, const STMultiFile *mf, - const ModelArchitecture *arch, - const TokenizerData *tokenizer, - OptimizerMode opt_mode, - const IMatrixData *imatrix, - int verbose) + const ModelArchitecture *arch, + const TokenizerData *tokenizer, + OptimizerMode opt_mode, + const IMatrixData *imatrix, + int verbose) { FILE *fp = fopen(output_path, "wb"); if (!fp) { @@ -4677,7 +5005,7 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, data_offset += tensor_sizes[i]; data_offset = (data_offset + GGUF_DEFAULT_ALIGNMENT - 1) & - ~(uint64_t)(GGUF_DEFAULT_ALIGNMENT - 1); + ~(uint64_t)(GGUF_DEFAULT_ALIGNMENT - 1); } /* ── Write header ── */ @@ -4721,24 +5049,24 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, if (has_tokenizer) { gguf_write_kv_string(fp, "tokenizer.ggml.model", tokenizer->model_type); gguf_write_kv_string_array(fp, "tokenizer.ggml.tokens", - (const char **)tokenizer->tokens, - (uint64_t)tokenizer->vocab_size); + (const char **)tokenizer->tokens, + (uint64_t)tokenizer->vocab_size); gguf_write_kv_float32_array(fp, "tokenizer.ggml.scores", - tokenizer->scores, - (uint64_t)tokenizer->vocab_size); - gguf_write_kv_int32_array(fp, "tokenizer.ggml.token_type", - tokenizer->token_types, + tokenizer->scores, (uint64_t)tokenizer->vocab_size); + gguf_write_kv_int32_array(fp, "tokenizer.ggml.token_type", + tokenizer->token_types, + (uint64_t)tokenizer->vocab_size); gguf_write_kv_uint32(fp, "tokenizer.ggml.bos_token_id", - (uint32_t)tokenizer->bos_id); + (uint32_t)tokenizer->bos_id); gguf_write_kv_uint32(fp, "tokenizer.ggml.eos_token_id", - (uint32_t)tokenizer->eos_id); + (uint32_t)tokenizer->eos_id); gguf_write_kv_uint32(fp, "tokenizer.ggml.unknown_token_id", - (uint32_t)tokenizer->unk_id); + (uint32_t)tokenizer->unk_id); if (tokenizer->n_merges > 0) { gguf_write_kv_string_array(fp, "tokenizer.ggml.merges", - (const char **)tokenizer->merges, - (uint64_t)tokenizer->n_merges); + (const char **)tokenizer->merges, + (uint64_t)tokenizer->n_merges); } printf(" Tokenizer metadata written (%d tokens, %d merges)\n\n", tokenizer->vocab_size, tokenizer->n_merges); @@ -4754,8 +5082,8 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, dims[d] = (uint64_t)ti->shape[nd - 1 - d]; } gguf_write_tensor_info(fp, gguf_names[i], - ti->n_dims, dims, - tensor_types[i], tensor_offsets[i]); + ti->n_dims, dims, + tensor_types[i], tensor_offsets[i]); } /* ── Alignment padding before data section ── */ @@ -4809,7 +5137,7 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, const float *imp = NULL; if (imatrix) { const IMatrixEntry *ime = imatrix_find_any(imatrix, - gguf_names[i], tensor_hf_names[i]); + gguf_names[i], tensor_hf_names[i]); if (ime && ime->n_values > 0) { imp = ime->normalized; if (verbose) @@ -4819,8 +5147,8 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, } quantize_tensor_q2k_hpc(f32_data, n_elements, - quant_data, &tensor_error, - opt_mode, imp, verbose); + quant_data, &tensor_error, + opt_mode, imp, verbose); fwrite(quant_data, sizeof(BlockQ2K), n_blocks, fp); @@ -4884,7 +5212,7 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, const float *imp = NULL; if (imatrix) { const IMatrixEntry *ime = imatrix_find_any(imatrix, - gguf_names[i], tensor_hf_names[i]); + gguf_names[i], tensor_hf_names[i]); if (ime && ime->n_values > 0) { imp = ime->normalized; if (verbose) @@ -4894,8 +5222,8 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, } quantize_tensor_q4_0_hpc(f32_data, n_elements, - q4_data, &tensor_error, - imp, verbose); + q4_data, &tensor_error, + imp, verbose); fwrite(q4_data, sizeof(BlockQ4_0), n_blocks_q4, fp); @@ -4986,14 +5314,14 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, original_f32_size += ti->n_elements * sizeof(float); } float compression_ratio = (original_f32_size > 0) ? - (float)original_f32_size / (float)final_size : 0.0f; + (float)original_f32_size / (float)final_size : 0.0f; float effective_bpw = (total_elements_quantized > 0) ? - 8.0f * (float)total_bytes_quantized / (float)total_elements_quantized : - 0.0f; + 8.0f * (float)total_bytes_quantized / (float)total_elements_quantized : + 0.0f; float total_rmse = (total_elements_quantized > 0) ? - sqrtf(total_error_sum / (float)total_elements_quantized) : 0.0f; + sqrtf(total_error_sum / (float)total_elements_quantized) : 0.0f; float mean_mse_per_tensor = (quant_count > 0) ? - total_error_sum / (float)quant_count : 0.0f; + total_error_sum / (float)quant_count : 0.0f; const char *overall_fid, *overall_icon; if (total_rmse <= 1.0e-04f) { overall_fid = "ULTRA (≤1e-04)"; overall_icon = "★★★★"; } @@ -5062,6 +5390,368 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, return 0; } +/* ═══════════════════════════════════════════════════════════════════════════ + * ANALOG Q2_K — OUTPUT-ERROR (HESSIAN) OBJECTIVE + * + * Every other path in this file minimises weight reconstruction error, + * L_w = sum_i w_i (x_i - deq_i)^2 + * and the bound analysis shows that objective is exhausted: the existing + * pipeline lands within 0.1-2.6% of the best RMSE the Q2_K format can + * physically achieve. There is nothing left there. + * + * The quantity that actually determines whether the MODEL works is the + * layer output error. For a linear layer y = Wx with calibration + * activations X and H = X^T X / n: + * + * L_out(row r) = (w_r - what_r)^T H (w_r - what_r) + * + * Crucially L_out DECOMPOSES EXACTLY over output rows -- row r of What + * only affects output r -- so each row is an independent subproblem and + * the whole tensor parallelises over rows with the FULL row Hessian, no + * approximation. (Measured: full-row H beats 256x256 block-local H by 11 + * points of output error, so restricting to block-local to fit the + * existing per-block loop would throw away a third of the gain.) + * + * Under this objective nearest-rounding is NOT optimal, and the codes that + * minimise L_out have visibly WORSE weight RMSE -- that is the trade being + * made deliberately, not a bug. Do not judge this path by RMSE. + * + * Three exact coordinate-descent moves, each accept-only-if-better, so + * L_out is monotonically non-increasing and the routine can never make a + * row worse than its nearest-rounding seed: + * + * A codes q_i in {0..3} best of 4 in closed form. With residual + * e = w - what and gradient cache g = H e, + * changing what_i by delta changes L_out by + * exactly -2*delta*g_i + delta^2*H_ii. + * B (ls_j, lm_j) 4-bit full 16x16 per sub-block; the 16 weights move + * together so the delta is a 16-vector and the + * quadratic term uses the 16x16 diagonal chunk. + * C (d, dmin) fp16 2x2 H-weighted least squares per super-block, + * rounded to fp16, accepted on strict decrease. + * + * COST: move A is O(n_in^2) per row per sweep (the g -= delta*H[j,:] + * update). That is ~HEX_HESS_SWEEPS times a single GPTQ pass. Set + * HEX_HESS_FLUSH > 1 to accumulate deltas across that many columns before + * refreshing g -- FLUSH = 16 (one sub-block) cuts move A by 16x and loses + * little, because within-sub-block coupling is handled exactly by move B. + * + * H is supplied by the caller as n_in*n_in row-major float32. It is NOT + * derivable from a standard llama.cpp imatrix, which stores only + * diag(X^T X); see hexstate_q2k_hessian_bytes() and the loader in the + * Python driver for the file format. + * ═══════════════════════════════════════════════════════════════════════════ */ + +#ifndef HEX_HESS_SWEEPS +#define HEX_HESS_SWEEPS 8 +#endif +#ifndef HEX_HESS_FLUSH +#define HEX_HESS_FLUSH 1 + /* 1 = exact; 16 = one sub-block per refresh */ +#endif + +/* Quadratic form e^T H e for one row. */ +static double hex_quad(const float *e, const float *H, int64_t n) +{ + double t = 0.0; + for (int64_t i = 0; i < n; i++) { + if (e[i] == 0.0f) continue; + const float *Hi = H + i * n; + double s = 0.0; + for (int64_t k = 0; k < n; k++) s += (double)Hi[k] * (double)e[k]; + t += (double)e[i] * s; + } + return t; +} + +/* Optimise one output row of length n_in (a whole number of Q2_K blocks). */ +static void hex_analog_row(const float *w, int64_t n_in, const float *H, + BlockQ2K *out, int n_sweeps, + float *E, float *Gd, uint8_t *q, + uint8_t *Ls, uint8_t *Lm, float *dv, float *mv) +{ + const int NS = QK_K / 16; /* 16 sub-blocks per block */ + const int64_t NB = n_in / QK_K; + + /* ── seed: the existing Phase-1 nearest-rounding construction ── */ + for (int64_t b = 0; b < NB; b++) { + const float *bx = w + b * QK_K; + float scales[16], mins[16], sw[16], wt[16]; + uint8_t Lq[QK_K], Laux[16]; + float sumx2 = 0.0f; + for (int i = 0; i < QK_K; i++) sumx2 += bx[i] * bx[i]; + float sigma2 = sumx2 / (float)QK_K; + for (int j = 0; j < NS; j++) { + const float *sx = bx + 16 * j; + sw[j] = 0.0f; + for (int l = 0; l < 16; l++) { + wt[l] = sqrtf(sigma2 + sx[l] * sx[l]); + sw[j] += wt[l]; + } + scales[j] = hpc_make_qkx2_quants(16, 3, sx, wt, Lq + 16 * j, + &mins[j], Laux); + } + uint8_t Lsb[16], Lmb[16]; + float dd = hpc_make_qp_quants(NS, 15, scales, Lsb, sw); + float mm = hpc_make_qp_quants(NS, 15, mins, Lmb, sw); + dd = gguf_fp16_to_fp32(gguf_fp32_to_fp16(dd)); + mm = gguf_fp16_to_fp32(gguf_fp32_to_fp16(mm)); + dv[b] = dd; mv[b] = mm; + for (int j = 0; j < NS; j++) { + Ls[b * NS + j] = Lsb[j]; + Lm[b * NS + j] = Lmb[j]; + float ds = dd * (float)Lsb[j], ms = mm * (float)Lmb[j]; + for (int k = 0; k < 16; k++) { + int idx = (int)(b * QK_K) + 16 * j + k; + int qq = 0; + if (ds >= 1e-15f) { + qq = gguf_nearest_int((w[idx] + ms) / ds); + if (qq < 0) qq = 0; if (qq > 3) qq = 3; + } + q[idx] = (uint8_t)qq; + } + } + } + + /* residual and gradient cache */ + for (int64_t i = 0; i < n_in; i++) { + int64_t b = i / QK_K, j = i / 16; + float ds = dv[b] * (float)Ls[j], ms = mv[b] * (float)Lm[j]; + E[i] = w[i] - (ds * (float)q[i] - ms); + } + for (int64_t i = 0; i < n_in; i++) { + const float *Hi = H + i * n_in; + double s = 0.0; + for (int64_t k = 0; k < n_in; k++) s += (double)Hi[k] * (double)E[k]; + Gd[i] = (float)s; + } + + for (int sweep = 0; sweep < n_sweeps; sweep++) { + int64_t changed = 0; + + /* ── move A: per-weight best of 4 ── */ + { + int64_t pend = 0; + float *acc = E; /* deltas applied directly to E; g refreshed below */ + (void)acc; + for (int64_t j = 0; j < n_in; j++) { + int64_t b = j / QK_K, sj = j / 16; + float ds = dv[b] * (float)Ls[sj]; + if (ds < 1e-15f) continue; + float Hjj = H[j * n_in + j]; + int qc = q[j]; + float bestd = 0.0f, bestg = 0.0f; int bestq = qc; + for (int qn = 0; qn < 4; qn++) { + float dq = ((float)qn - (float)qc) * ds; + float gain = 2.0f * dq * Gd[j] - dq * dq * Hjj; + if (gain > bestg + 1e-20f) { bestg = gain; bestd = dq; bestq = qn; } + } + if (bestd != 0.0f) { + q[j] = (uint8_t)bestq; + E[j] -= bestd; + changed++; + if (HEX_HESS_FLUSH <= 1) { + const float *Hj = H + j * n_in; + for (int64_t k = 0; k < n_in; k++) Gd[k] -= bestd * Hj[k]; + } else { + pend++; + } + } + if (HEX_HESS_FLUSH > 1 && pend && ((j + 1) % HEX_HESS_FLUSH == 0)) { + for (int64_t i = 0; i < n_in; i++) { + const float *Hi = H + i * n_in; + double s = 0.0; + for (int64_t k = 0; k < n_in; k++) s += (double)Hi[k] * (double)E[k]; + Gd[i] = (float)s; + } + pend = 0; + } + } + if (HEX_HESS_FLUSH > 1) { + for (int64_t i = 0; i < n_in; i++) { + const float *Hi = H + i * n_in; + double s = 0.0; + for (int64_t k = 0; k < n_in; k++) s += (double)Hi[k] * (double)E[k]; + Gd[i] = (float)s; + } + } + } + + /* ── move B: joint (ls, lm) per sub-block, full 16x16 ── */ + for (int64_t sj = 0; sj < n_in / 16; sj++) { + int64_t b = (sj * 16) / QK_K, c0 = sj * 16; + float dd = dv[b], mm = mv[b]; + float old[16], del[16], bestdel[16]; + for (int k = 0; k < 16; k++) + old[k] = dd * (float)Ls[sj] * (float)q[c0 + k] - mm * (float)Lm[sj]; + float bestg = 0.0f; int bl = Ls[sj], bm = Lm[sj]; int found = 0; + for (int tl = 0; tl <= 15; tl++) { + for (int tm = 0; tm <= 15; tm++) { + float lin = 0.0f, quad = 0.0f; + for (int k = 0; k < 16; k++) { + float nw = dd * (float)tl * (float)q[c0 + k] - mm * (float)tm; + del[k] = nw - old[k]; + lin += del[k] * Gd[c0 + k]; + } + for (int a = 0; a < 16; a++) { + if (del[a] == 0.0f) continue; + const float *Ha = H + (c0 + a) * n_in + c0; + float s = 0.0f; + for (int c = 0; c < 16; c++) s += Ha[c] * del[c]; + quad += del[a] * s; + } + float gain = 2.0f * lin - quad; + if (gain > bestg + 1e-20f) { + bestg = gain; bl = tl; bm = tm; found = 1; + memcpy(bestdel, del, sizeof(del)); + } + } + } + if (found) { + Ls[sj] = (uint8_t)bl; Lm[sj] = (uint8_t)bm; + for (int k = 0; k < 16; k++) E[c0 + k] -= bestdel[k]; + for (int64_t i = 0; i < n_in; i++) { + const float *Hi = H + i * n_in + c0; + float s = 0.0f; + for (int k = 0; k < 16; k++) s += Hi[k] * bestdel[k]; + Gd[i] -= s; + } + changed++; + } + } + + /* ── move C: (d, dmin) per super-block, 2x2 H-weighted LS ── */ + for (int64_t b = 0; b < NB; b++) { + int64_t c0 = b * QK_K; + double Saa=0, Sab=0, Sbb=0, Sat=0, Sbt=0; + /* a_i = ls*q, b_i = lm ; target t = w - (residual from other blocks) */ + static const int NSl = QK_K / 16; + float av[QK_K], bvv[QK_K], tv[QK_K]; + for (int j = 0; j < NSl; j++) + for (int k = 0; k < 16; k++) { + int i = 16 * j + k; + av[i] = (float)Ls[b * NSl + j] * (float)q[c0 + i]; + bvv[i] = (float)Lm[b * NSl + j]; + tv[i] = w[c0 + i]; + } + for (int i = 0; i < QK_K; i++) { + const float *Hi = H + (c0 + i) * n_in + c0; + double sa=0, sb=0, st=0; + for (int k = 0; k < QK_K; k++) { + sa += (double)Hi[k] * av[k]; + sb += (double)Hi[k] * bvv[k]; + st += (double)Hi[k] * tv[k]; + } + Saa += av[i]*sa; Sab += av[i]*sb; Sbb += bvv[i]*sb; + Sat += av[i]*st; Sbt += bvv[i]*st; + } + double det = Saa * Sbb - Sab * Sab; + if (fabs(det) < 1e-30) continue; + double dn = ( Sbb * Sat + Sab * Sbt) / det; + double mn = (-Sab * Sat - Saa * Sbt) / det; + if (dn <= 0.0) continue; + if (mn < 0.0) mn = 0.0; + float dtry = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)dn)); + float mtry = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)mn)); + /* exact accept test on the block's contribution */ + float eo[QK_K], en[QK_K]; + for (int i = 0; i < QK_K; i++) { + eo[i] = E[c0 + i]; + en[i] = w[c0 + i] - (dtry * av[i] - mtry * bvv[i]); + } + double Lo=0, Ln=0; + for (int i = 0; i < QK_K; i++) { + const float *Hi = H + (c0 + i) * n_in + c0; + double so=0, sn=0; + for (int k = 0; k < QK_K; k++) { so += Hi[k]*eo[k]; sn += Hi[k]*en[k]; } + Lo += eo[i]*so; Ln += en[i]*sn; + } + if (Ln < Lo) { + dv[b] = dtry; mv[b] = mtry; + for (int i = 0; i < QK_K; i++) { + float d0 = en[i] - E[c0 + i]; + E[c0 + i] = en[i]; + (void)d0; + } + for (int64_t i = 0; i < n_in; i++) { + const float *Hi = H + i * n_in; + double s = 0.0; + for (int64_t k = 0; k < n_in; k++) s += (double)Hi[k] * (double)E[k]; + Gd[i] = (float)s; + } + changed++; + } + } + + if (!changed) break; + } + + /* ── emit Q2_K blocks ── */ + for (int64_t b = 0; b < NB; b++) { + BlockQ2K *o = &out[b]; + o->d = gguf_fp32_to_fp16(dv[b]); + o->dmin = gguf_fp32_to_fp16(mv[b]); + for (int j = 0; j < NS; j++) + o->scales[j] = (uint8_t)(Ls[b * NS + j] | (Lm[b * NS + j] << 4)); + const uint8_t *qq = q + b * QK_K; + for (int j = 0; j < QK_K; j += 128) + for (int l = 0; l < 32; l++) + o->qs[j / 4 + l] = qq[j + l] + | (qq[j + l + 32] << 2) + | (qq[j + l + 64] << 4) + | (qq[j + l + 96] << 6); + } +} + +/* Tensor driver: parallel over output rows. */ +static void quantize_tensor_q2k_analog(const float *weights, int64_t n_elements, + int64_t n_in, const float *H, + BlockQ2K *output, float *out_total_error, + int n_sweeps, int verbose) +{ + int64_t n_rows = n_elements / n_in; + int64_t nb_row = n_in / QK_K; + double total_err = 0.0; + + #pragma omp parallel for schedule(dynamic, 1) reduction(+:total_err) + for (int64_t r = 0; r < n_rows; r++) { + float *E = (float *)malloc((size_t)n_in * sizeof(float)); + float *Gd = (float *)malloc((size_t)n_in * sizeof(float)); + uint8_t *q = (uint8_t *)malloc((size_t)n_in); + uint8_t *Ls = (uint8_t *)malloc((size_t)(n_in / 16)); + uint8_t *Lm = (uint8_t *)malloc((size_t)(n_in / 16)); + float *dv = (float *)malloc((size_t)nb_row * sizeof(float)); + float *mv = (float *)malloc((size_t)nb_row * sizeof(float)); + + const float *w = weights + r * n_in; + BlockQ2K *o = output + r * nb_row; + hex_analog_row(w, n_in, H, o, n_sweeps, E, Gd, q, Ls, Lm, dv, mv); + + for (int64_t b = 0; b < nb_row; b++) + total_err += (double)gguf_q2_k_block_error(w + b * QK_K, &o[b]); + + free(E); free(Gd); free(q); free(Ls); free(Lm); free(dv); free(mv); + } + + if (out_total_error) *out_total_error = (float)total_err; + + if (verbose) { + double sig2 = 0.0; + for (int64_t i = 0; i < n_elements; i++) + sig2 += (double)weights[i] * (double)weights[i]; + double sg = sqrt(sig2 / (double)n_elements); + double rmse = sqrt(total_err / (double)n_elements); + printf("\n +---- Analog Q2_K (output-error objective) ----------------------+\n"); + printf(" | Rows: %-10lld n_in: %-8lld sweeps: %-3d |\n", + (long long)n_rows, (long long)n_in, n_sweeps); + printf(" | Weight RMSE: %-12.4e RMSE/sigma: %-8.4f |\n", + rmse, (sg > 1e-15) ? rmse / sg : 0.0); + printf(" | NOTE: weight RMSE is EXPECTED to be higher than the RMSE- |\n"); + printf(" | optimal path. Judge this path by perplexity, not by RMSE. |\n"); + printf(" +---------------------------------------------------------------+\n"); + } +} + /* ═══════════════════════════════════════════════════════════════════════════ * LIBRARY API — Exported functions for Python ctypes integration * @@ -5093,30 +5783,66 @@ void hexstate_init(void) * verbose: 1 for per-block diagnostics */ void hexstate_quantize_tensor_q2k(const float *weights, int64_t n_elements, - void *output, float *out_error, - int opt_mode, int verbose) + void *output, float *out_error, + int opt_mode, int verbose) { hexstate_init(); quantize_tensor_q2k_hpc(weights, n_elements, - (BlockQ2K *)output, out_error, - (OptimizerMode)opt_mode, NULL, verbose); + (BlockQ2K *)output, out_error, + (OptimizerMode)opt_mode, NULL, verbose); } /* Same as above but with importance matrix weights */ void hexstate_quantize_tensor_q2k_imat(const float *weights, int64_t n_elements, - void *output, float *out_error, - int opt_mode, - const float *imat_importance, - int verbose) + void *output, float *out_error, + int opt_mode, + const float *imat_importance, + int verbose) { hexstate_init(); quantize_tensor_q2k_hpc(weights, n_elements, - (BlockQ2K *)output, out_error, - (OptimizerMode)opt_mode, imat_importance, verbose); + (BlockQ2K *)output, out_error, + (OptimizerMode)opt_mode, imat_importance, verbose); } /* Get the block size for Q2_K (84 bytes per 256 elements) */ int hexstate_q2k_block_bytes(void) { return sizeof(BlockQ2K); } + +/* Bytes of Hessian required for a tensor whose rows are n_in long. */ +int64_t hexstate_q2k_hessian_bytes(int64_t n_in) +{ + return n_in * n_in * (int64_t)sizeof(float); +} + +/* ANALOG Q2_K: assign fresh 2-bit codes minimising LAYER OUTPUT error. + * + * weights row-major, n_elements = n_rows * n_in + * n_in row length; MUST be a multiple of QK_K (256) + * H n_in*n_in row-major float32 = X^T X / n_samples, damped. + * NOT a standard imatrix (that is only diag(X^T X)). + * n_sweeps coordinate-descent sweeps; <= 0 uses HEX_HESS_SWEEPS + * + * Falls back to the RMSE pipeline when H is NULL or n_in is unusable, so a + * caller can pass H opportunistically per tensor. */ +void hexstate_quantize_tensor_q2k_analog(const float *weights, + int64_t n_elements, int64_t n_in, + const float *H, void *output, + float *out_error, int n_sweeps, + int verbose) +{ + if (!H || n_in <= 0 || (n_in % QK_K) != 0 || (n_elements % n_in) != 0) { + if (verbose) + printf(" [analog] no usable Hessian (n_in=%lld) -- " + "falling back to RMSE pipeline\n", (long long)n_in); + quantize_tensor_q2k_hpc(weights, n_elements, (BlockQ2K *)output, + out_error, OPT_HYBRID, NULL, verbose); + return; + } + quantize_tensor_q2k_analog(weights, n_elements, n_in, H, + (BlockQ2K *)output, out_error, + (n_sweeps > 0) ? n_sweeps : HEX_HESS_SWEEPS, + verbose); +} int hexstate_q2k_block_elements(void) { return QK_K; } /* HPC-optimized Q4_0 quantization for attention tensors. @@ -5129,15 +5855,15 @@ int hexstate_q2k_block_elements(void) { return QK_K; } * verbose: 1 for per-block diagnostics */ void hexstate_quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements, - void *output, float *out_error, - const float *imat_importance, - int verbose) + void *output, float *out_error, + const float *imat_importance, + int verbose) { hexstate_init(); float err = 0.0f; quantize_tensor_q4_0_hpc(weights, n_elements, - (BlockQ4_0 *)output, &err, - imat_importance, verbose); + (BlockQ4_0 *)output, &err, + imat_importance, verbose); if (out_error) *out_error = err; } @@ -5374,7 +6100,7 @@ int main(int argc, char **argv) /* ── Phase 3-5: Quantize and write GGUF ── */ printf(" Phase 3: HPC-Optimized Q2_K Quantization + GGUF Output...\n"); int result = write_gguf(output_path, mf, &arch, tokenizer, - opt_mode, imatrix, verbose); + opt_mode, imatrix, verbose); /* Wall-clock total: clock() sums CPU time over all OpenMP threads */ time_t t_end = time(NULL); @@ -5385,4 +6111,4 @@ int main(int argc, char **argv) st_multi_close(mf); return result; } -#endif /* HEXSTATE_LIBRARY */ \ No newline at end of file +#endif /* HEXSTATE_LIBRARY */