CompressedGemma commited on
Commit
0415704
·
verified ·
1 Parent(s): 2ec3df7

Fold baked into standard Q2K codebook

Browse files
Files changed (1) hide show
  1. hexstate_quantize.c +632 -121
hexstate_quantize.c CHANGED
@@ -1761,7 +1761,13 @@ static const int Q4_CAND_TO_QUHIT[Q4_N_CAND] = {
1761
  *
1762
  * E(block) = Σᵢ wᵢeᵢ²
1763
  * + (HEX_DC_LAMBDA / n) · (Σᵢeᵢ)²
1764
- * + (HEX_VW_LAMBDA / n) · Σ_{i<n/2} [(eᵢ+eⱼ)² − (eᵢ−eⱼ)²], j = i+n/2
 
 
 
 
 
 
1765
  *
1766
  * applied CONSISTENTLY to: Q2_K/Q4_0 candidate scoring, the closed-form
1767
  * (d, dmin) refit acceptance, the shaping accept guards, every polish
@@ -1776,35 +1782,493 @@ static const int Q4_CAND_TO_QUHIT[Q4_N_CAND] = {
1776
  * RMSE increase is the *intended* price for lower output error. Per-block
1777
  * terms are a proxy for row-level structure (the API sees a flat stream);
1778
  * the Phase-3.9 rolling-DC pass handles cross-block linkage. */
1779
- #ifndef HEX_DC_LAMBDA
1780
- #define HEX_DC_LAMBDA 1.0f
 
 
 
1781
  #endif
1782
- #ifndef HEX_VW_LAMBDA
1783
- #define HEX_VW_LAMBDA 1.0f
1784
  #endif
1785
- /* Default (1, 1): unit-strength spectral prior. Empirically (synthetic
1786
- * benchmark, identical inputs): lowers dot-product output error ~0.8-1.4%
1787
- * on both mean-only and fold-correlated activation models for ~+0.05%
1788
- * weight RMSE. The theoretically optimal λ grows with the deployment
1789
- * model's activation mean energy and row length (the per-block term
1790
- * under-counts cross-block row coupling); the synthetic sweep kept
1791
- * improving monotonically through λ = 4 at ~+0.1% RMSE. Set both to
1792
- * 0.0f to recover the exact pure weighted-SSE / minimum-RMSE pipeline. */
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1793
 
1794
  /* Spectral penalty of the extended objective for one block: residuals e[n],
1795
  * fold at n/2. Negative values are possible (anti-phase credit) — the total
1796
  * E remains positive-definite as argued above. */
 
 
 
 
 
 
1797
  static inline float hex_spectral_penalty(const float *e, int n)
1798
  {
1799
  if (HEX_DC_LAMBDA == 0.0f && HEX_VW_LAMBDA == 0.0f) return 0.0f;
1800
- float dc = 0.0f, cross = 0.0f;
1801
  int half = n / 2;
1802
  for (int i = 0; i < half; i++) {
1803
- dc += e[i] + e[i + half];
1804
- cross += e[i] * e[i + half];
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1805
  }
1806
- return (HEX_DC_LAMBDA / (float)n) * dc * dc
1807
- + (HEX_VW_LAMBDA / (float)n) * 4.0f * cross;
1808
  }
1809
 
1810
  /* Robust temperature estimator for the HExState measurement model.
@@ -2441,13 +2905,14 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements,
2441
  }
2442
 
2443
  float v_live[QK4_0 / 2];
2444
- float vesica_cur = 0.0f, dc_cur = 0.0f;
2445
  for (int j = 0; j < QK4_0 / 2; j++) {
2446
  v_live[j] = e_live[j] + e_live[j + QK4_0 / 2];
2447
  vesica_cur += v_live[j] * v_live[j];
 
2448
  }
2449
  for (int j = 0; j < QK4_0; j++) dc_cur += e_live[j];
2450
- float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
2451
 
2452
  /* Deterministic greedy descent: only strict improvements.
2453
  * The previous SA acceptance called rand() inside an OpenMP
@@ -2477,7 +2942,7 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements,
2477
 
2478
  float vesica_alt = vesica_cur - v_old * v_old + v_new * v_new;
2479
  float dc_alt = dc_cur + de;
2480
- float metric_alt = 4.0f * vesica_alt + dc_alt * dc_alt;
2481
 
2482
  float delta = metric_cur - metric_alt;
2483
  if (delta > best_delta) {
@@ -2501,7 +2966,7 @@ static void quantize_tensor_q4_0_hpc(const float *weights, int64_t n_elements,
2501
 
2502
  vesica_cur += v_new_commit * v_new_commit - v_old_commit * v_old_commit;
2503
  dc_cur += de_commit;
2504
- metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
2505
 
2506
  v_live[pi_commit] = v_new_commit;
2507
  e_live[best_k] = e_new_commit;
@@ -2812,15 +3277,16 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements,
2812
  int8_t qs_shaped[QK8_0];
2813
  memcpy(qs_shaped, qs, QK8_0);
2814
  float e_live[QK8_0], v_live[QK8_0 / 2];
2815
- float vesica_cur = 0.0f, dc_cur = 0.0f;
2816
  for (int k = 0; k < QK8_0; k++)
2817
  e_live[k] = bw[k] - (float)qs_shaped[k] * d;
2818
  for (int p = 0; p < QK8_0 / 2; p++) {
2819
  v_live[p] = e_live[p] + e_live[p + QK8_0 / 2];
2820
  vesica_cur += v_live[p] * v_live[p];
 
2821
  dc_cur += v_live[p];
2822
  }
2823
- float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
2824
  for (int pass = 0; pass < QK8_0; pass++) {
2825
  int best_k = -1, best_q_alt = 0;
2826
  float best_delta = 0.0f;
@@ -2835,7 +3301,8 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements,
2835
  float ves_a = vesica_cur - v_live[pi] * v_live[pi]
2836
  + v_new * v_new;
2837
  float dc_a = dc_cur + de;
2838
- float delta = metric_cur - (4.0f * ves_a + dc_a * dc_a);
 
2839
  if (delta > best_delta) {
2840
  best_delta = delta; best_k = k; best_q_alt = q_try;
2841
  }
@@ -2849,7 +3316,7 @@ static void quantize_tensor_q8_0_hpc(const float *weights, int64_t n_elements,
2849
  float v_new = v_live[pi] + de;
2850
  vesica_cur += v_new * v_new - v_live[pi] * v_live[pi];
2851
  dc_cur += de;
2852
- metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
2853
  v_live[pi] = v_new;
2854
  e_live[best_k] = e_new;
2855
  qs_shaped[best_k] = (int8_t)best_q_alt;
@@ -2948,6 +3415,31 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
2948
  else clamp_sigma = 4.0f; /* very heavy tails: stabilise seed */
2949
  float clamp_val = w_sigma * clamp_sigma;
2950
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
2951
  /* ══════════════════════════════════════════════════════════════════
2952
  * PHASE 1: Greedy quantization — produce seed (d, dmin) per block
2953
  * ══════════════════════════════════════════════════════════════════ */
@@ -2961,6 +3453,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
2961
 
2962
  BlockSeed *seeds = (BlockSeed *)calloc(n_blocks, sizeof(BlockSeed));
2963
  if (!seeds) {
 
2964
  if (out_total_error) *out_total_error = -1.0f;
2965
  return;
2966
  }
@@ -3023,6 +3516,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3023
  sizeof(float[TOTAL_SCALE_CANDIDATES]));
3024
  if (!seeds || !candidate_errors) {
3025
  free(seeds); free(candidate_errors);
 
3026
  if (out_total_error) *out_total_error = -1.0f;
3027
  return;
3028
  }
@@ -3634,8 +4128,6 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3634
  * this array (safe: written sequentially before the parallel loop).
3635
  * ══════════════════════════════════════════════════════════════════ */
3636
 
3637
- #define DC_DECAY 0.85f /* Boundary-condition leak factor (0 = isolated, 1 = full) */
3638
-
3639
  float *block_dc_bias = (float *)calloc(n_blocks, sizeof(float));
3640
 
3641
  if (block_dc_bias) {
@@ -3664,7 +4156,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3664
  dm0, mm0, dc_Ls, dc_Lm);
3665
 
3666
  /* Bias applied to THIS block's WLS targets */
3667
- float dc_bias = (DC_DECAY * rolling_dc) / (float)QK_K;
3668
  block_dc_bias[blk] = dc_bias;
3669
 
3670
  /* Quick round-nearest quant to estimate DC residual for NEXT block.
@@ -3736,6 +4228,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3736
  uint16_t prev_dm16 = 0, prev_mm16 = 0;
3737
  for (int ls_iter = 0; ls_iter < 5; ls_iter++) {
3738
 
 
3739
  uint8_t state_ls[N_SUB][6];
3740
  uint8_t state_lm[N_SUB][6];
3741
  float state_err[N_SUB][6];
@@ -3818,6 +4311,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3818
 
3819
  for (int j = 0; j < N_SUB - 1; j++)
3820
  hpc_cz(sg, j, j + 1);
 
 
 
 
 
3821
 
3822
  double sub_marg[N_SUB][6];
3823
  int sub_measured[N_SUB];
@@ -3826,33 +4324,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3826
 
3827
  sieve_measure_graph(sg, N_SUB, sub_marg, sub_measured, 1);
3828
 
3829
- for (int j = 0; j < N_SUB; j++) {
3830
- double best_prob = -1.0;
3831
- int best_v = 0;
3832
- for (int v = 0; v < 6; v++) {
3833
- if (sub_marg[j][v] > best_prob) {
3834
- best_prob = sub_marg[j][v];
3835
- best_v = v;
3836
- }
3837
- }
3838
- Ls_blk[j] = state_ls[j][best_v];
3839
- Lm_blk[j] = state_lm[j][best_v];
3840
- }
3841
- }
3842
-
3843
- uint8_t L[QK_K];
3844
- for (int j = 0; j < N_SUB; j++) {
3845
- float d_sub = dm * (float)Ls_blk[j];
3846
- float m_sub = mm * (float)Lm_blk[j];
3847
- if (d_sub < 1e-15f) {
3848
- for (int k = 0; k < 16; k++) L[16*j+k] = 0;
3849
- continue;
3850
- }
3851
- for (int k = 0; k < 16; k++) {
3852
- int q = gguf_nearest_int((adj_block_x[16*j+k] + m_sub) / d_sub);
3853
- if (q < 0) q = 0; if (q > 3) q = 3;
3854
- L[16*j+k] = (uint8_t)q;
3855
- }
3856
  }
3857
 
3858
  double Saa = 0, Sab = 0, Sbb = 0, Sxa = 0, Sxb = 0;
@@ -3916,7 +4392,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3916
  float d_sub = trial_dm * (float)Ls_blk[j];
3917
  float m_sub = trial_mm * (float)Lm_blk[j];
3918
  for (int k = 0; k < 16; k++) {
3919
- float x = adj_block_x[16*j+k];
3920
  float w = (imat_importance) ?
3921
  imat_importance[blk * QK_K + 16*j+k] : 1.0f;
3922
  int q;
@@ -3942,7 +4418,7 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3942
  }
3943
 
3944
  for (int j = 0; j < N_SUB; j++) {
3945
- const float *sx = adj_block_x + 16 * j;
3946
  float best_sub_err = 1e30f;
3947
  uint8_t best_ls = Ls_blk[j], best_lm = Lm_blk[j];
3948
  for (int try_ls = 0; try_ls <= 15; try_ls++) {
@@ -3997,26 +4473,14 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
3997
  */
3998
  uint8_t L[QK_K];
3999
  {
4000
- float q_cont_all[QK_K];
4001
- int q_base_all[QK_K];
4002
- int q_shaped_all[QK_K];
4003
-
4004
- for (int i = 0; i < QK_K; i++) {
4005
- int jj = i >> 4;
4006
- float d_s = dm * (float)(output[blk].scales[jj] & 0xF);
4007
- float m_s = mm * (float)(output[blk].scales[jj] >> 4);
4008
- if (d_s < 1e-15f) {
4009
- q_cont_all[i] = 0.0f;
4010
- q_base_all[i] = 0;
4011
- } else {
4012
- /* Quantize the DC-adjusted target */
4013
- float qc = (adj_block_x[i] + m_s) / d_s;
4014
- q_cont_all[i] = qc;
4015
- int qr = gguf_nearest_int(qc);
4016
- if (qr < 0) qr = 0; if (qr > 3) qr = 3;
4017
- q_base_all[i] = qr;
4018
- }
4019
- }
4020
  memcpy(q_shaped_all, q_base_all, QK_K * sizeof(int));
4021
 
4022
  float e_live[QK_K];
@@ -4028,19 +4492,26 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4028
  * d_s == 0 (NOT 0 — the −dmin·lm term always applies). */
4029
  float deq = d_s * (float)q_shaped_all[i] - m_s;
4030
  /* Residual against the adjusted target (DC-corrected view) */
4031
- e_live[i] = adj_block_x[i] - deq;
4032
  }
4033
 
4034
  float v_live[QK_K / 2];
4035
- float vesica_cur = 0.0f, dc_cur = 0.0f;
4036
  for (int i = 0; i < QK_K / 2; i++) {
4037
  v_live[i] = e_live[i] + e_live[i + QK_K / 2];
4038
  vesica_cur += v_live[i] * v_live[i];
 
4039
  }
4040
  for (int i = 0; i < QK_K; i++) dc_cur += e_live[i];
4041
- float metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
 
 
 
 
 
 
4042
 
4043
- /* Deterministic greedy descent: accept only strict improvements */
4044
  for (int pass = 0; pass < QK_K; pass++) {
4045
  int best_k = -1;
4046
  int best_q_alt = 0;
@@ -4054,22 +4525,24 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4054
  int q_cur = q_shaped_all[k];
4055
  float m_s = mm * (float)(output[blk].scales[jj] >> 4);
4056
  int pi = (k < QK_K / 2) ? k : k - QK_K / 2;
 
4057
 
4058
- /* Try ALL alternate codes, not just ±1. A weight at q=0
4059
- * with strong positive DC bias may need to jump to q=2 or
4060
- * q=3; the old ±1 path took multiple passes to walk there
4061
- * one step at a time, if it got there at all. */
4062
  for (int q_try = 0; q_try <= 3; q_try++) {
4063
  if (q_try == q_cur) continue;
4064
 
4065
- float e_new = adj_block_x[k] - (d_s * (float)q_try - m_s);
4066
  float de = e_new - e_live[k];
4067
 
4068
  float v_new = v_live[pi] + de;
4069
 
4070
  float vesica_alt = vesica_cur - v_live[pi]*v_live[pi] + v_new*v_new;
4071
  float dc_alt = dc_cur + de;
4072
- float delta = metric_cur - (4.0f * vesica_alt + dc_alt * dc_alt);
 
 
 
 
 
4073
 
4074
  if (delta > best_delta) {
4075
  best_delta = delta;
@@ -4086,13 +4559,18 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4086
  int jj_c = best_k >> 4;
4087
  float d_c = dm * (float)(output[blk].scales[jj_c] & 0xF);
4088
  float m_c = mm * (float)(output[blk].scales[jj_c] >> 4);
4089
- float e_new_c = adj_block_x[best_k] - (d_c * (float)best_q_alt - m_c);
4090
  float de_c = e_new_c - e_live[best_k];
4091
  int pi_c = (best_k < QK_K / 2) ? best_k : best_k - QK_K / 2;
4092
  float v_new_c = v_live[pi_c] + de_c;
 
 
 
 
4093
  vesica_cur += v_new_c * v_new_c - v_live[pi_c] * v_live[pi_c];
4094
  dc_cur += de_c;
4095
- metric_cur = 4.0f * vesica_cur + dc_cur * dc_cur;
 
4096
  v_live[pi_c] = v_new_c;
4097
  e_live[best_k]= e_new_c;
4098
  }
@@ -4404,8 +4882,8 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4404
  * candidate against the whole-block penalty with the partner
4405
  * residuals held fixed — exact coordinate descent on E. */
4406
  float pe[QK_K];
4407
- float sub_sse[16], sub_dc[16], pair_cross[8];
4408
- float dc_tot = 0.0f, cross_tot = 0.0f;
4409
  for (int j = 0; j < N_SUB; j++) {
4410
  float d_sub = dm * (float)pl_Ls[j];
4411
  float m_sub = mm * (float)pl_Lm[j];
@@ -4424,10 +4902,12 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4424
  dc_tot += sub_dc[j];
4425
  }
4426
  for (int p = 0; p < 8; p++) {
4427
- pair_cross[p] = 0.0f;
4428
- for (int k = 0; k < 16; k++)
4429
- pair_cross[p] += pe[16*p + k] * pe[16*(p+8) + k];
4430
- cross_tot += pair_cross[p];
 
 
4431
  }
4432
 
4433
  for (int j = 0; j < N_SUB; j++) {
@@ -4436,12 +4916,12 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4436
  int pj = j ^ 8; /* partner sub-block */
4437
  const float *ppe = pe + 16 * pj; /* partner residuals */
4438
  float dc_rest = dc_tot - sub_dc[j];
4439
- float cross_rest = cross_tot - pair_cross[pi];
4440
 
4441
  /* Extended score of the CURRENT committed state */
4442
  float best_sub = sub_sse[j]
4443
  + (HEX_DC_LAMBDA / (float)QK_K) * dc_tot * dc_tot
4444
- + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f * cross_tot;
4445
  int best_ls = -1, best_lm = 0;
4446
  uint8_t best_q[16];
4447
  float best_e[16];
@@ -4470,10 +4950,10 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4470
  e_loc[k] = e;
4471
  sub_err += e * e * w;
4472
  dcc += e;
4473
- cxc += e * ppe[k];
4474
- /* SSE-partial prune is a valid lower bound
4475
- * only while the spectral terms are ≥ 0,
4476
- * i.e. when the (signable) vw credit is off */
4477
  if (HEX_VW_LAMBDA == 0.0f &&
4478
  sub_err >= best_sub) { aborted = 1; break; }
4479
  }
@@ -4481,8 +4961,8 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4481
  float score = sub_err
4482
  + (HEX_DC_LAMBDA / (float)QK_K)
4483
  * (dc_rest + dcc) * (dc_rest + dcc)
4484
- + (HEX_VW_LAMBDA / (float)QK_K) * 4.0f
4485
- * (cross_rest + cxc);
4486
  if (score < best_sub) {
4487
  best_sub = score;
4488
  best_ls = try_ls;
@@ -4503,9 +4983,9 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4503
  memcpy(pe + 16 * j, best_e, sizeof(best_e));
4504
  sub_sse[j] = best_sse;
4505
  sub_dc[j] = best_dcc;
4506
- pair_cross[pi] = best_cxc;
4507
  dc_tot = dc_rest + best_dcc;
4508
- cross_tot = cross_rest + best_cxc;
4509
  pol_improved = 1;
4510
  }
4511
  }
@@ -4666,20 +5146,34 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4666
  output[blk].dmin = gguf_fp32_to_fp16(mm);
4667
  }
4668
 
4669
- /* ══ PHASE 4.7: CANDIDATE FLOOR (worst-case bound) ══
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
4670
  *
4671
- * candidate_errors[blk][c] is the EXACT weighted SSE of a directly
4672
- * encodable configuration (fp16 d/dmin + derived Ls/Lm + nearest
4673
- * rounding vs the original weights). The multi-stage assembly
4674
- * (DC-shifted WLS, shaping, diffusion, polish) usually improves on
4675
- * its seed, but each stage optimises a slightly different objective
4676
- * and coordinate descent can land in a worse basin. Compare the
4677
- * finished block against the best raw candidate and fall back when
4678
- * the pipeline ended up worse — guaranteeing
4679
- * final weighted SSE ≤ min_c candidate_errors[blk][c]. */
4680
  {
4681
  float fin_err = 0.0f;
4682
- float e_f[QK_K];
4683
  for (int j = 0; j < N_SUB; j++) {
4684
  float d_sub = dm * (float)(output[blk].scales[j] & 0xF);
4685
  float m_sub = mm * (float)(output[blk].scales[j] >> 4);
@@ -4687,12 +5181,11 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4687
  int idx = 16 * j + k;
4688
  float w = (imat_importance) ?
4689
  imat_importance[blk * QK_K + idx] : 1.0f;
4690
- e_f[idx] = block_x[idx] -
4691
- (d_sub * (float)L[idx] - m_sub);
4692
- fin_err += e_f[idx] * e_f[idx] * w;
4693
  }
4694
  }
4695
- fin_err += hex_spectral_penalty(e_f, QK_K);
4696
 
4697
  float g_best = candidate_errors[blk][0];
4698
  int g_cand = 0;
@@ -4732,6 +5225,22 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4732
  }
