FusionAtlasPaths={"~/projects/fusionatlas/code/package/","~/Documents/fusionatlas/code/package","~/Documents/math/fusionatlas/code/package"};
$Path=$Path~Join~FusionAtlasPaths;
<<FusionAtlas`

SetAttributes[Assert,HoldAll]

Assert[x_==y_]:=If[x==y,Null,Print[x," is not equal to ",y, "!"];$Failed]

Assert[1==1]

Assert[3+4==1+2]

TrainGirth[STrain[S_GeneratorData,deficits_List]]:=2GeneratorDepth[S]+2(Plus@@deficits)

GeneratorDimension[GeneratorData[d_,n_,ω_,r_,rcheck_]]:=d

GeneratorDepth[GeneratorData[d_,n_,ω_,r_,rcheck_]]:=n

GeneratorRotation[GeneratorData[d_,n_,ω_,r_,rcheck_]]:=ω

GeneratorRatio[GeneratorData[d_,n_,ω_,r_,rcheck_]]:=r

GeneratorDualRatio[GeneratorData[d_,n_,ω_,r_,rcheck_]]:=rcheck

TrainGirth[STrain[GeneratorData[d,4,ω,r,rcheck],{1}]]

TrainGirth[Rainbow[S_GeneratorData,k_,T_]]:=2k+TrainGirth[T]

Assert[TrainGirth[Rainbow[GeneratorData[d,4,ω,r,rcheck],1,STrain[GeneratorData[d,4,ω,r,rcheck],{}]]]==10]

cappablePositions[STrain[S_GeneratorData,deficits_List],d_]:=Union[Flatten[Outer[Plus,Table[i,{i,-d,-1}],GeneratorDepth[S]+Rest[FoldList[Plus,0,ListConvolve[{1,1},deficits,1,0]]]]]]

cappablePositions[STrain[S_GeneratorData,deficits_List]]:=cappablePositions[STrain[S,deficits],1]

cappablePositions[STrain[GeneratorData[d,4,ω,r,rcheck],{1,1}],1]

cappablePositions[STrain[GeneratorData[d,4,ω,r,rcheck],{1,1}],2]

cappablePositions[STrain[GeneratorData[d,8,ω,r,rcheck],{1,1}],3]

Subsets[cappablePositions[STrain[GeneratorData[d,8,ω,r,rcheck],{1,1}],3],{3}]

expandJWRight[product[X_,f[n_],Y_]]:=product[X,f[n-1]⊗f[1],Y]+Plus@@((((-1))^(#[[1]]+#[[2]]+1))((qi[#[[1]]]qi[#[[2]]])/(qi[n-1]qi[n]))product[X,cup[#[[1]]],f[n-2],cap[#[[2]]],Y]&/@Flatten[Outer[List,cappablePositions[X],cappablePositions[Y]],1])

explodeJW[product[X_,f[n_],Y_]]:=product[X,f[n-1]⊗f[1],Y]+Plus@@((((-1))^(#[[1]]+#[[2]]+1))((qi[#[[1]]]qi[#[[2]]])/(qi[n-1]qi[n]))product[X,cup[#[[1]]],f[n-2],cap[#[[2]]],Y]&/@Flatten[Outer[List,cappablePositions[X],cappablePositions[Y]],1])

GeneratorDepth[GeneratorData[d,n,ω,r,rcheck]]

expandJWRight[product[STrain[GeneratorData[d,n,ω,r,rcheck],{1}],f[2n+2],STrain[GeneratorData[d,n,ω,r,rcheck],{1}]]]

product[STrain[S_GeneratorData,deficits_List],cup[k_],Z___]:=Module[{p},
p=FindFirst[cappablePositions[STrain[S,deficits]],k];
If[p==-1,0,
product[STrain[S,deficits-UnitVector[Length[deficits],p]],Z]
]
]

product[Z___,cap[k_],STrain[S_GeneratorData,deficits_List]]:=Module[{p},
p=FindFirst[cappablePositions[STrain[S,deficits]],k];
If[p==-1,0,
product[Z,STrain[S,deficits-UnitVector[Length[deficits],p]]]
]
]

FindFirst[a_,b_]:=With[{p=Position[a,b,{1},1,Heads->False]},If[Length[p]==0,-1,p[[1,1]]]]

expandf[product[STrain[GeneratorData[d,n,ω,r,rcheck],{1}],f[2n+2],STrain[GeneratorData[d,n,ω,r,rcheck],{1}]]]

Clear[expandf]

expandf[product[X_,f[n_]⊗L_,Y_]]:=product[X,f[n-1]⊗f[1]⊗L,Y]+Plus@@((((-1))^(#[[1]]+#[[2]]+1))((qi[#[[1]]]qi[#[[2]]])/(qi[n-1]qi[n]))product[X,cup[#[[1]]],f[n-2]⊗L,cap[#[[2]]],Y]&/@Cases[Flatten[Outer[List,cappablePositions[X],cappablePositions[Y]],1],{a_,b_}/;a≤n-1∧b≤n-1])

product[STrain[S_GeneratorData,{x__,0,y__}], f[m_], z___]/;m==TrainGirth[STrain[S,{x,0,y}]]:=(1-GeneratorRatio[S])product[STrain[S,{x,y}], f[m], z]+ GeneratorRatio[S]*JWCoefficientFromDeficits[GeneratorDimension[S],GeneratorDepth[S],Last[{x}],First[{y}]]*product[STrain[S,{(Sequence@@Most[{x}]),Last[{x}]+First[{y}],(Sequence@@Rest[{y}])}],f[m],z]

product[STrain[S_GeneratorData,{0,y__}], f[m_], z___]/;m==TrainGirth[STrain[S,{0,y}]]:=(1-GeneratorRatio[S])product[STrain[S,{y}], f[m], z]+ GeneratorRatio[S]*product[STrain[S,Rest[{y}]],f[m],z]

product[STrain[S_GeneratorData,{x__,0}], f[m_], z___]/;m==TrainGirth[STrain[S,{x,0}]]:=(1-GeneratorRatio[S])product[STrain[S,{x}], f[m], z]+ GeneratorRatio[S]*product[STrain[S,Most[{x}]],f[m],z]

product[STrain[S_GeneratorData,{0}], f[m_], z___]/;m==TrainGirth[STrain[S,{0}]]:=(1-GeneratorRatio[S])product[STrain[S,{}], f[m], z]

product[STrain[S_GeneratorData,{}], f[m_]⊗CircleTimes[],, z___]

FullForm[f[1]⊗f[1]]

f[1]⊗f[1]⊗f[2]⊗f[3]/.CircleTimes[nf:(f[1]...),zf___f]:>{{nf},{zf}}

f[2]⊗f[3]⊗f[1]⊗f[1]/.CircleTimes[zf___f,nf:(f[1]...)]:>{{zf},{nf}}

f[1]⊗f[2]⊗f[3]/.{c:CircleTimes[x_f,y___f]:>{c,x,y}}

TrainGirth[STrain[GeneratorData[d,n,ω,r,rcheck],{1,0,1}]]

product[STrain[GeneratorData[d,n,ω,r,rcheck],{1,0,1}],f[2n+4]]

JWCoefficientFromDeficits[d_,n_,a_,b_]/;(n≥a+b):=
Module[{d1=Min[a,b],d3=Max[a,b],q=qFromd[d]},
RootReduce[((Product[qInteger[n-i][q],{i,d1,d1+d3-1}])/(Product[qInteger[n-i][q],{i,0,d3-1}]))]]

JWCoefficientFromDeficits[d_,n_,a_,b_]/;(n<a+b):=0

Assert[JWCoefficientFromDeficits[d,8,0,0]==1]
Assert[Simplify[JWCoefficientFromDeficits[2,8,1,1]]==(((8-1))/8)]
Assert[Simplify[JWCoefficientFromDeficits[2,8,1,2]]==(((8-2))/8)]
Assert[Simplify[JWCoefficientFromDeficits[2,2,1,2]]==0]