text:
# stuff done once per exucution--not depending on r
function InitWRE!(o::StrBootTest{T}) where T
	o.Repl.kZ>1 && (o.numer_b = Vector{T}(undef,nrows(o.Repl.RRpar)))

	iszero(o.granular) && (o.Repl.Zperp = o.DGP.Zperp = Matrix{T}(undef,0,0))  # drop this potentially large array

	o.liml && o.Repl.kZ==1 && o.Nw==1 && (o.As = o.β̈s = zeros(1, o.B+1))

	o.S✻ZperpU              = [Matrix{T}(undef, o.Repl.kZperp, o.N✻) for _ ∈ 0:o.Repl.kZ]
	o.invZperpZperpS✻ZperpU = [Matrix{T}(undef, o.Repl.kZperp, o.N✻) for _ ∈ 0:o.Repl.kZ]
	o.S✻YU                  = [Vector{T}(undef, o.N✻) for _ ∈ 0:o.Repl.kZ, _ ∈ 0:o.Repl.kZ]
	o.S✻XU                  = [Matrix{T}(undef, o.DGP.kX, o.N✻) for _ ∈ 0:o.Repl.kZ]
	o.invXXS✻XU             = [Matrix{T}(undef, o.DGP.kX, o.N✻) for _ ∈ 0:o.Repl.kZ]
	o.S✻UU                  = [Vector{T}(undef, o.N✻) for _ ∈ 0:o.Repl.kZ, _ ∈ 0:o.Repl.kZ]
	o.S✻⋂XU₂      = Array{T,3}(undef, o.Repl.kX, o.N✻⋂, o.kY₂)
	o.S✻⋂XU₂RparY = Array{T,3}(undef, o.Repl.kX, o.N✻⋂, o.Repl.kZ)
	o.invXXS✻XU₂  = Array{T,3}(undef, o.Repl.kX, o.N✻ , o.kY₂)
	o.invXXS✻XU₂RparY = Array{T,3}(undef, o.Repl.kX, o.N✻, o.Repl.kZ)
	o.S✻XU₂ = Array{T,3}(undef, o.Repl.kX, o.N✻, o.kY₂)
	o.S✻XU₂RparY = Array{T,3}(undef, o.Repl.kX, o.N✻, o.Repl.kZ)
	o.S✻ZperpU₂ = Array{T,3}(undef, o.Repl.kZperp, o.N✻, o.kY₂)
	o.S✻ZperpU₂RparY = Array{T,3}(undef, o.Repl.kZperp, o.N✻, o.Repl.kZ)
	o.invZperpZperpS✻ZperpU₂ = Array{T,3}(undef, o.Repl.kZperp, o.N✻, o.kY₂)
	o.invZperpZperpS✻ZperpU₂RparY = Array{T,3}(undef, o.Repl.kZperp, o.N✻, o.Repl.kZ)

	o.T1L = isone(o.Nw) ? [Matrix{T}(undef, o.Repl.kX, ncols(o.v))] :
												[Matrix{T}(undef, o.Repl.kX, length(o.WeightGrp[1])), Matrix{T}(undef, o.Repl.kX, length(o.WeightGrp[end]))]
	o.T1R = deepcopy(o.T1L)
	o.T2 = Matrix{T}(undef, o.N✻, o.N✻)
	o.robust && o.bootstrapt && iszero(o.granular) &&
		(o.negS✻UMZperpX = [Array{T,3}(undef, o.Repl.kX, o.N⋂, o.N✻) for _ in 0:o.Repl.kZ])

	if o.bootstrapt
		if o.robust && o.granular
			o.S✻UMZperp = [Matrix{T}(undef, o.Nobs, o.N✻) for _ ∈ 0:o.Repl.kZ]
			o.S✻UPX     = [Matrix{T}(undef, o.Nobs, o.N✻) for _ ∈ 0:o.Repl.kZ]
		end
		if o.liml || !o.robust
			o.YY✻_b   = zeros(o.Repl.kZ+1, o.Repl.kZ+1)
			o.YPXY✻_b = zeros(o.Repl.kZ+1, o.Repl.kZ+1)
		end
		o.NFE>0 && (o.bootstrapt || !isone(o.κ) || o.liml) && (o.CTFEU = Vector{Matrix{T}}(undef, o.Repl.kZ+1))
	end
	o.S✻⋂XY₂      = o.Repl.S✻⋂XY₂     - o.Repl.S✻⋂XZperp     * o.Repl.invZperpZperpZperpY₂  - o.Repl.invZperpZperpZperpX' * (o.Repl.S✻⋂ZperpY₂  - o.Repl.S✻⋂ZperpZperp * o.Repl.invZperpZperpZperpY₂ )
	o.S✻⋂XX       = o.Repl.S✻⋂XX      - o.Repl.S✻⋂XZperp     * o.Repl.invZperpZperpZperpX   - o.Repl.invZperpZperpZperpX' * (o.Repl.S✻⋂XZperp'  - o.Repl.S✻⋂ZperpZperp * o.Repl.invZperpZperpZperpX  )
	o.S✻⋂XDGPZ    = o.DGP.S✻⋂XZpar    - o.Repl.S✻⋂XZperp     * o.DGP.invZperpZperpZperpZpar - o.Repl.invZperpZperpZperpX' * (o.DGP.S✻⋂ZperpZpar - o.Repl.S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpZpar)
	o.S✻⋂Xy₁      = o.Repl.S✻⋂Xy₁     - o.Repl.S✻⋂XZperp     * o.Repl.invZperpZperpZperpy₁  - o.Repl.invZperpZperpZperpX' * (o.Repl.S✻⋂Zperpy₁  - o.Repl.S✻⋂ZperpZperp * o.Repl.invZperpZperpZperpy₁ )
	  S✻⋂ZperpX   = o.Repl.S✻⋂XZperp' - o.Repl.S✻⋂ZperpZperp * o.Repl.invZperpZperpZperpX
	o.DGP.restricted &&
		(o.S✻⋂X_DGPZR₁ = o.DGP.S✻⋂XZR₁     - o.Repl.S✻⋂XZperp     * o.DGP.invZperpZperpZperpZR₁  - o.Repl.invZperpZperpZperpX' * (o.DGP.S✻⋂ZperpZR₁  - o.Repl.S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpZR₁ ))

	o.invXXS✻⋂XY₂     = o.Repl.invXX * o.S✻⋂XY₂
	o.invXXS✻⋂XX      = o.Repl.invXX * o.S✻⋂XX
	o.invXXS✻⋂XDGPZ   = o.Repl.invXX * o.S✻⋂XDGPZ
	o.invXXS✻⋂Xy₁     = o.Repl.invXX * o.S✻⋂Xy₁
	o.DGP.restricted &&
		(o.invXXS✻⋂XDGPZR₁ = o.Repl.invXX * o.S✻⋂X_DGPZR₁)

	_S✻ZperpY₂      = @panelsum(o, o.Repl.S✻⋂ZperpY₂ , o.info✻_✻⋂)  # moments of variables _before_ FWL processing
	_S✻Zperpy₁      = dropdims(@panelsum(o, reshape(o.Repl.S✻⋂Zperpy₁,Val(3)), o.info✻_✻⋂); dims=3)
	_S✻ZperpDGPZpar = @panelsum(o, o.DGP.S✻⋂ZperpZpar, o.info✻_✻⋂)
	o.DGP.restricted &&
		(_S✻ZperpDGPZR₁  = @panelsum(o, o.DGP.S✻⋂ZperpZR₁, o.info✻_✻⋂))

	S✻ZperpZperp    = @panelsum(o, o.Repl.S✻⋂ZperpZperp, o.info✻_✻⋂)
	o.S✻XY₂         = @panelsum(o, o.S✻⋂XY₂   , o.info✻_✻⋂)
	o.S✻XX          = @panelsum(o, o.S✻⋂XX    , o.info✻_✻⋂)
	o.S✻XDGPZ       = @panelsum(o, o.S✻⋂XDGPZ, o.info✻_✻⋂)
	o.S✻Xy₁         = dropdims(@panelsum(o, reshape(o.S✻⋂Xy₁,Val(3)), o.info✻_✻⋂); dims=3)
	o.S✻ZperpX      = @panelsum(o, S✻⋂ZperpX, o.info✻_✻⋂)
	o.S✻ZperpY₂     = _S✻ZperpY₂ - S✻ZperpZperp * o.Repl.invZperpZperpZperpY₂
	o.S✻ZperpDGPZ   = _S✻ZperpDGPZpar - S✻ZperpZperp * o.DGP.invZperpZperpZperpZpar
	o.S✻Zperpy₁     = _S✻Zperpy₁ - S✻ZperpZperp * o.Repl.invZperpZperpZperpy₁
	if o.DGP.restricted
		o.S✻XZR₁        = @panelsum(o, o.S✻⋂X_DGPZR₁, o.info✻_✻⋂)
		o.S✻ZperpDGPZR₁ = @panelsum(o, o.DGP.S✻⋂ZperpZR₁ , o.info✻_✻⋂) - S✻ZperpZperp * o.DGP.invZperpZperpZperpZR₁
	end

	if o.NFE>0 && (o.liml || !isone(o.κ) || o.bootstrapt)
		  CT✻⋂FEX  = [crosstabFE(o, o.Repl.X₁, o.info✻⋂) crosstabFE(o, o.Repl.X₂, o.info✻⋂)]
		o.CT✻FEX   = @panelsum(o, CT✻⋂FEX, o.info✻_✻⋂)
		o.CT✻FEY₂  = crosstabFE(o, o.DGP.Y₂, o.info✻)
		o.CT✻FEZ   = crosstabFE(o, o.DGP.Z, o.info✻)
		o.CT✻FEy₁  = crosstabFE(o, o.DGP.y₁, o.info✻)
		o.DGP.restricted &&
			(o.CT✻FEZR₁ = crosstabFE(o, o.DGP.ZR₁, o.info✻))
	end

	if ((o.robust && o.bootstrapt) || o.liml || !o.robust || !isone(o.κ))
		S✻⋂ReplZX = (o.Repl.S✻⋂XZpar - o.Repl.S✻⋂XZperp * o.Repl.invZperpZperpZperpZpar - o.Repl.invZperpZperpZperpX' * (o.Repl.S✻⋂ZperpZpar - o.Repl.S✻⋂ZperpZperp * o.Repl.invZperpZperpZperpZpar))'
	end

	if o.bootstrapt & o.robust
		o.info⋂_✻⋂ = panelsetup(o.ID✻⋂, o.subcluster+1:o.NClustVar)

		o.S⋂ReplZX = @panelsum(o, S✻⋂ReplZX, o.info⋂_✻⋂)
		S⋂ZperpX  = @panelsum(o, S✻⋂ZperpX, o.info⋂_✻⋂)
		o.S⋂Xy₁  = dropdims(@panelsum(o, reshape(o.S✻⋂Xy₁,Val(3)), o.info⋂_✻⋂); dims=3)

		o.J⋂s = isone(o.Nw) ? [Array{T,3}(undef, o.N⋂, ncols(o.v), o.Repl.kZ)] :
		                      [Array{T,3}(undef, o.N⋂, length(o.WeightGrp[1]), o.Repl.kZ), Array{T,3}(undef, o.N⋂, length(o.WeightGrp[end]), o.Repl.kZ)]

		if o.granular
			o.crosstab✻ind = o.Nobs==o.N✻ ? Vector(diagind(FakeArray(o.N✻,o.N✻))) : LinearIndices(FakeArray(o.Nobs,o.N✻))[CartesianIndex.(1:o.Nobs, o.ID✻)]
			o.XinvXX = X₁₂B(o, o.Repl.X₁, o.Repl.X₂, o.Repl.invXX)
			o.PXZ    = X₁₂B(o, o.Repl.X₁, o.Repl.X₂, o.Repl.invXXXZ)
		else
			inds = o.subcluster>0 ?
							[CartesianIndex(j,i) for (j,v) ∈ enumerate(o.info⋂_✻⋂) for i ∈ v] :  # crosstab ∩,* is wide
							o.NClustVar == o.NBootClustVar ?
									[CartesianIndex(i,i) for i ∈ 1:o.N✻⋂] :  # crosstab *,∩ is square
									[CartesianIndex(i,j) for (j,v) ∈ enumerate(o.clust[o.BootClust].info) for i ∈ v]  # crosstab ∩,* is tall
			inds = [CartesianIndex(k,I) for I ∈ inds for k ∈ 1:o.Repl.kX]
			o.crosstab⋂✻ind = LinearIndices(FakeArray(Tuple(max(inds...))...))[inds]

			o.S⋂XZperpinvZperpZperp = S⋂ZperpX' * o.Repl.invZperpZperp

			o.Q    = Array{T,3}(undef, o.N✻, o.N⋂, o.N✻)
			o.NFE>0 && (o.CT⋂FEX = o.invFEwt .* @panelsum(o, CT✻⋂FEX, o.info⋂_✻⋂))

			o.β̈v = isone(o.Nw) ? [Matrix{T}(undef, o.N✻, ncols(o.v))] :
			                     [Matrix{T}(undef, o.N✻, length(o.WeightGrp[1])), Matrix{T}(undef, o.N✻, length(o.WeightGrp[end]))]

		end
	end

	if o.liml || !o.robust || !isone(o.κ)  # cluster-wise moments after FWL
		o.S✻Y₂Y₂     = o.Repl.S✻Y₂Y₂    - _S✻ZperpY₂'   * o.DGP.invZperpZperpZperpY₂   - o.DGP.invZperpZperpZperpY₂'   * o.S✻ZperpY₂
		o.S✻DGPZDGPZ = o.DGP.S✻ZparZpar - _S✻ZperpDGPZpar' * o.DGP.invZperpZperpZperpZpar - o.DGP.invZperpZperpZperpZpar' * o.S✻ZperpDGPZ
		o.S✻DGPZY₂   = o.DGP.S✻ZparY₂   - _S✻ZperpDGPZpar' * o.DGP.invZperpZperpZperpY₂   - o.DGP.invZperpZperpZperpZpar' * o.S✻ZperpY₂
		o.S✻DGPZy₁   = o.DGP.S✻Zpary₁   - _S✻ZperpDGPZpar' * o.DGP.invZperpZperpZperpy₁   - o.DGP.invZperpZperpZperpZpar' * o.S✻Zperpy₁
		o.S✻Y₂y₁     = o.Repl.S✻Y₂y₁    - _S✻ZperpY₂'   * o.DGP.invZperpZperpZperpy₁   - o.DGP.invZperpZperpZperpY₂'   * o.S✻Zperpy₁
		o.S✻y₁y₁     = o.Repl.S✻y₁y₁    - _S✻Zperpy₁'   * o.DGP.invZperpZperpZperpy₁   - o.S✻Zperpy₁' * o.DGP.invZperpZperpZperpy₁
		o.DGP.restricted &&
			(o.S✻DGPZR₁y₁ = o.DGP.S✻ZR₁y₁ - _S✻ZperpDGPZR₁' * o.DGP.invZperpZperpZperpy₁ - o.DGP.invZperpZperpZperpZR₁' * o.S✻Zperpy₁)

		if o.Repl.restricted
			_S✻ZperpReplZR₁ = @panelsum(o, o.Repl.S✻⋂ZperpZR₁, o.info✻_✻⋂)
			_S✻⋂XReplZR₁    = @panelsum(o, o.Repl.S✻⋂XZR₁    , o.info✻_✻⋂)

			o.r₁S✻ReplZR₁Y₂     = o.r₁' * (o.Repl.S✻ZR₁Y₂ - _S✻ZperpReplZR₁' * o.Repl.invZperpZperpZperpY₂ - o.Repl.invZperpZperpZperpZR₁' * o.S✻ZperpY₂)
			o.r₁S✻ReplZR₁X      = o.r₁' * (_S✻⋂XReplZR₁'  - _S✻ZperpReplZR₁' * o.Repl.invZperpZperpZperpX  - o.Repl.invZperpZperpZperpZR₁' * o.S✻ZperpX )
			o.r₁S✻ReplZR₁DGPZ   = o.r₁' * panelcross(o.Repl.ZR₁, o.DGP.Z, o.info✻)
			o.r₁S✻ReplZR₁y₁     = o.r₁' * (o.Repl.S✻ZR₁y₁ - _S✻ZperpReplZR₁' * o.Repl.invZperpZperpZperpy₁ - o.Repl.invZperpZperpZperpZR₁' * o.S✻Zperpy₁)
			o.DGP.restricted &&
				(o.r₁S✻ReplZR₁DGPZR₁ = o.r₁' * panelcross(o.Repl.ZR₁, o.DGP.ZR₁, o.info✻))
		end

		_S✻ZperpReplZpar = @panelsum(o, o.Repl.S✻⋂ZperpZpar, o.info✻_✻⋂)
		_S✻ReplXZ        = @panelsum(o, o.Repl.S✻⋂XZpar    , o.info✻_✻⋂)

		o.S✻ReplZY₂      = o.Repl.S✻ZparY₂ - _S✻ZperpReplZpar' * o.Repl.invZperpZperpZperpY₂ - o.Repl.invZperpZperpZperpZpar' * o.S✻ZperpY₂
		o.S✻ReplZX       = _S✻ReplXZ'      - _S✻ZperpReplZpar' * o.Repl.invZperpZperpZperpX  - o.Repl.invZperpZperpZperpZpar' * o.S✻ZperpX
		o.S✻ReplZDGPZ    = panelcross(o.Repl.Z, o.DGP.Z, o.info✻)
		o.S✻ReplZy₁      = o.Repl.S✻Zpary₁ - _S✻ZperpReplZpar' * o.Repl.invZperpZperpZperpy₁ - o.Repl.invZperpZperpZperpZpar' * o.S✻Zperpy₁

		if o.DGP.restricted
			_S✻⋂XDGPZR₁ = @panelsum(o, o.DGP.S✻⋂XZR₁, o.info✻_✻⋂)

			o.S✻ReplZDGPZR₁  = panelcross(o.Repl.Z, o.DGP.ZR₁, o.info✻)
			o.S✻DGPZR₁Y₂     = o.DGP.S✻ZR₁Y₂  - _S✻ZperpDGPZR₁' * o.Repl.invZperpZperpZperpY₂  - o.DGP.invZperpZperpZperpZR₁' * o.S✻ZperpY₂
			o.S✻DGPZR₁DGPZR₁ = o.DGP.S✻ZR₁ZR₁ - _S✻ZperpDGPZR₁' * o.DGP.invZperpZperpZperpZR₁  - o.DGP.invZperpZperpZperpZR₁' * o.S✻ZperpDGPZR₁
			o.S✻DGPZR₁DGPZ   = o.DGP.S✻ZR₁Z   - _S✻ZperpDGPZR₁' * o.DGP.invZperpZperpZperpZpar - o.DGP.invZperpZperpZperpZR₁' * o.S✻ZperpDGPZ
			o.S✻DGPZR₁X      = _S✻⋂XDGPZR₁'   - _S✻ZperpDGPZR₁' * o.Repl.invZperpZperpZperpX   - o.DGP.invZperpZperpZperpZR₁' * o.S✻ZperpX
		end
	end

	o.invXXS✻XY₂   = @panelsum(o, o.invXXS✻⋂XY₂  , o.info✻_✻⋂)
	o.invXXS✻XX    = @panelsum(o, o.invXXS✻⋂XX   , o.info✻_✻⋂)
	o.invXXS✻XDGPZ = @panelsum(o, o.invXXS✻⋂XDGPZ, o.info✻_✻⋂)
	o.invXXS✻Xy₁   = dropdims(@panelsum(o, reshape(o.invXXS✻⋂Xy₁,Val(3)), o.info✻_✻⋂); dims=3)
	o.invZperpZperpS✻ZperpY₂   = o.Repl.invZperpZperp * o.S✻ZperpY₂
	o.invZperpZperpS✻ZperpX    = o.Repl.invZperpZperp * o.S✻ZperpX
	o.invZperpZperpS✻Zperpy₁   = o.Repl.invZperpZperp * o.S✻Zperpy₁
	o.invZperpZperpS✻ZperpDGPZ = o.Repl.invZperpZperp * o.S✻ZperpDGPZ

	if o.DGP.restricted
		o.invXXS✻XDGPZR₁ = @panelsum(o, o.invXXS✻⋂XDGPZR₁, o.info✻_✻⋂)
		o.invZperpZperpS✻ZperpDGPZR₁ = o.Repl.invZperpZperp * o.S✻ZperpDGPZR₁
	end
end

function PrepWRE!(o::StrBootTest{T}) where T
	r₁ = o.null ? [o.r₁ ; o.r] : o.r₁
  EstimateIV!(o.DGP, o, r₁)
  MakeResidualsIV!(o.DGP, o)
  o.robust && o.bootstrapt && o.granular && (Ü₂par = view(o.DGP.Ü₂ * o.Repl.RparY,:,:))

	o.S✻⋂XU₂ .= o.S✻⋂XY₂ - o.S✻⋂XX * o.DGP.Π̂
	o.S✻⋂XU₂RparY .= o.S✻⋂XU₂ * o.Repl.RparY
	o.S✻XU₂ .= o.S✻XY₂ - o.S✻XX * o.DGP.Π̂
	o.S✻XU₂RparY .= o.S✻XU₂ * o.Repl.RparY
	o.S✻ZperpU₂ .= o.S✻ZperpY₂ - o.S✻ZperpX * o.DGP.Π̂
	o.S✻ZperpU₂RparY .= o.S✻ZperpU₂ * o.Repl.RparY
	o.invZperpZperpS✻ZperpU₂ .= o.invZperpZperpS✻ZperpY₂ - o.invZperpZperpS✻ZperpX * o.DGP.Π̂
	o.invZperpZperpS✻ZperpU₂RparY .= o.invZperpZperpS✻ZperpU₂ * o.Repl.RparY
	o.invXXS✻XU₂ .= o.Repl.invXX * o.S✻XU₂
	o.invXXS✻XU₂RparY .= o.invXXS✻XU₂ * o.Repl.RparY

	if o.liml || !o.robust || !isone(o.κ)
		S✻U₂y₁ = o.S✻Y₂y₁ - o.DGP.Π̂' * o.S✻Xy₁
		S✻U₂RparYy₁ = o.Repl.RparY' * S✻U₂y₁
		S✻ZU₂ = o.S✻ReplZY₂ - o.S✻ReplZX * o.DGP.Π̂
		S✻ZU₂RparY = S✻ZU₂ * o.Repl.RparY
		Π̂S✻XÜ₂γ̈ = o.DGP.Π̂' * o.S✻XU₂ * o.DGP.γ̈
		S✻Ü₂Y₂ = o.S✻Y₂Y₂ - o.DGP.Π̂' * o.S✻XY₂
		S✻Y₂Ü₂γ̈ = S✻Ü₂Y₂' * o.DGP.γ̈
		S✻Ü₂parÜ₂par = o.Repl.RparY' * (S✻Ü₂Y₂ - o.S✻XU₂' * o.DGP.Π̂) * o.Repl.RparY
		S✻UUterm = o.S✻Y₂y₁ - o.S✻DGPZY₂' * o.DGP.β̈ - o.DGP.Π̂' * (o.S✻Xy₁  - o.S✻XDGPZ * o.DGP.β̈)

		if o.Repl.restricted
			S✻ReplZR₁r₁U₂ = o.r₁S✻ReplZR₁Y₂ - o.r₁S✻ReplZR₁X * o.DGP.Π̂
			S✻ReplZR₁r₁U₂RparY = S✻ReplZR₁r₁U₂ * o.Repl.RparY
		end
		if o.DGP.restricted
			r₁S✻DGPZR₁y₁ = r₁' * o.S✻DGPZR₁y₁
			S✻ReplZDGPZR₁r₁ = o.S✻ReplZDGPZR₁ * r₁
			o.Repl.restricted && (S✻ReplZR₁r₁DGPZR₁r₁ = o.r₁S✻ReplZR₁DGPZR₁ * r₁)
		end
	end

	if (o.liml || o.bootstrapt || !isone(o.κ)) && o.NFE>0
		CT✻FEU = o.CT✻FEY₂ - o.CT✻FEX * o.DGP.Π̂
		CT✻FEURparY = CT✻FEU * o.Repl.RparY
	end

  @inbounds for i ∈ 0:o.Repl.kZ  # precompute various clusterwise sums
		# S_✻(u .* X), S_✻(u .* Zperp) for residuals u for each endog var; store transposed
		if iszero(i)
			o.S✻XU[1]      .= o.S✻Xy₁      - o.S✻XDGPZ      * o.DGP.β̈ + o.S✻XU₂      * o.DGP.γ̈
			o.invXXS✻XU[1] .= o.invXXS✻Xy₁ - o.invXXS✻XDGPZ * o.DGP.β̈ + o.invXXS✻XU₂ * o.DGP.γ̈
			if o.DGP.restricted
				o.S✻XU[1]      .-= o.S✻XZR₁         * r₁
				o.invXXS✻XU[1] .-= o.invXXS✻XDGPZR₁ * r₁
			end
		else
			o.S✻XU[i+1]      .= view(o.S✻XU₂RparY,:,:,i)
			o.invXXS✻XU[i+1] .= view(o.invXXS✻XU₂RparY,:,:,i)
		end

		if o.liml || !isone(o.κ) || o.bootstrapt
			if iszero(i)
				o.S✻ZperpU[1]              .= o.S✻Zperpy₁              - o.S✻ZperpDGPZ              * o.DGP.β̈ + o.S✻ZperpU₂              * o.DGP.γ̈
				o.invZperpZperpS✻ZperpU[1] .= o.invZperpZperpS✻Zperpy₁ - o.invZperpZperpS✻ZperpDGPZ * o.DGP.β̈ + o.invZperpZperpS✻ZperpU₂ * o.DGP.γ̈
				if o.DGP.restricted
					o.S✻ZperpU[1]              .-= o.S✻ZperpDGPZR₁              * r₁
					o.invZperpZperpS✻ZperpU[1] .-= o.invZperpZperpS✻ZperpDGPZR₁ * r₁
				end
			else
				o.S✻ZperpU[i+1]              .= view(o.S✻ZperpU₂RparY,:,:,i)
				o.invZperpZperpS✻ZperpU[i+1] .= view(o.invZperpZperpS✻ZperpU₂RparY,:,:,i)
			end

			if o.NFE>0
				if iszero(i)
					o.CTFEU[1] = o.CT✻FEy₁ - o.CT✻FEZ * o.DGP.β̈ + CT✻FEU * o.DGP.γ̈
					o.DGP.restricted &&
						(o.CTFEU[1] .-= o.CT✻FEZR₁ * r₁)
				else
					o.CTFEU[i+1] = @view CT✻FEURparY[:,:,i]
				end
			end
		end

		if o.liml || !isone(o.κ) || !o.robust
			if iszero(i)  # panelsum2(o, o.Repl.y₁par, o.Repl.Z, uwt, o.info✻)
				o.S✻YU[1,1] .= o.S✻y₁y₁ - o.S✻DGPZy₁'o.DGP.β̈ + S✻U₂y₁'o.DGP.γ̈
				o.DGP.restricted &&
					(o.S✻YU[1,1] .-= view(r₁S✻DGPZR₁y₁, 1,:))

				if o.Repl.restricted
					o.S✻YU[1,1] .-= view(o.r₁S✻ReplZR₁y₁ - o.r₁S✻ReplZR₁DGPZ * o.DGP.β̈ + S✻ReplZR₁r₁U₂ * o.DGP.γ̈, 1,:,1)
					o.DGP.restricted &&
						(o.S✻YU[1,1] .+= view(S✻ReplZR₁r₁DGPZR₁r₁, 1,:,1))
				end

				for j ∈ 1:o.Repl.kZ
					o.S✻YU[j+1,1] .= view(o.S✻ReplZy₁ - o.S✻ReplZDGPZ * o.DGP.β̈ + S✻ZU₂ * o.DGP.γ̈, j,:,1)
					o.DGP.restricted &&
						(o.S✻YU[j+1,1] .-= view(S✻ReplZDGPZR₁r₁, j,:,1))
				end
      else
				o.S✻YU[1,i+1] .= @view S✻U₂RparYy₁[i,:]
				o.Repl.restricted &&
					(o.S✻YU[1,i+1] .-= @view S✻ReplZR₁r₁U₂RparY[1,:,i])
				for j ∈ 1:o.Repl.kZ
					o.S✻YU[j+1,i+1] .= @view S✻ZU₂RparY[j,:,i]
				end
			end

			if iszero(i)
		 		o.S✻UU[1,i+1] .= o.S✻y₁y₁ - (2 * o.S✻DGPZy₁ - o.S✻DGPZDGPZ * o.DGP.β̈)'o.DGP.β̈ + (2 * S✻UUterm - Π̂S✻XÜ₂γ̈ + S✻Y₂Ü₂γ̈)'o.DGP.γ̈
				o.DGP.restricted
					(o.S✻UU[1,i+1] .+= dropdims(-2 * r₁S✻DGPZR₁y₁ + r₁' * (o.S✻DGPZR₁DGPZR₁ * r₁) + 2 * r₁' * (o.S✻DGPZR₁DGPZ * o.DGP.β̈ + (o.S✻DGPZR₁X * o.DGP.Π̂ - o.S✻DGPZR₁Y₂) * o.DGP.γ̈); dims=1))  # XXX dropdims vs view...
			else
				o.S✻UU[1,i+1] .= view(o.Repl.RparY' * (S✻UUterm + S✻Y₂Ü₂γ̈ - Π̂S✻XÜ₂γ̈), i,:)

				if o.DGP.restricted
					o.S✻UU[1,i+1] .-= view(r₁' * (o.S✻DGPZR₁Y₂ - o.S✻DGPZR₁X * o.DGP.Π̂) * o.Repl.RparY, 1,:,i)
				end
			end
			for j ∈ 1:i
				o.S✻UU[j+1,i+1] .= view(S✻Ü₂parÜ₂par, j,:,i)
			end
		end

		if o.robust && o.bootstrapt && o.granular
			i>0 && (o.S✻UPX[i+1] .= o.XinvXX * o.S✻XU[i+1])

			o.S✻UMZperp[i+1] .= o.Repl.Zperp * o.invZperpZperpS✻ZperpU[i+1]
			if iszero(i)  # subtract crosstab of observation by ∩-group of u
				o.S✻UMZperp[   1][o.crosstab✻ind] .-= o.DGP.u⃛₁
			else
				o.S✻UMZperp[i+1][o.crosstab✻ind] .-= view(Ü₂par,:,i)
			end

			o.NFE>0 &&
				(o.S✻UMZperp[i+1] .-= view(o.invFEwt .* o.CTFEU[i+1], o._FEID, :))  # CT_(*,FE) (U ̈_(parj) ) (S_FE S_FE^' )^(-1) S_FE
		end
  end

	if o.robust && o.bootstrapt && iszero(o.granular)
		@inbounds for j ∈ 0:o.Repl.kZ
			if o.Repl.Yendog[j+1]
				o.negS✻UMZperpX[j+1] = o.S⋂XZperpinvZperpZperp * o.S✻ZperpU[j+1]  # S_* diag(U ̈_(∥j) ) Z_⊥ (Z_⊥^' Z_⊥ )^(-1) Z_(⊥g)^' X_(∥g)
				if iszero(j)  # - S_*  diag(U ̈_(∥j) ) I_g^' X_(∥g)
					o.negS✻UMZperpX[j+1][o.crosstab⋂✻ind] .-= vec(o.S✻⋂Xy₁ - o.S✻⋂XDGPZ * o.DGP.β̈ + o.S✻⋂XU₂ * o.DGP.γ̈)
					o.DGP.restricted &&
						(o.negS✻UMZperpX[j+1][o.crosstab⋂✻ind] .+= vec(o.S✻⋂X_DGPZR₁ * r₁))
				else
					o.negS✻UMZperpX[j+1][o.crosstab⋂✻ind] .-= vec(view(o.S✻⋂XU₂RparY,:,:,j))
				end
				o.NFE>0 &&
					(o.negS✻UMZperpX[j+1] .+= o.CT⋂FEX'o.CTFEU[j+1])  # CT_(*,FE) (U ̈_(∥j) ) (S_FE S_FE^' )^(-1) S_FE
			end
		end
	end
	nothing
end

# For WRE, and with reference to Y = [y₁ Z], given 0-based columns indexes within it, i, j, return all bootstrap realizations of
# Y[:,i]'((1-κ)*M_Zperp-κ*M_Xpar)*Y[:,j] for κ constant across replications
# i can be a rowvector
# (only really the Hessian when we narrow Y to Z)
function HessianFixedkappa(o::StrBootTest{T}, is::Vector{S} where S<:Integer, j::Integer, κ::Number, w::Integer) where T
  dest = Matrix{T}(undef, length(is), ncols(o.v))
  @inbounds for i ∈ eachindex(is, axes(dest,1))
		_HessianFixedkappa!(o, view(dest,i:i,:), is[i], j, κ, w)
  end
  dest
end

function _HessianFixedkappa!(o::StrBootTest, dest::AbstractMatrix, i::Integer, j::Integer, κ::Number, w::Integer)
  if !(o.Repl.Yendog[i+1] || o.Repl.Yendog[j+1])  # if both vars exog, result = order-0 term only, same for all draws
		!iszero(κ) &&
			(dest .= dot(view(o.Repl.XZ,:,i), view(o.Repl.invXXXZ,:,j)))
		if !isone(κ)
			if iszero(κ)
				fill!(dest, o.Repl.YY[i+1,j+1])
			else
				dest .= κ .* dest .+ (1 - κ) .* o.Repl.YY[i+1,j+1]
			end
		end
	else
		if !iszero(κ)  # repititiveness in this section to maintain type stability
			if o.Repl.Yendog[i+1]
				T1L = o.T1L[isone(o.Nw) || w<o.Nw ? 1 : 2]  # preallocated destinations
				mul!(T1L, o.S✻XU[i+1], o.v)
				if iszero(i)
					T1L .+= o.Repl.Xy₁par
				else
					T1L .+= view(o.Repl.XZ,:,i)
				end
				if o.Repl.Yendog[j+1]
					T1R = o.T1R[isone(o.Nw) || w<o.Nw ? 1 : 2]  # preallocated destinations
					mul!(T1R, o.invXXS✻XU[j+1], o.v)
					if iszero(j)
						T1R .+=  o.Repl.invXXXy₁par
					else
						T1R .+= view(o.Repl.invXXXZ,:,j)
					end
					coldot!(o, dest, T1L, T1R)
				else
					dest .= view(o.Repl.invXXXZ,:,j)'T1L  # coldot!(o, dest, T1L, view(o.Repl.invXXXZ,:,j))
				end
			else
				if o.Repl.Yendog[j+1]
					T1R = o.T1R[isone(o.Nw) || w<o.Nw ? 1 : 2]  # use preallocated destinations
					mul!(T1R, o.invXXS✻XU[j+1], o.v)
					if iszero(j)
						T1R .+=  o.Repl.invXXXy₁par
					else
						T1R .+= view(o.Repl.invXXXZ,:,j)
					end
					dest .= view(o.Repl.XZ,:,i)'T1R
				else
					dest .= dot(view(o.Repl.XZ,:,i), view(o.Repl.invXXXZ,:,j))
				end
			end
		end
		if !isone(κ)
			if o.Repl.Yendog[i+1]
				o.T2 .= o.invZperpZperpS✻ZperpU[i+1]'o.S✻ZperpU[j+1]  # quadratic term
				o.T2[diagind(o.T2)] .-= i ≤ j ? o.S✻UU[i+1, j+1] : o.S✻UU[j+1, i+1]  # minus diagonal crosstab
				o.NFE>0 &&
					(o.T2 .+= o.CTFEU[i+1]' * (o.invFEwt .* o.CTFEU[j+1]))
				if iszero(κ)
					dest .= o.Repl.YY[i+1,j+1] .+ colquadformminus!(o, (                            o.S✻YU[i+1,j+1] .+ o.S✻YU[j+1,i+1])'o.v, o.T2, o.v)
				else
					dest .= κ .* dest .+ (1 - κ)     .* colquadformminus!(o, (o.Repl.YY[i+1,j+1] .+ o.S✻YU[i+1,j+1] .+ o.S✻YU[j+1,i+1])'o.v, o.T2, o.v)
				end
			elseif iszero(κ)
				dest .= o.Repl.YY[i+1,j+1]
			else
				dest .= κ .* dest .+ (1 - κ) .* o.Repl.YY[i+1,j+1]
			end
		end
  end
	nothing
end

# put threaded loops in functions to prevent type instability https://discourse.julialang.org/t/type-inference-with-threads/2004/3
function FillingLoop1!(o::StrBootTest{T}, dest::Matrix{T}, i::Integer, j::Integer, _β̈::AbstractMatrix{T}, β̈v::AbstractMatrix{T}) where T
	Threads.@threads for g ∈ 1:o.N⋂
		PXY✻ = [o.PXZ[g,i]]
		o.Repl.Yendog[i+1] && (PXY✻ = PXY✻ .+ view(o.S✻UPX[i+1],g,:)'o.v)

		if iszero(j)
			dest[g,:]   = dropdims(colsum(PXY✻ .* (o.Repl.y₁[g] .- view(o.S✻UMZperp[1],g,:)'o.v)); dims=1)
		elseif o.Repl.Yendog[j+1]
			dest[g,:] .-= dropdims(colsum(PXY✻ .* (o.Repl.Z[g,j] * _β̈ .- view(o.S✻UMZperp[j+1],g,:)'β̈v)); dims=1)
		else
			dest[g,:] .-= dropdims(colsum(PXY✻ .* (o.Repl.Z[g,j] * _β̈)); dims=1)
		end
	end
	nothing
end
function FillingLoop2!(o::StrBootTest{T}, dest::Matrix{T}, i::Integer, j::Integer, _β̈::AbstractMatrix{T}, β̈v::AbstractMatrix{T}) where T
	Threads.@threads for g ∈ 1:o.N⋂
		S = o.info⋂[g]
		PXY✻ = o.Repl.Yendog[i+1] ? view(o.PXZ,S,i) .+ view(o.S✻UPX[i+1],S,:) * o.v :
												 reshape(view(o.PXZ,S,i), :, 1)

		if iszero(j)
			dest[i,:]   = dropdims(colsum(PXY✻ .* (o.Repl.y₁[S] .- view(o.S✻UMZperp[1],S,:) * o.v)); dims=1)
		else
			dest[i,:] .-= dropdims(colsum(PXY✻ .* (o.Repl.Yendog[j+1] ? o.Repl.Z[S,j] * _β̈ .- view(o.S✻UMZperp[j+1],S,:) * β̈v :
																														       o.Repl.Z[S,j] * _β̈                                       )); dims=1)
		end
	end
	nothing
end

# Workhorse for WRE CRVE sandwich filling
# Given a zero-indexed column index, i>0, and a matrix β̈s of all the bootstrap estimates,
# return all bootstrap realizations of P_X * Z[:,i]_g ' û₁g^*b
# for all groups in the intersection of all error clusterings
# return value has one row per ⋂ cluster, one col per bootstrap replication
# that is, given i, β̈s = δ ̂_CRκ^(*), return, over all g, b (P_(X_∥ g) Z_(∥i)^(*b) )^' (M_(Z_⊥ ) y_(1∥)^(*b) )_g-(P_(X_∥ g) Z_(∥i)^(*b) )^' (M_(Z_⊥ ) Z_∥^(*b) )_g δ ̂_CRκ^(*b)
function Filling!(o::StrBootTest{T}, dest::AbstractMatrix{T}, i::Int64, β̈s::AbstractMatrix) where T
	if o.granular
		β̈v = _β̈ = Matrix{T}(undef,0,0)
   	if o.Nw == 1  # create or avoid NxB matrix?
			PXY✻ = reshape(view(o.PXZ,:,i), :, 1)
			o.Repl.Yendog[i+1] && (PXY✻ = PXY✻ .+ o.S✻UPX[i+1] * o.v)
			dest .= @panelsum(o, PXY✻ .* (o.Repl.y₁ .- o.S✻UMZperp[1] * o.v), o.info⋂)
			@inbounds for j ∈ 1:o.Repl.kZ
				_β̈ = view(β̈s,j,:)'
				dest .-= @panelsum(o, PXY✻ .* (o.Repl.Yendog[j+1] ? view(o.Repl.Z,:,j) * _β̈ .- o.S✻UMZperp[j+1] * (o.v .* _β̈) :
															                                 (view(o.Repl.Z,:,j) * _β̈)                                      ), o.info⋂)
			end
		else  # create pieces of each N x B matrix one at a time rather than whole thing at once
			@inbounds for j ∈ 0:o.Repl.kZ
				j>0 && (β̈v = o.v .* (_β̈ = view(β̈s,j,:)'))
				(o.purerobust ? FillingLoop1! : FillingLoop2!)(o, dest, i, j, _β̈, β̈v)
			end
    end
  else  # coarse error clustering
		# (P_(X_∥ g) Z_∥^* )^' (M_(Z_⊥ ) y_(1∥)^* )_g
		F1₀ = view(o.Repl.invXXXZ,:,i)
		F1₁ = o.invXXS✻XU[i+1]
		F2₀ = o.S⋂Xy₁
		F2₁ = o.negS✻UMZperpX[1]
		if o.Repl.Yendog[i+1]  # add terms that are zero only if Zpar[i] is exogenous, i.e. if a null refers only to exogenous variables
			dest .= reshape(F1₀'F2₀,:,1) .- (dropdims(F1₀'F2₁; dims=1) - F2₀'F1₁) * o.v  # 0th- & 1st-order terms
			o.Q .= F1₁'F2₁
			@inbounds for g ∈ 1:o.N⋂
				o.colquadformminus!(dest, g, o.v, o.Q[:,g,:], o.v)
			end
		else
			dest .= F2₀'F1₀ .- dropdims(F1₀'F2₁; dims=1) * o.v  # 0th- & 1st-order terms
		end

		# -(P_(X_∥ g) Z_∥^* )^' (M_(Z_⊥ ) Z_∥^(*b) )_g δ ̂_CRκ^*
		@inbounds for j ∈ 1:o.Repl.kZ
			F2₀ = view(o.S⋂ReplZX,j,:,:)'
			β̈v = o.β̈v[isone(o.Nw) || w<o.Nw ? 1 : 2]; β̈v .= o.v .* (_β̈ = -view(β̈s,j,:)')
			if o.Repl.Yendog[j+1]
				F2₁ = o.negS✻UMZperpX[j+1]
				dest .+= F2₀'F1₀ .* _β̈ .- (dropdims(F1₀'F2₁; dims=1) - F2₀'F1₁) * β̈v  # "-" because S✻UMZperpX is stored negated as F2₁=negS✻UMZperpX[j+1]
				o.Q .= F1₁'F2₁
				for g ∈ 1:o.N⋂
					o.colquadformminus!(dest, g, o.v, o.Q[:,g,:], β̈v)
				end
			elseif o.Repl.Yendog[i+1]
				dest .+= reshape(F1₀'F2₀,:,1) .* _β̈  .+ F2₀'F1₁ * β̈v
			else
				dest .+= reshape(F1₀'F2₀,:,1) .* _β̈
			end
		end
  end
  nothing
end

function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
  if isone(o.Repl.kZ)  # optimized code for 1 coefficient in bootstrap regression
		if o.liml
			YY₁₁   = HessianFixedkappa(o, [0], 0, zero(T), w)  # κ=0 => Y*MZperp*Y
			YY₁₂   = HessianFixedkappa(o, [0], 1, zero(T), w)
			YY₂₂   = HessianFixedkappa(o, [1], 1, zero(T), w)
			YPXY₁₁ = HessianFixedkappa(o, [0], 0, one(T) , w)  # κ=1 => Y*PXpar*Y
			YPXY₁₂ = HessianFixedkappa(o, [0], 1, one(T) , w)
			YPXY₂₂ = HessianFixedkappa(o, [1], 1, one(T) , w)
			YY₁₂YPXY₁₂ = YY₁₂ .* YPXY₁₂
			x₁₁ = YY₂₂ .* YPXY₁₁ .- YY₁₂YPXY₁₂      # elements of YY✻^-1 * YPXY✻ up to factor of det(YY✻)
			x₁₂ = YY₂₂ .* YPXY₁₂ .- YY₁₂ .* YPXY₂₂
			x₂₁ = YY₁₁ .* YPXY₁₂ .- YY₁₂ .* YPXY₁₁
			x₂₂ = YY₁₁ .* YPXY₂₂ .- YY₁₂YPXY₁₂
			κs = (x₁₁ .+ x₂₂)./2; κs .= 1 ./ (1 .- (κs .- sqrtNaN.(κs.^2 .- x₁₁ .* x₂₂ .+ x₁₂ .* x₂₁)) ./ (YY₁₁ .* YY₂₂ .- YY₁₂ .* YY₁₂))  # solve quadratic equation for smaller eignenvalue; last term is det(YY✻)
			!iszero(o.fuller) && (κs .-= o.fuller / (o._Nobs - o.kX))
			o.As = κs .* (YPXY₂₂ .- YY₂₂) .+ YY₂₂
			o.β̈s = (κs .* (YPXY₁₂ .- YY₁₂) .+ YY₁₂) ./ o.As
		else
			o.As = HessianFixedkappa(o, [1], 1, o.κ, w)
			o.β̈s = HessianFixedkappa(o, [1], 0, o.κ, w) ./ o.As
		end

		if o.null
			o.numerw = o.β̈s .+ (o.Repl.Rt₁ - o.r) / o.Repl.RRpar
		else
			o.numerw = o.β̈s .- o.DGP.β̈
			isone(w) && (o.numerw[1] = o.β̈s[1] + (o.Repl.Rt₁[1] - o.r[1]) / o.Repl.RRpar[1])
		end

		@storeWtGrpResults!(o.numer, o.numerw)
		if o.bootstrapt
			if o.robust
				J⋂s1 = dropdims(o.J⋂s[isone(o.Nw) || w<o.Nw ? 1 : 2]; dims=3)
				Filling!(o, J⋂s1, 1, o.β̈s);
				J⋂s1 ./= o.As
				@inbounds for c ∈ 1:o.NErrClustCombs  # sum sandwich over error clusterings
					nrows(o.clust[c].order)>0 &&
						(J⋂s1 .= J⋂s1[o.clust[c].order,:])
					@clustAccum!(denom, c, coldot(o, @panelsum(o, J⋂s1, o.clust[c].info)))
				end
			else
				denom = (HessianFixedkappa(o, [0], 0, zero(T), w) .- 2 .* o.β̈s .* HessianFixedkappa(o, [0], 1, zero(T), w) .+ o.β̈s.^2 .* HessianFixedkappa(o, [1], 1, zero(T), w)) ./ o._Nobs ./ o.As  # classical error variance
			end
			@storeWtGrpResults!(o.dist, o.sqrt ? o.numerw ./ sqrtNaN.(denom) : o.numerw .^ 2 ./ denom)
			denom *= o.Repl.RRpar[1]^2
		end
		w==1 && o.bootstrapt && (o.statDenom = fill(denom[1],1,1))  # original-sample denominator

  else  # WRE bootstrap for more than 1 retained coefficient in bootstrap regression

		β̈s = zeros(T, o.Repl.kZ, ncols(o.v))
		A = Vector{Matrix{T}}(undef, ncols(o.v))

		if o.liml
			YY✻   = [HessianFixedkappa(o, collect(0:i), i, zero(T), w) for i ∈ 0:o.Repl.kZ] # κ=0 => Y*MZperp*Y
			YPXY✻ = [HessianFixedkappa(o, collect(0:i), i,  one(T), w) for i ∈ 0:o.Repl.kZ] # κ=1 => Y*PXpar*Y

			@inbounds for b ∈ axes(o.v,2)
				for i ∈ 0:o.Repl.kZ
					o.YY✻_b[1:i+1,i+1]   = YY✻[i+1][:,b]  # fill uppper triangles, which is all that invsym() looks at
					o.YPXY✻_b[1:i+1,i+1] = YPXY✻[i+1][:,b]
				end
				o.κ = 1/(1 - real(eigvals(invsym(o.YY✻_b) * Symmetric(o.YPXY✻_b))[1]))
				!iszero(o.fuller) && (o.κ -= o.fuller / (o._Nobs - o.kX))
				β̈s[:,b] = (A[b] = invsym(o.κ*o.YPXY✻_b[2:end,2:end] + (1-o.κ)*o.YY✻_b[2:end,2:end])) * (o.κ*o.YPXY✻_b[1,2:end] + (1-o.κ)*o.YY✻_b[1,2:end])
			end
		else
			δnumer =  HessianFixedkappa(o, collect(1:o.Repl.kZ), 0, o.κ, w)
			δdenom = [HessianFixedkappa(o, collect(1:i), i, o.κ, w) for i ∈ 1:o.Repl.kZ]

			nt = Threads.nthreads()
  		cs = [round(Int, size(o.v,2)/nt*i) for i ∈ 0:nt]
			@inbounds Threads.@threads for t ∈ 1:nt
				δdenom_b = zeros(T, o.Repl.kZ, o.Repl.kZ)  # thread-safe scratch pad
				@inbounds for b ∈ cs[t]+1:cs[t+1]
					for i ∈ 1:o.Repl.kZ
						δdenom_b[1:i,i] = view(δdenom[i],:,b)  # fill uppper triangle
					end
					β̈s[:,b] = (A[b] = invsym(δdenom_b)) * view(δnumer,:,b)
				end
			end
		end

		if o.bootstrapt
			if o.robust
				J⋂s = o.J⋂s[isone(o.Nw) || w<o.Nw ? 1 : 2]
				@inbounds for i ∈ 1:o.Repl.kZ  # avoid list comprehension construction for compiler-perceived type stability
					Filling!(o, view(J⋂s,:,:,i), i, β̈s)
				end
			else
				YY✻ = [HessianFixedkappa(o, collect(i:o.Repl.kZ), i, zero(T), w) for i ∈ 0:o.Repl.kZ]  # κ=0 => Y*MZperp*Y
			end
		end

		@inbounds for b ∈ reverse(axes(o.v,2))
			o.numer_b .= o.null || w==1 && b==1 ? o.Repl.RRpar * view(β̈s,:,b) + o.Repl.Rt₁ - o.r : o.Repl.RRpar * (view(β̈s,:,b) - o.DGP.β̈₀)
			if o.bootstrapt
				if o.robust  # Compute denominator for this WRE test stat
					J⋂ = view(J⋂s,:,b,:) * (A[b] * o.Repl.RRpar')
					for c ∈ 1:o.NErrClustCombs
						(!isone(o.NClustVar) && nrows(o.clust[c].order)>0) &&
							(J⋂ = J⋂[o.clust[c].order,:])
						J_b = @panelsum(o, J⋂, o.clust[c].info)
						@clustAccum!(denom, c, J_b'J_b)
					end
				else  # non-robust
					for i ∈ 0:o.Repl.kZ
						o.YY✻_b[i+1,i+1:o.Repl.kZ+1] = view(YY✻[i+1],:,b)  # fill upper triangle
					end
					denom = (o.Repl.RRpar * A[b] * o.Repl.RRpar') * [-one(T) ; β̈s[:,b]]'Symmetric(o.YY✻_b) * [-one(T) ; β̈s[:,b]] / o._Nobs  # 2nd half is sig2 of errors
				end
				o.dist[b+first(o.WeightGrp[w])-1] = o.sqrt ? o.numer_b[1] / sqrtNaN(denom[1]) : o.numer_b'invsym(denom)*o.numer_b  # hand-code for 2-dimensional?
			end
			o.numer[:,b+first(o.WeightGrp[w])-1] = o.numer_b  # slight inefficiency: in usual bootstrap-t case, only need to save numerators in numer if getdist("numer") is coming because of svmat(numer)
		end
		w==1 && o.bootstrapt && (o.statDenom = denom)  # original-sample denominator
	end
	nothing
end


text:
export iauFk5hip
"""
FK5 to Hipparcos rotation and spin.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Returned:
   r5h   double[3][3]  r-matrix: FK5 rotation wrt Hipparcos (Note 2)
   s5h   double[3]     r-vector: FK5 spin wrt Hipparcos (Note 3)

Notes:

   1. This function models the FK5 to Hipparcos transformation as a
   pure rotation and spin;  zonal errors in the FK5 catalogue are
   not taken into account.

   2. The r-matrix r5h operates in the sense:

         P_Hipparcos = r5h x P_FK5

   where P_FK5 is a p-vector in the FK5 frame, and P_Hipparcos is
   the equivalent Hipparcos p-vector.

   3. The r-vector s5h represents the time derivative of the FK5 to
   Hipparcos rotation.  The units are radians per year (Julian,
   TDB).

Called:
   iauRv2m      r-vector to r-matrix

Reference:

   F.Mignard & M.Froeschle, Astron.Astrophys., 354, 732-739 (2000).

This revision:  2017 October 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauFk5hip(double r5h[3][3], double s5h[3])

function iauFk5hip()
   r5h = zeros(Float64, 3, 3)
   s5h = zeros(Float64, 3)

   ccall((:iauFk5hip, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
          r5h, s5h)

   return SMatrix{3,3}(r5h'), SVector{3}(s5h)
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 67477abe-0b56-49a6-818d-3fccc669c6ce
begin
	# using Pkg; Pkg.add("Primes")
	using Primes
end

# ╔═╡ 727bf063-5fef-4126-a615-8a483aab9c99
begin
	# using Pkg; Pkg.add("Plots")
	using Plots
end

# ╔═╡ 80fc74b2-e038-11eb-39c8-fbd84caf7338
md"
# Fun with Arithmetic Derivatives

In this notebook I use Julia and Pluto to entertain the concept of the [Arithmetic Derivative](https://en.wikipedia.org/wiki/Arithmetic_derivative), also known as **number derivative**. It is named due to its axiomatic definition, which shares some of the properties of [real derivatives](https://en.wikipedia.org/wiki/Derivative), except it is a function acting on numbers.

**Definition (Arithmetic Derivative).** _Let $n\in \mathbb{N}$. Then $D:\mathbb{N}\rightarrow\mathbb{N}$ is a function with the following properties:_

* _$D(p)=1$ for any prime $p$, and_

* _$D(pq)=D(p)q+D(q)p$ for any $p,q\in\mathbb{N}$._

Thus we can generalize the arithmetic derivative to the set of natural numbers, ie. if we take $p=q=0$, we obtain that $D(0)=0$, and if we take $p=q=1$ we also obtain that $D(1)=0$. Thus we can calculate the arithmetic derivative for any non-prime natural by means of prime factorization, using the \"power rule\"

$$D(p^n)=np^{n-1}D(p)$$

for $p,n\in\mathbb{N}$.

**Proof (Power rule).** _The \"power rule\" can be proven in a straight-forward manner by means of mathematical induction. We take $n=0$ and see easily that the equation holds, for $D(1)=0$. Next we assume that the equation holds for a given $n$ and apply the induction step $n\rightarrow n+1$:_

$\begin{align}
D(p^{n+1})&=D(p^np)\\
&=D(p^n)p+p^nD(p)\\
&=np^{n-1}D(p)p + p^nD(p)\\
&=np^{n}D(p) + p^nD(p)\\
&=(n+1)p^{n}D(p)
\end{align}$

_Thus the proof is concluded._

Therefore we conclude that the arithmetic derivatie of any natural number $q$ greater than $1$ can be determined using its prime factorization, as in

$$D(q)=D\left(\prod_{j=1}^Np_j^{m_j}\right),$$

where $p_1,...,p_N$ are the prime factors of $q$ and $m_1,...,m_N$ their multiplicities.

Now we have all the building blocks to implement the arithmetic derivative in Julia. To this end, we make use of the Julia package `Primes.jl`, which provides optimized routines to calculate prime numbers and perform prime factorization of `Integer` values.
"

# ╔═╡ 8b83eccb-c629-45ff-b8cb-4b4879664ace
md"
Next we implement $D(p)$ as the Julia function `arithmetic_derivative(::Integer)`. We also implement the generalization upon the set of whole numbers $\mathbb{Z}$, which will be useful later. This generalization rests upon the fact that we consider $D(-q)=-D(q)$. This statement can be proven, but we will simply pretend that the proof is trivial, as in the prime factorization of a negative integer simply yields a factor of $-1$ multiplied with the other (positive) prime factors.
"

# ╔═╡ c0d8d582-364c-4b06-8b08-cf902a359897
function arithmetic_derivative(n::Integer)
	if (n < 0)
		return -arithmetic_derivative(-n)
	end

	if (n < 2)
		return 0
	elseif (isprime(n))
		return 1
	end

	fac = factor(n)
	sum = 0

	map(collect(fac)) do pair
		m = pair.second
		p = pair.first
		sum += m * p^(m-1) * (n ÷ p^m)
	end
	sum
end

# ╔═╡ fa3e00d3-7f8e-4158-98ee-4e22285952d7
md"
Thus we have implemented the basic algorithm that determines the arithmetic derivative of an integer number. Indeed, the arithmetic derivative is sometimes generalized to the set of rational numbers $\mathbb{Q}$, and in this generalization the special case of a negative integer is indeed treated that way (see later).

To test our implementation, we let Julia compute `arithmetic_derivative(n)` for $0\le n \le 25$. We make sure that the values calculated are correct by calculating them by hand ([see this handwritten document](https://cdn.calza.xyz/pics/notebooks/maths/arithmetic_derivatives/Primes_Test.pdf) to check the results). While the standard output is not shown in Pluto, the reader is welcome to download the notebook and run it locally to double-check the test output with the results calculated by hand.
"

# ╔═╡ 49eea834-b3bf-426b-993b-d2ea01da742d
md"
At this point it would be interesting to show the values in a `scatter` plot. Indeed, the arithmetic derivative is a very interesting construction from a number theory standpoint since it is a reasonably simble and yet fun axiomatic construction, however it is difficult to visualize. A plot will hopefully shed more light on the nature of the function $D(p)$.

To this end, we plot the first $80$ values for $D(p)$ against their corresponding values $p$:
"

# ╔═╡ 5fbd5946-9e1e-444b-9ab6-03bfed725620
md"
This is perhaps less useful than we had hoped. We can try to connect the dots to visualize the path of the function $D(n)$:
"

# ╔═╡ 5166fe7d-50fc-4e0e-b839-c0db8bf73b01
md"
This makes a moderate amount of sense because we expect some kind of oscillation as the derivative is bound to rise as $n$ grows, but will keep falling to $1$ as primes are encountered every now and again in order to rise again. The \"amplitude\" of the oscillation (as in, the maximum) is also expected to rise as larger natural numbers are likelier to have a longer array of prime factors and higher multiplicities, which means that the sum due to the power rule can reach higher values.

Julia is more than capable of doing better, though, such that we can venture into higher bounds. Let's try larger values of $N$ (note that the code cells containing the calls to `plot_derivatives(n, connect_dots = true)` have been hidden for aesthetic purposes):
"

# ╔═╡ a3a14843-7b62-4f6f-8e81-ed772e4a27f3
md"
These plots share in the character described earlier, but fail to provide further insights. Interestingly enough, [a generalization over the rational numbers $\mathbb{Q}$ is given](https://en.wikipedia.org/wiki/Arithmetic_derivative#Extensions_beyond_natural_numbers) by the well-known quotient rule

$$D\left(\frac pq\right)=\frac{D(p)q-D(q)p}{q^2}$$

for any rational number expressed as the fraction $\frac pq$. We implement the function `arithmetic_derivative(::Rational)` which calls `arithmetic_derivative(::Integer)`:
"

# ╔═╡ 8af8ac5e-616d-407b-8001-6fa7529c5729
begin
	function arithmetic_derivative(r::Rational)
		p = numerator(r)
		q = denominator(r)
		(arithmetic_derivative(p) * q - arithmetic_derivative(q) * p) / q^2
	end
end

# ╔═╡ 3c9f9459-9828-47ec-bdfc-8bd8a0187553
for n in 0:25
	println("n=$n\tD(n)=$(arithmetic_derivative(n))")
end

# ╔═╡ 1f2295eb-188f-4498-9d8c-1c0f303b2396
function plot_derivatives(N; connect_dots = false, plot_label = nothing, new_plot = true, dots_color = :lightblue)
	mode_func = nothing
	if (connect_dots)
		mode_func = new_plot ? plot : plot!
	else
		mode_func = new_plot ? scatter : scatter!
	end

	if (plot_label === nothing)
		mode_func(0:N, arithmetic_derivative.(0:N), legend=false, color = dots_color)
	else
		mode_func(0:N, arithmetic_derivative.(0:N), label=plot_label, legend=:topleft, color = dots_color)
	end

	xlabel!("Natural number n")
	ylabel!("Number derivative D(n)")
	title!("Plot of the first $N derivatives")
end

# ╔═╡ 9ebc2997-8db4-4aeb-a2e1-4614302d4f56
plot_derivatives(80)

# ╔═╡ f9b2636f-231c-44f4-95a3-013ba47df516
plot_derivatives(80, connect_dots = true)

# ╔═╡ 652bf146-8059-4c52-8f9a-63c3b5a22487
plot_derivatives(100, connect_dots = true)

# ╔═╡ 9cbc71f0-d902-4276-8a3d-b025355dddef
plot_derivatives(200, connect_dots = true)

# ╔═╡ 9c06d571-24e9-4678-9c0f-a7ac01f7d8f6
plot_derivatives(500, connect_dots = true)

# ╔═╡ 0374c7ce-c485-4712-a95b-55e2f3c30a67
plot_derivatives(3000, connect_dots = true)

# ╔═╡ c3248ade-6106-47df-8e6f-38847aee52ce
md"
We can plot this function and overlay it with the scatter plot for $N=100$:
"

# ╔═╡ 2ae73301-1458-4d49-80f4-249706d2eb09
begin
	# Calculate array of rationals to use:

	N_max = 101
	rationals = zeros(Rational, N_max^2 + 1)
	for i in 0:N_max, j in 1:N_max
		rationals[i * (N_max-1) + j] = i//j
	end
	sort!(rationals)

	# Plot the rational derivatives:

	plot(rationals, arithmetic_derivative.(rationals), label = "Rational generalization", color = :lightblue)

	# ----------------------------
	# Overlay integer derivatives:

	plot_derivatives(100, connect_dots = false, plot_label = "Integer derivatives", new_plot = false, dots_color = :lightgreen)
end

# ╔═╡ 7f491dc8-d921-4637-acc9-09ccc3d2b9d6
md"
With this interesting plot we conclude our small exploration.

About this very last plot a few questions could be asked, for example why does the line belonging to the rational generalization become less dense after $n=50$? Could it be due to the way the `rationals` array is generated? If so, are the rationals in the array all equidistant from one another? Or does it depend on the nature of the generalized function $D(p/q)$?
"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"

[compat]
Plots = "~1.18.1"
Primes = "~0.5.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.6+5"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.0+6"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
git-tree-sha1 = "ed268efe58512df8c7e224d2e170afd76dd6a417"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.13.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "3cc57ad0a213808473eafef4845a74766242e05f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.3.1+4"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.1+14"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.1+5"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+0"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "b83e3125048a9c3158cbb7ca423790c7b1b57bea"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.57.5"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e14907859a1d3aee73a019e7b3c98e9e7b8b5b3e"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.57.3+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "15ff9a14b9e1218958d3530cc288cf31465d9ae2"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.3.13"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "47ce50b742921377301e15005c96e979574e130b"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.1+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.12"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[LibVPX_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "12ee7e23fa4d18361e7c2cde8f8337d4c3101bc7"
uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a"
version = "1.10.0+0"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.6"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "501c20a63a34ac1d015d5304da0e645f42d91c9f"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.11"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "b93181645c1209d912d5632ba2d0094bc00703ad"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.18.1"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Primes]]
git-tree-sha1 = "afccf037da52fa596223e5a0e331ff752e0e845c"
uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
version = "0.5.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "2a7a2469ed5d94a98dea0e85c46fa653d76be0cd"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.3.4"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "a43a7b58a6e7dc933b2fa2e0ca653ccf8bb8fd0e"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.6"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.8"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "000e168f5cc9aded17b6999a560b7c11dda69095"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.14.0+4"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "0.1.6+4"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d713c1ce4deac133e3334ee12f4adff07f81778f"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2020.7.14+2"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "487da2f8f2f0c8ee0e83f39d13037d6bbf0a45ab"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.0.0+3"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─80fc74b2-e038-11eb-39c8-fbd84caf7338
# ╠═67477abe-0b56-49a6-818d-3fccc669c6ce
# ╟─8b83eccb-c629-45ff-b8cb-4b4879664ace
# ╠═c0d8d582-364c-4b06-8b08-cf902a359897
# ╟─fa3e00d3-7f8e-4158-98ee-4e22285952d7
# ╠═3c9f9459-9828-47ec-bdfc-8bd8a0187553
# ╟─49eea834-b3bf-426b-993b-d2ea01da742d
# ╠═727bf063-5fef-4126-a615-8a483aab9c99
# ╠═1f2295eb-188f-4498-9d8c-1c0f303b2396
# ╠═9ebc2997-8db4-4aeb-a2e1-4614302d4f56
# ╟─5fbd5946-9e1e-444b-9ab6-03bfed725620
# ╠═f9b2636f-231c-44f4-95a3-013ba47df516
# ╟─5166fe7d-50fc-4e0e-b839-c0db8bf73b01
# ╟─652bf146-8059-4c52-8f9a-63c3b5a22487
# ╟─9cbc71f0-d902-4276-8a3d-b025355dddef
# ╟─9c06d571-24e9-4678-9c0f-a7ac01f7d8f6
# ╟─0374c7ce-c485-4712-a95b-55e2f3c30a67
# ╟─a3a14843-7b62-4f6f-8e81-ed772e4a27f3
# ╠═8af8ac5e-616d-407b-8001-6fa7529c5729
# ╟─c3248ade-6106-47df-8e6f-38847aee52ce
# ╠═2ae73301-1458-4d49-80f4-249706d2eb09
# ╟─7f491dc8-d921-4637-acc9-09ccc3d2b9d6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using PkgSRA, Plots, LinearAlgebra, Statistics
using BSON: @load
pyplot()

#####
##### Load variables
#####
include("paper_settings.jl")
this_dat_name = DAT_FOLDERNAME*"TMP_dat_flowchart_inkscape_"

# Uncontrolled data
fname = this_dat_name*"uncontrolled_lorenz.bson"
@load fname dat_raw numerical_grad_raw

# Original ODE and controller variables
fname = this_dat_name*"ode_vars.bson";
@load fname dat dat_grad true_grad U_true

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_sindy.bson";
@load fname sindy_dat_unctr #sindy_unctr

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_bayes.bson";
@load fname dat_unctr sample_trajectory_noise residual accepted_ind ctr_guess

# SINDY variables 2: after control
fname = this_dat_name*"ctr_vars_sindy.bson";
@load fname sindy_dat_ctr sindy_grad_ctr #sindy_sub

# Bayesian variables 2: after control
fname = this_dat_name*"ctr_vars_turing.bson";
@load fname dat_ctr sample_trajectory_noise_ctr sample_gradients_ctr turing_dat_ctr ctr_guess2

#####
##### Plot everything
#####
this_fig_name = FIGURE_FOLDERNAME*"TMP_fig_flowchart_inkscape_"
plot_opt = Dict(:titlefontsize=>48,
        :xticks=>false, :yticks=>false, :zticks=>false,
        :legend=>false, :fontfamily=>:serif)

## 0: Raw, uncontrolled system
plot_raw = plot3d(dat_raw[1, :], dat_raw[2, :], dat_raw[3, :],
        color=COLOR_DICT["data_uncontrolled"], lw=4; plot_opt...)
        title!("Unforced Lorenz")
fname = this_fig_name * "uncontrolled.png";
savefig(plot_raw, fname)

## 1:  3d example plot: data
plot_attractor = plot3d(dat[1, :], dat[2, :], dat[3, :],
        color=COLOR_DICT["true"], lw=3; plot_opt...)
        title!("Perturbed Data")
fname = this_fig_name * "data.png";
savefig(plot_attractor, fname)

## 2: Two 1d example plots: naive model
uncontrolled_ind = 1:500;
    ind = uncontrolled_ind;
let d = sindy_dat_unctr
        global plot_uncontrolled1 = plot(d[1, ind],
                color=COLOR_DICT["model_uncontrolled"],
                alpha=0.8, lw=4;
                plot_opt...)
end
    plot!(dat[1, ind], color=COLOR_DICT["true"], lw=4)
    ylabel!("X", guidefontsize=34)
    title!("Naive Model")
let d = sindy_dat_unctr
    global plot_uncontrolled2 = plot(d[2, ind],
            color=COLOR_DICT["model_uncontrolled"],
            alpha=0.8, lw=4;
            plot_opt...)
end
    plot!(dat[2, ind], color=COLOR_DICT["true"], lw=4)
    ylabel!("dY/dt", guidefontsize=34)

fname = this_fig_name * "uncontrolledX.png";
savefig(plot_uncontrolled1, fname)
fname = this_fig_name * "uncontrolledY.png";
savefig(plot_uncontrolled2, fname)

## 3: Residual with noise envelope (x and y coordinates)
plot_coordinate = 1;
    residual_ind = 1:500;
    ind = residual_ind;
plot_residual1 = plot(residual[plot_coordinate,ind],
                    ribbon=mean(sample_trajectory_noise),
                    fillalpha=0.5,
                    color=COLOR_DICT["residual"], lw=4;
                    plot_opt...)
    ylabel!("dX/dt", guidefontsize=34)
    title!("Residual")
plot_coordinate = 2;
plot_residual2 = plot(residual[plot_coordinate,ind],
                ribbon=mean(sample_trajectory_noise),
                fillalpha=0.5,
                color=COLOR_DICT["residual"], lw=4;
                plot_opt...)
    ylabel!("dY/dt", guidefontsize=24)

fname = this_fig_name * "residualX.png";
savefig(plot_residual1, fname)
fname = this_fig_name * "residualY.png";
savefig(plot_residual2, fname)

## 4: Subsampling
plot_coordinate = 1;
    sub_ind = 1:500;
    ind = sub_ind;
plot_subsample = plot(dat[plot_coordinate,ind],
                    color=COLOR_DICT["true"], lw=4;
                    plot_opt...)
    ylabel!("X", guidefontsize=34)
    title!("Subsample")
    subsample_ind = intersect(ind, accepted_ind)
    scatter!(ind[subsample_ind], dat[plot_coordinate,subsample_ind],
            color=COLOR_DICT["intrinsic"],
            markershape=:o, markersize=20;
            plot_opt...)
fname = this_fig_name * "subsample.png";
savefig(plot_subsample, fname)

## 5: Partial model (same ind as above)
plot_coordinate = 1;
    ind = 1:500
plot_partial = plot(sindy_grad_ctr[plot_coordinate,ind],
                    color=COLOR_DICT["model_partial"], lw=8;
                    plot_opt...)
    plot!(dat_grad[plot_coordinate,ind],
        color=COLOR_DICT["true"], lw=4;
        plot_opt...)
    ylabel!("dX/dt", guidefontsize=34)
    title!("Partial Model")
fname = this_fig_name * "partial.png";
savefig(plot_partial, fname)

## 6: Control signal guess (x coordinate)
# Same ind as above
plot_ctr_guess = plot(ctr_guess[plot_coordinate,ind],
                    color=COLOR_DICT["control_time"], lw=6;
                    plot_opt...)
    plot!(U_true[plot_coordinate,ind],
            color=COLOR_DICT["true"], lw=3;
            plot_opt...)
    title!("Control Signal")
fname = this_fig_name * "control_guess.png";
savefig(plot_ctr_guess, fname)

## 7: "Revealed underlying system"
controlled_ind = 1:5000;
    ind = controlled_ind;
let d = sindy_dat_ctr
    global plot_controlled = plot3d(
                d[1, ind], d[2, ind], d[3, ind],
                color=COLOR_DICT["model_controlled"], alpha=0.8, lw=3;
                plot_opt...)
end
    # plot3d!(dat[1, ind], dat[2, ind], dat[3, ind],
    #         color=COLOR_DICT["true"], lw=2)
    title!("Learned System")
fname = this_fig_name * "controlled.png";
savefig(plot_controlled, fname)


text:
module FileWatcher

using FileWatching

const watchedFolders = Set{String}()

basePath = "."

function getFilesRecursively(path::String)
    files = String[]
    folders = String[]

    if isdir(path)
        for file in readdir(path)
            filename = joinpath(path, file)

            if isdir(filename)
                subFiles, subFolders = getFilesRecursively(filename)

                push!(folders, filename)
                append!(folders, subFolders)
                append!(files, subFiles)

            elseif isfile(filename)
                push!(files, filename)
            end
        end
    end

    return files, folders
end

function attemptToWatchFolder(path::String, timeout::Number=0.0)
    if isdir(path) && !(path in watchedFolders)
        try
            watch_folder(path, timeout)
            push!(watchedFolders, path)

        catch e
            println(Base.stderr, "Failed to watch $path")

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println()
            end
        end
    end
end

# Returns all added/updated filenames since last call
# This is good enough for us, removed files don't need to be handled
function processWatchEvents(basePath::String, recursive::Bool=true, timeout::Number=0.0)
    attemptToWatchFolder(basePath, timeout)

    addedFiles = Set{String}()
    removedFolders = Set{String}()

    for path in watchedFolders
        try
            # TODO - Check for better way for unwatching
            if !ispath(path)
                unwatch_folder(path)
                push!(removedFolders, path)

                continue
            end

        catch
            # Do nothing, check again next time
        end

        while true
            relative, event = watch_folder(path, timeout)

            if event.timedout
                # No changes since last check

                break
            end

            filename = joinpath(path, relative)

            if event.changed
                # Check if file was removed or added
                # Folders add all sub files recursively
                # Watch folder as well if recursive

                if isdir(filename)
                    files, folders = getFilesRecursively(filename)

                    if recursive
                        watch_folder(filename, timeout)
                        push!(watchedFolders, filename)

                        for folder in folders
                            attemptToWatchFolder(folder, timeout)
                        end
                    end

                    union!(addedFiles, files)

                elseif isfile(filename)
                    push!(addedFiles, filename)
                end
            end
        end

        setdiff!(watchedFolders, removedFolders)
    end

    return addedFiles
end

function watchAllFoldersRecursively(basePath::String, timeout::Number=0.0)
    files, folders = getFilesRecursively(basePath)

    for folder in folders
        watch_folder(folder, timeout)
        push!(watchedFolders, folder)
    end
end

function initFileWatcher(targetPath::String)
    global basePath = targetPath

    # Call once to start watching
    watchAllFoldersRecursively(basePath)
    processWatchEvents(basePath)
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# set true if using scip initialization, false if using IPOPT initialization
const scip_init = true

if(scip_init)
	const initialObjValue = 48.27745824188737
	user_initial_guess = Float64[0.0 0.0 0.0 0.0 43.52912917750612 0.0 0.0 0.0 6.345178399087325e-16 0.0 0.0 25.078984684205707 0.0 0.0 0.0 0.0 15.524911271438539 0.0 0.0 0.0]
else
	const initialObjValue = 49.31927629671476
	user_initial_guess = Float64[8.00743603082653e-8 7.353973701611809e-8 9.394197568065582e-8 8.708630221162697e-8 43.082467219608795 8.523171541639657e-8 6.581006027607289e-8 6.202145514365337e-8 12.44034130710884 6.677873136922174e-8 6.433461920978504e-8 24.843886953698174 9.062880795235767e-8 8.015355067589583e-8 9.655886628481626 7.96120708055623e-8 1.0557601429767548e-7 1.1308307032517315e-7 6.715543849744971e-8 7.926991356692931e-8]
end

const objGap = 0.005*initialObjValue

const riskLowerBound = 5E-04

user_initial_guess = user_initial_guess'

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using CSDP
const CSDPOPT = with_optimizer(CSDP.Optimizer, printlevel = 0)

@testset "CSDP SDP triangle Problems" begin
    list_of_sdp_triang_problems = [
        # sdpt1_test, # CSDP is returning SLOW_PROGRESS
        # sdpt2_test, # CSDP is returning SLOW_PROGRESS
        sdpt3_test,
        sdpt4_test
    ]
    test_strong_duality(list_of_sdp_triang_problems, CSDPOPT)
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractHyperPlane end

abstract type HyperPlaneType end
abstract type OptimalityCut <: HyperPlaneType end
abstract type FeasibilityCut <: HyperPlaneType end
abstract type LinearConstraint <: HyperPlaneType end
abstract type Infeasible <: HyperPlaneType end
abstract type Unbounded <: HyperPlaneType end

struct HyperPlane{H <: HyperPlaneType, T <: AbstractFloat, A <: AbstractVector} <: AbstractHyperPlane
    δQ::A
    q::T
    id::Int

    function HyperPlane(δQ::AbstractVector, q::AbstractFloat, id::Int, ::Type{H}) where H <: HyperPlaneType
        T = promote_type(eltype(δQ), Float32)
        δQ_ = convert(AbstractVector{T}, δQ)
        new{H, T, typeof(δQ_)}(δQ_, q, id)
    end
end
OptimalityCut(δQ::AbstractVector, q::AbstractFloat, id::Int) = HyperPlane(δQ, q, id, OptimalityCut)
FeasibilityCut(δQ::AbstractVector, q::AbstractFloat, id::Int) = HyperPlane(δQ, q, id, FeasibilityCut)
LinearConstraint(δQ::AbstractVector, q::AbstractFloat, id::Int) = HyperPlane(δQ, q, id, LinearConstraint)
Infeasible(Q::AbstractFloat, id::Int) = HyperPlane(sparsevec(Float64[]), Q, id, Infeasible)
Unbounded(Q::AbstractFloat, id::Int) = HyperPlane(sparsevec(Float64[]), Q, id, Unbounded)

struct AggregatedOptimalityCut{T <: AbstractFloat, A <: AbstractVector} <: AbstractHyperPlane
    δQ::A
    q::T
    ids::Vector{Int}

    function AggregatedOptimalityCut(δQ::AbstractVector, q::AbstractFloat, ids::Vector{Int})
        T = promote_type(eltype(δQ), Float32)
        δQ_ = convert(AbstractVector{T}, δQ)
        new{T, typeof(δQ_)}(δQ_, q, ids)
    end
end

SparseFeasibilityCut{T <: AbstractFloat} = HyperPlane{FeasibilityCut, T, SparseVector{T,Int64}}
SparseLinearConstraint{T <: AbstractFloat} = HyperPlane{LinearConstraint, T, SparseVector{T,Int64}}
AnyOptimalityCut{T, A} = Union{HyperPlane{OptimalityCut, T, A}, AggregatedOptimalityCut{T, A}}
SparseOptimalityCut{T <: AbstractFloat} = HyperPlane{OptimalityCut, T, SparseVector{T,Int64}}
SparseAggregatedOptimalityCut{T <: AbstractFloat} = AggregatedOptimalityCut{T, SparseVector{T,Int64}}
AnySparseOptimalityCut{T <: AbstractFloat} = Union{SparseOptimalityCut, SparseAggregatedOptimalityCut{T}}
SparseHyperPlane{T <: AbstractFloat} = Union{HyperPlane{<:HyperPlaneType, T, SparseVector{T,Int64}}, SparseAggregatedOptimalityCut{T}}
QCut{T} = Tuple{Int,SparseHyperPlane{T}}
CutQueue{T} = RemoteChannel{Channel{QCut{T}}}

function (hyperplane::HyperPlane{FeasibilityCut})(x::AbstractVector)
    return Inf
end
function (cut::HyperPlane{OptimalityCut})(x::AbstractVector)
    if length(cut.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match", length(cut.δQ), length(x))))
    end
    return cut.q-cut.δQ⋅x
end
function (cut::AggregatedOptimalityCut)(x::AbstractVector)
    if length(cut.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match", length(cut.δQ), length(x))))
    end
    return cut.q-cut.δQ⋅x
end
function (hyperplane::HyperPlane{Infeasible})(x::AbstractVector)
    return hyperplane.q
end
function (hyperplane::HyperPlane{Unbounded})(x::AbstractVector)
    return hyperplane.q
end

infeasible(hyperplane::AbstractHyperPlane) = false
infeasible(hyperplane::HyperPlane{Infeasible}) = true
bounded(hyperplane::AbstractHyperPlane) = true
bounded(hyperplane::HyperPlane{Unbounded}) = false
function optimal(cut::AnyOptimalityCut, x::AbstractVector, θ::AbstractFloat, τ::AbstractFloat)
    Q = cut(x)
    return θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
end
function active(hyperplane::AbstractHyperPlane, x::AbstractVector, τ::AbstractFloat)
    return abs(gap(hyperplane,x)) <= τ
end
function satisfied(hyperplane::AbstractHyperPlane, x::AbstractVector, τ::AbstractFloat)
    return gap(hyperplane,x) >= -τ
end
function satisfied(cut::HyperPlane{OptimalityCut}, x::AbstractVector, θ::AbstractFloat, τ::AbstractFloat)
    Q = cut(x)
    return θ > -Inf && θ >= Q - τ
end
function gap(hyperplane::AbstractHyperPlane, x::AbstractVector)
    if length(hyperplane.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match", length(hyperplane.δQ), length(x))))
    end
    return hyperplane.δQ⋅x-hyperplane.q
end
function gap(cut::AnyOptimalityCut, x::AbstractVector, θ::AbstractFloat)
    if θ > -Inf
        return θ-cut(x)
    else
        return Inf
    end
end
function num_subproblems(::AbstractHyperPlane)
    return 1
end
function num_subproblems(cut::AggregatedOptimalityCut)
    return length(cut.ids)
end

moi_constraint(hyperplane::AbstractHyperPlane, master_variables::Vector{MOI.VariableIndex}) = moi_constraint(hyperplane, master_variables, 1.)
function moi_constraint(hyperplane::HyperPlane{H,T,SparseVector{T,Int}}, ::Vector{MOI.VariableIndex}, scaling::T) where {H <: HyperPlaneType, T <: AbstractFloat}
    terms = map(zip(hyperplane.δQ.nzval, hyperplane.δQ.nzind)) do (coeff, idx)
        MOI.ScalarAffineTerm(scaling * coeff, MOI.VariableIndex(idx))
    end
    f = MOI.ScalarAffineFunction{Float64}(terms, 0.0)
    set = MOI.GreaterThan{Float64}(scaling * hyperplane.q)
    return f, set
end
function moi_constraint(cut::HyperPlane{FeasibilityCut,T,SparseVector{T,Int}}, ::Vector{MOI.VariableIndex}, scaling::T) where T <: AbstractFloat
    terms = map(zip(cut.δQ.nzval, cut.δQ.nzind)) do (coeff, idx)
        MOI.ScalarAffineTerm(scaling * coeff, MOI.VariableIndex(idx))
    end
    f = AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, zero(T)),
                               MOI.ScalarAffineFunction(terms, zero(T)))
    set = MOI.GreaterThan{Float64}(scaling * cut.q)
    return f, set
end
function moi_constraint(cut::HyperPlane{OptimalityCut,T,SparseVector{T,Int}}, master_variables::Vector{MOI.VariableIndex}, scaling::T) where T <: AbstractFloat
    terms = map(zip(cut.δQ.nzval, cut.δQ.nzind)) do (coeff, idx)
        MOI.ScalarAffineTerm(scaling * coeff, MOI.VariableIndex(idx))
    end
    # Add model term
    model_terms = [MOI.ScalarAffineTerm(scaling, master_variables[cut.id])]
    f = AffineDecisionFunction(MOI.ScalarAffineFunction(model_terms, zero(T)),
                               MOI.ScalarAffineFunction(terms, zero(T)))
    set = MOI.GreaterThan{Float64}(scaling * cut.q)
    return f, set
end
function moi_constraint(cut::SparseAggregatedOptimalityCut, master_variables::Vector{MOI.VariableIndex}, scaling::T) where T <: AbstractFloat
    terms = map(zip(cut.δQ.nzval, cut.δQ.nzind)) do (coeff, idx)
        MOI.ScalarAffineTerm(scaling * coeff, MOI.VariableIndex(idx))
    end
    # Add model terms
    model_terms = map(cut.ids) do idx
        MOI.ScalarAffineTerm(scaling, master_variables[idx])
    end
    f = AffineDecisionFunction(MOI.ScalarAffineFunction(model_terms, zero(T)),
                               MOI.ScalarAffineFunction(terms, zero(T)))
    set = MOI.GreaterThan{Float64}(scaling * cut.q)
    return f, set
end

function aggregate(cuts::Vector{HyperPlane{OptimalityCut,T,A}}, id::Integer) where {T <: AbstractFloat, A <: AbstractVector}
    return HyperPlane(sum([cut.δQ for cut in cuts]),
                      sum([cut.q for cut in cuts]),
                      id,
                      OptimalityCut)
end

function aggregate(cuts::Vector{<:AbstractHyperPlane})
    return sum(cuts)
end

struct ZeroVector{T <: AbstractFloat} <: AbstractVector{T} end
length(::ZeroVector) = 0
size(::ZeroVector) = 0
ZeroOptimalityCut{T <: AbstractFloat} = AggregatedOptimalityCut{T, ZeroVector{T}}

function Base.zero(::Type{AggregatedOptimalityCut{T}}) where T <: AbstractFloat
    return AggregatedOptimalityCut(ZeroVector{T}(), zero(T), Int[])
end
Base.iszero(::AggregatedOptimalityCut) = false
Base.iszero(::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat = true

function +(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return c₁
end
function +(c₁::HyperPlane{OptimalityCut, T, A}, c₂::HyperPlane{OptimalityCut, T, A}) where {T <: AbstractFloat, A <: AbstractVector}
    return AggregatedOptimalityCut(c₁.δQ + c₂.δQ, c₁.q + c₂.q, vcat(c₁.id, c₂.id))
end
function +(c₁::AggregatedOptimalityCut{T, A}, c₂::AggregatedOptimalityCut{T, A}) where {T <: AbstractFloat, A <: AbstractVector}
    return AggregatedOptimalityCut(c₁.δQ + c₂.δQ, c₁.q + c₂.q, vcat(c₁.ids, c₂.ids))
end
function +(c₁::AggregatedOptimalityCut{T, A}, c₂::HyperPlane{OptimalityCut, T, A}) where {T <: AbstractFloat, A <: AbstractVector}
    return AggregatedOptimalityCut(c₁.δQ + c₂.δQ, c₁.q + c₂.q, vcat(c₁.ids, c₂.id))
end
function +(c₁::HyperPlane{OptimalityCut, T, A}, c₂::AggregatedOptimalityCut{T, A}) where {T <: AbstractFloat, A <: AbstractVector}
    return AggregatedOptimalityCut(c₁.δQ + c₂.δQ, c₁.q + c₂.q, vcat(c₁.id, c₂.ids))
end
function +(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::AggregatedOptimalityCut{T}) where T <: AbstractFloat
    return c₂
end
function +(c₁::AggregatedOptimalityCut{T}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return c₁
end
function +(c₁::AggregatedOptimalityCut{T,ZeroVector{T}}, c₂::HyperPlane{OptimalityCut,T}) where T <: AbstractFloat
    return AggregatedOptimalityCut(c₂.δQ, c₂.q, [c₂.id])
end
function +(c₁::HyperPlane{OptimalityCut,T}, c₂::AggregatedOptimalityCut{T,ZeroVector{T}}) where T <: AbstractFloat
    return AggregatedOptimalityCut(c₁.δQ, c₁.q, [c₁.id])
end


text:
function solve(prob::Problem, settings::Settings)
    stats = Stats(settings)
    prob, equil = equilibrate(prob, max_iters=settings.prob_equil_iters)
    term_cache = TermCache()
    scaling = CompositeConeScaling(prob.C)
    kkt = KktCache(prob, settings, scaling)
    iter = get_init_iter(prob, kkt, settings)
    best_iter = copy(iter)
    print_first_line(stats, settings)
    for stats.iters in 1:settings.max_iters
        term = iterate!(iter, best_iter, term_cache, kkt, prob, scaling, settings, stats)
        if term
            break
        end
        print_stats_iter(stats, settings)
    end
    sol = unpack_sol(best_iter, prob, stats, settings)
    unscale!(sol, equil)
    truncate_stats!(stats)
    return sol, stats
end

function iterate!(iter::Iterate, best_iter::Iterate, term_cache::TermCache,
                  kkt::KktCache, prob::Problem,
                  scaling::CompositeConeScaling, settings::Settings, stats::Stats)

    # Step 1.
    resid = get_residual(iter, prob, settings)
    if check_term_conds!(iter, best_iter, term_cache, resid, prob, settings, stats)
        return true
    end

    # Step 2.
    update!(scaling, iter.s, iter.z)
    λ = get_scaled_var(scaling)
    update!(kkt, iter, prob, scaling, settings)
    cache_first_solve!(kkt, prob, settings, stats)
    aff_dir = get_aff_dir!(iter, resid, λ, prob, kkt, scaling, settings, stats)

    # Step 3.
    α1 = get_max_step(iter, aff_dir, λ, prob, scaling)
    σ = get_center_param(α1)
    σ = max(σ, settings.min_centering)

    # Step 4.
    comb_dir = get_comb_dir!(iter, resid, λ, prob, kkt, scaling, aff_dir, settings, σ, stats)

    # Step 5.
    α2 = settings.step_param*get_max_step(iter, comb_dir, λ, prob, scaling)
    α2 = min(α2, settings.max_step)
    if α2 > settings.min_step && is_valid(iter, prob) && is_valid(comb_dir) && isfinite(α2)
        take_step!(iter, comb_dir, α2, prob)
    end

    update_stats!(stats, iter, σ, α1, α2)
    return false
end

function update_stats!(stats::Stats{T}, iter::Iterate{T}, σ::T, α1::T, α2::T) where T
    i = stats.iters
    stats.σ[i] = σ
    stats.α_aff[i] = α1
    stats.α_comb[i] = α2
    stats.κ[i] = iter.κ
    stats.τ[i] = iter.τ
end

function take_step!(iter::Iterate{T}, dir::StepDirection{T},
                    α::T, prob::Problem{T}) where T <: AbstractFloat
    if iter.s + α * dir.Δs ∈ prob.C && iter.z + α * dir.Δz ∈ prob.C
        iter.x = iter.x + α * dir.Δx
        iter.y = iter.y + α * dir.Δy
        iter.z = iter.z + α * dir.Δz
        iter.s = iter.s + α * dir.Δs
        iter.τ = iter.τ + α * dir.Δτ
        iter.κ = iter.κ + α * dir.Δκ
    end
end

get_center_param(α::AbstractFloat) = (1 - α)^3

function get_max_step(iter::Iterate{T}, dir::StepDirection{T}, λ::Vector{T},
                      prob::Problem{T}, scaling::CompositeConeScaling{T}) where T
    Δs_tilde = scaling \ dir.Δs
    Δz_tilde = scaling * dir.Δz
    α = min(get_max_step(λ, Δs_tilde, Δz_tilde, prob.C),
            get_max_step([iter.τ], [dir.Δτ], NonnegCone(1)),
            get_max_step([iter.κ], [dir.Δκ], NonnegCone(1)))
    return α
end

function get_residual(iter::Iterate{T}, prob::Problem{T}, settings::Settings{T}) where T
    rx = -prob.A' * iter.y - prob.G' * iter.z - prob.c * iter.τ
    ry = prob.A * iter.x - prob.b * iter.τ
    rz = iter.s + prob.G * iter.x - prob.h * iter.τ
    rτ = iter.κ + prob.c' * iter.x + prob.b' * iter.y + prob.h' * iter.z
    return Residual{T}(rx, ry, rz, rτ)
end

function get_init_iter(prob::Problem{T}, kkt::KktCache{T}, settings::Settings{T}) where T
    x_hat, s_hat = get_init_xs(prob, kkt::KktCache)
    y_hat, z_hat = get_init_yz(prob, kkt::KktCache)
    #if s_hat ∈ prob.C && z_hat ∈ prob.C
    if false
        iter = Iterate(x_hat, y_hat, z_hat, s_hat)
    else
        iter = Iterate(prob)
    end
    return iter
end

function get_init_xs(prob::Problem{T}, kkt::KktCache{T}) where T <: AbstractFloat
    n = get_num_vars(prob)
    (x_hat, _, s_tilde), _, _, _, _ = solve_and_refine(kkt, [zeros(T, n), prob.b, prob.h])
    e = convert(Vector{T}, get_id(prob.C))
    αp = dist_along_id_ray(s_tilde, prob.C)
    s_hat = s_tilde + (1 + 1.1*αp).*e
    return x_hat, s_hat
end

function get_init_yz(prob::Problem{T}, kkt::KktCache{T}) where T <: AbstractFloat
    _, m, p = get_prob_sizes(prob)
    (_, y_hat, z_tilde), _, _, _, _ = solve_and_refine(kkt, [-prob.c, zeros(T, m), zeros(T, p)])
    e = convert(Vector{T}, get_id(prob.C))
    αd = dist_along_id_ray(z_tilde, prob.C)
    z_hat = z_tilde + (1 + αd).*e
    return y_hat, z_hat
end


text:
# NEW MACROS: @planar and @plansor
macro planar(ex::Expr)
    return esc(planar_parser(ex))
end

function planar_parser(ex::Expr)
    parser = TO.TensorParser()

    braidingtensors = Vector{Any}()

    pop!(parser.preprocessors) # remove TO.extracttensorobjects
    push!(parser.preprocessors, _conj_to_adjoint)
    push!(parser.preprocessors, _extract_tensormap_objects)
    push!(parser.preprocessors, _construct_braidingtensors)
    treebuilder = parser.contractiontreebuilder
    treesorter = parser.contractiontreesorter
    push!(parser.preprocessors, ex->TO.processcontractions(ex, treebuilder, treesorter))
    push!(parser.preprocessors, ex->_check_planarity(ex))
    temporaries = Vector{Symbol}()
    push!(parser.preprocessors, ex->_decompose_planar_contractions(ex, temporaries))

    pop!(parser.postprocessors) # remove TO.addtensoroperations
    push!(parser.postprocessors, ex->_update_temporaries(ex, temporaries))
    push!(parser.postprocessors, ex->_annotate_temporaries(ex, temporaries))
    push!(parser.postprocessors, _add_modules)
    return parser(ex)
end

macro plansor(ex::Expr)
    return esc(plansor_parser(ex))
end

function plansor_parser(ex)
    inputtensors = TO.getinputtensorobjects(ex)
    newtensors = TO.getnewtensorobjects(ex)

    # find the first non-braiding tensor to determine the braidingstyle
    targetobj = inputtensors[findfirst(x->x != :τ, inputtensors)]
    targetsym = gensym()

    ex = TO.replacetensorobjects(ex) do obj, leftind, rightind
        obj == targetobj ? targetsym : obj
    end

    defaultparser = TO.TensorParser()
    insert!(defaultparser.preprocessors, 3, _remove_braidingtensors)
    defaultex = defaultparser(ex)
    planarex = planar_parser(ex)

    ex = Expr(:block)
    push!(ex.args, Expr(:(=), targetsym, targetobj))
    push!(ex.args, :(if BraidingStyle(sectortype($targetsym)) isa Bosonic
                $(defaultex)
            else
                $(planarex)
            end))
    if targetobj in newtensors
        push!(ex.args, Expr(:(=), targetobj, targetsym))
        push!(ex.args, newtensors[end])
    end
    return ex
end