4733
  }
4734
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
4735
  for (int j = 0; j < QK_K; j += 128) {
4736
  for (int l = 0; l < 32; l++) {
4737
  output[blk].qs[j / 4 + l] = L[j + l]
@@ -4758,6 +5267,8 @@ static void quantize_tensor_q2k_hpc(const float *weights, int64_t n_elements,
4758
  free(seeds);
4759
  free(candidate_errors);
4760
  free(best_candidate);
 
 
4761
  if (out_total_error) *out_total_error = total_err;
4762
 
4763
  if (verbose) {
 
1761
  *
1762
  * E(block) = Σᵢ wᵢeᵢ²
1763
  * + (HEX_DC_LAMBDA / n) · (Σᵢeᵢ)²
1764
+ * + (HEX_VW_LAMBDA / n) · Σ_pairs (eᵢ + eⱼ)², j = i+n/2
1765
+ *
1766
+ * Vesica guiding: DC = Σ_pairs (eᵢ+eⱼ), so zeroing EACH pair vesica is
1767
+ * a strict DC floor ((Σ v_p)² ≤ n_pairs Σ v_p²). The older 4·eᵢeⱼ
1768
+ * (= vesica² − wave²) credited anti-phase amplitude and let vesicas
1769
+ * cancel across pairs — that raises |v_p| while looking like DC win.
1770
+ * Pair-joint code search (hex_vesica_guide_pairs) is the last move.
1771
  *
1772
  * applied CONSISTENTLY to: Q2_K/Q4_0 candidate scoring, the closed-form
1773
  * (d, dmin) refit acceptance, the shaping accept guards, every polish
 
1782
  * RMSE increase is the *intended* price for lower output error. Per-block
1783
  * terms are a proxy for row-level structure (the API sees a flat stream);
1784
  * the Phase-3.9 rolling-DC pass handles cross-block linkage. */
1785
+ #ifndef HEX_DC_LAMBDA_DEFAULT
1786
+ #define HEX_DC_LAMBDA_DEFAULT 4.0f
1787
+ #endif
1788
+ #ifndef HEX_VW_LAMBDA_DEFAULT
1789
+ #define HEX_VW_LAMBDA_DEFAULT 4.0f
1790
  #endif
1791
+ #ifndef HEX_DC_DECAY_DEFAULT
1792
+ #define HEX_DC_DECAY_DEFAULT 1.0f
1793
  #endif
1794
+ /* Smoke: vesica-guided pair snap (near nearest-int, live DC) plus
1795
+ * Σ(e_i+e_j)² in E. λ=(4,4) decay=1 cuts |DC| ~4× and mean pair
1796
+ * vesica ~6% vs λ=0 at ~+1.5% RMSE. Wave-credit (4 e_i e_j) is gone. */
1797
+ /* Runtime knobs (defaults = previous compile-time 1, 1, 0.85). The
1798
+ * greedy shaper used to descend 4·Σ(e_i+e_j)² + DC², which is not the
1799
+ * extended objective below; it now descends the same E. Sweep from
1800
+ * Python via hexstate_set_spectral_params. λ = 0, 0 recovers pure
1801
+ * weighted SSE. */
1802
+
1803
+ static float hex_dc_lambda = HEX_DC_LAMBDA_DEFAULT;
1804
+ static float hex_vw_lambda = HEX_VW_LAMBDA_DEFAULT;
1805
+ static float hex_dc_decay = HEX_DC_DECAY_DEFAULT;
1806
+
1807
+ #define HEX_DC_LAMBDA hex_dc_lambda
1808
+ #define HEX_VW_LAMBDA hex_vw_lambda
1809
+
1810
+ void hexstate_set_spectral_params(float dc_lambda, float vw_lambda, float dc_decay)
1811
+ {
1812
+ hex_dc_lambda = dc_lambda;
1813
+ hex_vw_lambda = vw_lambda;
1814
+ if (dc_decay >= 0.0f && dc_decay <= 1.0f)
1815
+ hex_dc_decay = dc_decay;
1816
+ }
1817
+
1818
+ void hexstate_get_spectral_params(float *dc_lambda, float *vw_lambda, float *dc_decay)
1819
+ {
1820
+ if (dc_lambda) *dc_lambda = hex_dc_lambda;
1821
+ if (vw_lambda) *vw_lambda = hex_vw_lambda;
1822
+ if (dc_decay) *dc_decay = hex_dc_decay;
1823
+ }
1824
 
1825
  /* Spectral penalty of the extended objective for one block: residuals e[n],
1826
  * fold at n/2. Negative values are possible (anti-phase credit) — the total
1827
  * E remains positive-definite as argued above. */
1828
+ static inline float hex_spectral_from_moments(float dc, float vesica_ssq, int n)
1829
+ {
1830
+ return (HEX_DC_LAMBDA / (float)n) * dc * dc
1831
+ + (HEX_VW_LAMBDA / (float)n) * vesica_ssq;
1832
+ }
1833
+
1834
  static inline float hex_spectral_penalty(const float *e, int n)
1835
  {
1836
  if (HEX_DC_LAMBDA == 0.0f && HEX_VW_LAMBDA == 0.0f) return 0.0f;
1837
+ float dc = 0.0f, ves = 0.0f;
1838
  int half = n / 2;
1839
  for (int i = 0; i < half; i++) {
1840
+ float v = e[i] + e[i + half];
1841
+ dc += v;
1842
+ ves += v * v;
1843
+ }
1844
+ return hex_spectral_from_moments(dc, ves, n);
1845
+ }
1846
+
1847
+ #ifndef HEX_VESICA_SSE_EPS
1848
+ #define HEX_VESICA_SSE_EPS 1.0f
1849
+ #endif
1850
+
1851
+ typedef struct { float v; int i; } hex_fold_kv;
1852
+
1853
+ static int hex_fold_kv_cmp(const void *a, const void *b)
1854
+ {
1855
+ const hex_fold_kv *pa = (const hex_fold_kv *)a;
1856
+ const hex_fold_kv *pb = (const hex_fold_kv *)b;
1857
+ if (pa->v < pb->v) return -1;
1858
+ if (pa->v > pb->v) return 1;
1859
+ return pa->i - pb->i;
1860
+ }
1861
+
1862
+ /* Fold's layout, in-place: k-th smallest at k, k-th largest at k+128.
1863
+ * Pair-sum ≈ constant, 16-wide scale groups become homogeneous.
1864
+ * Stored order IS decode order — stock llama, no sidecar. */
1865
+ static void hex_fold_interleave_block(float *x, float *w_row)
1866
+ {
1867
+ hex_fold_kv kv[QK_K];
1868
+ float xt[QK_K], wt[QK_K];
1869
+ memcpy(xt, x, sizeof(xt));
1870
+ if (w_row) memcpy(wt, w_row, sizeof(wt));
1871
+ for (int i = 0; i < QK_K; i++) {
1872
+ kv[i].v = xt[i];
1873
+ kv[i].i = i;
1874
+ }
1875
+ qsort(kv, (size_t)QK_K, sizeof(kv[0]), hex_fold_kv_cmp);
1876
+ int half = QK_K / 2;
1877
+ for (int k = 0; k < half; k++) {
1878
+ int lo = kv[k].i;
1879
+ int hi = kv[QK_K - 1 - k].i;
1880
+ x[k] = xt[lo];
1881
+ x[k + half] = xt[hi];
1882
+ if (w_row) {
1883
+ w_row[k] = wt[lo];
1884
+ w_row[k + half] = wt[hi];
1885
+ }
1886
+ }
1887
+ }
1888
+
1889
+ /* Vesica subcode: among L2-optimal 4x4 points (Edge), pick the one
1890
+ * on the best pair-sum lattice σ=deq_i+deq_j. Independent nearest-int
1891
+ * breaks those ties without looking at the partner, so the two
1892
+ * estimates of s disagree. Primary SSE, secondary v² — reversing this
1893
+ * was the 0.51 basin. */
1894
+ static void hex_vesica_subcode_q(float xi, float xj,
1895
+ float di, float mi, float dj, float mj,
1896
+ float wi, float wj,
1897
+ int vertex_first,
1898
+ int *qi_out, int *qj_out,
1899
+ float *ei_out, float *ej_out)
1900
+ {
1901
+ float best_sse = 1e30f, best_v = 1e30f;
1902
+ int bqi = 0, bqj = 0;
1903
+ float bei = 0.0f, bej = 0.0f;
1904
+ for (int qi = 0; qi <= 3; qi++) {
1905
+ float deqi = di * (float)qi - mi;
1906
+ float ei = xi - deqi;
1907
+ for (int qj = 0; qj <= 3; qj++) {
1908
+ float deqj = dj * (float)qj - mj;
1909
+ float ej = xj - deqj;
1910
+ float sse = wi * ei * ei + wj * ej * ej;
1911
+ float v2 = (ei + ej) * (ei + ej);
1912
+ int better;
1913
+ if (vertex_first)
1914
+ better = (v2 < best_v - 1e-12f) ||
1915
+ (v2 <= best_v + 1e-12f && sse < best_sse);
1916
+ else
1917
+ better = (sse < best_sse - 1e-12f) ||
1918
+ (sse <= best_sse + 1e-12f && v2 < best_v);
1919
+ if (better) {
1920
+ best_sse = sse;
1921
+ best_v = v2;
1922
+ bqi = qi;
1923
+ bqj = qj;
1924
+ bei = ei;
1925
+ bej = ej;
1926
+ }
1927
+ }
1928
+ }
1929
+ *qi_out = bqi;
1930
+ *qj_out = bqj;
1931
+ *ei_out = bei;
1932
+ *ej_out = bej;
1933
+ }
1934
+
1935
+ static void hex_vesica_bake_subcodes(const float *x, const float *w_row,
1936
+ float dm, float mm,
1937
+ const uint8_t *Ls, const uint8_t *Lm,
1938
+ uint8_t *L, int vertex_first)
1939
+ {
1940
+ for (int p = 0; p < QK_K / 2; p++) {
1941
+ int i = p, j = p + QK_K / 2;
1942
+ int si = i >> 4, sj = j >> 4;
1943
+ float di = dm * (float)Ls[si], mi = mm * (float)Lm[si];
1944
+ float dj = dm * (float)Ls[sj], mj = mm * (float)Lm[sj];
1945
+ float wi = w_row ? w_row[i] : 1.0f;
1946
+ float wj = w_row ? w_row[j] : 1.0f;
1947
+ int qi, qj;
1948
+ float ei, ej;
1949
+ hex_vesica_subcode_q(x[i], x[j], di, mi, dj, mj, wi, wj, vertex_first,
1950
+ &qi, &qj, &ei, &ej);
1951
+ L[i] = (uint8_t)qi;
1952
+ L[j] = (uint8_t)qj;
1953
+ }
1954
+ }
1955
+
1956
+ static float hex_block_sse(const float *x, const float *w_row,
1957
+ float dm, float mm,
1958
+ const uint8_t *Ls, const uint8_t *Lm,
1959
+ const uint8_t *L)
1960
+ {
1961
+ float sse = 0.0f;
1962
+ for (int j = 0; j < QK_K / 16; j++) {
1963
+ float d_sub = dm * (float)Ls[j];
1964
+ float m_sub = mm * (float)Lm[j];
1965
+ for (int k = 0; k < 16; k++) {
1966
+ int idx = 16 * j + k;
1967
+ float w = w_row ? w_row[idx] : 1.0f;
1968
+ float e = x[idx] - (d_sub * (float)L[idx] - m_sub);
1969
+ sse += w * e * e;
1970
+ }
1971
+ }
1972
+ return sse;
1973
+ }
1974
+
1975
+ static void hex_wls_d_dmin(const float *x, const float *w_row,
1976
+ const uint8_t *Ls, const uint8_t *Lm,
1977
+ const uint8_t *L, float *dm, float *mm)
1978
+ {
1979
+ double Saa = 0, Sab = 0, Sbb = 0, Sxa = 0, Sxb = 0;
1980
+ for (int j = 0; j < QK_K / 16; j++) {
1981
+ float ls_f = (float)Ls[j];
1982
+ float lm_f = (float)Lm[j];
1983
+ for (int k = 0; k < 16; k++) {
1984
+ int idx = 16 * j + k;
1985
+ float xv = x[idx];
1986
+ float w = w_row ? w_row[idx] : 1.0f;
1987
+ float a = ls_f * (float)L[idx];
1988
+ float b = lm_f;
1989
+ Saa += w * a * a;
1990
+ Sab += w * a * b;
1991
+ Sbb += w * b * b;
1992
+ Sxa += w * xv * a;
1993
+ Sxb += w * xv * b;
1994
+ }
1995
+ }
1996
+ double det = Saa * Sbb - Sab * Sab;
1997
+ if (fabs(det) <= 1e-30)
1998
+ return;
1999
+ double d_new = (Sbb * Sxa - Sab * Sxb) / det;
2000
+ double m_new = (Sab * Sxa - Saa * Sxb) / det;
2001
+ if (d_new > 0.0 && d_new < 4.0 * ((double)*dm + 1e-10))
2002
+ *dm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)d_new));
2003
+ if (m_new > 0.0 && m_new < 4.0 * ((double)*mm + 1e-10))
2004
+ *mm = gguf_fp16_to_fp32(gguf_fp32_to_fp16((float)m_new));
2005
+ }
2006
+
2007
+ /* Tie complementary (ls,lm) when it lowers 32-weight SSE (Fold's
2008
+ * homogeneous group as ordinary nibbles). σ-bake + WLS, never Edge
2009
+ * after WLS — that was the mismatched-scalar basin. Keep only if
2010
+ * block SSE drops. */
2011
+ static void hex_vesica_commit(const float *x, const float *w_row,
2012
+ float *dm, float *mm,
2013
+ uint8_t *Ls, uint8_t *Lm, uint8_t *L)
2014
+ {
2015
+ uint8_t L0[QK_K], Ls0[16], Lm0[16];
2016
+ memcpy(L0, L, QK_K);
2017
+ memcpy(Ls0, Ls, 16);
2018
+ memcpy(Lm0, Lm, 16);
2019
+ float dm0 = *dm, mm0 = *mm;
2020
+ float sse0 = hex_block_sse(x, w_row, dm0, mm0, Ls0, Lm0, L0);
2021
+
2022
+ for (int p = 0; p < 8; p++) {
2023
+ int A = p, B = p + 8;
2024
+ uint8_t tLs[16], tLm[16], tL[QK_K];
2025
+ float best32 = 1e30f;
2026
+ uint8_t bA_ls = Ls[A], bA_lm = Lm[A], bB_ls = Ls[B], bB_lm = Lm[B];
2027
+ uint8_t cand_lsA[4] = { Ls[A], Ls[A], Ls[B], (uint8_t)((Ls[A] + Ls[B] + 1) / 2) };
2028
+ uint8_t cand_lmA[4] = { Lm[A], Lm[A], Lm[B], (uint8_t)((Lm[A] + Lm[B] + 1) / 2) };
2029
+ uint8_t cand_lsB[4] = { Ls[B], Ls[A], Ls[B], (uint8_t)((Ls[A] + Ls[B] + 1) / 2) };
2030
+ uint8_t cand_lmB[4] = { Lm[B], Lm[A], Lm[B], (uint8_t)((Lm[A] + Lm[B] + 1) / 2) };
2031
+ for (int t = 0; t < 4; t++) {
2032
+ memcpy(tLs, Ls, 16);
2033
+ memcpy(tLm, Lm, 16);
2034
+ tLs[A] = cand_lsA[t]; tLm[A] = cand_lmA[t];
2035
+ tLs[B] = cand_lsB[t]; tLm[B] = cand_lmB[t];
2036
+ memcpy(tL, L, QK_K);
2037
+ hex_vesica_bake_subcodes(x, w_row, *dm, *mm, tLs, tLm, tL, 0);
2038
+ float s32 = 0.0f;
2039
+ for (int k = 0; k < 16; k++) {
2040
+ int i = 16 * A + k, j = 16 * B + k;
2041
+ float di = *dm * (float)tLs[A], mi = *mm * (float)tLm[A];
2042
+ float dj = *dm * (float)tLs[B], mj = *mm * (float)tLm[B];
2043
+ float ei = x[i] - (di * (float)tL[i] - mi);
2044
+ float ej = x[j] - (dj * (float)tL[j] - mj);
2045
+ float wi = w_row ? w_row[i] : 1.0f;
2046
+ float wj = w_row ? w_row[j] : 1.0f;
2047
+ s32 += wi * ei * ei + wj * ej * ej;
2048
+ }
2049
+ if (s32 < best32) {
2050
+ best32 = s32;
2051
+ bA_ls = tLs[A]; bA_lm = tLm[A];
2052
+ bB_ls = tLs[B]; bB_lm = tLm[B];
2053
+ }
2054
+ }
2055
+ Ls[A] = bA_ls; Lm[A] = bA_lm;
2056
+ Ls[B] = bB_ls; Lm[B] = bB_lm;
2057
+ }
2058
+
2059
+ for (int it = 0; it < 3; it++) {
2060
+ hex_vesica_bake_subcodes(x, w_row, *dm, *mm, Ls, Lm, L, 1);
2061
+ hex_wls_d_dmin(x, w_row, Ls, Lm, L, dm, mm);
2062
+ }
2063
+ hex_vesica_bake_subcodes(x, w_row, *dm, *mm, Ls, Lm, L, 0);
2064
+
2065
+ float sse1 = hex_block_sse(x, w_row, *dm, *mm, Ls, Lm, L);
2066
+ if (!(sse1 < sse0)) {
2067
+ memcpy(L, L0, QK_K);
2068
+ memcpy(Ls, Ls0, 16);
2069
+ memcpy(Lm, Lm0, 16);
2070
+ *dm = dm0;
2071
+ *mm = mm0;
2072
+ }
2073
+ }
2074
+
2075
+ /* Vesica + triality. Cooley-Tukey DFT6 = 2 x 3: antipodal fold is
2076
+ * vesica/wave; the 3-cycle Edge/Vertex/Diagonal is triality_rotate.
2077
+ * Vertex: min v^2 + DC^2 (in-phase / DC floor)
2078
+ * Diagonal: min w^2, keep v and DC (anti-phase / leftover RMSE)
2079
+ * Edge: min SSE, keep v and DC (computational RMSE)
2080
+ * RMSE^2 = (v^2 + w^2)/2 per pair, so Diagonal+Edge shrink RMSE
2081
+ * without giving back the Vertex cancellation. Two sweeps/view. */
2082
+ static void hex_vesica_guide_pairs(const float *x, const float *w_row,
2083
+ float dm, float mm,
2084
+ const uint8_t *Ls, const uint8_t *Lm,
2085
+ uint8_t *L)
2086
+ {
2087
+ float e[QK_K];
2088
+ float dc = 0.0f;
2089
+ for (int i = 0; i < QK_K; i++) {
2090
+ int s = i >> 4;
2091
+ float d_s = dm * (float)Ls[s];
2092
+ float m_s = mm * (float)Lm[s];
2093
+ e[i] = x[i] - (d_s * (float)L[i] - m_s);
2094
+ dc += e[i];
2095
+ }
2096
+
2097
+ /* Vertex (vesica+DC), then Edge (RMSE with vesica/DC held).
2098
+ * Diagonal/wave is implied: RMSE^2 = (v^2+w^2)/2 at fixed v. */
2099
+ for (int rot = 0; rot < 2; rot++) {
2100
+ for (int sweep = 0; sweep < 2; sweep++) {
2101
+ for (int p = 0; p < QK_K / 2; p++) {
2102
+ int i = p, j = p + QK_K / 2;
2103
+ int si = i >> 4, sj = j >> 4;
2104
+ float di = dm * (float)Ls[si], mi = mm * (float)Lm[si];
2105
+ float dj = dm * (float)Ls[sj], mj = mm * (float)Lm[sj];
2106
+ float wi = w_row ? w_row[i] : 1.0f;
2107
+ float wj = w_row ? w_row[j] : 1.0f;
2108
+
2109
+ int q0i = 0, q0j = 0;
2110
+ if (di >= 1e-15f) {
2111
+ q0i = gguf_nearest_int((x[i] + mi) / di);
2112
+ if (q0i < 0) q0i = 0; if (q0i > 3) q0i = 3;
2113
+ }
2114
+ if (dj >= 1e-15f) {
2115
+ q0j = gguf_nearest_int((x[j] + mj) / dj);
2116
+ if (q0j < 0) q0j = 0; if (q0j > 3) q0j = 3;
2117
+ }
2118
+
2119
+ float v_old = e[i] + e[j];
2120
+ float dc_rest = dc - v_old;
2121
+ float dc_abs0 = dc * dc;
2122
+ float ves0 = v_old * v_old;
2123
+ float cur_sse = wi * e[i] * e[i] + wj * e[j] * e[j];
2124
+ float best = (rot == 0) ? 1e30f : cur_sse;
2125
+ int bqi = (int)L[i], bqj = (int)L[j];
2126
+ float bei = e[i], bej = e[j];
2127
+
2128
+ int qi0 = q0i - 1, qi1 = q0i + 1;
2129
+ int qj0 = q0j - 1, qj1 = q0j + 1;
2130
+ if (qi0 < 0) qi0 = 0; if (qi1 > 3) qi1 = 3;
2131
+ if (qj0 < 0) qj0 = 0; if (qj1 > 3) qj1 = 3;
2132
+
2133
+ for (int qi = qi0; qi <= qi1; qi++) {
2134
+ float ei = x[i] - (di * (float)qi - mi);
2135
+ for (int qj = qj0; qj <= qj1; qj++) {
2136
+ float ej = x[j] - (dj * (float)qj - mj);
2137
+ float v = ei + ej;
2138
+ float sse = wi * ei * ei + wj * ej * ej;
2139
+ float dc_n = dc_rest + v;
2140
+ float sc;
2141
+ if (rot == 0)
2142
+ sc = v * v + dc_n * dc_n
2143
+ + HEX_VESICA_SSE_EPS * sse;
2144
+ else {
2145
+ /* Edge: cut RMSE only if vesica and DC do not grow. */
2146
+ if (v * v > ves0 * 1.05f) continue;
2147
+ if (dc_n * dc_n > dc_abs0 * 1.05f) continue;
2148
+ sc = sse;
2149
+ }
2150
+ if (sc < best) {
2151
+ best = sc;
2152
+ bqi = qi; bqj = qj;
2153
+ bei = ei; bej = ej;
2154
+ }
2155
+ }
2156
+ }
2157
+ L[i] = (uint8_t)bqi;
2158
+ L[j] = (uint8_t)bqj;
2159
+ dc = dc_rest + bei + bej;
2160
+ e[i] = bei;
2161
+ e[j] = bej;
2162
+ }
2163
+ }
2164
+ }
2165
+ }
2166
+
2167
+ /* Joint HPC-graph decode of a vesica sub-block pair (A, A^8).
2168
+ * Edge 4x4 is L2-optimal codes for a fixed scale pair and is separable.
2169
+ * Vertex (Σv²) and DC² on the 32-weight group are not: they couple the
2170
+ * two HPC 6-states. Fold's RMSE win is that pair-sum; scoring
2171
+ * Edge SSE + λ_vw Vertex + λ_dc DC²
2172
+ * picks complementary (ls,lm) so WLS (d, dmin) sees one s, not two.
2173
+ * Stock layout, no permutation. */
2174
+ static void hex_joint_fold_subpairs(const float *x, const float *w_row,
2175
+ float dm, float mm,
2176
+ uint8_t state_ls[][6], uint8_t state_lm[][6],
2177
+ uint8_t *Ls, uint8_t *Lm, uint8_t *L)
2178
+ {
2179
+ for (int p = 0; p < 8; p++) {
2180
+ int A = p, B = p + 8;
2181
+ float best = 1e30f;
2182
+ int bva = 0, bvb = 0;
2183
+ uint8_t bqA[16], bqB[16];
2184
+
2185
+ for (int va = 0; va < 6; va++) {
2186
+ float dA = dm * (float)state_ls[A][va];
2187
+ float mA = mm * (float)state_lm[A][va];
2188
+ for (int vb = 0; vb < 6; vb++) {
2189
+ float dB = dm * (float)state_ls[B][vb];
2190
+ float mB = mm * (float)state_lm[B][vb];
2191
+ float sse = 0.0f, ves = 0.0f, dcg = 0.0f;
2192
+ uint8_t qA[16], qB[16];
2193
+ for (int k = 0; k < 16; k++) {
2194
+ int i = 16 * A + k;
2195
+ int j = 16 * B + k;
2196
+ float wi = w_row ? w_row[i] : 1.0f;
2197
+ float wj = w_row ? w_row[j] : 1.0f;
2198
+ float bestp = 1e30f;
2199
+ int bqi = 0, bqj = 0;
2200
+ float bei = 0.0f, bej = 0.0f;
2201
+ for (int qi = 0; qi <= 3; qi++) {
2202
+ float ei = x[i] - (dA * (float)qi - mA);
2203
+ for (int qj = 0; qj <= 3; qj++) {
2204
+ float ej = x[j] - (dB * (float)qj - mB);
2205
+ float sc = wi * ei * ei + wj * ej * ej;
2206
+ if (sc < bestp) {
2207
+ bestp = sc;
2208
+ bqi = qi;
2209
+ bqj = qj;
2210
+ bei = ei;
2211
+ bej = ej;
2212
+ }
2213
+ }
2214
+ }
2215
+ qA[k] = (uint8_t)bqi;
2216
+ qB[k] = (uint8_t)bqj;
2217
+ sse += bestp;
2218
+ ves += (bei + bej) * (bei + bej);
2219
+ dcg += bei + bej;
2220
+ }
2221
+ float sc = HEX_VESICA_SSE_EPS * sse
2222
+ + HEX_VW_LAMBDA * ves
2223
+ + HEX_DC_LAMBDA * dcg * dcg;
2224
+ if (sc < best) {
2225
+ best = sc;
2226
+ bva = va;
2227
+ bvb = vb;
2228
+ memcpy(bqA, qA, 16);
2229
+ memcpy(bqB, qB, 16);
2230
+ }
2231
+ }
2232
+ }
2233
+ Ls[A] = state_ls[A][bva];
2234
+ Lm[A] = state_lm[A][bva];
2235
+ Ls[B] = state_ls[B][bvb];
2236
+ Lm[B] = state_lm[B][bvb];
2237
+ memcpy(L + 16 * A, bqA, 16);
2238
+ memcpy(L + 16 * B, bqB, 16);
2239
+ }
2240
+ }
2241
+
2242
+ /* Edge 4x4 on every fold pair for committed (d, dmin, Ls, Lm). */
2243
+ static void hex_pair_edge_codes(const float *x, const float *w_row,
2244
+ float dm, float mm,
2245
+ const uint8_t *Ls, const uint8_t *Lm,
2246
+ uint8_t *L)
2247
+ {
2248
+ for (int p = 0; p < QK_K / 2; p++) {
2249
+ int i = p, j = p + QK_K / 2;
2250
+ int si = i >> 4, sj = j >> 4;
2251
+ float di = dm * (float)Ls[si], mi = mm * (float)Lm[si];
2252
+ float dj = dm * (float)Ls[sj], mj = mm * (float)Lm[sj];
2253
+ float wi = w_row ? w_row[i] : 1.0f;
2254
+ float wj = w_row ? w_row[j] : 1.0f;
2255
+ float best = 1e30f;
2256
+ int bqi = 0, bqj = 0;
2257
+ for (int qi = 0; qi <= 3; qi++) {
2258
+ float ei = x[i] - (di * (float)qi - mi);
2259
+ for (int qj = 0; qj <= 3; qj++) {
2260
+ float ej = x[j] - (dj * (float)qj - mj);
2261
+ float sc = wi * ei * ei + wj * ej * ej;
2262
+ if (sc < best) {
2263
+ best = sc;
2264
+ bqi = qi;
2265
+ bqj = qj;
2266
+ }
2267
+ }
2268
+ }
2269
+ L[i] = (uint8_t)bqi;
2270
+ L[j] = (uint8_t)bqj;
2271
  }
 
 
2272
  }
2273
 
2274
  /* Robust temperature estimator for the HExState measurement model.
 
2905
  }
2906
 
2907
  float v_live[QK4_0 / 2];
2908
+ float vesica_cur = 0.0f, dc_cur = 0.0f, cross_cur = 0.0f;
2909
  for (int j = 0; j < QK4_0 / 2; j++) {
2910
  v_live[j] = e_live[j] + e_live[j + QK4_0 / 2];
2911
  vesica_cur += v_live[j] * v_live[j];
2912
+ cross_cur += e_live[j] * e_live[j + QK4_0 / 2];
2913
  }
2914
  for (int j = 0; j < QK4_0; j++) dc_cur += e_live[j];
2915
+ float metric_cur = hex_spectral_from_moments(dc_cur, vesica_cur, QK4_0);
2916
 
2917
  /* Deterministic greedy descent: only strict improvements.
2918
  * The previous SA acceptance called rand() inside an OpenMP
 
2942
 
2943
  float vesica_alt = vesica_cur - v_old * v_old + v_new * v_new;
2944
  float dc_alt = dc_cur + de;
2945
+ float metric_alt = hex_spectral_from_moments(dc_alt, vesica_alt, QK4_0);
2946
 
2947
  float delta = metric_cur - metric_alt;
2948
  if (delta > best_delta) {
 
2966
 
2967
  vesica_cur += v_new_commit * v_new_commit - v_old_commit * v_old_commit;
2968
  dc_cur += de_commit;
2969
+ metric_cur = hex_spectral_from_moments(dc_cur, vesica_cur, QK4_0);
2970
 
2971
  v_live[pi_commit] = v_new_commit;
2972
  e_live[best_k] = e_new_commit;
 
3277
  int8_t qs_shaped[QK8_0];
3278
  memcpy(qs_shaped, qs, QK8_0);
3279
  float e_live[QK8_0], v_live[QK8_0 / 2];
3280
+ float vesica_cur = 0.0f, dc_cur = 0.0f, cross_cur = 0.0f;
3281
  for (int k = 0; k < QK8_0; k++)
3282
  e_live[k] = bw[k] - (float)qs_shaped[k] * d;
3283
  for (int p = 0; p < QK8_0 / 2; p++) {
3284
  v_live[p] = e_live[p] + e_live[p + QK8_0 / 2];
3285
  vesica_cur += v_live[p] * v_live[p];
3286
+ cross_cur += e_live[p] * e_live[p + QK8_0 / 2];
3287
  dc_cur += v_live[p];
3288
  }
3289
+ float metric_cur = hex_spectral_from_moments(dc_cur, vesica_cur, QK8_0);
3290
  for (int pass = 0; pass < QK8_0; pass++) {
3291
  int best_k = -1, best_q_alt = 0;
3292
  float best_delta = 0.0f;
 
3301
  float ves_a = vesica_cur - v_live[pi] * v_live[pi]
3302
  + v_new * v_new;
3303
  float dc_a = dc_cur + de;
3304
+ float delta = metric_cur
3305
+ - hex_spectral_from_moments(dc_a, ves_a, QK8_0);
3306
  if (delta > best_delta) {
3307
  best_delta = delta; best_k = k; best_q_alt = q_try;
3308
  }
 
3316
  float v_new = v_live[pi] + de;
3317
  vesica_cur += v_new * v_new - v_live[pi] * v_live[pi];
3318
  dc_cur += de;
3319
+ metric_cur = hex_spectral_from_moments(dc_cur, vesica_cur, QK8_0);
3320
  v_live[pi] = v_new;
3321
  e_live[best_k] = e_new;
3322
  qs_shaped[best_k] = (int8_t)best_q_alt;
 
3415
  else clamp_sigma = 4.0f; /* very heavy tails: stabilise seed */
3416
  float clamp_val = w_sigma * clamp_sigma;
3417
 
3418
+ /* Fold layout: order-statistic pair i ↔ i+128 so pair-sum is ~constant
3419
+ * and each 16-wide scale group is magnitude-homogeneous. The GGUF
3420
+ * stores this order; llama.cpp dequants it as ordinary Q2_K. */
3421
+ float *fold_w = (float *)malloc((size_t)n_elements * sizeof(float));
3422
+ float *fold_i = NULL;
3423
+ if (!fold_w) {
3424
+ if (out_total_error) *out_total_error = -1.0f;
3425
+ return;
3426
+ }
3427
+ memcpy(fold_w, weights, (size_t)n_elements * sizeof(float));
3428
+ if (imat_importance) {
3429
+ fold_i = (float *)malloc((size_t)n_elements * sizeof(float));
3430
+ if (!fold_i) {
3431
+ free(fold_w);
3432
+ if (out_total_error) *out_total_error = -1.0f;
3433
+ return;
3434
+ }
3435
+ memcpy(fold_i, imat_importance, (size_t)n_elements * sizeof(float));
3436
+ }
3437
+ for (int64_t b = 0; b < n_blocks; b++)
3438
+ hex_fold_interleave_block(fold_w + b * QK_K,
3439
+ fold_i ? fold_i + b * QK_K : NULL);
3440
+ weights = fold_w;
3441
+ imat_importance = fold_i;
3442
+
3443
  /* ══════════════════════════════════════════════════════════════════
3444
  * PHASE 1: Greedy quantization — produce seed (d, dmin) per block
3445
  * ══════════════════════════════════════════════════════════════════ */
 
3453
 
3454
  BlockSeed *seeds = (BlockSeed *)calloc(n_blocks, sizeof(BlockSeed));
3455
  if (!seeds) {
3456
+ free(fold_w); free(fold_i);
3457
  if (out_total_error) *out_total_error = -1.0f;
3458
  return;
3459
  }
 
3516
  sizeof(float[TOTAL_SCALE_CANDIDATES]));
