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,47 @@ 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, "gemma4") == 0 || + strcmp(cfg.model_type, "gemma4_text") == 0 || + strcmp(cfg.model_type, "gemma4_unified") == 0 || + strcmp(cfg.model_type, "gemma4_unified_text") == 0) { + strcpy(arch->architecture, "gemma4"); + } else if (strcmp(cfg.model_type, "gemma") == 0 || + strcmp(cfg.model_type, "gemma2") == 0 || + strcmp(cfg.model_type, "gemma3") == 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); } @@ -322,12 +328,17 @@ static void detect_architecture(const STMultiFile *mf, ModelArchitecture *arch, /* Architecture-specific detection */ int has_qkv_proj = count_tensors_with_prefix(mf, "model.layers.0.self_attn.qkv_proj"); int has_kv_a_proj = count_tensors_with_prefix(mf, "model.layers.0.self_attn.kv_a_proj_with_mqa"); + int has_gemma4 = count_tensors_with_prefix(mf, "model.layers.0.inp_gate") || + count_tensors_with_prefix(mf, "model.layers.0.layer_output_scale") || + count_tensors_with_prefix(mf, "model.layers.0.post_ffw_norm_1"); int has_final_norm = (st_multi_find_tensor(mf, "model.final_norm.weight") >= 0); if (has_qkv_proj > 0 && !cfg.valid) { strcpy(arch->architecture, "phi3"); } else if (has_kv_a_proj > 0 && !cfg.valid) { strcpy(arch->architecture, "llama"); /* DeepSeek uses llama arch */ + } else if (has_gemma4 && !cfg.valid) { + strcpy(arch->architecture, "gemma4"); } else if (has_final_norm && !cfg.valid) { strcpy(arch->architecture, "gemma"); } @@ -563,6 +574,21 @@ static void map_tensor_name(const char *hf_name, char *gguf_name, int buflen) * - Keep F32: norms, biases, embeddings, 1D tensors * ═══════════════════════════════════════════════════════════════════════════ */ +static inline uint64_t gguf_row_width(const STTensorInfo *ti) +{ + return (ti && ti->n_dims > 0) ? (uint64_t)ti->shape[ti->n_dims - 1] : 0; +} + +static inline int q2k_row_compatible(const STTensorInfo *ti) +{ + return ti && ti->n_dims >= 2 && gguf_row_width(ti) % QK_K == 0; +} + +static inline int q4_row_compatible(const STTensorInfo *ti) +{ + return ti && ti->n_dims >= 2 && gguf_row_width(ti) % QK4_0 == 0; +} + static int should_quantize(const STTensorInfo *ti, const char *gguf_name) { /* Never quantize 1D tensors (norms, biases) */ @@ -586,11 +612,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. @@ -645,6 +671,26 @@ static int is_attention_tensor(const char *gguf_name) #define N_CAND_M 24 /* dmin multiplier candidates (expanded) */ #define TOTAL_SCALE_CANDIDATES (N_CAND_D * N_CAND_M) +static const float HEX_NEIGHBOR_MULTS_D[N_CAND_D] = { + 0.780f, 0.835f, 0.880f, 0.915f, 0.943f, 0.963f, + 0.978f, 0.988f, 0.994f, 0.997f, 0.999f, 1.000f, + 1.002f, 1.005f, 1.011f, 1.021f, 1.035f, 1.054f, + 1.080f, 1.115f, 1.160f, 1.215f, 1.275f, 1.340f +}; +static const float HEX_NEIGHBOR_MULTS_M[N_CAND_M] = { + 0.750f, 0.800f, 0.840f, 0.870f, 0.900f, 0.920f, + 0.940f, 0.955f, 0.970f, 0.985f, 0.995f, 1.000f, + 1.005f, 1.015f, 1.030f, 1.045f, 1.060f, 1.080f, + 1.100f, 1.130f, 1.160f, 1.200f, 1.250f, 1.300f +}; +static inline void hex_candidate_pair(float base_d, float base_m, int cidx, + uint16_t *d16, uint16_t *m16) +{ + int di = cidx / N_CAND_M, mi = cidx % N_CAND_M; + *d16 = gguf_fp32_to_fp16(base_d * HEX_NEIGHBOR_MULTS_D[di]); + *m16 = gguf_fp32_to_fp16(base_m * HEX_NEIGHBOR_MULTS_M[mi]); +} + /* ════════════════════════════════════════════════════════════════════════ * EXPERIMENTAL / CURRENTLY-UNUSED CODE PATHS * @@ -727,8 +773,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 +813,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 +844,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 +860,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 +879,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 +961,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,30 +1075,13 @@ 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 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. */ + * 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. */ 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]; @@ -1071,20 +1100,19 @@ 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_sse = 0; + float best_mad = 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_sse += w[i] * diff * diff; + best_mad += w[i] * fabsf(diff); } - /* 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); + /* 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); float sl = 0, sl2 = 0, sxl = 0; for (int i = 0; i < n; i++) { int l = gguf_nearest_int(try_iscale * (x[i] - xmin)); @@ -1103,14 +1131,14 @@ static float hpc_make_qkx2_quants(int n, int nmax, const float *x, this_min = 0; this_scale = sxl / sl2; } - float sse = 0; + float mad = 0; for (int i = 0; i < n; i++) { float diff = this_scale * (float)Laux[i] + this_min - x[i]; - sse += w[i] * diff * diff; + mad += w[i] * fabsf(diff); } - if (sse < best_sse) { + if (mad < best_mad) { for (int i = 0; i < n; i++) L[i] = Laux[i]; - best_sse = sse; + best_mad = mad; scale = this_scale; xmin = this_min; } @@ -1123,16 +1151,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; @@ -1253,7 +1281,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++) { @@ -1286,7 +1314,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]; @@ -1306,8 +1334,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--) { @@ -1338,8 +1366,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++) { @@ -1411,7 +1439,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) { @@ -1503,8 +1531,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; @@ -1513,7 +1541,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); @@ -1531,7 +1559,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++) { @@ -1558,7 +1586,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) { @@ -1581,13 +1609,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) { @@ -1609,7 +1637,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; @@ -1692,53 +1720,25 @@ static const int Q4_CAND_TO_QUHIT[Q4_N_CAND] = { * one being optimised mis-ranks them; aligning the two strictly lowers the * final weighted RMSE (and is bit-identical when no imatrix is used). */ -/* ── Cross-block prior override ratio ── - * Q2_K and Q4_0 blocks are decoded INDEPENDENTLY by every GGUF runtime: - * there is no cross-block coupling in the dequantizer, so a smoothness - * prior that keeps a block on a worse candidate can only raise the true - * reconstruction RMSE. With 1.00f the per-block argmin over the candidate - * grid always wins (provably optimal seed for the assembly phase); the HPC - * graph/Viterbi/Born machinery still shapes ties and seeds the search. - * Set to e.g. 0.95f to restore the old 5%-hysteresis smoothness prior. */ +/* ── HExState preservation threshold ───────────────────────────────────── + * The HExState/Viterbi path is the primary global optimizer. A conventional + * per-block argmin is allowed to override it only when it is materially + * better on the actual candidate objective. With 0.995f the local winner + * must improve the current HExState proposal by at least 0.5%%. This prevents + * the old 1.00f setting from collapsing the exotic search back into a + * conventional independent-block quantizer. + * + * Set to 1.0f for a pure local-candidate floor; lower values give HExState + * progressively more authority over the final candidate field. */ #ifndef HEX_GREEDY_OVERRIDE_RATIO -#define HEX_GREEDY_OVERRIDE_RATIO 1.00f +#define HEX_GREEDY_OVERRIDE_RATIO 0.995f #endif /* fp16-ULP radius of the monotone (d, dmin) micro-search in the Phase-4.6 * 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 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 +#define HEX_POLISH_ULP 4 #endif /* ── DC + vesica/wave extended objective (dot-product error cancellation) ── @@ -1777,29 +1777,19 @@ 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 0.0f +#define HEX_DC_LAMBDA 1.0f #endif #ifndef HEX_VW_LAMBDA -#define HEX_VW_LAMBDA 0.0f +#define HEX_VW_LAMBDA 1.0f #endif -/* DEFAULT CHANGED (1,1) → (0,0): PURE-RMSE MODE. - * - * (1, 1) is the unit-strength spectral prior. Empirically (synthetic +/* Default (1, 1): 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. - * - * λ = 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. */ + * improving monotonically through λ = 4 at ~+0.1% RMSE. Set both to + * 0.0f to recover the exact pure weighted-SSE / minimum-RMSE pipeline. */ /* 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 @@ -1814,16 +1804,80 @@ 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; +} + +/* Robust temperature estimator for the HExState measurement model. + * + * The old path estimated T from the mean of each block's MAXIMUM candidate + * error. That measures the width of the worst tail, not the local energy + * landscape near the optimum, and can make exp(-(E-Emin)/(2T)) almost flat. + * We instead use the median excess energy E-Emin across a small sample of + * blocks. The resulting temperature tracks the actual candidate basin and + * keeps the six-state marginals informative without becoming brittle. */ +static float hex_q2k_robust_temperature(const float *candidate_errors, + int64_t n_blocks) +{ + enum { SAMPLE_BLOCKS = 128, MAX_CAND = TOTAL_SCALE_CANDIDATES }; + float sample[MAX_CAND]; + float medians[SAMPLE_BLOCKS]; + int ns = 0; + + if (!candidate_errors || n_blocks <= 0) return 1e-10f; + + int64_t step = n_blocks / SAMPLE_BLOCKS; + if (step < 1) step = 1; + + for (int64_t b = 0; b < n_blocks && ns < SAMPLE_BLOCKS; b += step) { + const float *row = candidate_errors + b * TOTAL_SCALE_CANDIDATES; + float mn = row[0]; + for (int c = 1; c < TOTAL_SCALE_CANDIDATES; c++) + if (row[c] < mn) mn = row[c]; + + for (int c = 0; c < TOTAL_SCALE_CANDIDATES; c++) + sample[c] = fmaxf(0.0f, row[c] - mn); + + /* insertion sort is cheap here (<=128 × 576) and avoids a qsort + * callback in the hot tensor loop. */ + for (int i = 1; i < TOTAL_SCALE_CANDIDATES; i++) { + float v = sample[i]; + int j = i - 1; + while (j >= 0 && sample[j] > v) { + sample[j + 1] = sample[j]; + --j; + } + sample[j + 1] = v; + } + medians[ns++] = sample[TOTAL_SCALE_CANDIDATES / 2]; + } + + if (ns == 0) return 1e-10f; + + for (int i = 1; i < ns; i++) { + float v = medians[i]; + int j = i - 1; + while (j >= 0 && medians[j] > v) { + medians[j + 1] = medians[j]; + --j; + } + medians[j + 1] = v; + } + + float t = medians[ns / 2]; + if (!(t > 1e-10f) || !isfinite(t)) t = 1e-10f; + return t; } 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) { + if (!weights || !output || n_elements <= 0 || n_elements % QK4_0 != 0) { + if (out_total_error) *out_total_error = -1.0f; + return; + } int64_t n_blocks = n_elements / QK4_0; - /* double: same float-saturation hazard as the Q2_K accumulator. */ - double total_err = 0.0; + float total_err = 0.0f; (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)); @@ -1843,9 +1897,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++) { @@ -1864,7 +1918,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]; @@ -1883,7 +1937,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 ── */ @@ -1920,572 +1974,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 */ - - 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]); + 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 */ - /* 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; - } - } + 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]); - /* 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]; + /* 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++) - agg_errors[c] = 0.0f; + 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; + } + } - 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; + /* 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; - 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; - } + 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; - float min_err = 1e30f; + 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++) - if (agg_errors[c] < min_err) - min_err = agg_errors[c]; + agg_errors[c] *= inv_gs; + } - 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; - } + float min_err = 1e30f; + for (int c = 0; c < Q4_N_CAND; c++) + if (agg_errors[c] < min_err) + min_err = agg_errors[c]; - 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]); + 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; } - /* Neighbor edges */ - for (int64_t i = 0; i < graph_blocks - 1; i++) - hpc_cz(graph, i, i + 1); + 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]); + } - /* ── 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); + } + + 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++; } + } - 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; - } + 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; } + 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; + 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 hist_idx = i * Q4_N_BEAMS + k; - history[hist_idx].cand_idx = cand; - history[hist_idx].parent_idx = beams[sb].history_idx; + int hist_idx = i * Q4_N_BEAMS + k; + history[hist_idx].cand_idx = cand; + history[hist_idx].parent_idx = beams[sb].history_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; + 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 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; - } + 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; - } + 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]]; - } + } + free(history); - /* 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; } - } + /* ══════════════════════════════════════════════════════════════ + * 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]]; + } - /* 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)); - shot_sparse_q4[gi] = best_bin_cand; - shot_err += cand_errors[blk][best_bin_cand]; + 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; } } - /* 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; + /* 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]; } - free(shot_sparse_q4); + /* 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; + } } - /* 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; - } + 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; } - 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 - * ══════════════════════════════════════════════════════════════════ */ + /* ══════════════════════════════════════════════════════════════════ + * 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]; + #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]; - /* Start from the grid-selected scale (the "assembled frequency") */ - float d_current = gguf_fp16_to_fp32(cand_d16[blk][cidx]); + /* Start from the grid-selected scale (the "assembled frequency") */ + float d_current = gguf_fp16_to_fp32(cand_d16[blk][cidx]); - /* 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; + /* 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; - 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; - } + 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; + } - 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); - } + 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 ── */ + /* ── 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; + } { - 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; + 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)base_d16 + delta; + int cand16 = (int)neg_d16 + delta; if (cand16 >= 0 && cand16 <= 0x7BFF) ulp_candidates[n_ulp++] = (uint16_t)cand16; } - { - 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; - } - } + } - 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]; - } + 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; - 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)); - /* ── D₆ Hadamard Error Shaping with Simulated Annealing ── */ - int q_base[QK4_0], q_shaped[QK4_0]; - float q_cont[QK4_0]; + { + float e_live[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; + float deq = ((float)q_shaped[j] - 8.0f) * actual_d; + e_live[j] = bw[j] - deq; } - memcpy(q_shaped, q_base, QK4_0 * sizeof(int)); - { - 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; - } + 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 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; - } + /* 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 = (float)total_err; - free(greedy_d); - free(cand_errors); - free(cand_d16); - free(best_candidate); + *out_total_error = total_err; + free(greedy_d); + free(cand_errors); + free(cand_d16); + free(best_candidate); } /* ════════════════════════════════════════════════════════════════════════ @@ -2544,15 +2598,18 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements, float *out_total_error, const float *imat_importance, int verbose) { + if (!weights || !output || n_elements <= 0 || n_elements % QK8_0 != 0) { + if (out_total_error) *out_total_error = -1.0f; + return; + } int64_t n_blocks = n_elements / QK8_0; - /* double: same float-saturation hazard as the Q2_K accumulator. */ - double total_err = 0.0; + float total_err = 0.0f; (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); @@ -2650,12 +2707,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); @@ -2720,7 +2777,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; } } @@ -2769,14 +2826,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) { @@ -2788,7 +2845,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; @@ -2824,7 +2881,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 = (float)total_err; + if (out_total_error) *out_total_error = total_err; } @@ -2849,365 +2906,18 @@ static inline void hex_derive_subscales(const float *scales, const float *mins, } } -/* 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) -{ - const int N_SUB = QK_K / 16; - 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; - } - } - 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) + BlockQ2K *output, float *out_total_error, + OptimizerMode opt_mode, + const float *imat_importance, + int verbose) { + if (!weights || !output || n_elements <= 0 || n_elements % QK_K != 0) { + if (out_total_error) *out_total_error = -1.0f; + return; + } 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; + float total_err = 0.0f; const int N_SUB = QK_K / 16; /* ── Outlier Clamping for WLS Seeds ── @@ -3233,9 +2943,9 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t 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; + 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 @@ -3243,12 +2953,16 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, typedef struct { float dm, mm; - uint16_t d_fp16, dmin_fp16; + float base_dm, base_mm; uint8_t Ls[16], Lm[16]; float scales[16], mins[16], sw[16]; } BlockSeed; BlockSeed *seeds = (BlockSeed *)calloc(n_blocks, sizeof(BlockSeed)); + if (!seeds) { + if (out_total_error) *out_total_error = -1.0f; + return; + } #pragma omp parallel for schedule(dynamic, 64) for (int64_t blk = 0; blk < n_blocks; blk++) { @@ -3278,58 +2992,39 @@ 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].d_fp16 = gguf_fp32_to_fp16(seeds[blk].dm); - seeds[blk].dmin_fp16 = gguf_fp32_to_fp16(seeds[blk].mm); + seeds[blk].Lm, seeds[blk].sw); + seeds[blk].base_dm = seeds[blk].dm; + seeds[blk].base_mm = seeds[blk].mm; } /* ══════════════════════════════════════════════════════════════════ * PHASE 2: WLS-Optimal Candidate Generation * ══════════════════════════════════════════════════════════════════ */ - /* Expanded neighborhood around WLS optimum: ±30% with 24 candidates */ - /* d is the sensitive axis, so concentrate resolution near 1.0 while - * keeping wide tails for blocks whose WLS seed is off. 1.000 stays at - * index 11 so the neutral-candidate fallback/init remains valid. */ - static const float NEIGHBOR_MULTS_D[N_CAND_D] = { - 0.780f, 0.835f, 0.880f, 0.915f, 0.943f, 0.963f, - 0.978f, 0.988f, 0.994f, 0.997f, 0.999f, 1.000f, - 1.002f, 1.005f, 1.011f, 1.021f, 1.035f, 1.054f, - 1.080f, 1.115f, 1.160f, 1.215f, 1.275f, 1.340f - }; - static const float NEIGHBOR_MULTS_M[N_CAND_M] = { - 0.750f, 0.800f, 0.840f, 0.870f, 0.900f, 0.920f, - 0.940f, 0.955f, 0.970f, 0.985f, 0.995f, 1.000f, - 1.005f, 1.015f, 1.030f, 1.045f, 1.060f, 1.080f, - 1.100f, 1.130f, 1.160f, 1.200f, 1.250f, 1.300f - }; - /* Map 24 candidates → 6 quhit states for BP encoding */ static const int CAND_TO_QUHIT[24] = { 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5 }; float (*candidate_errors)[TOTAL_SCALE_CANDIDATES] = NULL; - uint16_t (*candidate_d)[TOTAL_SCALE_CANDIDATES] = NULL; - uint16_t (*candidate_dmin)[TOTAL_SCALE_CANDIDATES] = NULL; - candidate_errors = (float (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(float[TOTAL_SCALE_CANDIDATES])); - candidate_d = (uint16_t (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); - candidate_dmin = (uint16_t (*)[TOTAL_SCALE_CANDIDATES])calloc(n_blocks, - sizeof(uint16_t[TOTAL_SCALE_CANDIDATES])); + sizeof(float[TOTAL_SCALE_CANDIDATES])); + if (!seeds || !candidate_errors) { + free(seeds); free(candidate_errors); + if (out_total_error) *out_total_error = -1.0f; + return; + } /* 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 @@ -3379,7 +3074,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; @@ -3415,20 +3110,19 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, } /* ── Step 2b: Generate Candidates ── */ + seeds[blk].base_dm = wls_dm; + seeds[blk].base_mm = wls_mm; for (int di = 0; di < N_CAND_D; di++) { - float trial_dm = wls_dm * NEIGHBOR_MULTS_D[di]; + float trial_dm = wls_dm * HEX_NEIGHBOR_MULTS_D[di]; uint16_t trial_d16 = gguf_fp32_to_fp16(trial_dm); float actual_dm = gguf_fp16_to_fp32(trial_d16); for (int mi = 0; mi < N_CAND_M; mi++) { int cidx = di * N_CAND_M + mi; - float trial_mm = wls_mm * NEIGHBOR_MULTS_M[mi]; + float trial_mm = wls_mm * HEX_NEIGHBOR_MULTS_M[mi]; uint16_t trial_dmin16 = gguf_fp32_to_fp16(trial_mm); float actual_mm = gguf_fp16_to_fp32(trial_dmin16); - candidate_d[blk][cidx] = trial_d16; - candidate_dmin[blk][cidx] = trial_dmin16; - uint8_t trial_Ls[16], trial_Lm[16]; for (int j = 0; j < N_SUB; j++) { if (actual_dm > 1e-15f) { @@ -3467,7 +3161,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); } } } @@ -3480,821 +3174,740 @@ 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; - - HPCGraph *graph = hpc_create(n_sites); - if (graph) { - for (int64_t i = 0; i < n_sites; i++) - triality_dft(&graph->locals[i]); + if (opt_mode != OPT_MSE && n_blocks >= 2) { + /* Give HExState a materially larger contiguous neighborhood before + * falling back to stride-group aggregation. 8192 is also the cap + * used by the experimental graph builder and keeps graph state bounded. */ + int64_t graph_blocks = (n_blocks > 8192) ? 8192 : n_blocks; + int64_t stride = n_blocks / graph_blocks; + float temperature = 0.5f; + int64_t n_sites = graph_blocks * QUHITS_PER_BLOCK; - { - 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; - } - } + HPCGraph *graph = hpc_create(n_sites); + if (graph) { + for (int64_t i = 0; i < n_sites; i++) + triality_dft(&graph->locals[i]); - 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; + /* Robust local-landscape temperature. This replaces the old + * mean-of-max-errors heuristic, which could flatten the Born + * amplitudes and make HExState effectively irrelevant. */ + temperature = hex_q2k_robust_temperature(&candidate_errors[0][0], + n_blocks); - 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 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 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; - } + 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; - float min_err = 1e30f; + for (int64_t b = blk_start; b < blk_end; b++) { 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)); - } - } - 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; - } + 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; + } - double fine_re[6]; - double fine_norm = 0.0; - for (int qi = 0; qi < 6; qi++) fine_re[qi] = 0.0; + 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 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; + 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; + } - 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; + 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)); } - 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 (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 (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); - } + 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]); + } + + 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); - { - int64_t vit_gi, vit_b; - int vit_s, vit_sp; + /* ══════════════════════════════════════════════════════════════ + * 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_energy[i][s] = soft-min/free-energy of 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_energy[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; - /* 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; - } - } - } + /* Soft-min/free-energy bin reduction. + * Instead of throwing away every candidate except the + * single minimum, retain the width of the local basin: + * + * F_bin = Emin - 2T log(mean(exp(-(E-Emin)/(2T)))) + * + * A broad cluster of nearly-equivalent candidates therefore + * survives into the HExState/Viterbi decision. The actual + * code emitted later is still the exact minimum candidate + * inside the chosen bin, so this changes the SEARCH prior, + * not GGUF decoding semantics. */ + double bin_sum[VIT_N_STATES]; + int bin_count[VIT_N_STATES]; + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + bin_sum[vit_s] = 0.0; + bin_count[vit_s] = 0; + } - /* 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)); + float bin_min[VIT_N_STATES]; + float rep_bin_min[VIT_N_STATES]; + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + bin_min[vit_s] = 1e30f; + rep_bin_min[vit_s] = 1e30f; + } + + for (vit_b = vit_gi * stride; + vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; + vit_b++) { + for (int 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 < bin_min[vit_s]) + bin_min[vit_s] = e; + if (vit_b == vit_gi * stride && e < rep_bin_min[vit_s]) + rep_bin_min[vit_s] = e; } } - /* ── 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++; - } + /* The emitted representative candidate must be chosen + * from the representative block itself, not whichever + * block happened to define the group's global bin minimum. */ + for (int 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_gi * stride][vit_c]; + if (e <= rep_bin_min[vit_s]) { + vit_bin_cand[vit_gi][vit_s] = vit_c; } } - 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 ── */ + for (vit_b = vit_gi * stride; + vit_b < (vit_gi + 1) * stride && vit_b < n_blocks; + vit_b++) { + for (int 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]; + float ex = e - bin_min[vit_s]; + double a = (temperature > 1e-30f) + ? exp(-(double)ex / (2.0 * (double)temperature)) + : (ex <= 0.0f ? 1.0 : 0.0); + bin_sum[vit_s] += a; + bin_count[vit_s]++; - /* 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] - : 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; - - 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; + for (vit_s = 0; vit_s < VIT_N_STATES; vit_s++) { + if (bin_count[vit_s] > 0) { + double mean_exp = bin_sum[vit_s] / (double)bin_count[vit_s]; + if (mean_exp < 1e-300) mean_exp = 1e-300; + vit_bin_err[vit_gi][vit_s] = + bin_min[vit_s] - + 2.0f * temperature * (float)log(mean_exp); } } - /* ── 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; - } + /* 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++; } - 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]]; } + } + 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 ── */ + + /* 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] + : 1e30f; + vit_dp [0][vit_s] = local; + vit_back[0][vit_s] = -1; + } - /* ── Step E: Map Viterbi path → best_candidate[] ── */ - for (vit_gi = 0; vit_gi < graph_blocks; vit_gi++) { - vit_s = vit_path[vit_gi]; + /* 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; - 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; - } + + 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; } - best_candidate[vit_b] = best_c; + } + vit_dp [vit_gi][vit_s] = (best_pred < 1e29f) + ? best_pred + local : 1e30f; + vit_back[vit_gi][vit_s] = best_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_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 F: 5 % greedy override (pure MSE safety net) ── */ - for (vit_b = 0; vit_b < n_blocks; vit_b++) { + /* ── 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 cur_err = candidate_errors[vit_b][best_candidate[vit_b]]; - float g_best = cur_err; - int g_cand = best_candidate[vit_b]; + float best_e = 1e30f; + int best_c = best_candidate[blk_rep]; 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 (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 (g_best < cur_err * HEX_GREEDY_OVERRIDE_RATIO) - best_candidate[vit_b] = g_cand; + best_candidate[vit_b] = best_c; } + } - free(vit_path); - free(vit_dp); - free(vit_back); - free(vit_bin_err); - free(vit_bin_cand); - free(vit_log_pri); + /* ── Step F: 0.5 % local-MSE rescue (HExState remains primary) ── */ + 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(coarse_marg); - free(fine_marg); - hpc_destroy(graph); + free(vit_path); + free(vit_dp); + free(vit_back); + free(vit_bin_err); + free(vit_bin_cand); + free(vit_log_pri); } - } 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; - } + + 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; } - 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). - * ══════════════════════════════════════════════════════════════════ */ + /* ══════════════════════════════════════════════════════════════════ + * 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). + * ══════════════════════════════════════════════════════════════════ */ - /* 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; + #define DC_DECAY 0.85f /* Boundary-condition leak factor (0 = isolated, 1 = full) */ + + float *block_dc_bias = (float *)calloc(n_blocks, sizeof(float)); + + 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]; + uint16_t c_d16, c_m16; + hex_candidate_pair(seeds[blk].base_dm, seeds[blk].base_mm, cidx, &c_d16, &c_m16); + float dm0 = gguf_fp16_to_fp32(c_d16); + float mm0 = gguf_fp16_to_fp32(c_m16); + + 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; } + rolling_dc = dc_res; } + } + + /* ══════════════════════════════════════════════════════════════════ + * PHASE 4: Assemble blocks via least-squares (d, dmin) extraction + * ══════════════════════════════════════════════════════════════════ */ - /* ══════════════════════════════════════════════════════════════════ - * 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) + 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); + + #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; - const float *blk_imat = imat_importance ? imat_importance + blk * QK_K - : NULL; + const float *block_x = weights + blk * QK_K; + int cidx = best_candidate[blk]; + uint8_t Ls_blk[16], Lm_blk[16]; - /* ══════════════════════════════════════════════════════════════════ - * 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; + /* ── 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; + + /* 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]; { - /* 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]; - } - } + int _i; + for (_i = 0; _i < QK_K; _i++) + adj_block_x[_i] = block_x[_i] - dc_adj; } - /* 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; - - for (int kstart = 0; kstart < n_topk; kstart++) { + uint16_t base_c_d16, base_c_m16; + hex_candidate_pair(seeds[blk].base_dm, seeds[blk].base_mm, cidx, &base_c_d16, &base_c_m16); + float dm = gguf_fp16_to_fp32(base_c_d16); + float mm = gguf_fp16_to_fp32(base_c_m16); - /* 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; + hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, + dm, mm, Ls_blk, Lm_blk); - int cidx = topk_c[kstart]; - uint8_t Ls_blk[16], Lm_blk[16]; - - /* ── 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; - - /* 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; - } - - float dm = gguf_fp16_to_fp32(candidate_d[blk][cidx]); - float mm = gguf_fp16_to_fp32(candidate_dmin[blk][cidx]); - - hex_derive_subscales(seeds[blk].scales, seeds[blk].mins, - dm, mm, Ls_blk, Lm_blk); - - uint16_t prev_dm16 = 0, prev_mm16 = 0; - for (int ls_iter = 0; ls_iter < 5; ls_iter++) { + uint16_t prev_dm16 = 0, prev_mm16 = 0; + for (int ls_iter = 0; ls_iter < 5; ls_iter++) { - uint8_t state_ls[N_SUB][6]; - uint8_t state_lm[N_SUB][6]; - float state_err[N_SUB][6]; + uint8_t state_ls[N_SUB][6]; + uint8_t state_lm[N_SUB][6]; + float state_err[N_SUB][6]; - 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 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 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 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; + 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; + } - 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; + 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; } } } } + } - 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; - } - 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; 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]); + } - 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; - } - 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; - } + 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; } - - 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; - } + 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 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 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; } - - 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); - - 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); + 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; + uint16_t seed_d16, seed_m16; + hex_candidate_pair(seeds[blk].base_dm, seeds[blk].base_mm, cidx, &seed_d16, &seed_m16); + float d_seed = gguf_fp16_to_fp32(seed_d16); + float m_seed = gguf_fp16_to_fp32(seed_m16); + 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)); } - 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; + 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); + + 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 = sx[k]; + float x = adj_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; int q; if (d_sub < 1e-15f) { q = 0; } else { @@ -4303,491 +3916,732 @@ 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; - sub_err += diff * diff * w; + 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; + } + 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; } + 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); + 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); + 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]; + /* ── 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; - } + 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; + } + 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 */ + + for (int k = 0; k < QK_K; k++) { + int jj = k >> 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; - } + if (d_s < 1e-15f) continue; - 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; - } - } + 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; - if (best_k < 0) break; /* converged — no further improvement */ + 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]; - 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; + 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; } } - /* 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); + if (best_k < 0) break; /* converged — no further improvement */ + + q_shaped_all[best_k] = best_q_alt; { - 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]); + 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; } + } - } else { - /* lambda = 0 path: commit nearest rounding directly. */ + /* 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); for (int i = 0; i < QK_K; i++) - L[i] = (uint8_t)q_base_all[i]; - } + L[i] = (uint8_t)(use_shaped ? q_shaped_all[i] : 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. - * - * 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; + /* ── 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; + } - 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; - } + /* 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; + } + } + } + } - /* 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; + /* ── 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; } } + 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); + } - /* ── 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 */ + /* ══ 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; 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); + 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 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; + 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; } + dc_tot += sub_dc[j]; } - /* 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 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]; } - 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 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++) { - 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; + 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; } } - 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; + 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; + } } - 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; + /* ── (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; + } + } } } - 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); + 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 = c_dm * (float)c_Ls[j]; - float m_sub = c_mm * (float)c_Lm[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; - 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; + 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; } - L[idx] = (uint8_t)q; } } - 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 (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; + } } - } - /* ── 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; - } + if (!pol_improved) break; /* converged on true objective */ } - } /* ── end multi-start restart loop ── */ + /* 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); + } - /* ── Commit the winning restart ── */ + /* ══ 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]. */ { - output[blk].d = ms_best_d16; - output[blk].dmin = ms_best_m16; - memcpy(output[blk].scales, ms_best_scales, N_SUB); - + float fin_err = 0.0f; + float e_f[QK_K]; for (int j = 0; j < N_SUB; j++) { - 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++; + 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; + } } - 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); + if (g_best < fin_err) { + /* Rebuild the block exactly as the candidate was scored */ + uint16_t c_d16, c_m16; + hex_candidate_pair(seeds[blk].base_dm, seeds[blk].base_mm, g_cand, &c_d16, &c_m16); + float c_dm = gguf_fp16_to_fp32(c_d16); + float c_mm = gguf_fp16_to_fp32(c_m16); + 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); } + dm = c_dm; mm = c_mm; + output[blk].d = c_d16; + output[blk].dmin = c_m16; } + } - 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); + 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); } - 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++) @@ -4797,14 +4651,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements, free(block_dc_bias); free(seeds); free(candidate_errors); - free(candidate_d); - free(candidate_dmin); free(best_candidate); - if (out_total_error) *out_total_error = (float)total_err; - + if (out_total_error) *out_total_error = total_err; + if (verbose) { - /* Divide in double, narrow once at the end. */ - float rmse = (float)sqrt(total_err / (double)n_elements); + float rmse = sqrtf(total_err / (float)n_elements); double w_sum2 = 0.0; for (int64_t i = 0; i < n_elements; i++) @@ -4833,30 +4684,14 @@ 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.6g │\n", total_err); + printf(" │ Total MSE: %-12.6f │\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 + %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(" │ Engine: Sieve sequential (log-sieve + parity) │\n"); printf(" └─────────────────────────────────────────────────────────────────┘\n"); } - (void)sat_ls_hi; (void)sat_lm_hi; (void)sat_ls_lo; (void)sat_subs; } @@ -4865,7 +4700,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; @@ -4895,11 +4730,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) { @@ -4913,6 +4748,7 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, /* ── Determine which tensors to include ── */ int *include_list = (int *)calloc(mf->n_tensors, sizeof(int)); + if (!include_list) { fclose(fp); remove(output_path); return -1; } int n_include = 0; for (int i = 0; i < mf->n_tensors; i++) { if (!should_skip_tensor(mf->tensor_map[i].name)) { @@ -4971,6 +4807,11 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, uint64_t *tensor_offsets = calloc(total_tensors, sizeof(uint64_t)); int *tensor_src_idx = calloc(total_tensors, sizeof(int)); /* map to unified ST index */ char (*tensor_hf_names)[ST_MAX_NAME_LEN] = calloc(total_tensors, ST_MAX_NAME_LEN); + if (!gguf_names || !tensor_types || !tensor_sizes || !tensor_offsets || !tensor_src_idx || !tensor_hf_names) { + free(include_list); free(gguf_names); free(tensor_types); free(tensor_sizes); + free(tensor_offsets); free(tensor_src_idx); free(tensor_hf_names); + fclose(fp); remove(output_path); return -1; + } GGMLType quant_type = GGML_TYPE_Q2_K; @@ -4982,16 +4823,18 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, tensor_src_idx[i] = src; if (should_quantize(ti, gguf_names[i])) { - if (is_attention_tensor(gguf_names[i])) { + if (is_attention_tensor(gguf_names[i]) && q4_row_compatible(ti)) { tensor_types[i] = GGML_TYPE_Q4_0; - int64_t n_blocks_q4 = (ti->n_elements + QK4_0 - 1) / QK4_0; - tensor_sizes[i] = n_blocks_q4 * sizeof(BlockQ4_0); - if (verbose) - printf(" [ATTN→Q4_0] %s (%ld elements)\n", - gguf_names[i], (long)ti->n_elements); - } else { + tensor_sizes[i] = (ti->n_elements / QK4_0) * sizeof(BlockQ4_0); + } else if (q2k_row_compatible(ti)) { tensor_types[i] = quant_type; tensor_sizes[i] = ggml_type_size(quant_type, ti->n_elements); + } else if (q4_row_compatible(ti)) { + tensor_types[i] = GGML_TYPE_Q4_0; + tensor_sizes[i] = (ti->n_elements / QK4_0) * sizeof(BlockQ4_0); + } else { + tensor_types[i] = GGML_TYPE_F16; + tensor_sizes[i] = ti->n_elements * (int64_t)sizeof(uint16_t); } } else if (ti->n_dims >= 2) { tensor_types[i] = GGML_TYPE_F16; @@ -5005,7 +4848,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 ── */ @@ -5049,24 +4892,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); + 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); + 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); @@ -5082,8 +4925,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 ── */ @@ -5110,34 +4953,24 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, if (!f32_data) { fprintf(stderr, "\n ERROR: Failed to convert tensor '%s' to F32\n", ti->name); - continue; + goto write_fail; } int64_t n_elements = ti->n_elements; float tensor_error = 0.0f; - int64_t padded = (n_elements + QK_K - 1) / QK_K * QK_K; - if (padded > n_elements) { - float *grown = realloc(f32_data, padded * sizeof(float)); - if (!grown) { - fprintf(stderr, "\n ERROR: Out of memory padding '%s'\n", - ti->name); - free(f32_data); - continue; - } - f32_data = grown; - for (int64_t j = n_elements; j < padded; j++) - f32_data[j] = 0.0f; - n_elements = padded; + if (!q2k_row_compatible(ti) || n_elements % QK_K != 0) { + fprintf(stderr, "\n ERROR: Q2_K row alignment violation for '%s'\n", ti->name); + free(f32_data); goto write_fail; } - int64_t n_blocks = n_elements / QK_K; BlockQ2K *quant_data = calloc(n_blocks, sizeof(BlockQ2K)); + if (!quant_data) { free(f32_data); goto write_fail; } 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) @@ -5147,10 +4980,13 @@ 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); + if (tensor_error < 0.0f) { + free(quant_data); free(f32_data); goto write_fail; + } - fwrite(quant_data, sizeof(BlockQ2K), n_blocks, fp); + if (fwrite(quant_data, sizeof(BlockQ2K), n_blocks, fp) != (size_t)n_blocks) { free(quant_data); free(f32_data); goto write_fail; } float rmse = sqrtf(tensor_error / (float)ti->n_elements); @@ -5185,34 +5021,24 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, if (!f32_data) { fprintf(stderr, "\n ERROR: Failed to convert tensor '%s' to F32\n", ti->name); - continue; + goto write_fail; } int64_t n_elements = ti->n_elements; - int64_t padded = (n_elements + QK4_0 - 1) / QK4_0 * QK4_0; - if (padded > n_elements) { - float *grown = realloc(f32_data, padded * sizeof(float)); - if (!grown) { - fprintf(stderr, "\n ERROR: Out of memory padding '%s'\n", - ti->name); - free(f32_data); - continue; - } - f32_data = grown; - for (int64_t j = n_elements; j < padded; j++) - f32_data[j] = 0.0f; - n_elements = padded; + if (!q4_row_compatible(ti) || n_elements % QK4_0 != 0) { + fprintf(stderr, "\n ERROR: Q4_0 row alignment violation for '%s'\n", ti->name); + free(f32_data); goto write_fail; } - int64_t n_blocks_q4 = n_elements / QK4_0; BlockQ4_0 *q4_data = calloc(n_blocks_q4, sizeof(BlockQ4_0)); + if (!q4_data) { free(f32_data); goto write_fail; } float tensor_error = 0.0f; 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) @@ -5222,10 +5048,13 @@ 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); + if (tensor_error < 0.0f) { + free(q4_data); free(f32_data); goto write_fail; + } - fwrite(q4_data, sizeof(BlockQ4_0), n_blocks_q4, fp); + if (fwrite(q4_data, sizeof(BlockQ4_0), n_blocks_q4, fp) != (size_t)n_blocks_q4) { free(q4_data); free(f32_data); goto write_fail; } float rmse = sqrtf(tensor_error / (float)ti->n_elements); @@ -5314,14 +5143,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 = "★★★★"; } @@ -5388,425 +5217,18 @@ static int write_gguf(const char *output_path, const STMultiFile *mf, free(tensor_hf_names); 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, - int h_diag, 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); - } - if (h_diag) { - for (int64_t i = 0; i < n_in; i++) Gd[i] = H[i * n_in + i] * E[i]; - } else { - 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 (h_diag) { - /* diagonal H: only Gd[j] depends on E[j]. O(1). */ - Gd[j] = Hjj * E[j]; - } else { - const float *Hj = H + j * n_in; - for (int64_t k = 0; k < n_in; k++) Gd[k] -= bestd * Hj[k]; - } - } - } - (void)pend; - } - - /* ── 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]; - /* del(tl,tm) is AFFINE in the candidate: - * del = a*u - b*e, a = tl-Ls, b = tm-Lm - * u[k] = d*q[k], e[k] = dmin (constant vector) - * so the gain expands with five precomputed scalars - * gain = 2(a*Lu - b*Le) - (a^2*Quu - 2ab*Que + b^2*Qee) - * and each of the 256 candidates costs O(1) instead of the - * O(256) 16x16 quadratic form. This inner loop dominated the - * whole routine: ~65k ops per sub-block became ~3k. */ - float u[16], del[16], bestdel[16]; - for (int k = 0; k < 16; k++) u[k] = dd * (float)q[c0 + k]; - const float ev = mm; - - double Quu = 0, Que = 0, Qee = 0, Lu = 0, Le = 0; - if (h_diag) { - for (int k = 0; k < 16; k++) { - double h = H[(c0 + k) * n_in + c0 + k]; - Quu += h * u[k] * u[k]; - Que += h * u[k] * ev; - Qee += h * ev * ev; - } - } else { - for (int a = 0; a < 16; a++) { - const float *Ha = H + (c0 + a) * n_in + c0; - double su = 0, se = 0; - for (int c = 0; c < 16; c++) { su += Ha[c] * u[c]; se += Ha[c] * ev; } - Quu += u[a] * su; Que += u[a] * se; Qee += ev * se; - } - } - for (int k = 0; k < 16; k++) { Lu += u[k] * Gd[c0 + k]; Le += ev * Gd[c0 + k]; } - - float bestg = 0.0f; int bl = Ls[sj], bm = Lm[sj]; int found = 0; - double ba = 0.0, bb2 = 0.0; - for (int tl = 0; tl <= 15; tl++) { - double a = (double)tl - (double)Ls[sj]; - for (int tm = 0; tm <= 15; tm++) { - double b = (double)tm - (double)Lm[sj]; - double gain = 2.0 * (a * Lu - b * Le) - - (a*a*Quu - 2.0*a*b*Que + b*b*Qee); - if (gain > bestg + 1e-20) { - bestg = (float)gain; bl = tl; bm = tm; found = 1; - ba = a; bb2 = b; - } - } - } - if (found) - for (int k = 0; k < 16; k++) - bestdel[k] = (float)(ba * u[k] - bb2 * ev); - (void)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]; - if (h_diag) { - for (int k = 0; k < 16; k++) - Gd[c0 + k] = H[(c0 + k) * n_in + c0 + k] * E[c0 + k]; - } else { - 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]; - } - if (h_diag) { - for (int i = 0; i < QK_K; i++) { - double h = H[(c0 + i) * n_in + c0 + i]; - Saa += h*av[i]*av[i]; Sab += h*av[i]*bvv[i]; - Sbb += h*bvv[i]*bvv[i]; - Sat += h*av[i]*tv[i]; Sbt += h*bvv[i]*tv[i]; - } - } else { - 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)); - /* Accept test. E changes only inside this block, by dl = en - eo, - * so the exact objective change is - * dL = 2*dl^T Gd + dl^T H_blk dl - * -- the SAME gain form used by moves A and B. The previous - * version compared eo^T H_blk eo against en^T H_blk en, which - * silently dropped the cross-block coupling carried in Gd and so - * could accept a move that raised the true objective. */ - float eo[QK_K], en[QK_K], dl[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]); - dl[i] = en[i] - eo[i]; - } - double lin = 0.0, quad = 0.0; - for (int i = 0; i < QK_K; i++) lin += (double)dl[i] * Gd[c0 + i]; - if (h_diag) { - for (int i = 0; i < QK_K; i++) - quad += (double)dl[i]*dl[i]*H[(c0 + i) * n_in + c0 + i]; - } else { - for (int i = 0; i < QK_K; i++) { - if (dl[i] == 0.0f) continue; - const float *Hi = H + (c0 + i) * n_in + c0; - double sn = 0.0; - for (int k = 0; k < QK_K; k++) sn += (double)Hi[k] * dl[k]; - quad += (double)dl[i] * sn; - } - } - if (2.0 * lin + quad < 0.0) { /* strict decrease */ - dv[b] = dtry; mv[b] = mtry; - for (int i = 0; i < QK_K; i++) E[c0 + i] = en[i]; - /* Banded gradient update: Gd += H[:, blk] @ dl. - * Was a full O(n_in^2) recompute -- 10x waste at n_in=2560, - * 24x at 6144, and it dominated the whole runtime. */ - if (h_diag) { - for (int i = 0; i < QK_K; i++) - Gd[c0 + i] = H[(c0 + i) * n_in + c0 + i] * E[c0 + i]; - } else { - for (int64_t i = 0; i < n_in; i++) { - const float *Hi = H + i * n_in + c0; - double sm2 = 0.0; - for (int k = 0; k < QK_K; k++) sm2 += (double)Hi[k] * dl[k]; - Gd[i] += (float)sm2; - } - } - 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; - - /* Detect a diagonal Hessian once. An imatrix-derived H is diagonal, and - * then Gd = diag(H) .* E has NO coupling: every gradient update collapses - * from O(n_in) or O(n_in^2) to O(1). Same output, orders of magnitude - * less work -- worth the one-time O(n_in^2) scan. */ - int h_diag = 1; - for (int64_t i = 0; i < n_in && h_diag; i++) { - const float *Hi = H + i * n_in; - for (int64_t k = 0; k < n_in; k++) - if (k != i && Hi[k] != 0.0f) { h_diag = 0; break; } - } - if (verbose) - printf(" [analog] Hessian is %s -> %s path\n", - h_diag ? "DIAGONAL" : "dense", h_diag ? "fast" : "full"); - - #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, h_diag, 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"); - } +write_fail: + fclose(fp); + remove(output_path); + free(include_list); + free(gguf_names); + free(tensor_types); + free(tensor_sizes); + free(tensor_offsets); + free(tensor_src_idx); + free(tensor_hf_names); + return -1; } /* ═══════════════════════════════════════════════════════════════════════════ @@ -5840,66 +5262,30 @@ 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. @@ -5912,15 +5298,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; } @@ -6157,7 +5543,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); @@ -6168,4 +5554,4 @@ int main(int argc, char **argv) st_multi_close(mf); return result; } -#endif /* HEXSTATE_LIBRARY */ +#endif /* HEXSTATE_LIBRARY */ \ No newline at end of file