import sympy as sp, json
from pathlib import Path
h,b,c,d,e,f,g,r,q,u=sp.symbols('h b c d e f g r q u')
mut=[h,b,c,d,e,f,g,r,q,u]
a=sp.symbols('a0:18'); variables=mut+list(a)
chi=[b*d*a[4]*a[17]/(q*a[3]*a[16]),h*c*r*a[13]/(g*u*a[14]),q*a[16]/(h*b*a[15]),f*a[4]/r,d*a[1]*a[6]/(e*a[3]*a[16]),f*a[6]*a[11]/(c*e*a[13]),a[9]*a[14]/(c*a[13]),u/(h*e*a[12]),u/(h*f*a[8]),b*e*g*a[14]/(f*q*a[6]),h*q*a[5]/(b*d*g*a[4]),f*a[10]/(c*a[13]),b*f*a[7]/(c*e*a[13]),r*a[3]/(d*g*a[4]),b*d*g*a[0]/(q*a[3]*a[16]),b*e*g*a[15]/(f*q*a[5]),d*e*g*a[12]/(f*q*a[5]),d*a[2]/(a[3]*a[16]),g*a[8]/(c*a[5]),f*q*a[5]/(e*r*a[13])]
def exponent(mon):
 pd=sp.sympify(mon).as_powers_dict(); return sp.Matrix([pd.get(v,0) for v in variables])
C=sp.Matrix.vstack(*[exponent(x).T for x in chi+list(a[:8])]); inv=C.T.inv()
assert abs(C.det())==1
A0=g*(b*(a[0]*f+a[1]*a[14]+a[15]*h)+a[11]*a[16]*h)+a[17]*b*(a[4]*f+r)
A1=g*(b*(a[0]*f+a[1]*a[14]+a[15]*h)+h*(a[11]*a[16]+a[12]*d))+a[17]*a[4]*b*f
A2=a[3]*f+d*g
A3=a[1]*a[6]*c+a[2]*c*e+a[4]*a[9]*h+a[5]*c*h+a[8]*g*h
Fh=a[1]*b*c+u; Fb=a[11]*h+q; Fg=a[4]*f+r
Fc=a[9]*h*Fg+a[8]*f*g*h+g*u
Fd=A0+a[16]*g*q
Fe=a[6]*q*Fh+h*(a[6]*a[11]+a[7]*b)*(b*r+u)+b*h*q*(a[4]*a[9]+a[5]*c+a[8]*g)
Ff=a[14]*g*Fh+h*(c*g*(a[11]*a[16]+a[12]*d+a[15]*b)+r*(a[13]*c+a[14]*a[9]))
Fr=c*A2*(e*A1+d*g*u)+d*g*q*(f*A3+g*u)+a[16]*a[3]*c*e*f*g*q
Fq=e*A1+b*r*(a[17]*e+d*g)+d*g*u
Fu=Fr.subs(u,0)+r*(d*g*h*(f*(a[10]*e+a[11]*a[6]+a[7]*b)+e*(a[13]*c+a[14]*a[9]))+a[17]*a[3]*b*c*e*f)
Fs=[Fh,Fb,Fc,Fd,Fe,Ff,Fg,Fr,Fq,Fu]; norms=[1,1,g,b*g,h*b,h,1,1,g,d*g]
leads=[]; supports=[]; H=[]
for i,F in enumerate(Fs):
 terms=sp.Poly(F,*variables).terms(); exps=[sp.Matrix(t[0]) for t in terms]; ce=[inv*v for v in exps]
 idx=[j for j,v in enumerate(ce) if all(all(v[t]<=w[t] for t in range(20)) for w in ce)]
 assert len(idx)==1,(i,len(idx))
 L=exps[idx[0]]; leads.append(L)
 sup=[v-ce[idx[0]] for v in ce]; supports.append(sup)
 H.append(L-exponent(norms[i])-sp.eye(28)[:,i])
 print(mut[i], 'terms',len(terms),'leading',sp.prod(v**int(w) for v,w in zip(variables,L)),flush=True)