3517
  if (!seeds || !candidate_errors) {
3518
  free(seeds); free(candidate_errors);
3519
+ free(fold_w); free(fold_i);
3520
  if (out_total_error) *out_total_error = -1.0f;
3521
  return;
3522
  }
 
4128
  * this array (safe: written sequentially before the parallel loop).
4129
  * ══════════════════════════════════════════════════════════════════ */
4130
 
 
 
4131
  float *block_dc_bias = (float *)calloc(n_blocks, sizeof(float));
4132
 
4133
  if (block_dc_bias) {
 
4156
  dm0, mm0, dc_Ls, dc_Lm);
4157
 
4158
  /* Bias applied to THIS block's WLS targets */
4159
+ float dc_bias = (hex_dc_decay * rolling_dc) / (float)QK_K;
4160
  block_dc_bias[blk] = dc_bias;
4161
 
4162
  /* Quick round-nearest quant to estimate DC residual for NEXT block.
 
4228
  uint16_t prev_dm16 = 0, prev_mm16 = 0;
4229
  for (int ls_iter = 0; ls_iter < 5; ls_iter++) {
4230
 
4231
+ uint8_t L[QK_K];
4232
  uint8_t state_ls[N_SUB][6];
4233
  uint8_t state_lm[N_SUB][6];
4234
  float state_err[N_SUB][6];
 
4311
 
4312
  for (int j = 0; j < N_SUB - 1; j++)
4313
  hpc_cz(sg, j, j + 1);
4314
+ /* Vesica partners (sub j <-> j^8) share fold pairs i <-> i+128.
4315
+ * Entangle them so sieve back-action diversifies complementary
4316
+ * bins (triality pair: each view of one is a role of the other). */
4317
+ for (int j = 0; j < 8; j++)
4318
+ hpc_cz(sg, j, j + 8);
4319
 
