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]));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')