# prescribed target binomial exchanges
v=sp.symbols('v1:21')
I1=[1,5,15,18];I2=[6,7,12,13,20];I3=[11,14];I4=[16,17,19]
ex={}
for I,F in [(I1,v[10]+v[2]),(I2,v[1]+v[7]*v[8]),(I3,v[0]*v[3]+v[18]*v[2]),(I4,v[13]*v[1]+v[19]*v[3])]:
 for i in I: ex[i]=(F,sp.prod(v[j-1] for j in I))
for i,F in {2:v[9]*v[15]*v[16]*v[8]+v[11]*v[12]*v[19]*v[3]*v[5],3:v[0]*v[14]*v[16]*v[17]*v[4]+v[10]*v[13],4:v[9]*v[14]*v[15]*v[16]*v[18]+v[10]*v[13]*v[1],8:v[12]*v[19]+v[16],9:v[11]*v[12]*v[18]*v[5]*v[6]+v[1],10:v[11]*v[3]+v[1]*v[4]}.items():ex[i]=(F,v[i-1])
Bmap={}
for group,inds in [(I1,[0,9]),(I2,[9]),(I3,[6,7,8]),(I4,[7]),([2],[5]),([3],[3]),([4],[6]),([8],[1,8]),([9],[2]),([10],[4,5,9])]:
 for i in group:
  beta=sp.zeros(28,1)
  for j in inds:beta[j]=1
  Bmap[i]=beta
rows=[]
for i in range(1,21):
 beta=Bmap[i];delta=C*beta
 F,den=ex[i];pows=[]
 for term in sp.expand(F).as_ordered_terms():
  P=(term/den).as_powers_dict();pows.append(sp.Matrix([P.get(z,0) for z in v]+[0]*8))
 if pows[1]-pows[0]==delta:mu=pows[0]
 elif pows[0]-pows[1]==delta:mu=pows[1]
 else:raise Exception(('cap difference',i))
 source_m=C.inv()*mu
 choices=[]
 for j in range(10):
  dP=min((alpha.dot(mu) for alpha in supports[j]))
  hd=H[j].dot(beta)
  if sum(beta)==1 and hd==-1 and dP==-1:
   # choose principal torus function (x_j')^N prod(other initial x)
   other=sum(beta[t] for t in range(10) if t!=j)
   N=sp.floor(other/(-hd))+1
   delta_f=N*hd+other
   assert delta_f<0
   choices.append((j,int(N),int(hd),int(dP),int(delta_f)))
 rows.append({'target':i,'beta':[int(x) for x in beta[:10]],'label':[int(x) for x in mu[:20]],'choices':choices})
 print(i,'beta',[(str(mut[j]),int(beta[j])) for j in range(10) if beta[j]],'choices',[(str(mut[j]),N,hd,defect) for j,N,hd,dP,defect in choices],flush=True)
assert [r['target'] for r in rows if r['choices']]==[2,3,4,6,7,9,12,13,16,17,19,20]
source_data=[]
for i,F in enumerate(Fs):
 source_data.append({'variable':str(mut[i]),'exchange_terms':[[list(m),str(c0)] for m,c0 in sp.Poly(F,*variables).terms()],
   'normalizer_exponents':[int(x) for x in exponent(norms[i])],
   'least_exponent':[int(x) for x in leads[i]],
   'normalized_support':[[int(x) for x in alpha] for alpha in supports[i]],
   'first_mutant_leading_exponent':[int(x) for x in H[i]]})
Path(__file__).with_name('e7_calibration.json').write_text(json.dumps({'determinant':int(C.det()),'rows':rows,'C':[[int(x) for x in row] for row in C.tolist()],
  'source_variables':[str(x) for x in variables], 'source_exchanges':source_data,
  'scope':'Twelve single-contact directions only. The other eight retain their contact equations.'},indent=2))
print('PASS: 28-dimensional integral marking and twelve single-contact normalizations.')