cat > /workspace/scratch/fitpchip.py <<'PY'
import json,numpy as np
from scipy.interpolate import PchipInterpolator
from scipy.optimize import least_squares
mods=json.load(open('/workspace/scratch/colourmodels.json'));ref=np.array(json.load(open('/workspace/input/chart_reference.json'))['srgb8'])/255
lin=lambda v:np.where(v<=.04045,v/12.92,((v+.055)/1.055)**2.4)
enc=lambda v:np.where(v<=.0031308,12.92*v,1.055*np.maximum(v,0)**(1/2.4)-.055)
for s,m in mods.items():
 obs=np.array(m['observed']);a=lin(obs);t=lin(ref);curves=[];z=[]
 for c in range(3):
  xx=a[18:,c][::-1];yy=t[18:,c][::-1]
  # fit power + offset for useful extrapolation
  fp=least_squares(lambda p:p[0]*xx**p[1]+p[2]-yy,[1,1,0],bounds=([.01,.2,-.05],[20,4,.05]))
  grid=np.linspace(0,1,65536)
  vals=fp.x[0]*grid**fp.x[1]+fp.x[2]
  interp=PchipInterpolator(xx,yy,extrapolate=True)
  vals[(grid>=xx[0])&(grid<=xx[-1])]=interp(grid[(grid>=xx[0])&(grid<=xx[-1])])
  # positive branch below the lowest step is a power anchored on its slope
  slope=interp.derivative()(xx[0]);g=max(.4,min(2.5,slope*xx[0]/yy[0]));lo=grid<xx[0];vals[lo]=yy[0]*(grid[lo]/xx[0])**g
  vals=np.maximum(vals,0);curves.append(vals.tolist())
  z.append(np.interp(a[:,c],grid,vals))
 z=np.array(z).T
 M=np.linalg.lstsq(np.r_[z,np.ones((1,3))*2],np.r_[t,np.ones((1,3))*2],rcond=None)[0]
 err=enc(z@M)-ref
 print(s,round(abs(err).mean()*255,3),round(abs(err).max()*255,3))
 np.savez('/workspace/scratch/'+s+'_curve.npz',curves=np.array(curves),matrix=M,anchor=m['anchor_linear'])
PY
OPENBLAS_NUM_THREADS=1 python /workspace/scratch/fitpchip.py