---
title: "High-precision refinement of the Italian opera time weights"
subtitle: "Same data and default penalty · resolving weakly regularized directions"
format:
  html:
    embed-resources: true
    toc: true
    code-fold: true
jupyter: python3
execute:
  warning: false
---

## Why refine these cases

Three Italian opera specifications have different double-precision QP weights
from different initializers despite very small objective gaps. At the default
penalty, the objective is almost flat along some weight directions. All eleven
nondegenerate opera blocks are therefore solved at **80 decimal digits**, rather
than refining only the three most conspicuous cases. Default tuning is retained.
The representative example selected before fitting is unchanged.

This is a primal active-set solution of the same nonnegative, unit-sum quadratic
program. The intercept is profiled out by centering across donors. The intended
least-squares objective is evaluated in high precision using the same stored
outcome values; the default noise-scale recipe is unchanged. The tiny rounding
difference in its evaluation is checked against the original saved penalty.
Every solution is checked for feasibility and the simplex dual gap. The three
initializer-sensitive cases are solved again from the default weights.

```{python}
#| label: high-precision-solver
import csv, json, math, time
from pathlib import Path
import mpmath as mp
import numpy as np
mp.mp.dps=80
root=Path.cwd().parents[1]
inputs=root/'analyses/sdid-time-weights-default/inputs'
def read(path):
    with open(path) as f:return list(csv.DictReader(f))
manifest={r['fit_id']:r for r in read(inputs/'manifest.csv')}
base=read('fits.csv'); vectors=read('time-weights.csv')
original={r['fit_id']:r for r in read('../sdid-time-weights-default/fits.csv')}
initial=read('initialization-sensitivity.csv')
sensitive={r['fit_id'] for r in initial if float(r['l1_difference'])>1e-3}
lookup={}
for r in vectors:lookup.setdefault((r['fit_id'],r['method']),[]).append(float(r['lambda']))

def solve(H,g,start):
    n=len(g);tol=mp.mpf('1e-55')
    x=mp.matrix([mp.mpf(float(v)) for v in start]);x=x/sum(x)
    S=[i for i in range(n) if x[i]>0]
    for it in range(1000):
        k=len(S);K=mp.matrix(k+1);rhs=mp.matrix(k+1,1)
        for i,ii in enumerate(S):
            rhs[i]=g[ii];K[i,k]=1;K[k,i]=1
            for j,jj in enumerate(S):K[i,j]=H[ii,jj]
        rhs[k]=1
        ans=mp.lu_solve(K,rhs)
        y=mp.matrix(n,1)
        for i,ii in enumerate(S):y[ii]=ans[i]
        negative=[i for i in S if y[i]<-tol]
        if negative:
            alpha=min(x[i]/(x[i]-y[i]) for i in negative)
            x=x+alpha*(y-x)
            for i in S:
                if abs(x[i])<tol:x[i]=0
            S=[i for i in S if x[i]>0]
            continue
        x=y
        gradient=2*(H*x-g)
        mu=sum(gradient[i] for i in S)/len(S)
        outside=[i for i in range(n) if i not in S]
        if not outside or min(gradient[i]-mu for i in outside)>=-tol:
            gap=sum(x[i]*gradient[i] for i in range(n))-min(gradient)
            assert min(x)>=-tol and abs(sum(x)-1)<tol and abs(gap)<mp.mpf('1e-50')
            return x,max(mp.mpf(0),gap),it+1
        S.append(min(outside,key=lambda i:gradient[i]))
    raise RuntimeError('Active-set iteration limit')

refined={};receipts=[];export=[]
for r in [r for r in base if r['method']=='qp_uniform' and r['paper_id']=='study-08']:
    tic=time.time();fid=r['fit_id'];m=manifest[fid];n=int(m['T0'])
    units=read(inputs/m['units_file'])
    with (inputs/m['matrix_file']).open() as f:all_y=list(csv.reader(f))
    # Convert via binary64 to use exactly the same stored outcome values as R.
    Y=[[mp.mpf(float(v)) for v in row] for row,u in zip(all_y,units) if int(u['first_treated_column_zero_based'])==-1]
    N=len(Y);P=len(Y[0])-n
    diffs=[row[t]-row[t-1] for row in Y for t in range(1,n)]
    dm=sum(diffs)/len(diffs);sigma=mp.sqrt(sum((v-dm)**2 for v in diffs)/(len(diffs)-1))
    zeta=mp.mpf('1e-6')*sigma
    assert abs(float(zeta)/float(original[fid]['zeta_lambda'])-1)<1e-12
    means=[sum(Y[i][t] for i in range(N))/N for t in range(n)]
    target=[sum(row[n:])/P for row in Y];target_mean=sum(target)/N
    A=mp.matrix([[(row[t]-means[t])/sigma for t in range(n)] for row in Y])
    b=mp.matrix([(v-target_mean)/sigma for v in target])
    H=A.T*A/N+mp.eye(n)*mp.mpf('1e-12');g=A.T*b/N
    x,gap,it=solve(H,g,[1/n]*n)
    init_diff=mp.mpf(0)
    if fid in sensitive:
        xx,gg,ii=solve(H,g,lookup[fid,'default'])
        init_diff=sum(abs(x[i]-xx[i]) for i in range(n))
        assert init_diff<mp.mpf('1e-45')
    xf=np.asarray([float(v) for v in x]);xf[np.abs(xf)<1e-50]=0;xf/=xf.sum()
    Af=np.asarray(A.tolist(),dtype=float);bf=np.asarray(list(b),dtype=float)
    gradient=2*Af.T@(Af@xf-bf)/N+2e-12*xf
    double_gap=max(0.,float(xf@gradient-gradient.min()))
    default=np.asarray(lookup[fid,'default']);qpu=np.asarray(lookup[fid,'qp_uniform'])
    obj=float(mp.fsum(v*v for v in A*x-b)/N+mp.mpf('1e-12')*mp.fsum(v*v for v in x))
    rr=dict(r);rr.update(method='qp_reference',objective=obj,gap=double_gap,accepted='TRUE',iterations=it,
        iteration_cap='FALSE',solver_status='80-digit active-set KKT verified',raw_violation=0,repair_l1=0,
        seconds=time.time()-tic,l1_from_default=float(np.abs(xf-default).sum()),
        max_change_from_default=float(np.abs(xf-default).max()),
        donor_prediction_rms_change=float(np.sqrt(np.mean((Af@(xf-default))**2))),
        support_1e6=int((xf>1e-6).sum()),effective_periods=float(1/(xf@xf)),
        mass_recent=float(xf[-math.ceil(.2*n):].sum()),mean_position=float(xf@np.linspace(0,1,n)))
    refined[fid]=rr
    template=[w for w in vectors if w['fit_id']==fid and w['method']=='qp_uniform']
    for j,w in enumerate(template):export.append(dict(w,method='qp_reference',lambda_value=float(xf[j]),raw_lambda=float(xf[j])))
    receipts.append({'fit_id':fid,'digits':80,'active_set_steps':it,'mp_gap':mp.nstr(gap,30),
      'initializer_check':fid in sensitive,'mp_initializer_l1':mp.nstr(init_diff,30),
      'l1_from_double_qp':float(np.abs(xf-qpu).sum()),'max_from_double_qp':float(np.abs(xf-qpu).max()),
      'rounded_double_gap':double_gap,'support_1e6':int((xf>1e-6).sum())})
    print(fid,'support',int((xf>1e-6).sum()),'L1 from double QP',float(np.abs(xf-qpu).sum()),flush=True)

def write(name,rows):
    with open(name,'w') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]),lineterminator="\n");w.writeheader();w.writerows(rows)
# The reference uses high precision for opera and the independently checked
# double QP for the remaining papers. The original five variants stay unchanged.
reference=[];reference_weights=[]
for r in [r for r in base if r['method']=='qp_uniform']:
    fid=r['fit_id'];rr=refined.get(fid,dict(r,method='qp_reference'))
    rr['precision_digits']=80 if fid in refined else 16
    reference.append(rr)
for w in vectors:
    if w['method']=='qp_uniform' and w['fit_id'] not in refined:reference_weights.append(dict(w,method='qp_reference'))
for w in export:
    w['lambda']=w.pop('lambda_value');reference_weights.append(w)
write('reference-fits.csv',reference);write('reference-weights.csv',reference_weights)
write('precision-checks.csv',receipts)
Path('precision-verification.json').write_text(json.dumps({'refined_fits':len(refined),'digits':80,
    'initializer_checks':len(sensitive),'all_mp_gaps_below':'1e-50',
    'all_checked_initializer_l1_below':'1e-45','mpmath_version':mp.__version__},indent=2)+'\n')
```

## Results

```{python}
#| echo: false
import pandas as pd
pd.DataFrame(receipts)
```

The high-precision gap concerns the stated numerical objective, not data accuracy
or statistical uncertainty. The near-zero coefficient makes its minimizer very
sensitive to the precise objective. Large differences from default sparsified
weights can coexist with essentially unchanged fitted donor predictions.

- [High-precision checks](precision-checks.csv)
- [Full optimized reference metadata](reference-fits.csv)
- [Full optimized reference weights](reference-weights.csv)
- [Executable precision-refinement source](precision.qmd)