4320
  double sub_marg[N_SUB][6];
4321
  int sub_measured[N_SUB];
 
4324
 
4325
  sieve_measure_graph(sg, N_SUB, sub_marg, sub_measured, 1);
4326
 
4327
+ hex_joint_fold_subpairs(
4328
+ adj_block_x,
4329
+ imat_importance ? (imat_importance + blk * QK_K) : NULL,
4330
+ dm, mm, state_ls, state_lm,
4331
+ Ls_blk, Lm_blk, L);
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
4332
  }
4333
 
4334
  double Saa = 0, Sab = 0, Sbb = 0, Sxa = 0, Sxb = 0;
 
4392
  float d_sub = trial_dm * (float)Ls_blk[j];
4393
  float m_sub = trial_mm * (float)Lm_blk[j];
4394
  for (int k = 0; k < 16; k++) {
4395
+ float x = block_x[16*j+k];
4396
  float w = (imat_importance) ?
4397
  imat_importance[blk * QK_K + 16*j+k] : 1.0f;
4398
  int q;
 
4418
  }
4419
 
4420
  for (int j = 0; j < N_SUB; j++) {
4421
+ const float *sx = block_x + 16 * j;
4422
  float best_sub_err = 1e30f;
4423
  uint8_t best_ls = Ls_blk[j], best_lm = Lm_blk[j];
4424
  for (int try_ls = 0; try_ls <= 15; try_ls++) {
 
4473
  */
4474
  uint8_t L[QK_K];
4475
  {
4476
+ hex_vesica_bake_subcodes(
4477
+ block_x,
4478
+ imat_importance ? (imat_importance + blk * QK_K) : NULL,
4479
+ dm, mm, Ls_blk, Lm_blk, L, 0);
4480
+ int q_base_all[QK_K];
4481
+ int q_shaped_all[QK_K];
4482
+ for (int i = 0; i < QK_K; i++)
4483
+ q_base_all[i] = (int)L[i];
 
 
 
 
 
 
 
 
 
 
 
 
4484
  memcpy(q_shaped_all, q_base_all, QK_K * sizeof(int));
4485
 
4486
  float e_live[QK_K];
 
4492
  * d_s == 0 (NOT 0 — the −dmin·lm term always applies). */
4493
  float deq = d_s * (float)q_shaped_all[i] - m_s;
4494
  /* Residual against the adjusted target (DC-corrected view) */
4495
+ e_live[i] = block_x[i] - deq;
4496
  }
4497
 
4498
  float v_live[QK_K / 2];
4499
+ float vesica_cur = 0.0f, dc_cur = 0.0f, cross_cur = 0.0f;
4500
  for (int i = 0; i < QK_K / 2; i++) {
4501
  v_live[i] = e_live[i] + e_live[i + QK_K / 2];
4502
  vesica_cur += v_live[i] * v_live[i];
4503
+ cross_cur += e_live[i] * e_live[i + QK_K / 2];
4504
  }
4505
  for (int i = 0; i < QK_K; i++) dc_cur += e_live[i];
4506
+ float sse_cur = 0.0f;
4507
+ for (int i = 0; i < QK_K; i++) {
4508
+ float w = (imat_importance) ? imat_importance[blk * QK_K + i] : 1.0f;
4509
+ sse_cur += w * e_live[i] * e_live[i];
4510
+ }
4511
+ float metric_cur = HEX_VESICA_SSE_EPS * sse_cur
4512
+ + hex_spectral_from_moments(dc_cur, vesica_cur, QK_K);
4513
 
4514
+ /* Deterministic greedy: Edge SSE + Vertex/DC in one metric. */
4515
  for (int pass = 0; pass < QK_K; pass++) {
4516
  int best_k = -1;
4517
  int best_q_alt = 0;
 
4525
  int q_cur = q_shaped_all[k];
4526
  float m_s = mm * (float)(output[blk].scales[jj] >> 4);
4527
  int pi = (k < QK_K / 2) ? k : k - QK_K / 2;
4528
+ float wk = (imat_importance) ? imat_importance[blk * QK_K + k] : 1.0f;
4529
 
 
 
 
 
4530
  for (int q_try = 0; q_try <= 3; q_try++) {
4531
  if (q_try == q_cur) continue;
4532
 
4533
+ float e_new = block_x[k] - (d_s * (float)q_try - m_s);
4534
  float de = e_new - e_live[k];
4535
 
4536
  float v_new = v_live[pi] + de;
4537
 
4538
  float vesica_alt = vesica_cur - v_live[pi]*v_live[pi] + v_new*v_new;
4539
  float dc_alt = dc_cur + de;
4540
+ float sse_alt = sse_cur
4541
+ - wk * e_live[k] * e_live[k]
4542
+ + wk * e_new * e_new;
4543
+ float met_alt = HEX_VESICA_SSE_EPS * sse_alt
4544
+ + hex_spectral_from_moments(dc_alt, vesica_alt, QK_K);
4545
+ float delta = metric_cur - met_alt;
4546
 
4547
  if (delta > best_delta) {
4548
  best_delta = delta;
 
4559
  int jj_c = best_k >> 4;
4560
  float d_c = dm * (float)(output[blk].scales[jj_c] & 0xF);
4561
  float m_c = mm * (float)(output[blk].scales[jj_c] >> 4);
4562
+ float e_new_c = block_x[best_k] - (d_c * (float)best_q_alt - m_c);
4563
  float de_c = e_new_c - e_live[best_k];
4564
  int pi_c = (best_k < QK_K / 2) ? best_k : best_k - QK_K / 2;
4565
  float v_new_c = v_live[pi_c] + de_c;
4566
+ float wk_c = (imat_importance)
4567
+ ? imat_importance[blk * QK_K + best_k] : 1.0f;
4568
+ sse_cur += wk_c * (e_new_c * e_new_c
4569
+ - e_live[best_k] * e_live[best_k]);
4570
  vesica_cur += v_new_c * v_new_c - v_live[pi_c] * v_live[pi_c];
4571
  dc_cur += de_c;
4572
+ metric_cur = HEX_VESICA_SSE_EPS * sse_cur
4573
+ + hex_spectral_from_moments(dc_cur, vesica_cur, QK_K);
4574
  v_live[pi_c] = v_new_c;
4575
  e_live[best_k]= e_new_c;
4576
  }
 
4882
  * candidate against the whole-block penalty with the partner
4883
  * residuals held fixed — exact coordinate descent on E. */
4884
  float pe[QK_K];
4885
+ float sub_sse[16], sub_dc[16], pair_ves[8];
4886
+ float dc_tot = 0.0f, ves_tot = 0.0f;
4887
  for (int j = 0; j < N_SUB; j++) {
4888
  float d_sub = dm * (float)pl_Ls[j];
4889
  float m_sub = mm * (float)pl_Lm[j];
 
4902
  dc_tot += sub_dc[j];
4903
  }
4904
  for (int p = 0; p < 8; p++) {
4905
+ pair_ves[p] = 0.0f;
4906
+ for (int k = 0; k < 16; k++) {
4907
+ float v = pe[16*p + k] + pe[16*(p+8) + k];
4908
+ pair_ves[p] += v * v;
4909
+ }
4910
+ ves_tot += pair_ves[p];
4911
  }
4912
 
4913
  for (int j = 0; j < N_SUB; j++) {
 
4916
  int pj = j ^ 8; /* partner sub-block */
4917
  const float *ppe = pe + 16 * pj; /* partner residuals */
4918
  float dc_rest = dc_tot - sub_dc[j];
4919
+ float ves_rest = ves_tot - pair_ves[pi];
4920
 
4921
  /* Extended score of the CURRENT committed state */
4922
  float best_sub = sub_sse[j]
4923
  + (HEX_DC_LAMBDA / (float)QK_K) * dc_tot * dc_tot
4924
+ + (HEX_VW_LAMBDA / (float)QK_K) * ves_tot;
4925
  int best_ls = -1, best_lm = 0;
4926
  uint8_t best_q[16];
4927
  float best_e[16];
 
4950
  e_loc[k] = e;
4951
  sub_err += e * e * w;
4952
  dcc += e;
4953
+ {
4954
+ float vpk = e + ppe[k];
4955
+ cxc += vpk * vpk;
4956
+ }
4957
  if (HEX_VW_LAMBDA == 0.0f &&
4958
  sub_err >= best_sub) { aborted = 1; break; }
4959
  }
 
4961
  float score = sub_err
4962
  + (HEX_DC_LAMBDA / (float)QK_K)
4963
  * (dc_rest + dcc) * (dc_rest + dcc)
4964
+ + (HEX_VW_LAMBDA / (float)QK_K)
4965
+ * (ves_rest + cxc);
4966
  if (score < best_sub) {
4967
  best_sub = score;
4968
  best_ls = try_ls;
 
4983
  memcpy(pe + 16 * j, best_e, sizeof(best_e));
4984
  sub_sse[j] = best_sse;
4985
  sub_dc[j] = best_dcc;
4986
+ pair_ves[pi] = best_cxc;
4987
  dc_tot = dc_rest + best_dcc;
4988
+ ves_tot = ves_rest + best_cxc;
4989
  pol_improved = 1;
4990
  }
4991
  }
 
5146
  output[blk].dmin = gguf_fp32_to_fp16(mm);
5147
  }
5148
 
5149
+ /* Bake vesica subcodes into stock 2-bit + nibble fields.
5150
+ * Accept only if block SSE drops; 4.7 is the backstop. */
5151
+ {
5152
+ uint8_t gLs[16], gLm[16];
5153
+ const float *wrow = imat_importance
5154
+ ? (imat_importance + blk * QK_K) : NULL;
5155
+ for (int sj = 0; sj < N_SUB; sj++) {
5156
+ gLs[sj] = output[blk].scales[sj] & 0xF;
5157
+ gLm[sj] = output[blk].scales[sj] >> 4;
5158
+ }
5159
+ hex_vesica_commit(block_x, wrow, &dm, &mm, gLs, gLm, L);
5160
+ for (int sj = 0; sj < N_SUB; sj++)
5161
+ output[blk].scales[sj] = gLs[sj] | (gLm[sj] << 4);
5162
+ output[blk].d = gguf_fp32_to_fp16(dm);
5163
+ output[blk].dmin = gguf_fp32_to_fp16(mm);
5164
+ memcpy(Ls_blk, gLs, 16);
5165
+ memcpy(Lm_blk, gLm, 16);
5166
+ }
5167
+
5168
+ /* ══ PHASE 4.7: CANDIDATE FLOOR ══
5169
  *
5170
+ * candidate_errors is pure nearest-int SSE (the "theoretical floor"
5171
+ * Fold beat by changing coordinates). Compare that to pipeline SSE
5172
+ * only — never SSE+spectral. Adding the vesica/DC penalty here made
5173
+ * fin_err > g_best whenever λ>0, so every block snapped back to
5174
+ * independent nearest-int and threw away HPC+triality+vesica. */
 
 
 
 
5175
  {
5176
  float fin_err = 0.0f;
 
5177
  for (int j = 0; j < N_SUB; j++) {
5178
  float d_sub = dm * (float)(output[blk].scales[j] & 0xF);
5179
  float m_sub = mm * (float)(output[blk].scales[j] >> 4);
 
5181
  int idx = 16 * j + k;
5182
  float w = (imat_importance) ?
5183
  imat_importance[blk * QK_K + idx] : 1.0f;
5184
+ float e = block_x[idx] -
5185
+ (d_sub * (float)L[idx] - m_sub);
5186
+ fin_err += e * e * w;
5187
  }
5188
  }
 
5189
 
5190
  float g_best = candidate_errors[blk][0];
5191
  int g_cand = 0;
 
5225
  }
5226
  }
