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-tranche2/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 base if r['method']=='default'}
initial=read('initialization-sensitivity.csv')
sensitive={r['fit_id'] for r in initial if float(r['l1_difference'])>1e-5}
lookup={}
for r in vectors:lookup.setdefault((r['fit_id'],r['method']),[]).append(float(r['lambda']))
def solve(A,b,start):
N,n=A.rows,A.cols;rho=mp.mpf('1e-12');tol=mp.mpf('1e-55')
z=np.asarray(start,dtype=float).copy();z[z<1e-9]=0
x=mp.matrix([mp.mpf(float(v)) for v in z]);x=x/sum(x)
S=[i for i in range(n) if x[i]>0];gram={}
def hij(i,j):
key=(min(i,j),max(i,j))
if key not in gram:gram[key]=mp.fsum(A[r,i]*A[r,j] for r in range(N))/N+(rho if i==j else 0)
return gram[key]
g=A.T*b/N
for it in range(1000):
k=len(S);y=mp.matrix(n,1)
if k>N:
# Woodbury solve on the donor dimension, preserving tiny ridge.
AS=mp.matrix([[A[r,j] for j in S] for r in range(N)])
K=AS*AS.T+mp.eye(N)*(N*rho)
v=AS.T*mp.lu_solve(K,b)
one=mp.matrix([1]*k)
h=(one-AS.T*mp.lu_solve(K,AS*one))/rho
mu=(sum(v)-1)/sum(h);ans=v-mu*h
for i,j in enumerate(S):y[j]=ans[i]
else:
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]=hij(ii,jj)
rhs[k]=1;ans=mp.lu_solve(K,rhs)
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*(A.T*(A*x-b)/N+rho*x)
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 (int(r['N0'])<=20 or r['fit_id'] in sensitive)]:
tic=time.time();fid=r['fit_id'];m=manifest[fid];n=int(m['T0'])
source=np.memmap(inputs/m['matrix_file'],dtype='<f8',mode='r',shape=(int(m['source_N']),int(m['source_T'])))
selections=read(inputs/m['selection_file'])
ci=[int(z['matrix_row_1based'])-1 for z in selections if z['role']=='control']
yf=np.asarray(source[ci,:n+int(m['T1'])])
Y=[[mp.mpf(float(v)) for v in row] for row in yf]
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)/(1e-6*float(original[fid]['noise_scale']))-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])
x,gap,it=solve(A,b,lookup[fid,'qp_uniform'])
init_diff=mp.mpf(0)
if True:
xx,gg,ii=solve(A,b,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':True,'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)
# Refine low-donor and initializer-sensitive blocks; retain checked QP elsewhere.
# 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(refined),'all_mp_gaps_below':'1e-50',
'all_checked_initializer_l1_below':'1e-45','mpmath_version':mp.__version__},indent=2)+'\n')