5227
 
5228
+ if (HEX_VW_LAMBDA != 0.0f || HEX_DC_LAMBDA != 0.0f) {
5229
+ uint8_t gLs[16], gLm[16];
5230
+ const float *wrow = imat_importance
5231
+ ? (imat_importance + blk * QK_K) : NULL;
5232
+ for (int sj = 0; sj < N_SUB; sj++) {
5233
+ gLs[sj] = output[blk].scales[sj] & 0xF;
5234
+ gLm[sj] = output[blk].scales[sj] >> 4;
5235
+ }
5236
+ float sse_pre = hex_block_sse(block_x, wrow, dm, mm, gLs, gLm, L);
5237
+ uint8_t Lpre[QK_K];
5238
+ memcpy(Lpre, L, QK_K);
5239
+ hex_vesica_guide_pairs(block_x, wrow, dm, mm, gLs, gLm, L);
5240
+ if (hex_block_sse(block_x, wrow, dm, mm, gLs, gLm, L) > sse_pre)
5241
+ memcpy(L, Lpre, QK_K);
5242
+ }
5243
+
5244
  for (int j = 0; j < QK_K; j += 128) {
5245
  for (int l = 0; l < 32; l++) {
5246
  output[blk].qs[j / 4 + l] = L[j + l]
 
5267
  free(seeds);
5268
  free(candidate_errors);
5269
  free(best_candidate);
5270
+ free(fold_w);
5271
+ free(fold_i);
5272
  if (out_total_error) *out_total_error = total_err;
5273
 
5274
  if (verbose) {