text:
# OLS, ML (score bootstrap), and Anderson-Rubin tests

# Construct stuff that depends linearly or quadratically on r, possibly by interpolation
function MakeInterpolables!(o::StrBootTest{T}) where T
	if o.interpolable
		if iszero(nrows(o.anchor))  # first call? save current r as anchor for interpolation
			o.anchor = o.r
			_MakeInterpolables!(o, o.anchor)
			o.numer₀ = o.numer
			o.interpolate_u && (o.ü₀ = copy(o.ü))
			o.robust && (o.Jcd₀ = deepcopy(o.Jcd))
			return
		end

		if iszero(nrows(o.poles))  # second call: from anchor make set of orthogonal poles, which equal anchor except in one dimension
			o.poles = o.r - o.anchor
			o.robust && (o.denom₀ = deepcopy(o.denom))  # grab quadratic denominator from *previous* (1st) evaluation
			newPole = trues(o.q)  # all poles new
		else  # been here at least twice? interpolate unless current r stretches range > 2X in some dimension(s)
			newPole = abs.(o.r .- o.anchor) .> T(2) .* abs.(o.poles)
		end

		if any(newPole)  # prep interpolation
			@inbounds for h₁ ∈ 1:o.q
				if newPole[h₁]
					o.poles[h₁] = o.r[h₁] - o.anchor[h₁]
					thisr = copy(o.anchor); thisr[h₁] = o.r[h₁]  # if q>1 this creates anchor points that are not graphed, an inefficiency. But simpler to make the deviations from 1st point orthogonal
					_MakeInterpolables!(o, thisr)  # calculate linear stuff at new anchor

					o.∂numer∂r[h₁] = (o.numer .- o.numer₀) ./ o.poles[h₁]
					o.interpolate_u && (o.∂u∂r[h₁] = (o.ü .- o.ü₀) ./ o.poles[h₁])
					if o.robust  # dof > 1 for an ARubin test with >1 instruments.
						for d₁ ∈ 1:o.dof
							for c ∈ 1:o.NErrClustCombs
								o.∂Jcd∂r[h₁][c,d₁] = (o.Jcd[c,d₁] .- o.Jcd₀[c,d₁]) ./ o.poles[h₁]
								for d₂ ∈ 1:d₁
									tmp = coldot(o, o.Jcd₀[c,d₁], o.∂Jcd∂r[h₁][c,d₂])
									d₁ ≠ d₂ && (o.coldotplus!(tmp, o.Jcd₀[c,d₂], o.∂Jcd∂r[h₁][c,d₁]))  # for diagonal items, faster to just double after the c loop
									@clustAccum!(o.∂denom∂r[h₁][d₁,d₂], c, tmp)
								end
							end
							o.∂denom∂r[h₁][d₁,d₁] .*= T(2)  # double diagonal terms
						end
					end
				end
			end
			if o.robust  # quadratic interaction terms
				@inbounds for h₁ ∈ 1:o.q, h₂ ∈ 1:h₁
					if newPole[h₁] || newPole[h₂]
						for d₁ ∈ 1:o.dof, d₂ ∈ 1:d₁, c ∈ 1:o.NErrClustCombs
							@clustAccum!(o.∂²denom∂r²[h₁,h₂][d₁,d₂], c, coldot(o, o.∂Jcd∂r[h₁][c,d₁], o.∂Jcd∂r[h₂][c,d₂]))
						end
					end
				end
			end
			Δ = o.poles
			o.interpolating = true

			if o.q==2  # in this case we haven't yet actually computed interpolables at r, so interpolate them
				o.numerw .= o.numer₀ .+ o.∂numer∂r[1] .* Δ[1] .+ o.∂numer∂r[2] .* Δ[2]
				if o.interpolate_u
					o.ü = o.ü₀ .+ o.∂u∂r[1] .* Δ[1] .+ o.∂u∂r[2] .* Δ[2]
				end
			end

		else  # routine linear interpolation if the anchors not moved
			Δ = o.r - o.anchor
			o.numerw .= o.numer₀ .+ o.∂numer∂r[1] .* Δ[1]
			o.q > 1 && (o.numerw .+= o.∂numer∂r[2] .* Δ[2])
			if o.interpolate_u
				o.ü .= o.ü₀ .+ o.∂u∂r .* Δ[1]
				o.q > 1 && (o.ü .+= o.∂u∂r[2] .* Δ[2])
			end
		end

		if o.robust  # even if an anchor was just moved, and linear components just computed from scratch, do the quadratic interpolation now, from the updated linear factors
			if isone(o.q)
				@inbounds for d₁ ∈ 1:o.dof, d₂ ∈ 1:d₁
					o.denom[d₁,d₂] .= o.denom₀[d₁,d₂] .+ o.∂denom∂r[d₁,d₂][1,1] .* Δ .+ o.∂²denom∂r²[d₁,d₂][1,1] .* Δ.^2
				end
			else  # q==2
				@inbounds for d₁ ∈ 1:o.dof, d₂ ∈ 1:d₁
					o.denom[d₁,d₂] .= o.denom₀[d₁,d₂] .+
														o.∂denom∂r[1][d₁,d₂] .* Δ[1] .+
														o.∂denom∂r[2][d₁,d₂] .* Δ[2] .+
														o.∂²denom∂r²[1,1][d₁,d₂] .* (Δ[1] .^ 2) .+
														o.∂²denom∂r²[2,1][d₁,d₂] .* (Δ[1] .* Δ[2]) .+
														o.∂²denom∂r²[2,2][d₁,d₂] .* (Δ[2] .^ 2)
				end
			end
		end
	else  # non-interpolable cases
		_MakeInterpolables!(o, o.r)
	end
	nothing
end

# Construct stuff that depends linearly or quadratically on r and doesn't depend on v. No interpolation.
function _MakeInterpolables!(o::StrBootTest{T}, thisr::AbstractVector) where T
	if o.ML
		o.uXAR = o.sc * (o.AR = o.A * o.R')
	else
		if o.ARubin
			EstimateARubin!(o.DGP, o, thisr)
			MakeResidualsOLSARubin!(o.DGP, o)
		elseif iszero(o.κ)  # regular OLS
			EstimateOLS!(o.DGP, o.null ? [o.r₁ ; thisr] : o.r₁)
			MakeResidualsOLSARubin!(o.DGP, o)
		elseif o.null  # in score bootstrap for IV/GMM, if imposing null, then DGP constraints, κ, Hessian, etc. do vary with r and must be set now
			EstimateIV!(o.DGP, o, [o.r₁ ; thisr])
			InitTestDenoms!(o.DGP, o)
			MakeResidualsIV!(o.DGP, o)
		end

		o.ü = o.DGP.ü₁

		(o.scorebs || (o.robust && o.granular < o.NErrClustCombs)) &&
			(o.uXAR = o.DGP.ü₁ .* o.M.XAR)
	end

	o.SuwtXA = o.scorebs ?
				o.B>0 ?
					 o.NClustVar ?
				          @panelsum(o, o.uXAR, o.wt, o.info✻) :
					      vHadw(o.uXAR, o.wt)                    :
				        o.haswt ? reshape(o.uXAR'o.wt,1,:) : sum(o.uXAR,dims=1)  :
			  o.DGP.A * panelsum2(o, o.X₁, o.X₂, vHadw(o.ü, o.wt), o.info✻)'  # same calc as in score BS but broken apart to grab intermediate stuff, and assuming residuals defined; X₂ empty except in Anderson-Rubin

	if o.robust && o.bootstrapt && o.granular < o.NErrClustCombs
		u✻XAR = @panelsum(o, o.uXAR, o.wt, o.info✻⋂)  # collapse data to all-boot && error-cluster-var intersections. If no collapsing needed, panelsum() will still fold in any weights
		if o.B>0
			if o.scorebs
				K = [zeros(T, o.clust[1].N, o.N✻) for _ in o.dof]::Vector{Matrix{T}}  # inefficient, but not optimizing for the score bootstrap
			else
				K = [panelsum2(o, o.X₁, o.X₂, vHadw(view(o.DGP.XAR,:,d), o.wt), o.info⋂) * o.SuwtXA for d ∈ 1:o.dof]::Vector{Matrix{T}}
			end

			o.NFE>0 && !o.FEboot && (o.CT_WE = crosstabFE(o, vHadw(o.ü, o.wt), o.info✻))

			if o.NFE>0 && !o.FEboot
				tmp = o.invFEwt .* o.CT_WE
				@inbounds for d ∈ 1:o.dof
					K[d] .+= o.M.CT_XAR[d] * tmp
				end
			end
			@inbounds for d ∈ 1:o.dof
				K[d][o.crosstab⋂✻ind] .-= view(u✻XAR,:,d)  # subtract crosstab of u✻XAR wrt bootstrapping cluster and all-cluster-var intersections from M
				o.scorebs && (K[d] .-= o.ClustShare * colsum(K[d]))  # recenter
			end

			@inbounds for c ∈ 1+o.granular:o.NErrClustCombs
				for d ∈ 1:o.dof
					nrows(o.clust[c].order)>0 &&
						(K[d] = K[d][o.clust[c].order,:])
					o.Kcd[c,d] = @panelsum(o, K[d], o.clust[c].info)
				end
			end
		else  # B = 0. In this case, only 1st term of (64) is non-zero after multiplying by v* (= all 1's), and it is then a one-way sum by c
			o.scorebs &&
				(u✻XAR .-= o.ClustShare * colsum(u✻XAR))  # recenter if OLS
			@inbounds for c ∈ 1:o.NErrClustCombs
				nrows(o.clust[c].order)>0 &&
					(u✻XAR = u✻XAR[o.clust[c].order,:])
				tmp = @panelsum(o, u✻XAR, o.clust[c].info)
				for d ∈ 1:o.dof
					o.Kcd[c,d] = reshape(view(tmp,:,d),:,1)
				end
			end
		end
	end
	MakeNumerAndJ!(o, 1, thisr)  # compute J = κ * v; if Nw > 1, then this is for 1st group; if interpolating, it is only group, and may be needed now to prep interpolation
	nothing
end

# compute stuff depending linearly on v, needed to prep for interpolation
function MakeNumerAndJ!(o::StrBootTest{T}, w::Integer, r::AbstractVector=Vector{T}(undef,0)) where T  # called to *prepare* interpolation, or when w>1, in which case there is no interpolation
	o.numerw = o.scorebs ?
			   (o.B>0 ?
				 	o.SuwtXA'o.v :
				 	o.SuwtXA * o.v_sd    ) :
			   (!o.robust || o.granular || o.purerobust ?
				  	 o.R * (o.β̈dev = o.SuwtXA * o.v) :
				 		(o.R * o.SuwtXA) * o.v)

	if isone(w)
		if o.ARubin
			o.numerw[:,1] = o.v_sd * o.DGP.β̈[o.kX₁+1:end]  # coefficients on excluded instruments in ARubin OLS
		elseif !o.null
			o.numerw[:,1] = o.v_sd * (o.R * (o.ML ? o.β̈ : o.M.β̈) - r)  # Analytical Wald numerator; if imposing null then numer[:,1] already equals this. If not, then it's 0 before this.
		end
	end

	@storeWtGrpResults!(o.numer, o.numerw)

	if o.B>0 && o.robust && o.bootstrapt
		if o.granular || o.purerobust  # optimized treatment when bootstrapping by many/small groups
			if o.purerobust
				o.u✻ = o.ü .* o.v
				o.NFE>0 && partialFE!(o, o.u✻)
				minusX₁₂B(o.u✻, o.X₁, o.X₂, o.β̈dev)
			else  # clusters small but not all singletons
				if o.NFE>0 && !o.FEboot
					o.u✻ = o.ü .* view(o.v, o.ID✻, :)
					partialFE!(o, o.u✻)
					@inbounds for d ∈ 1:o.dof
						o.Jcd[1,d] = @panelsum(o, o.u✻, view(o.M.WXAR,:,d), o.info⋂)                                - panelsum2(o, o.X₁, o.X₂, view(o.M.WXAR,:,d), o.info⋂) * o.β̈dev
					end
				else
					_v = view(o.v,o.ID✻_✻⋂,:)
					@inbounds for d ∈ 1:o.dof
						o.Jcd[1,d] = panelsum(o, panelsum(o, o.ü, view(o.M.WXAR,:,d), o.info✻⋂) .* _v, o.info⋂_✻⋂) - panelsum2(o, o.X₁, o.X₂, view(o.M.WXAR,:,d), o.info⋂) * o.β̈dev
					end
				end
			end
		end
		@inbounds	for c ∈ o.granular+1:o.NErrClustCombs, d ∈ eachindex(axes(o.Jcd, 2), axes(o.Kcd, 2))
			o.Jcd[c,d] = o.Kcd[c,d] * o.v
		end
	end
	nothing
end

function MakeNonWRELoop1!(o::StrBootTest, tmp::Matrix, w::Integer)
	@inbounds Threads.@threads for k ∈ 1:ncols(o.v)
		for i ∈ 1:o.dof
			for j ∈ 1:i
				tmp[j,i] = o.denom[i,j][k]  # fill upper triangle, which is all invsym() looks at
			end
		end
		numer_l = view(o.numerw,:,k)
		o.dist[k+first(o.WeightGrp[w])-1] = numer_l'invsym(tmp)*numer_l  # in degenerate cases, cross() would turn cross(.,.) into 0
	end
	nothing
end

function MakeNonWREStats!(o::StrBootTest{T}, w::Integer) where T
	w > 1 && MakeNumerAndJ!(o, w)
	!o.bootstrapt && return

	if o.robust
    if !o.interpolating  # these quadratic computation needed to *prepare* for interpolation but are superseded by interpolation once it is going
    	o.purerobust && (u✻2 = o.u✻ .^ 2)
    	@inbounds for i ∈ 1:o.dof, j ∈ 1:i
    		o.purerobust &&
  	   		(o.denom[i,j] = cross(view(o.M.WXAR,:,i), view(o.M.WXAR,:,j), u✻2) * (o.clust[1].even ? o.clust[1].multiplier : -o.clust[1].multiplier))
				for c ∈ o.purerobust+1:o.NErrClustCombs
					@clustAccum!(o.denom[i,j], c, j==i ? coldot(o, o.Jcd[c,i]) : coldot(o, o.Jcd[c,i],o.Jcd[c,j]))
				end
  		end
   	end

		if isone(o.dof)
			@storeWtGrpResults!(o.dist, o.numerw ./ sqrtNaN.(o.denom[1,1]))
			isone(w) &&
				(o.statDenom = hcat(o.denom[1,1][1]))  # original-sample denominator
		elseif o.dof==2  # hand-code 2D numer'inv(denom)*numer
			t1 = view(o.numerw,1,:)'; t2 = view(o.numerw,2,:)'; t12 = t1.*t2
			@storeWtGrpResults!(o.dist, (t1.^2 .* o.denom[2,2] .- 2 .* t12 .* o.denom[2,1] .+ t2.^2 .* o.denom[1,1]) ./ (o.denom[1,1].*o.denom[2,2] .- o.denom[2,1].^2))
			isone(w) &&
				(o.statDenom = [o.denom[1,1][1] o.denom[2,1][1] ; o.denom[2,1][1] o.denom[2,2][1]])  # original-sample denominator
		else  # build each replication's denominator from vectors that hold values for each position in denominator, all replications
			tmp = Matrix{T}(undef, o.dof, o.dof)
			MakeNonWRELoop1!(o,tmp,w)
			isone(w) && (o.statDenom = tmp)  # original-sample denominator
		end
	else  # non-robust
		AR = o.ML ? o.AR : o.M.AR
		if isone(o.dof)  # optimize for one null constraint
			o.denom[1,1] = o.R * AR
			if !o.ML
				o.u✻ = o.B>0 ? o.v .* o.ü : reshape(o.ü,:,1)  # reshape for type stability
				if o.scorebs
					if o.haswt  # Center variance if interpolated
						o.u✻ .-= o.ClustShare'o.u✻
					else
						o.u✻ .-= colsum(o.u✻) * o.ClustShare  # Center variance if interpolated
					end
				else
					minusX₁₂B(o.u✻, o.X₁, o.X₂, o.β̈dev)  # residuals of wild bootstrap regression are the wildized residuals after partialling out X (or XS) (Kline && Santos eq (11))
				end
				if o.haswt
					o.denom[1,1] .*= o.wt'(o.u✻ .^ 2)
				else
					o.denom[1,1] .*= coldot(o,o.u✻)
				end
			end

			@storeWtGrpResults!(o.dist, o.numerw ./ sqrtNaN.(o.denom[1,1]))
			isone(w) && (o.statDenom = o.denom[1,1])  # original-sample denominator
		else
			o.denom[1,1] = o.R * AR
			if o.ML
				for k ∈ 1:ncols(o.v)
					numer_l = view(o.numerw,:,k)
					o.dist[k+first(o.WeightGrp[w])-1] = numer_l'invsym(o.denom[1,1])*numer_l
				end
				isone(w) && (o.statDenom = o.denom[1,1])  # original-sample denominator
			else
				invdenom = invsym(o.denom[1,1])
				if o.B>0
					for k ∈ 1:ncols(o.v)
						numer_l = view(o.numerw,:,k)
						o.dist[k+first(o.WeightGrp[w])-1] = o.numer_l'invdenom*numer_l
						o.u✻ = view(o.v,:,k) .* o.ü
						if o.scorebs
							if o.haswt  # Center variance if interpolated
								o.u✻ .-= o.wt'o.u✻ * o.ClustShare
							else
								o.u✻ .-= colsum(o.u✻) * o.ClustShare
							end
						else
							minusX₁₂B(o.u✻, o.X₁, o.X₂, view(o.β̈dev,:,k))  # residuals of wild bootstrap regression are the wildized residuals after partialling out X (or XS) (Kline && Santos eq (11))
						end
						o.dist[k+first(o.WeightGrp[w])-1] ./= (tmp = symcross(o.u✻, o.wt))
					end
				else
					o.dist[1] = o.numerw'invdenom*o.numerw
					o.u✻ = o.ü
					if o.scorebs
						if o.haswt  # Center variance if interpolated
							o.u✻ .-= o.wt'o.u✻ * o.ClustShare
						else
							o.u✻ .-= colsum(o.u✻) * o.ClustShare
						end
					else
						minusX₁₂B(o.u✻, o.X₁, o.X₂, o.β̈dev)  # residuals of wild bootstrap regression are the wildized residuals after partialling out X (or XS) (Kline && Santos eq (11))
					end
					o.dist[1] /= (tmp = symcross(o.u✻, o.wt))
				end
				isone(w) && (o.statDenom = o.denom[1,1] * tmp)  # original-sample denominator
			end
		end
	end
	nothing
end


text:
function medianvalue(min, max, step)
    rg = min:step:max
    rg[div(length(rg), 2)+1]
end

function slider(; min = 0, max = 100, step = 1, value = nothing, type = "", kwargs...)
    (min isa AbstractObservable) || (min = Observable(min))
    (max isa AbstractObservable) || (max = Observable(max))
    (step isa AbstractObservable) || (step = Observable(step))
    value === nothing && (value = medianvalue(min[], max[], step[]))
    (value isa AbstractObservable) || (value = Observable(value))

    oninput = js"function() {data.value = parseFloat(this.value);}"
    m = input(; min = min, max = max, step = step, value = value, type = "range",
              oninput = oninput, selector = "input.slider", kwargs...)
    MithrilWidget{:slider}(component(m), value)
end

ext: jl


text:
using Distributions
using Cytof5
using JLD2, FileIO
using PyPlot

# Directory containing CB results
data_dir = "/scratchdata/alui2/cytof/results/cb/"

# Path to mm1 and best output
path_to_mm0_output = "$(data_dir)/best/output.jld2"
path_to_mm1_output = "$(data_dir)/mm1/output.jld2"
path_to_mm2_output = "$(data_dir)/mm2/output.jld2";

# Get output for miss-mech-0 and miss-mech-1
mm0 = load(path_to_mm0_output)
mm1 = load(path_to_mm1_output)
mm2 = load(path_to_mm2_output)

# Get K from output
K = mm0["c"].K

# Zs
Zs_mm0 = [m[:Z] for m in mm0["out"][1]]
Zs_mm1 = [m[:Z] for m in mm1["out"][1]]
Zs_mm2 = [m[:Z] for m in mm2["out"][1]]

# Last Z
Zhat_mm0 = Zs_mm0[end]
Zhat_mm1 = Zs_mm1[end]
Zhat_mm2 = Zs_mm2[end]

# Ws
Ws_mm0 = cat([m[:W] for m in mm0["out"][1]]..., dims=3)
Ws_mm1 = cat([m[:W] for m in mm1["out"][1]]..., dims=3)
Ws_mm2 = cat([m[:W] for m in mm2["out"][1]]..., dims=3)

What_mm0 = mean(Ws_mm0, dims=3)[:, :, 1]
What_mm1 = mean(Ws_mm1, dims=3)[:, :, 1]
What_mm2 = mean(Ws_mm2, dims=3)[:, :, 1]

function countmatches(Z1, Z2; discrepancy=0)
  K1 = size(Z1, 2)
  K2 = size(Z2, 2)
  k1s = Int8[]
  k2s = Int8[]

  for k1 in 1:K1
    for k2 in 1:K2
      if sum(abs.(Z1[:, k1] - Z2[:, k2])) <= discrepancy
        append!(k1s, k1)
        append!(k2s, k2)
      end
    end
  end

  return unique(k1s), unique(k2s)
end

function visualize(Z1, Z2, W1, i; discrepancy=0, plot=true)
  @assert size(Z1, 2) == size(Z2, 2)
  k1s, k2s = countmatches(Z1, Z2, discrepancy=discrepancy)
  K = size(Z1, 2)
  matched_features = k1s
  unmatched_features = setdiff(1:K, k1s)

  if ndims(W1) == 3
    if plot
      PyPlot.boxplot(W1[i, [matched_features; unmatched_features], :]',
                     showfliers=false)
    end
    W1_hat = mean(W1, dims=3)[:, :, 1]
  elseif ndims(W1) == 2
    if plot
      PyPlot.scatter(1:K, W1[i, [matched_features; unmatched_features]])
    end
    W1_hat = W1 .+ 0
  else
    throw("ndims(W1) needs to be 2 or 3!")
  end

  num_matches_1 = length(unique(k1s))
  num_matches_2 = length(unique(k2s))

  if plot
    PyPlot.axvline(num_matches_1 + .5, ls="--", color="grey")
    PyPlot.axhline(0, ls="--", color="grey")
    PyPlot.xticks(1:K, [matched_features; unmatched_features])
  end

  println("number of matched features: $((num_matches_1, num_matches_2))")
  println("mean W for matched features: $(mean(W1_hat[i, matched_features]))")
  println("mean W for unmatched features: $(mean(W1_hat[i, unmatched_features]))")
  println("sum W for matched features: $(sum(W1_hat[i, matched_features]))")
  println("sum W for unmatched features: $(sum(W1_hat[i, unmatched_features]))")
end

# MM0 vs MM1
for d in [0, 1, 2, 3]
  for i in 1:3
    println("i: $(i), d: $(d)")
    visualize(Zhat_mm0, Zhat_mm1, Ws_mm0, i, discrepancy=d, plot=false)
  end
  println()
end

# MM0 vs MM2
for d in [0, 1, 2, 3]
  for i in 1:3
    println("i: $(i), d: $(d)")
    visualize(Zhat_mm0, Zhat_mm2, Ws_mm0, i, discrepancy=d, plot=false)
  end
  println()
end


text:
using Plots, Revise, IJulia, Distributions, GLM, DataFrames, DifferentialEquations

plot(rand(Normal(),5),rand(Normal(),5))
ols = lm(@formula(X ~ Y), DataFrame(X=[1,2],Y=[4,5]))
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

max_stars_repo_name: olynch/dotfiles

max_issues_repo_name: olynch/dotfiles

max_forks_repo_name: olynch/dotfiles


text:
using Test
using ParametricSurfaces
using ParametricSurfaces: WavePropBase.Geometry

@testset "Disk" begin
    disk = ParametricSurfaces.Circle() # abstract entity
    Γ = boundary(disk) |> Domain
    M = meshgen(Γ,(10,))
    # plot(M,Γ)
    @test entities(Γ) == boundary(disk)
    @test geometric_dimension(disk) == 2
end

@testset "Ball" begin
    ball = ParametricSurfaces.Sphere() # abstract entity
    Γ = boundary(ball) |> Domain
    M = meshgen(Γ,(1,1))
    # plot(M,Γ)
    @test entities(Γ) == boundary(ball)
    @test geometric_dimension(ball) == 3
end

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
__precompile__()
module IndirectLikelihood

import Base: size, show

using ArgCheck: @argcheck
using ContinuousTransformations: transform, inverse
using DocStringExtensions: SIGNATURES
using ForwardDiff: jacobian
using MacroTools
using Parameters: @unpack
using StatsBase: mean_and_cov, AbstractWeights

"""
The default random number generator for methods in this package, used when not
specified in the first argument.
"""
const RNG = Base.Random.GLOBAL_RNG

include("utilities.jl")
include("general.jl")
include("mvnormal.jl")
include("ols.jl")

end # module


text:
module TestLevels
using Compat
using Compat.Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, levels!

@testset "CategoricalPool{Int} updates levels/index/order correctly" begin
    pool = CategoricalPool([2, 1, 3])

    @test isa(levels(pool), Vector{Int})
    @test length(levels(pool)) === 3
    @test levels(pool) == pool.index == [2, 1, 3]
    @test pool.invindex == Dict(1=>2, 2=>1, 3=>3)
    @test pool.order == [1, 2, 3]
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:3]

    for rep in 1:3
        push!(pool, 4)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 4
        @test pool.index == [2, 1, 3, 4]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4)
        @test pool.order == [1, 2, 3, 4]
        @test pool.levels == [2, 1, 3, 4]
        @test get(pool, 4) === DefaultRefType(4)
        @test pool[4] === CategoricalArrays.catvalue(4, pool)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:4]
    end

    for rep in 1:3
        push!(pool, 0)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 5
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5)
        @test pool.order == [1, 2, 3, 4, 5]
        @test pool.levels == [2, 1, 3, 4, 0]
        @test get(pool, 0) === DefaultRefType(5)
        @test pool[5] === CategoricalArrays.catvalue(5, pool)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:5]
    end

    for rep in 1:3
        push!(pool, 10, 11)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 7
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11]
        @test get(pool, 10) === DefaultRefType(6)
        @test get(pool, 11) === DefaultRefType(7)
        @test pool[6] === CategoricalArrays.catvalue(6, pool)
        @test pool[7] === CategoricalArrays.catvalue(7, pool)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:7]
    end

    for rep in 1:3
        push!(pool, 12, 13)

        @test isa(pool.index, Vector{Int})
        @test length(pool) === 9
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11, 12, 13]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7, 12=>8, 13=>9)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7, 8, 9]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11, 12, 13]
        @test get(pool, 12) === DefaultRefType(8)
        @test get(pool, 13) === DefaultRefType(9)
        @test pool[8] === CategoricalArrays.catvalue(8, pool)
        @test pool[9] === CategoricalArrays.catvalue(9, pool)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:9]
    end

    for rep in 1:3
        delete!(pool, 13)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 8
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10, 11, 12]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6, 11=>7, 12=>8)
        @test pool.order == [1, 2, 3, 4, 5, 6, 7, 8]
        @test pool.levels == [2, 1, 3, 4, 0, 10, 11, 12]
        @test_throws KeyError get(pool, 13)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:8]
    end

    for rep in 1:3
        delete!(pool, 12, 11)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 6
        @test levels(pool) == pool.index == [2, 1, 3, 4, 0, 10]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 4=>4, 0=>5, 10=>6)
        @test pool.order == [1, 2, 3, 4, 5, 6]
        @test pool.levels == [2, 1, 3, 4, 0, 10]
        @test_throws KeyError get(pool, 11)
        @test_throws KeyError get(pool, 12)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:6]
    end

    for rep in 1:3
        delete!(pool, 4)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 5
        @test levels(pool) == pool.index == [2, 1, 3, 0, 10]
        @test pool.invindex == Dict(1=>2, 2=>1, 3=>3, 0=>4, 10=>5)
        @test pool.order == [1, 2, 3, 4, 5]
        @test pool.levels == [2, 1, 3, 0, 10]
        @test_throws KeyError get(pool, 4)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:5]
    end

    @test levels!(pool, [1, 2, 3]) === pool
    @test levels(pool) == [1, 2, 3]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [1, 2, 3]
    @test pool.invindex == Dict(1=>1, 2=>2, 3=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [1, 2, 3]
    @test get(pool, 1) === DefaultRefType(1)
    @test_throws KeyError get(pool, 0)
    @test_throws KeyError get(pool, 10)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:3]

    @test levels!(pool, [1, 2, 4]) === pool
    @test levels(pool) == [1, 2, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [1, 2, 4]
    @test pool.invindex == Dict(1=>1, 2=>2, 4=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [1, 2, 4]
    @test get(pool, 1) === DefaultRefType(1)
    @test_throws KeyError get(pool, 3)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:3]

    @test levels!(pool, [6, 5, 4]) === pool
    @test levels(pool) == [6, 5, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == pool.index == [6, 5, 4]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3)
    @test pool.order == [1, 2, 3]
    @test pool.levels == [6, 5, 4]
    @test get(pool, 5) === DefaultRefType(2)
    @test_throws KeyError get(pool, 3)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:3]

    # Changing order while preserving existing levels
    @test levels!(pool, [5, 6, 4]) === pool
    @test levels(pool) == [5, 6, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 3
    @test length(pool.valindex) == 3
    @test levels(pool) == [5, 6, 4]
    @test pool.index == [6, 5, 4]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3)
    @test pool.order == [2, 1, 3]
    @test pool.levels == [5, 6, 4]
    @test get(pool, 5) === DefaultRefType(2)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:3]

    # Adding levels while preserving existing ones
    @test levels!(pool, [5, 2, 3, 6, 4]) === pool
    @test levels(pool) == [5, 2, 3, 6, 4]

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 5
    @test length(pool.valindex) == 5
    @test levels(pool) == [5, 2, 3, 6, 4]
    @test pool.index == [6, 5, 4, 2, 3]
    @test pool.invindex == Dict(6=>1, 5=>2, 4=>3, 2=>4, 3=>5)
    @test pool.order == [4, 1, 5, 2, 3]
    @test pool.levels == [5, 2, 3, 6, 4]
    @test get(pool, 2) === DefaultRefType(4)
    @test get(pool, 3) === DefaultRefType(5)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:5]

    for rep in 1:3
        delete!(pool, 6)

        @test isa(pool.index, Vector{Int})
        @test length(pool) == 4
        @test length(pool.valindex) == 4
        @test levels(pool) == [5, 2, 3, 4]
        @test pool.index == [5, 4, 2, 3]
        @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4)
        @test pool.order == [1, 4, 2, 3]
        @test pool.levels == [5, 2, 3, 4]
        @test get(pool, 4) === DefaultRefType(2)
        @test_throws KeyError get(pool, 6)
        @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:4]
    end

    # get!
    ordered!(pool, true)
    @test_throws OrderedLevelsException get!(pool, 10)
    ordered!(pool, false)

    @test get!(pool, 10) === DefaultRefType(5)

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 5
    @test length(pool.valindex) == 5
    @test levels(pool) == [5, 2, 3, 4, 10]
    @test pool.index == [5, 4, 2, 3, 10]
    @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4, 10=>5)
    @test pool.order == [1, 4, 2, 3, 5]
    @test pool.levels == [5, 2, 3, 4, 10]
    @test get(pool, 10) === DefaultRefType(5)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:5]

    # get! with CategoricalValue adding new levels
    v = CategoricalArrays.catvalue(2, CategoricalPool([100, 99, 5, 2]))

    ordered!(pool, true)
    @test_throws OrderedLevelsException get!(pool, v)
    ordered!(pool, false)

    @test get!(pool, v) === DefaultRefType(7)

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 7
    @test length(pool.valindex) == 7
    @test levels(pool) == [100, 99, 5, 2, 3, 4, 10]
    @test pool.index == [5, 4, 2, 3, 10, 100, 99]
    @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4, 10=>5, 100=>6, 99=>7)
    @test pool.order == [3, 6, 4, 5, 7, 1, 2]
    @test pool.levels == [100, 99, 5, 2, 3, 4, 10]
    @test get(pool, 99) === DefaultRefType(7)
    @test get(pool, 100) === DefaultRefType(6)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:7]

    # get! with CategoricalValue not adding new levels
    v = CategoricalArrays.catvalue(1, CategoricalPool([100, 2]))
    @test get!(pool, v) === DefaultRefType(6)

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 7
    @test length(pool.valindex) == 7
    @test levels(pool) == [100, 99, 5, 2, 3, 4, 10]
    @test pool.index == [5, 4, 2, 3, 10, 100, 99]
    @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4, 10=>5, 100=>6, 99=>7)
    @test pool.order == [3, 6, 4, 5, 7, 1, 2]
    @test pool.levels == [100, 99, 5, 2, 3, 4, 10]
    @test get(pool, 99) === DefaultRefType(7)
    @test get(pool, 100) === DefaultRefType(6)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:7]

    # get! with CategoricalValue from same pool
    @test get!(pool, pool[1]) === DefaultRefType(1)

    @test isa(pool.index, Vector{Int})
    @test length(pool) == 7
    @test length(pool.valindex) == 7
    @test levels(pool) == [100, 99, 5, 2, 3, 4, 10]
    @test pool.index == [5, 4, 2, 3, 10, 100, 99]
    @test pool.invindex == Dict(5=>1, 4=>2, 2=>3, 3=>4, 10=>5, 100=>6, 99=>7)
    @test pool.order == [3, 6, 4, 5, 7, 1, 2]
    @test pool.levels == [100, 99, 5, 2, 3, 4, 10]
    @test get(pool, 99) === DefaultRefType(7)
    @test get(pool, 100) === DefaultRefType(6)
    @test pool.valindex == [CategoricalArrays.catvalue(i, pool) for i in 1:7]

    # get! with CategoricalValue not adding new levels
    v = CategoricalArrays.catvalue(1, CategoricalPool(["a", "b"]))
    @test_throws MethodError get!(pool, v)

    # get! with CategoricalValue{Any} (#220)
    p1 = CategoricalPool(Any['a', 'b', 'c'])
    p2 = CategoricalPool(Any['a', 'b', 'x'])
    @test get!(p1, p2[1]) === UInt32(1)
    @test get!(p1, p2[3]) === UInt32(4)

    # get! with ordered CategoricalValue marks unordered empty pool as ordered
    p1 = CategoricalPool(['b', 'c', 'a'])
    ordered!(p1, true)
    p2 = CategoricalPool(Char[])
    @test get!(p2, p1[1]) === UInt32(1)
    @test isordered(p2)
    # But push! does not
    p2 = CategoricalPool(Char[])
    @test push!(p2, p1[1]) === p2
    @test !isordered(p2)
end

@testset "overflow of reftype is detected and doesn't corrupt levels" begin
    res = @test_throws LevelsException{Int, UInt8} CategoricalPool{Int, UInt8}(collect(256:-1:1))
    @test res.value.levels == [1]
    @test sprint(showerror, res.value) == "cannot store level(s) 1 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."

    pool = CategoricalPool(collect(30:288))
    res = @test_throws LevelsException{Int, UInt8} convert(CategoricalPool{Int, UInt8}, pool)
    @test res.value.levels == collect(285:288)
    @test sprint(showerror, res.value) == "cannot store level(s) 285, 286, 287 and 288 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."

    pool = CategoricalPool{String, UInt8}(string.(318:-1:65))
    res = @test_throws LevelsException{String, UInt8} levels!(pool, vcat("az", levels(pool), "bz", "cz"))
    @test res.value.levels == ["bz", "cz"]
    @test sprint(showerror, res.value) == "cannot store level(s) \"bz\" and \"cz\" since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."
    lev = copy(levels(pool))
    levels!(pool, vcat(lev, "az"))
    @test levels(pool) == vcat(lev, "az")
end

end


text:
function checkPW(word)
	if word[end] in ['i','l','o']
		# println(bytestring(word), " contains ilo")
		return false
	end
	sequence = false
	doubletCount = 0
	current = ""
	for w in 1:length(word)-1
		if word[w] in ['i','l','o']
			# println(bytestring(word), " contains ilo")
			return false
		end
		# println(w, " ", word[w], " ", word[w+1], " ", bytestring(word))
		if word[w] == word[w+1]
			if current != word[w]
				doubletCount += 1
				current = word[w]
			end
		else
			current = ""
		end
		if w < length(word)-1 && word[w]+2 == word[w+1]+1 == word[w+2]
			sequence = true
		end
	end
	if doubletCount < 2
		# println(bytestring(word), " contains not enough duplicates: ", doubletCount)
		return false
	end
	if !sequence
		# println(bytestring(word), " contains no sequence")
		return false
	end
	return true
end

function makeNew(pw)
	last = length(pw)
	pw[last] += 1
	while last > 1 && pw[last] > 'z'
		pw[last] = 'a'
		last -= 1
		pw[last] += 1
	end
	if bytestring(pw) == join(['z' for i in 1:length(pw)])
		return false
	end
	return true
end

passwd = b"vzbxxzzz"
while !checkPW(passwd)
	# println(bytestring(passwd))
	if !makeNew(passwd)
		break
	end
end
println(bytestring(passwd))

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
# BUG: The current result does not align with R's psych package
using Distributions

function kappa(x, α = 0.05)
    tot = sum(x)
    x = x ./ tot        # convert to probabilities
    rs = sum(x, 2)      # row sum
    cs = sum(x, 1)      # column sum

    prob = rs * cs
    po = trace(x)
    pc = trace(prob)
    kappa = (po - pc)/(1 - pc)

    var_kappa = (1 / (tot*(1-pc)^4)) * (trace(x * (I * (1-pc) - (rs .+ (cs))*(1-po))^2 ) + (1-po)^2 * (sum(x * (cs .+ (rs))^2) - trace(x * (cs .+ (rs))^2))  -(po*pc - 2*pc +po)^2)

    lower_ci = kappa + quantile(Normal(0, 1), (α/2)) * sqrt(var_kappa)
    upper_ci = kappa - quantile(Normal(0, 1), (α/2)) * sqrt(var_kappa)

    return (kappa, var_kappa, lower_ci, upper_ci)
end

# kappa(cm1.matrix)

# a = cm1.matrix


text:
module TianJun

end


text:
#
# Heap Sort
#

function heapSort(A)

  buildHeap(A)

  n = length(A)
  i = n
  while i > 0

    tmp = A[1]
    A[1] = A[i]
    A[i] = tmp

    heapify(A, 1, i)

    i = i - 1
  end

end

function buildHeap(A)

  n = length(A)
  i = Int(n/2)
  while i > 0

    heapify(A, i, n)

    i = i - 1
  end

end

function heapify(A, idx, idx_max)
  largest = idx
  left = 2 * (idx-1) + 1
  right = 2 * (idx-1) + 2

  if left < idx_max && A[left] > A[idx]
    largest = left
  end
  if right < idx_max && A[right] > A[largest]
    largest = right
  end
  if largest != idx
    tmp = A[idx]
    A[idx] = A[largest]
    A[largest] = tmp
    heapify(A,  largest,  idx_max)
  end
end

function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  A = shuffle(rng, Vector(1:100))
  Aorigin = copy(A)

  # sorting
  @time heapSort(A)

  print("The original list: ")
  println(Aorigin)

  print("The sorted list: ")
  println(A)

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
using ChainRulesCore
using ForwardDiff
using FiniteDifferences
using Transits: compute_grad, compute_gn, compute_gn_jac

# @testset "PolynomialLimbDark" begin
#     u = [0.4, 0.26, 0.1, 0.1]
#     # test_frule(PolynomialLimbDark, u)
#     # test_rrule(PolynomialLimbDark, u)

#     ld = PolynomialLimbDark(u)

#     test_frule(compute, ld, 0.0, 0.1; fdm=forward_fdm(10, 1))
#     test_rrule(compute, ld, 0.0, 0.1; fdm=forward_fdm(10, 1))
# end

# @testset "QuadLimbDark" begin
#     u = [0.4, 0.26]
#     test_frule(QuadLimbDark, u)
#     test_rrule(QuadLimbDark, u)

#     ld = QuadLimbDark(u)

#     test_frule(compute, ld, 0.0, 0.1)
#     test_rrule(compute, ld, 0.0, 0.1)
# end

function finite_diff(b_, r_, u_n_, law=PolynomialLimbDark; diff=big(1e-18))
    b = big(b_)
    r = big(r_)
    u_n = big.(u_n_)
    ld = law(u_n)
    f = ld(b, r)

    # modulate r
    f₊ = ld(b, r + diff)
    if r == 1 && length(u_n) > 2
        dfdr = (f₊ - f) / diff
    else
        f₋ = ld(b, r - diff)
        dfdr = (f₊ - f₋) / (2 * diff)
    end

    # modulate b
    f₊ = ld(b + diff, r)
    if b > diff
        f₋ = ld(b - diff, r)
        dfdb = (f₊ - f₋) / (2 * diff)
    elseif iszero(b)
        dfdb = zero(b)
    else
        dfdb = (f₊ - f) / diff
    end

    return Float64.([dfdb, dfdr])
end

function test_compute_grad(b, r, u_n, law=PolynomialLimbDark)
    ld = law(u_n)
    _, _, dfdb, dfdr = @inferred compute_grad(ld, b, r)
    grad_analytical = [dfdb, dfdr]

    atol = r != 1 ? 1e-7 : 3e-6

    @testset "compute_grad - b=$b, r=$r" begin
        # test against BigFloat for numerical accuracy
        _, _, dfdb_big, dfdr_big = compute_grad(ld, big(b), big(r))
        grad_big = [dfdb_big, dfdr_big]
        @test grad_analytical ≈ grad_big atol = atol

        # test against finite diff for correctness
        grad_numerical = finite_diff(b, r, u_n, law)
        @test grad_analytical ≈ grad_numerical atol = atol
    end
end

logspace(start, stop, N) = 10 .^range(log10(start), log10(stop), length=N)

test_coeffs = [[0.0], [1.0], [2.0, -1.0], [3.0, -3.0, 1.0], fill(0.1, 4), fill(0.1, 5), fill(0.1, 10)]
test_names = ["uniform", "linear", "quadratic", "cubic", "quartic", "quintic", "10th order"]

@testset "`compute` grads - $name" for (name, u_n) in zip(test_names, test_coeffs)
    # b + r < 1
    rs = [0.01, 0.1, 0.5, 0.999999]
    nb = 50
    ϵ, δ = 1e-9, 1e-3
    for r in rs
        bs = unique(abs.([
            range(0, ϵ, length=nb);
            range(ϵ, δ, length=nb);
            range(δ, r - δ, length=nb);
            -logspace(δ, ϵ, nb) .+ r;
            range(r - ϵ, r + ϵ, length=nb);
            logspace(ϵ, δ, nb) .+ r;
            range(r + δ, 1 - r - δ, length=nb);
            -logspace(δ, ϵ, nb) .+ (1 - r);
            range(1 - r - ϵ, 1 - r + ϵ, length=nb);
            logspace(ϵ, δ, nb) .+ (1 - r);
            range(1 - r + δ, 1 + r - δ, length=nb);
            -logspace(δ, ϵ, nb) .+ (1 + r);
            range(1 + r - ϵ, 1 + r - 1e-13, length=nb)
        ]))
        test_compute_grad.(bs, r, (u_n,), PolynomialLimbDark)
        if length(u_n) < 3
            test_compute_grad.(bs, r, (u_n,), QuadLimbDark)
        end
    end

    rs = [1, 1.000001, 2, 10, 100]
    for r in rs
        bs = unique(abs.([
            r - 1 + 1e-13;
            logspace(ϵ, δ, nb) .+ (r - 1);
            range(r - 1 + δ, r - δ, length=nb);
            -logspace(δ, ϵ, nb) .+ r;
            range(r - ϵ, r + ϵ, length=nb);
            logspace(ϵ, δ, nb) .+ r;
            range(r + δ, r + 1 - δ, length=nb);
            -logspace(δ, ϵ, nb) .+ (r + 1);
            r + 1 - 1e-13
        ]))
        test_compute_grad.(bs, r, (u_n,), PolynomialLimbDark)
        if length(u_n) < 3
            test_compute_grad.(bs, r, (u_n,), QuadLimbDark)
        end
    end
end

@testset "`compute_gn` jacs - $name" for (name, u_n) in zip(test_names, test_coeffs)
    un_full = [-1, u_n...]
    gn_primal = @inferred compute_gn(un_full)
    gn, jac_analytical = @inferred compute_gn_jac(un_full)

    @test gn_primal == gn

    jac_numerical = FiniteDifferences.jacobian(central_fdm(5, 1), compute_gn, un_full)[1]'

    @test jac_analytical ≈ jac_numerical atol=1e-8
end

@testset "frule, rrule comparison - $name" for (name, u_n) in zip(test_names, test_coeffs)
    @testset "PolynomialLimbDark" begin
        function primal(X)
            ld = PolynomialLimbDark(X[3:end])
            return compute(ld, X[1], X[2])
        end

        function gradr(X)
            ld, ld_pullback = rrule(PolynomialLimbDark, X[3:end])
            f, f_pullback = rrule(compute, ld, X[1], X[2])

            f̄ = one(eltype(X))
            _, l̄d, b̄, r̄ = f_pullback(f̄)
            _, ū_n = ld_pullback(l̄d)
            return [b̄, r̄, ū_n...]
        end

        function gradf(X)
            N = length(X) - 2
            ∂us = zeros(N)
            for i in eachindex(∂us)
                Δus = zeros(N)
                Δus[i] = 1
                ld, ∂ld = frule((0, Δus), PolynomialLimbDark, X[3:end])
                _, ∂us[i] = frule((0, ∂ld, 0, 0), compute, ld, X[1], X[2])
            end
            ld, ∂ld3 = frule((0, zeros(N)), PolynomialLimbDark, X[3:end])
            _, ∂b = frule((0, ∂ld3, 1, 0), compute, ld, X[1], X[2])
            _, ∂r = frule((0, ∂ld3, 0, 1), compute, ld, X[1], X[2])
            return [∂b, ∂r, ∂us...]
        end

        function test_grads(X)
            grad_rev = @inferred gradr(X)
            grad_for = @inferred gradf(X)
            grad_FD = ForwardDiff.gradient(primal, X)
            @test grad_rev ≈ grad_for atol=1e-7
            @test grad_rev ≈ grad_FD atol=1e-7
        end

        test_grads([0.3, 0.2, u_n...])
        test_grads([1.1, 0.3, u_n...])
    end

    if length(u_n) < 3
        @testset "QuadLimbDark" begin
            function primal(X)
                ld = QuadLimbDark(X[3:end])
                return compute(ld, X[1], X[2])
            end

            function gradr(X)
                ld, ld_pullback = rrule(QuadLimbDark, X[3:end])
                f, f_pullback = rrule(compute, ld, X[1], X[2])

                f̄ = one(eltype(X))
                _, l̄d, b̄, r̄ = f_pullback(f̄)
                _, ū_n = ld_pullback(l̄d)
                return [b̄, r̄, ū_n...]
            end

            function gradf(X)
                N = length(X) - 2
                ∂us = zeros(N)
                for i in eachindex(∂us)
                    Δus = zeros(N)
                    Δus[i] = 1
                    ld, ∂ld = frule((0, Δus), QuadLimbDark, X[3:end])
                    _, ∂us[i] = frule((0, ∂ld, 0, 0), compute, ld, X[1], X[2])
                end
                ld, ∂ld3 = frule((0, zeros(N)), QuadLimbDark, X[3:end])
                _, ∂b = frule((0, ∂ld3, 1, 0), compute, ld, X[1], X[2])
                _, ∂r = frule((0, ∂ld3, 0, 1), compute, ld, X[1], X[2])
                return [∂b, ∂r, ∂us...]
            end

            function test_grads(X)
                grad_rev = @inferred gradr(X)
                grad_for = @inferred gradf(X)
                grad_FD = ForwardDiff.gradient(primal, X)
                @test grad_rev ≈ grad_for atol=1e-7
                @test grad_rev ≈ grad_FD atol=1e-7
            end

            test_grads([0.3, 0.2, u_n...])
            test_grads([1.1, 0.3, u_n...])
        end
    end
end


text:
include("basic_setup.jl")

using DataFrames

function prepare(mean_s)
    t_s = 100
    t_x = 20

    κ = mean_s / t_s
    λ = 1 / t_s
    ρ = 1 / t_x
    μ = 1 / t_x
    mean_x = mean_s * ρ / μ

    # see Tostevin, ten Wolde, eq. 27
    sigma_squared_ss = mean_s
    sigma_squared_sx = ρ * mean_s / (λ + μ)
    sigma_squared_xx = mean_x * (1 + ρ / (λ + μ))

    joint_stationary = MvNormal([mean_s, mean_x], [sigma_squared_ss sigma_squared_sx; sigma_squared_sx sigma_squared_xx])
    signal_stationary = Poisson(mean_s)

    gen = Trajectories.configuration_generator(sn, rn, [κ, λ], [ρ, μ], signal_stationary, joint_stationary)

    (system, initial) = Trajectories.generate_configuration(gen; duration=100.0)
end

function get_results(algorithm, num_evals, system, initial)
    results = zeros(num_evals)
    Threads.@threads for i in 1:num_evals
        r = Trajectories.simulate(algorithm, initial, system)
        results[i] = Trajectories.log_marginal(r)
    end
    DataFrame(log_marginal=results, κ=system.sparams[1], λ=system.sparams[2], estimate=Trajectories.name(algorithm))
end

an = AnnealingEstimate(15, 50, 100)
ti = TIEstimate(1024, 6, 2^14)
direct = DirectMCEstimate(50_000)

function compute_data()
    df = DataFrame()
    for mean_s in [20, 50, 100, 200]
        system, initial = prepare(mean_s)
        df = vcat(df, get_results(an, 128, system, initial))
        df = vcat(df, get_results(ti, 128, system, initial))
        df = vcat(df, get_results(direct, 128, system, initial))
    end
    df
end

df.mean_s

# df = compute_data()
# using CSV
# CSV.write(projectdir("data", "estimator_variance", "estimator_comparison.csv"), df)

df["mean_s"] = df.κ ./ df.λ
df["mean_s_cat"] = string.(df.mean_s)

df2 = transform(groupby(df, [:mean_s]), :log_marginal => (x-> x .- mean(x[x.estimate == "AIS"])) => :deviation)
show(df2)

using Plots
using StatsPlots
using Statistics

@df df2 dotplot(:mean_s_cat, :deviation, group=(:estimate))

@show combine(groupby(df, [:estimate, :mean_s]), :log_marginal => mean, :log_marginal => std)

jitter_ti = randn(size(ti_samples)...)
jitter_an = randn(size(an_samples)...)
jitter_di = randn(size(di_samples)...)
p = scatter(jitter_ti, ti_samples, label="TI", xlim=(-15, 15))
scatter!(jitter_an, an_samples, label="AIS")
scatter!(jitter_di, di_samples, label="Direct MC")
plot!(p, [-3,3], mean(ti_samples) .* ones(2))
plot!(p, [-3,3], mean(an_samples) .* ones(2))
plot!(p, [-3,3], mean(di_samples) .* ones(2))

density([ti_samples, an_samples, di_samples], label=hcat("TI", "AIS", "Direct MC"))


text:
export BrickletAccelerometerV2Acceleration
struct BrickletAccelerometerV2Acceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickletAccelerometerV2Configuration
struct BrickletAccelerometerV2Configuration
    data_rate::Integer
    full_scale::Integer
end

export BrickletAccelerometerV2AccelerationCallbackConfiguration
struct BrickletAccelerometerV2AccelerationCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletAccelerometerV2ContinuousAccelerationConfiguration
struct BrickletAccelerometerV2ContinuousAccelerationConfiguration
    enable_x::Bool
    enable_y::Bool
    enable_z::Bool
    resolution::Integer
end

export BrickletAccelerometerV2FilterConfiguration
struct BrickletAccelerometerV2FilterConfiguration
    iir_bypass::Integer
    low_pass_filter::Integer
end

export BrickletAccelerometerV2SPITFPErrorCount
struct BrickletAccelerometerV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletAccelerometerV2Identity
struct BrickletAccelerometerV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAccelerometerV2
"""
Measures acceleration in three axis
"""
mutable struct BrickletAccelerometerV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAccelerometerV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_accelerometer_v2")
        deviceInternal = package.BrickletAccelerometerV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the acceleration in x, y and z direction. The values
are given in gₙ/10000 (1gₙ = 9.80665m/s²). The range is
configured with :func:`Set Configuration`.

If you want to get the acceleration periodically, it is recommended
to use the :cb:`Acceleration` callback and set the period with
:func:`Set Acceleration Callback Configuration`.
"""
function get_acceleration(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_acceleration()
end

export set_configuration
"""
    $(SIGNATURES)

Configures the data rate and full scale range.
Possible values are:

* Data rate of 0.781Hz to 25600Hz.
* Full scale range of ±2g up to ±8g.

Decreasing data rate or full scale range will also decrease the noise on
the data.
"""
function set_configuration(device::BrickletAccelerometerV2, data_rate, full_scale)
    device.deviceInternal.set_configuration(data_rate, full_scale)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_configuration()
end

export set_acceleration_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Acceleration`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

If this callback is enabled, the :cb:`Continuous Acceleration 16 Bit` callback
and :cb:`Continuous Acceleration 8 Bit` callback will automatically be disabled.
"""
function set_acceleration_callback_configuration(device::BrickletAccelerometerV2, period, value_has_to_change)
    device.deviceInternal.set_acceleration_callback_configuration(period, value_has_to_change)
end

export get_acceleration_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set Acceleration Callback Configuration`.
"""
function get_acceleration_callback_configuration(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_acceleration_callback_configuration()
end

export set_info_led_config
"""
    $(SIGNATURES)

Configures the info LED (marked as "Force" on the Bricklet) to be either turned off,
turned on, or blink in heartbeat mode.
"""
function set_info_led_config(device::BrickletAccelerometerV2, config)
    device.deviceInternal.set_info_led_config(config)
end

export get_info_led_config
"""
    $(SIGNATURES)

Returns the LED configuration as set by :func:`Set Info LED Config`
"""
function get_info_led_config(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_info_led_config()
end

export set_continuous_acceleration_configuration
"""
    $(SIGNATURES)

For high throughput of acceleration data (> 1000Hz) you have to use the
:cb:`Continuous Acceleration 16 Bit` or :cb:`Continuous Acceleration 8 Bit`
callbacks.

You can enable the callback for each axis (x, y, z) individually and choose a
resolution of 8 bit or 16 bit.

If at least one of the axis is enabled and the resolution is set to 8 bit,
the :cb:`Continuous Acceleration 8 Bit` callback is activated. If at least
one of the axis is enabled and the resolution is set to 16 bit,
the :cb:`Continuous Acceleration 16 Bit` callback is activated.

The returned values are raw ADC data. If you want to put this data into
a FFT to determine the occurrences of specific frequencies we recommend
that you use the data as is. It has all of the ADC noise in it. This noise
looks like pure noise at first glance, but it might still have some frequnecy
information in it that can be utilized by the FFT.

Otherwise you have to use the following formulas that depend on the configured
resolution (8/16 bit) and the full scale range (see :func:`Set Configuration`) to calculate
the data in gₙ/10000 (same unit that is returned by :func:`Get Acceleration`):

* 16 bit, full scale 2g: acceleration = value * 625 / 1024
* 16 bit, full scale 4g: acceleration = value * 1250 / 1024
* 16 bit, full scale 8g: acceleration = value * 2500 / 1024

If a resolution of 8 bit is used, only the 8 most significant bits will be
transferred, so you can use the following formulas:

* 8 bit, full scale 2g: acceleration = value * 256 * 625 / 1024
* 8 bit, full scale 4g: acceleration = value * 256 * 1250 / 1024
* 8 bit, full scale 8g: acceleration = value * 256 * 2500 / 1024

If no axis is enabled, both callbacks are disabled. If one of the continuous
callbacks is enabled, the :cb:`Acceleration` callback is disabled.

The maximum throughput depends on the exact configuration:

.. csv-table::
 :header: "Number of axis enabled", "Throughput 8 bit", "Throughout 16 bit"
 :widths: 20, 20, 20

 "1", "25600Hz", "25600Hz"
 "2", "25600Hz", "15000Hz"
 "3", "20000Hz", "10000Hz"
"""
function set_continuous_acceleration_configuration(device::BrickletAccelerometerV2, enable_x, enable_y, enable_z, resolution)
    device.deviceInternal.set_continuous_acceleration_configuration(enable_x, enable_y, enable_z, resolution)
end

export get_continuous_acceleration_configuration
"""
    $(SIGNATURES)

Returns the continuous acceleration configuration as set by
:func:`Set Continuous Acceleration Configuration`.
"""
function get_continuous_acceleration_configuration(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_continuous_acceleration_configuration()
end

export set_filter_configuration
"""
    $(SIGNATURES)

Configures IIR Bypass filter mode and low pass filter roll off corner frequency.

The filter can be applied or bypassed and the corner frequency can be
half or a ninth of the output data rate.

.. image:: /Images/Bricklets/bricklet_accelerometer_v2_filter.png
   :scale: 100 %
   :alt: Accelerometer filter
   :align: center
   :target: ../../_images/Bricklets/bricklet_accelerometer_v2_filter.png

.. versionadded:: 2.0.2\$nbsp;(Plugin)
"""
function set_filter_configuration(device::BrickletAccelerometerV2, iir_bypass, low_pass_filter)
    device.deviceInternal.set_filter_configuration(iir_bypass, low_pass_filter)
end

export get_filter_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Filter Configuration`.

.. versionadded:: 2.0.2\$nbsp;(Plugin)
"""
function get_filter_configuration(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_filter_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletAccelerometerV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletAccelerometerV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletAccelerometerV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletAccelerometerV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletAccelerometerV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletAccelerometerV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletAccelerometerV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAccelerometerV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAccelerometerV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
## ======================================================================================================
## Remove below before v0.6
## ======================================================================================================

@deprecate setPointsManiPartial!( Mdest::AbstractManifold,
                                  dest,
                                  Msrc::AbstractManifold,
                                  src,
                                  partial::AbstractVector{<:Integer},
                                  asPartial::Bool=true ) setPointPartial!( Mdest, dest, Msrc, src, partial, asPartial )

export productbelief

"""
    $SIGNATURES

Take product of `dens` (including optional partials beliefs) as proposals to be multiplied together.

Notes
-----
- Return points of full dimension, even if only partial dimensions in proposals.
  - 'Other' dimensions left unchanged from incoming `denspts`
- `d` dimensional product approximation
- Incorporate ApproxManifoldProducts to process variables in individual batches.

DevNotes
- TODO Consolidate with [`AMP.manifoldProduct`](@ref), especially concerning partials.
"""
function productbelief( denspts::AbstractVector{P},
                        manifold::MB.AbstractManifold,
                        dens::Vector{<:ManifoldKernelDensity},
                        N::Int;
                        asPartial::Bool=false,
                        dbg::Bool=false,
                        logger=ConsoleLogger()  ) where P
  #

  @warn "productbelief is being deprecated, use manifoldProduct together with getPoints instead."
  mkd = AMP.manifoldProduct(dens, manifold, Niter=1, oldPoints=denspts, logger=logger)
  pGM = getPoints(mkd, asPartial)

  return pGM
end

# function calcMean(mkd::ManifoldKernelDensity{M}) where {M <: ManifoldsBase.AbstractManifold}
#   data = getPoints(mkd)
#   # Returns the mean point on manifold for consitency
#   mean(mkd.manifold, data)
# end

@deprecate calcVariableCovarianceBasic(M::AbstractManifold, vecP::Vector{P}) where P calcCovarianceBasic(M, vecP)

## ======================================================================================================
## Remove below before v0.5
## ======================================================================================================

#
# """
#     $SIGGNATURES
#
# Assemble oplus and ominus operations from given manifols.
#
# Related
#
# buildHybridManifoldCallbacks
# """
# function getManifoldOperations(manis::T) where {T <: Tuple}
#
# end
#

uncoords(::Type{<:MB.AbstractManifold}, w...; kw...) = error("uncoords(::Type{<:MB.AbstractManifold}, w...; kw...) is obsolete, use makePointFromCoords(M::MB.AbstractManifold, w...) instead.")

coords(::Type{<:MB.AbstractManifold}, w...;   kw...) = error("coords(::Type{<:typeof(SpecialEuclidean(2))}, w...; kw...) is obsolete, use makeCoordsFromPoint(M::MB.AbstractManifold, w...) instead.")

# coords(::Type{<:typeof(SpecialEuclidean(2))}, p::ProductRepr) = [p.parts[1][1], p.parts[1][2], atan(p.parts[2][2,1],p.parts[2][1,1])]

# function coords(::Type{<:typeof(SpecialEuclidean(3))}, p::ProductRepr)
#   wELo = TU.convert(Euler, SO3(p.parts[2]))
#   [p.parts[1][1:3]; wELo.R; wELo.P; wELo.Y]
# end

# function uncoords(::Type{<:typeof(SpecialEuclidean(2))}, p::AbstractVector{<:Real}, static::Bool=true)
#   α = p[3]
#   ArrConst = static ? SA : eltype(α)
#   return ProductRepr((ArrConst[p[1], p[2]]), ArrConst[cos(α) -sin(α); sin(α) cos(α)])
# end
# # function uncoords(::Type{<:typeof(SpecialEuclidean(2))}, p::AbstractVector{<:Real})
# #   α = p[3]
# #   return ProductRepr(([p[1], p[2]]), [cos(α) -sin(α); sin(α) cos(α)])
# # end

# function uncoords(::Type{<:typeof(SpecialEuclidean(3))}, p::AbstractVector{<:Real})
#   # α = p[3]
#   wRo = TU.convert(SO3, Euler(p[4:6]...))
#   return ProductRepr(([p[1], p[2], p[3]]), wRo.R)
# end

function getPointsManifold(mkd::ManifoldKernelDensity{M}) where {M <: Euclidean}
  @warn "getPointsManifold is being deprecated, use getPoints(::MKD)::Vector{P} instead"
  data_ = getPoints(mkd.belief)
  TensorCast.@cast data[i][j] := data_[j,i]
  return data
end

function getPointsManifold(mkd::ManifoldKernelDensity{M}) where {M <: Circle}
  @warn "getPointsManifold is being deprecated, use getPoints(::MKD)::Vector{P} instead"
  data_ = getPoints(mkd.belief)
  return data_[:]
end

function getPointsManifold(mkd::ManifoldKernelDensity{M}) where {M <: SpecialEuclidean}
  @warn "getPointsManifold is being deprecated, use getPoints(::MKD)::Vector{P} instead"
  data_ = getPoints(mkd.belief)
  [uncoords(M, view(data_, :, i)) for i in 1:size(data_,2)]
end

@deprecate mmd!(v::AbstractVector{<:Real}, a::AbstractArray,b::AbstractArray,MF::MB.AbstractManifold, w...; kw...) mmd!(MF, v, a, b, w...; kw...)
@deprecate mmd(a::AbstractArray,b::AbstractArray,MF::MB.AbstractManifold, w...; kw...) mmd(MF, a, b, w...; kw...)

# function ker( ::typeof(Euclidean(1)),
#               x::AbstractVector{P1},
#               y::AbstractVector{P2},
#               dx::Vector{<:Real},
#               i::Int,
#               j::Int;
#               sigma::Real=0.001 ) where {P1<:AbstractVector, P2<:AbstractVector}
#   #
#   dx[1] = x[i][1]
#   dx[1] -= y[j][1]
#   dx[1] *= dx[1]
#   dx[1] *= -sigma
#   exp( dx[1] )
# end

# function ker( ::typeof(Euclidean(2)),
#               x::AbstractVector{P1},
#               y::AbstractVector{P2},
#               dx::Vector{<:Real},
#               i::Int,
#               j::Int;
#               sigma::Real=0.001 ) where {P1<:AbstractVector, P2<:AbstractVector}
#   #
#   dx[1] = x[i][1]
#   dx[2] = x[i][2]
#   dx[1] -= y[j][1]
#   dx[2] -= y[j][2]
#   dx .^= 2
#   dx[1] += dx[2]
#   dx[1] *= -sigma
#   exp( dx[1] )
# end

# function ker( ::typeof(SE2_Manifold),
#               x::AbstractVector{P1},
#               y::AbstractVector{P2},
#               dx::Vector{<:Real},
#               i::Int,
#               j::Int;
#               sigma::Real=0.001  ) where {P1<:AbstractVector, P2<:AbstractVector}
#   #
#   innov = se2vee(SE2(x[i][:])\SE2(y[j][:]))
#   exp( -sigma*(  innov'*innov  ) )
# end

# # This functin is still very slow, needs speedup
# # Obviously want to get away from the Euler angles throughout
# function ker( ::typeof(SE3_Manifold),
#               x::AbstractVector{P1},
#               y::AbstractVector{P2},
#               dx::Vector{<:Real},
#               i::Int,
#               j::Int;
#               sigma::Real=0.001  )  where {P1<:AbstractVector, P2<:AbstractVector}
#   #
#   innov = veeEuler(SE3(x[i][1:3],Euler((x[i][4:6])...))\SE3(y[j][1:3],Euler((y[j][4:6])...)))
#   exp( -sigma*(  innov'*innov  ) )
# end

# use Circle1 instead
# struct Circular <: MB.AbstractManifold{MB.ℝ}
#   dof::Int
#   addop::Function
#   diffop::Function
#   getMu
#   getLambda
#   domain::Tuple{Float64, Float64}
# end

# Circular() = Circular(1,
#                       addtheta,
#                       difftheta,
#                       getCircMu,
#                       getCircLambda,
#                       (-pi+0.0,pi-1e-15))

# struct SO2Manifold <: MB.AbstractManifold
# end
#
#
# # should not be defined in AMP, since we want IIF indepent of manifolds
# function *(PP::Vector{MKD{SO2Manifold,B}}) where B
#   @info "taking manifold product of $(length(PP)) terms"
#   @warn "SO2Manifold: work in progress"
# end
#
# mbr1 = ManifoldKernelDensity(SO2Manifold, 0.0)
# mbr2 = ManifoldKernelDensity(SO2Manifold, 0.0)
#
# *([mbr1;mbr2])

# # take the full pGM in, but only update the coordinate dimensions that are actually affected by new information.
# function _partialProducts!( pGM::AbstractVector{P},
#                             partials::Dict{Any, <:AbstractVector{<:ManifoldKernelDensity}},
#                             manifold::MB.AbstractManifold;
#                             inclFull::Bool=true  ) where P <: AbstractVector
#   #
#   # manis = convert(Tuple, manifold)
#   keepold = inclFull ? deepcopy(pGM) : typeof(pGM)()

#   # TODO remove requirement for P <: AbstractVector
#   allPartDimsMask = 0 .== zeros(Int, length(pGM[1]))
#   # FIXME, remove temporary Tuple manifolds method
#   for (dimnum,pp) in partials
#     # mark dimensions getting partial information
#     for d in dimnum
#       allPartDimsMask[d] = true
#     end
#     # change to vector
#     dimv = [dimnum...]
#     # include previous calcs (if full density products were done before partials)
#     partialMani = _buildManifoldPartial(manifold, dimv)
#     # take product of this partial's subset of dimensions
#     partial_GM = AMP.manifoldProduct(pp, partialMani, Niter=1) |> getPoints
#     # partial_GM = AMP.manifoldProduct(pp, (manis[dimv]...,), Niter=1) |> getPoints

#     for i in 1:length(pGM)
#       pGM[i][dimv] = partial_GM[i]
#     end
#   end

#   # multiply together previous full dim and new product of various partials
#   if inclFull
#     partialPts = [pGM[i][dimv] for i in 1:length(pGM)]
#     push!( pp, AMP.manikde!(partialPts, partialMani) )
#   end

#   nothing
# end

@deprecate manikde!( vecP::AbstractVector, M::MB.AbstractManifold ) manikde!(M, vecP)
@deprecate manikde!( vecP::AbstractVector, bw::AbstractVector{<:Real}, M::MB.AbstractManifold ) manikde!(M, vecP, bw)

# function ManifoldKernelDensity( M::MB.AbstractManifold,
#                                 ptsArr::AbstractVector{P} ) where P
#   #
#   # FIXME obsolete
#   arr = Matrix{Float64}(undef, length(ptsArr[1]), length(ptsArr))
#   @cast arr[i,j] = ptsArr[j][i]
#   manis = convert(Tuple, M)
#   bw = getKDEManifoldBandwidths(arr, manis )
#   addopT, diffopT, _, _ = buildHybridManifoldCallbacks(manis)
#   bel = KernelDensityEstimate.kde!(arr, bw, addopT, diffopT)
#   return ManifoldKernelDensity(M, bel)
# end

export ensurePeriodicDomains!
function ensurePeriodicDomains!( pts::AA, manif::T1 ) where {AA <: AbstractArray{Float64,2}, T1 <: Tuple}
  @warn "ensurePeriodicDomains! is being deprecated without replacement."
  i = 0
  for mn in manif
    i += 1
    if manif[i] == :Circular
      pts[i,:] = TUs.wrapRad.(pts[i,:])
    end
  end

  nothing
end

function manifoldProduct( ff::Union{Vector{BallTreeDensity},<:Vector{<:ManifoldKernelDensity}},
  manis::Tuple;
  kwargs... )
#
error("Obsolete, use manifoldProduct(::Vector{MKD}, <:AbstractManifold) instead.\n`::NTuple{Symbol}` for manifolds is outdated, use `getManifold(...)` and `ManfoldsBase.AbstractManifold` instead.")
end

@deprecate getManifolds(::Type{<:T}) where {T <: ManifoldsBase.AbstractManifold} convert(Tuple, T)
getManifolds(::T) where {T <: ManifoldsBase.AbstractManifold} = getManifolds(T)

# struct Euclid2 <: MB.AbstractManifold{MB.ℝ}
#   dof::Int
#   addop::Function
#   diffop::Function
#   getMu
#   getLambda
#   domain::Tuple{Tuple{Float64,Float64},Tuple{Float64, Float64}}
# end

# Euclid2() = Euclid2(2, +, -, KDE.getEuclidMu, KDE.getEuclidLambda, ((-Inf,Inf),(-Inf,Inf)))

# function *(PP::AbstractVector{<:MKD{typeof(EuclideanManifold),BallTreeDensity}})
#   bds = (p->p.belief).(PP)
#   *(bds)
# end

# """
#     $(SIGNATURES)

# Legacy extension of KDE.kde! function to approximate smooth functions based on samples, using likelihood cross validation for bandwidth selection.  This method allows approximation over hybrid manifolds.
# """
# function manikde!(pts::AbstractMatrix{<:Real},
#                   bws::AbstractVector{<:Real},
#                   manifolds::T  ) where {T <: Tuple}
#   #
#   error("NTuple{Symbol} for manifolds definition is now obsolete")
#   addopT, diffopT, _, _ = buildHybridManifoldCallbacks(manifolds)
#   @show size(pts), typeof(pts)
#   KernelDensityEstimate.kde!(pts, bws, addopT, diffopT)
# end

# function manikde!(ptsArr::AbstractVector{P},
#                   bws::AbstractVector{<:Real},
#                   manifolds::T  ) where {P <: AbstractVector, T <: Tuple}
#   #
#   error("NTuple{Symbol} for manifolds definition is now obsolete")
#   pts = Array{Float64,2}(undef, length(ptsArr[1]),length(ptsArr))
#   @cast pts[i,j] = ptsArr[j][i]
#   manikde!(pts, bws, manifolds)
# end

# function manikde!(pts::AbstractMatrix{<:Real},
#                   manifolds::T  ) where {T <: Tuple}
#   #
#   error("NTuple{Symbol} for manifolds definition is now obsolete")
#   bws = getKDEManifoldBandwidths(pts, manifolds)
#   ensurePeriodicDomains!(pts, manifolds)
#   AMP.manikde!(pts, bws, manifolds)
# end

# function manikde!(ptsArr::AbstractVector{P},
#                   manifolds::T  ) where {P <: AbstractVector, T <: Tuple}
#   #
#   error("NTuple{Symbol} for manifolds definition is now obsolete")
#   pts = Array{Float64,2}(undef, length(ptsArr[1]),length(ptsArr))
#   @cast pts[i,j] = ptsArr[j][i]
#   manikde!(pts, manifolds)
# end

@deprecate ManifoldBelief(w...;kw...) ManifoldKernelDensity(w...;kw...)
function ManifoldBelief(::M, mkd::ManifoldKernelDensity{M,T}) where {M <: MB.AbstractManifold{MB.ℝ}, T}
  @warn "ManifoldBelief is deprecated, use ManifoldKernelDensity instead"
  return mkd
end


text:
######################################################################
# cerc_trajectory.jl
# CErC (Constrained Ergodic Control)
# Does NOT assume linear dynamics
# written 2/07/2017 for Dubins stuff
######################################################################

using Convex, SCS

function cerc_trajectory(em::ErgodicManager, tm::TrajectoryManager; verbose::Bool=true, logging::Bool=false, max_iters::Int=30, es_crit::Float64=0.003,right::Bool=true)
	xd0, ud0 = initialize(tm.initializer, em, tm)
	cerc_trajectory(em, tm, xd0, ud0; verbose=verbose, logging=logging, max_iters=max_iters, es_crit=es_crit,right=right)
end

function cerc_trajectory(em::ErgodicManager, tm::TrajectoryManager, xd0::VVF, ud0::VVF; verbose::Bool=true, logging::Bool=false, max_iters::Int=30, es_crit::Float64=0.003, right::Bool=true)

	# let's not overwrite the initial trajectories
	xd = deepcopy(xd0)
	ud = deepcopy(ud0)
	N = tm.N
	start_idx = right ? 1 : 0

	# matrices for gradients
	ad = zeros(tm.dynamics.n, N)
	bd = zeros(tm.dynamics.m, N)

	# prepare variables for convex optimization
	z = Variable(tm.dynamics.n, N+1)
	v = Variable(tm.dynamics.m, N)
	c = Array(Constraint, 0)
	push!(c, z[:,1] == 0)
	#for n = 0:(N-1)
	#	push!(c, z[:,n+1+1] == tm.A*z[:,n+1] + tm.B*v[:,n+1])
	#end

	# prepare for logging if need be
	if logging
		outfile = open("temp.csv", "w")
		save(outfile, xd)
	end

	if verbose; print_header(); end
	i = 1
	not_finished = true
	es = 0.; cs = 0.; ts = 0.; dd = 0.; step_size = 0.
	while not_finished

		# Find gradients, descent step, and step size. Then descend!
		gradients!(ad, bd, em, tm, xd, ud, start_idx)

		# make constraints
		c = Array(Constraint, 0)
		push!(c, z[:,1] == 0)
		for n = 0:(N-1)
			temp_A = eye(3)
			temp_A[1,3] = -tm.h * sin(xd[n+1][3])
			temp_A[2,3] = tm.h * cos(xd[n+1][3])
			A,B = linearize(tm, xd[n+1], ud[n+1])
			push!(c, z[:,n+1+1] == A*z[:,n+1] + B*v[:,n+1])
		end

		zd, vd = convex_descent(ad, bd, N, z, v, c, tm.R[1,1], start_idx)
		step_size = get_step_size(tm.descender, i)
		descend!(xd, ud, zd, vd, step_size, N)

		# compute statistics and report
		es, cs, ts = all_scores(em, tm, xd, ud, start_idx)
		# only temporarily commented out 2/07/2017
		#dd = directional_derivative(ad, bd, zd, vd)
		dd = 0.05
		#sdd = scaled_dd(ad, bd, zd, vd)
		if verbose; step_report(i, es, cs, ts, dd, step_size); end
		if logging; save(outfile, xd); end

		# check convergence
		i += 1
		not_finished = check_convergence(es,es_crit,i,max_iters,verbose)
	end

	# now that we are done, print a special finished report
	if verbose
		print_header()
		if verbose; step_report(i-1, es, cs, ts, dd, step_size); end
	end

	if logging; close(outfile); end

	return xd, ud
end

#function check_convergence(es::Float64, es_crit::Float64, i::Int, max_iters::Int, verbose::Bool)
#	not_finished = true
#	if es < es_crit
#		not_finished = false
#		if verbose
#			println("reached ergodic criterion...")
#		end
#	end
#	if i > max_iters
#		not_finished = false
#		if verbose
#			println("max iterations reached...")
#		end
#	end
#	return not_finished
#end
#
#function step_report(i::Int, es::Float64, cs::Float64, ts::Float64, dd::Float64, step_size::Float64)
#	@printf " %-7i" i
#	@printf " %-14.7f" es
#	@printf " %-14.7f" cs
#	@printf " %-12.7f" ts
#	@printf " %-12.7f" dd
#	@printf " %-7.5f" step_size
#	println()
#end
#
#function print_header()
#	print_dashes()
#	println(" iter  |ergodic score |control score |total score |direc deriv |step size")
#	print_dashes()
#end
#function print_dashes()
#	println("--------------------------------------------------------------------------")
#end
#
#function gradients!(ad::Matrix{Float64}, bd::Matrix{Float64}, em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF, start_idx::Int)
#	ni =  1
#	for n = 0:(tm.N-1)
#		an_x, an_y = compute_ans(em, xd, tm, n, start_idx)
#		ad[1,ni] = an_x
#		ad[2,ni] = an_y
#		# assume that only the first two state variables matter
#		# in fact, this isn't necessary if it's been initalized to zeros...
#		for i = 3:tm.n
#			ad[i,ni] = 0.0
#		end
#
#		#bd[n+1] = tm.h * tm.R * ud[n+1]
#		# TODO: allow for different tm.m
#		bd[1,ni] = tm.h * (tm.R[1,1]*ud[ni][1] + tm.R[1,2]*ud[ni][2])
#		bd[2,ni] = tm.h * (tm.R[2,1]*ud[ni][1] + tm.R[2,2]*ud[ni][2])
#
#		ni += 1
#	end
#end
#
#function compute_ans(em::ErgodicManager, xd::VVF, tm::TrajectoryManager, n::Int, start_idx::Int)
#	xnx = xd[n + start_idx + 1][1]
#	xny = xd[n + start_idx + 1][2]
#	L = em.L
#
#	an_x = 0.0
#	an_y = 0.0
#
#	for k1 = 0:em.K
#		for k2 = 0:em.K
#			hk = em.hk[k1+1,k2+1]
#
#			dFk_dxn1 = -k1*pi*sin(k1*pi*xnx/L)*cos(k2*pi*xny/L) / (hk*L)
#			dFk_dxn2 = -k2*pi*cos(k1*pi*xnx/L)*sin(k2*pi*xny/L) / (hk*L)
#
#			fk = 0.0
#			for i in 0:(tm.N-1)
#				x = xd[i + start_idx + 1]
#				fk += cos(k1*pi*x[1]/L) * cos(k2*pi*x[2]/L) / hk
#			end
#			c = em.Lambdak[k1+1,k2+1] * (tm.h*fk/tm.T - em.phik[k1+1,k2+1])
#			an_x += c*dFk_dxn1
#			an_y += c*dFk_dxn2
#		end
#	end
#	an_x *= 2.0*tm.h
#	an_y *= 2.0*tm.h
#	return an_x, an_y
#end
#
## descent direction using convex optimization
#function convex_descent(a::Matrix{Float64}, b::Matrix{Float64}, N::Int, z::Variable, v::Variable, c::Vector{Constraint}, r::Float64, start_idx::Int)
#
#	# create the problem
#	N_range = (0:N-1) + start_idx + 1
#	problem = minimize(vecdot(a,z[:,N_range]) + vecdot(b,v) + sumsquares(z) + r*sumsquares(v), c)
#
#	# solve the problem
#	solve!(problem, SCSSolver(verbose=0,max_iters=100000))
#
#	return z.value, v.value
#end
#
#
## modifies (xd,ud) by moving step_size in direction (zd,vd)
#function descend!(xd::VVF, ud::VVF, zd::Matrix{Float64}, vd::Matrix{Float64}, step_size::Float64, N::Int)
#	for i = 0:(N-1)
#		ud[i+1][1] += step_size*vd[1,i+1]
#		ud[i+1][2] += step_size*vd[2,i+1]
#		xd[i+1][1] += step_size*zd[1,i+1]
#		xd[i+1][2] += step_size*zd[2,i+1]
#	end
#	xd[N+1][1] += step_size*zd[1,N+1]
#	xd[N+1][2] += step_size*zd[2,N+1]
#end

# called if logging, not meant for general use
#function save(outfile::IOStream, xd::VVF)
#	n = length(xd[1])
#	for xi in xd
#		for i = 1:(n-1)
#			wi = xi[i]
#			write(outfile,"$(xi[i]),")
#		end
#		write(outfile,"$(xi[n])\n")
#	end
#end


text:
struct SBitSet{N,T<:Unsigned}
    pieces::NTuple{N,T}
end

@inline function SBitSet(s::SBitSet{N,T}) where {N,T<:Unsigned}
    SBitSet{N,T}(s.pieces)
end

@inline function Base.length(s::SBitSet{N,T}) where {N,T<:Unsigned}
    sum(count_ones.(s.pieces))
end

@inline function Base.isempty(s::SBitSet{N,T}) where {N,T<:Unsigned}
    all(==(0), s.pieces)
end

@inline function Base.:(==)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    s1.pieces == s2.pieces
end

@inline function Base.:(&)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    SBitSet{N,T}((&).(s1.pieces,s2.pieces))
end
Base.:(∩)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned} = s1 & s2

@inline function Base.:(|)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    SBitSet{N,T}((|).(s1.pieces,s2.pieces))
end
Base.:(∪)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned} = s1 | s2

@inline function Base.:(~)(s1::SBitSet{N,T}) where {N,T<:Unsigned}
    SBitSet{N,T}((~).(s1.pieces))
end
Base.:(~)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned} = s1 ∩ ~s2

@inline function Base.:(⊆)(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    @inbounds @simd for i in 1:N
        !(s1.pieces[i]&s2.pieces[i] == s1.pieces[i]) && return false
    end
    return true
end

@inline function Base.xor(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    SBitSet{N,T}(xor.(s1.pieces,s2.pieces))
end

"""
    _divrem(n,N,)

For ``1≤n≤N*W``, return ``d,r`` with ``1≤r≤W`` and ``W*(d-1)+r=n``.

"""
@inline function _divrem(n::Integer,N::Integer,W::Integer)
   @assert 1 ≤ n ≤ N*W
   (d,r) = divrem(n,W)
    if r == 0
        (d,r) = (d,W)
    else
        (d,r) = (d+1,r)
    end
    return (d,r)
end

@inline function Base.in(x::Integer,s::SBitSet{N,T}) where {N,T<:Unsigned}
    (d,r) = _divrem(x,N,8*sizeof(T))
    Bool((s.pieces[d] >> (r-1)) & T(1))
end

@inline function Base.isdisjoint(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    @inbounds @simd for i in 1:N
        (s1.pieces[i]&s2.pieces[i]) ≠ 0 && return false
    end
    return true
end
⟂(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned} = isdisjoint(s1,s2) # WARNING: ⟂ is \perp, not \bot

@inline function singleton_intersection(s1::SBitSet{N,T},s2::SBitSet{N,T}) where {N,T<:Unsigned}
    ones_tot = 0
    @inbounds for i in 1:N
        ones_tot += count_ones(s1.pieces[i]&s2.pieces[i])
        ones_tot > 1 && return false
    end
    return Bool(ones_tot)
end

@inline function SBitSet{N,T}() where {N,T<:Unsigned}
    SBitSet{N,T}(ntuple(x->T(0),N))
end

@inline function SBitSet{N,T}(n::Integer) where {N,T<:Unsigned}
    (d,r) = _divrem(n,N,8*sizeof(T))
    SBitSet{N,T}(ntuple(x-> x == d ? (T(1) << (r-1)) : T(0),N))
end

@inline function SBitSet{N,T}(nn::Integer...) where {N,T<:Unsigned}
    s = SBitSet{N,T}()
    for i in nn
        s = push(s,i)
    end
    return s
end

@inline function SBitSet{N,T}(v::Vector) where {N,T<:Unsigned}
    s = SBitSet{N,T}()
    for i in v
        s = push(s,i)
    end
    return s
end

@inline function push(s::SBitSet{N,T},n::Integer) where {N,T<:Unsigned}
    (d,r) = _divrem(n,N,8*sizeof(T))
    SBitSet{N,T}(ntuple(x-> x == d ? (@inbounds s.pieces[x]) | (T(1) << (r-1)) : (@inbounds s.pieces[x]),N))
end

@inline function pop(s::SBitSet{N,T},n::Integer) where {N,T<:Unsigned}
    (d,r) = _divrem(n,N,8*sizeof(T))
    SBitSet{N,T}(ntuple(x-> x == d ? (@inbounds s.pieces[x]) & ~(T(1) << (r-1)) : (@inbounds s.pieces[x]),N))
end

@inline function Base.iterate(s::SBitSet{N,T}) where {N,T<:Unsigned}
    Base.iterate(s,(1,1,1, @inbounds s.pieces[1]))
end

@inline function Base.iterate(s::SBitSet{N,T},(d,r,current,current_block)) where {N,T<:Unsigned}
    while d ≤ N
        if current_block ≠ T(0)
            tz = trailing_zeros(current_block)
            return (current+tz, (d,r+tz+1,current+tz+1,current_block>>(tz+1)))
        end
        d += 1
        r = 1
        current = (8*sizeof(T))*(d-1)+1
        d > N && return nothing
        @inbounds current_block = s.pieces[d]
    end

    return nothing
end

Base.eltype(::Type{SBitSet{N,T}}) where {N,T<:Unsigned} = T

function Base.show(io::IO, s::SBitSet{N,T}) where {N,T<:Unsigned}
    print(io,"SBitSet{$N,$T}{")
    for n in s
        print(io,n,",")
    end
    print(io,"}")
end


text:
struct ROCSAlignmentResult{D,S,T,F<:AbstractAffineMap,X<:AbstractGMM{D,S},Y<:AbstractGMM{D,T}} <: AlignmentResults
    x::X
    y::Y
    minimum::T
    tform::F
end

"""
    com = center_of_mass(gmm)

Returns the center of mass of `gmm`, where its first order moments are equal to 0.
"""
function center_of_mass(positions::AbstractMatrix{<:Real}, weights=ones(eltype(positions),size(positions,2)))
    return SVector{size(positions,1)}([sum(hcat([weights[i]*positions[:,i] for i=1:size(positions,2)]...), dims=2)...]/sum(weights))
end

center_of_mass(gaussians::AbstractVector{<:AbstractIsotropicGaussian}) =
    return center_of_mass(hcat([g.μ for g in gaussians]...), [g.ϕ for g in gaussians])

center_of_mass(gmm::AbstractIsotropicGMM) = center_of_mass(gmm.gaussians)

center_of_mass(mgmm::AbstractMultiGMM) = center_of_mass(collect(Iterators.flatten([gmm.gaussians for (k,gmm) in mgmm])))

"""
    m = second_moment(gmm, center, dim1, dim2)

Returns the second order moment of `gmm`
"""
function mass_matrix(positions::AbstractMatrix{<:Real},
                     widths=zeros(eltype(positions),size(positions,2)),
                     weights=ones(eltype(positions),size(positions,2)),
                     center=center_of_mass(positions,weights))
    t = eltype(positions)
    npts = size(positions,2)
    M = fill(zero(t), 3, 3)
    dists = [positions[:,i].-center for i=1:npts]
    # diagonal terms
    for i=1:3
        for j=1:npts
            M[i,i] += weights[j] * (dists[j][i]^2 + widths[j]^2)
        end
    end
    # off-diagonal terms
    for i=1:3
        for j=i+1:3
            for k=1:npts
                inc = weights[k] * dists[k][i] * dists[k][j]
                M[i,j] += inc
                M[j,i] += inc
            end
        end
    end
    return M
end

"""
    tforms = inertial_transforms(gmm)

Returns 10 transformations to put `gmm` in an inertial frame. That is, the mass matrix of the GMM
is made diagonal, and the GMM center of mass is made the origin.
"""
function inertial_transforms(positions::AbstractMatrix{<:Real},
                             widths=zeros(eltype(positions),size(positions,2)),
                             weights=ones(eltype(positions),size(positions,2));
                             invert = false)
    com = center_of_mass(positions, weights)
    massmat = mass_matrix(positions, widths, weights, com)
    evecs = eigvecs(massmat)

    T = eltype(com)
    N = length(com)

    # obtain all rotations that align the mass matrix to the coordinate system
    tforms = AffineMap{SMatrix{N,N,T,N^2},SVector{N,T}}[]

    # first, align the the eigenvectors to the coordinate system axes
    # make sure that a reflection is not performed
    if det(evecs) < 0
        evecs[:,end] = -evecs[:,end]
    end

    push!(tforms, inv(LinearMap(SMatrix{length(com),length(com)}(evecs))) ∘ Translation(-com))

    # then rotate the first transformation about each coordinate axis by -π/2, π/2, and π
    for i=1:3
        axis = zeros(N)
        axis[i] = one(T)
        for angle in (π)
            push!(tforms, LinearMap(AngleAxis(angle, axis...)) ∘  tforms[1])
        end
    end
    if invert
        return [inv(tform) for tform in tforms]
    else
        return tforms
    end
end

inertial_transforms(gaussians; kwargs...
    ) = inertial_transforms(hcat([g.μ for g in gaussians]...), [g.σ for g in gaussians], [g.ϕ for g in gaussians]; kwargs...)

inertial_transforms(gmm::AbstractIsotropicGMM; kwargs...
    ) = inertial_transforms(gmm.gaussians; kwargs...)

inertial_transforms(mgmm::AbstractMultiGMM; kwargs...
    ) = inertial_transforms(collect(Iterators.flatten([gmm.gaussians for (k,gmm) in mgmm])); kwargs...)

"""
    score, tform, nevals = rocs_align_gmms(gmmfixed, gmmmoving; maxevals=1000)

Finds the optimal alignment between the two supplied GMMs using steric multipoles,
based on the [ROCS alignment algorithm.](https://docs.eyesopen.com/applications/rocs/theory/shape_shape.html)
"""
function rocs_align(gmmmoving::AbstractGMM, gmmfixed::AbstractGMM; kwargs...)
    # align both GMMs to their inertial frames
    tformfixed = inertial_transforms(gmmfixed)[1]    # Only need one alignment for the fixed GMM
    tformsmoving = inertial_transforms(gmmmoving)    # Take all 4 rotations for the GMM to be aligned

    # perform local alignment starting at each inertial rotation for the "moving" GMM
    results = Tuple{Float64, NTuple{6,Float64}}[]
    for tformm in tformsmoving
        push!(results, local_align(tformm(gmmmoving), tformfixed(gmmfixed); kwargs...))
    end

    # combine the inertial transform with the subsequent alignment transform for the best result
    minoverlap, mindex = findmin([r[1] for r in results])
    tformmoving = AffineMap(results[mindex][2]...) ∘ tformsmoving[mindex]

    # Apply the inverse of `tformfixed` to the optimized transformation
    alignment_tform = inv(tformfixed) ∘ tformmoving

    # return the result
    return ROCSAlignmentResult(gmmmoving, gmmfixed, minoverlap, alignment_tform)
end


text:
using DataFrames, CSV
data1 = CSV.read("../data/purchaseData.csv", copycols=true)
println("Original value: ", data1.Name[1],"\n")

data2 = data1
data2.Name[1] = "EMILY"
@show data1.Name[1]

data1 = CSV.read("../data/purchaseData.csv", copycols=true)
data2 = copy(data1)
data2.Name[1] = "EMILY"
@show data1.Name[1]
println()

data1 = DataFrame()
data1.X = [[0,1],[100,101]]
data2 = copy(data1)
data2.X[1][1] = -1
@show data1.X[1][1]

data1 = DataFrame(X = [[0,1],[100,101]])
data2 = deepcopy(data1)
data2.X[1][1] = -1
@show data1.X[1][1];

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# ---
# title: 1475. Final Prices With a Special Discount in a Shop
# id: problem1475
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/final-prices-with-a-special-discount-in-a-shop/description/>
# hidden: true
# ---
#
# Given the array `prices` where `prices[i]` is the price of the `ith` item in a
# shop. There is a special discount for items in the shop, if you buy the `ith`
# item, then you will receive a discount equivalent to `prices[j]` where `j` is
# the **minimum**  index such that `j > i` and `prices[j] <= prices[i]`,
# otherwise, you will not receive any discount at all.
#
# _Return an array where the`ith` element is the final price you will pay for
# the `ith` item of the shop considering the special discount._
#
#
#
# **Example 1:**
#
#
#
#     Input: prices = [8,4,6,2,3]
#     Output: [4,2,4,2,3]
#     Explanation:
#     For item 0 with price[0]=8 you will receive a discount equivalent to prices[1]=4, therefore, the final price you will pay is 8 - 4 = 4.
#     For item 1 with price[1]=4 you will receive a discount equivalent to prices[3]=2, therefore, the final price you will pay is 4 - 2 = 2.
#     For item 2 with price[2]=6 you will receive a discount equivalent to prices[3]=2, therefore, the final price you will pay is 6 - 2 = 4.
#     For items 3 and 4 you will not receive any discount at all.
#
#
# **Example 2:**
#
#
#
#     Input: prices = [1,2,3,4,5]
#     Output: [1,2,3,4,5]
#     Explanation: In this case, for all items, you will not receive any discount at all.
#
#
# **Example 3:**
#
#
#
#     Input: prices = [10,1,1,6]
#     Output: [9,0,1,6]
#
#
#
#
# **Constraints:**
#
#   * `1 <= prices.length <= 500`
#   * `1 <= prices[i] <= 10^3`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@adjoint (::Type{T})(args...) where T<:Array = T(args...), Δ -> nothing

@nograd size, length, eachindex, Colon(), findfirst, rand, randn

@adjoint Base.vect(xs...) = Base.vect(xs...), Δ -> (Δ...,)

Base.zero(xs::AbstractArray{Any}) = fill!(similar(xs), nothing)

@adjoint function getindex(xs::Array, i...)
  xs[i...], function (Δ)
    Δ′ = zero(xs)
    Δ′[i...] = Δ
    (Δ′, map(_ -> nothing, i)...)
  end
end

@adjoint! setindex!(xs::AbstractArray, x...) = setindex!(xs, x...),
  _ -> error("Mutating arrays is not supported")

# General

@adjoint collect(x) = collect(x), Δ -> (Δ,)

@adjoint permutedims(xs, dims) = permutedims(xs, dims),
  Δ -> (permutedims(Δ, invperm(dims)), nothing)

@adjoint reshape(xs, dims...) = reshape(xs, dims...),
  Δ -> (reshape(Δ, size(xs)),map(_->nothing,dims)...)

@adjoint function hvcat(rows::Tuple{Vararg{Int}}, xs::T...) where T<:Number
  hvcat(rows, xs...), ȳ -> (nothing, ȳ...)
end

@adjoint function repeat(xs; inner=ntuple(_->1, ndims(xs)), outer=ntuple(_->1, ndims(xs)))
  repeat(xs, inner = inner, outer = outer), function (Δ)
    Δ′ = zero(xs)
    S = size(xs)

    # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
    for (dest_idx, val) in pairs(IndexCartesian(), Δ)
        # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
        # wrap around based on original size S.
        src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
        Δ′[src_idx...] += val
    end
    return (Δ′,)
  end
end

# Reductions

@adjoint function sum(xs::AbstractArray; dims = :)
  if dims === (:)
    sum(xs), Δ -> (FillArray(Δ, size(xs)),)
  else
    sum(xs, dims = dims), Δ -> (similar(xs) .= Δ,)
  end
end

function _forward(cx::Context, ::typeof(sum), f, xs::AbstractArray)
  y, back = forward(cx, (xs -> sum(f.(xs))), xs)
  y, ȳ -> (nothing, nothing, back(ȳ)...)
end

@adjoint function prod(xs; dims = :)
  if dims === (:)
    prod(xs), Δ -> (prod(xs) ./ xs .* Δ,)
  else
    prod(xs, dims = dims),
      Δ -> (reshape(.*(circshift.([reshape(xs, length(xs))], 1:length(xs)-1)...), size(xs)) .* Δ,)
  end
end

@adjoint function maximum(xs; dims = :)
  max, i = findmax(xs, dims = dims)
  max, function (Δ)
    Δ isa Real && Δ <= sqrt(eps(float(Δ))) && return nothing
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

@adjoint function minimum(xs; dims = :)
  min, i = findmin(xs, dims = dims)
  min, function (Δ)
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

# LinAlg

@adjoint a::AbstractVecOrMat * b::AbstractVecOrMat = a * b,
  Δ -> (Δ * transpose(b), transpose(a) * Δ)

@adjoint transpose(x) = transpose(x), Δ -> (transpose(Δ),)
@adjoint Base.adjoint(x) = x', Δ -> (Δ',)
@adjoint parent(x::LinearAlgebra.Adjoint) = parent(x), ȳ -> (LinearAlgebra.Adjoint(ȳ),)

function _kron(mat1::AbstractMatrix,mat2::AbstractMatrix)
    m1, n1 = size(mat1)
    mat1_rsh = reshape(mat1,(1,m1,1,n1))

    m2, n2 = size(mat2)
    mat2_rsh = reshape(mat2,(m2,1,n2,1))

    return reshape(mat1_rsh.*mat2_rsh, (m1*m2,n1*n2))
end

@adjoint kron(a::AbstractMatrix, b::AbstractMatrix) = forward(_kron, a, b)

@adjoint iterate(r::UnitRange, i...) = iterate(r, i...), _ -> nothing


text:
#!/usr/bin/env julia

# Required packages
# -----------------------------------------------------------------------------

using JSON

# Function definitions
# -----------------------------------------------------------------------------

# Helper to read and parse a JSON
function read_json(path::String)::Array{Any,1}
  return JSON.parse(read(path, String))
end

# Helper to return SMILES codes, or nothing
function get_smiles(props)
  local subset = filter(p -> p["urn"]["label"] === "SMILES", props)
  if length(subset) === 0
    return nothing
  end
  return map(
    s -> Dict{String,String}(
      "smiles" => s["value"]["sval"],
      "type" => lowercase(s["urn"]["name"])
    ),
    subset
  )
end

# Helper to get props data for the (k, v) pair in the target info
function get_props_data(props, k, v)
  local subset = filter(p -> p["urn"]["label"] === v, props)
  if length(subset) === 0
    return nothing
  end
  return map(s -> s["value"]["sval"], subset)
end

# Abstracted map function to sift and arrange the PubChem data
function sift_data(d, info)
  # Initialise output and shorthand variables
  local out = Dict{String,Any}()
  local props = d["props"]
  local pcid = [string(d["id"]["id"]["cid"])]
  # Set the SMILES code and the PubChem ID
  out["smiles"] = get_smiles(props)
  out["id"] = Dict{String,Array{String,1}}("pubchem" => pcid)
  # Loop over the info key-value pairs
  for (k, v) in info[:main]
    out[k] = get_props_data(props, k, v)
  end
  for (k, v) in info[:id]
    out["id"][k] = get_props_data(props, k, v)
  end
  # If there are any aliases, get sorted and unique only
  if out["alias"] !== nothing
    sort!(unique!(out["alias"]))
  end
  return out
end

# Operations
# -----------------------------------------------------------------------------

# Define paths to use
path = Dict{Symbol,String}(
  :in => "../data/outCollect04Pubchem_rb.json",
  :out => "../data/outArrange04Pubchem.json"
)

# Set the info to look for, and under what name
info = Dict{Symbol,Dict{String,String}}(
  :main => Dict{String,String}(
    "molecular_formula" => "Molecular Formula",
    "alias" => "IUPAC Name",
  ),
  :id => Dict{String,String}(
    "inchi" => "InChI",
    "inchi_key" => "InChIKey"
  )
)

# Get the data to target
data = reduce(
  (a, c) -> push!(a, c["data"]["PC_Compounds"][1]),
  read_json(path[:in])["found"];
  init = Array{Any,1}()
)

# Initialise sifted output
sift = map(d -> sift_data(d, info), data)

# Write out of the park
open(path[:out], "w") do file
  write(file, JSON.json(sift, 2))
end


text:
module Advection

export plankton_advection!
export plankton_diffusion!

export get_xf_index, get_yf_index, get_zf_index

using KernelAbstractions
using CUDA
using Random

using PlanktonIndividuals.Architectures: device, Architecture, rng_type
using PlanktonIndividuals.Grids

include("interpolation.jl")
include("plankton_advection_kernels.jl")
include("plankton_advection.jl")
include("plankton_diffusion.jl")

end


text:
include("FunField/DegreeLocalization.jl")
include("FunField/Round2.jl")
include("FunField/Factor.jl")
include("FunField/HessQR.jl")
include("FunField/IntClsZx.jl")
include("FunField/FakeFracFldMatrix.jl")
include("FunField/FfOrdIdl.jl")
include("FunField/FfOrdFracIdl.jl")


text:
function run_pandamodels_vstab(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)

    result = _run_vstab(
        pm,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
        ext = extract_params!(pm),
    )
    return result
end

function run_pandamodels_multi_vstab(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)
    mn = set_pq_values_from_timeseries(pm)

    result = _run_multi_vstab(
        mn,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
        ext = extract_params!(pm),
    )
    return result
end

function run_pandamodels_qflex(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)

    result = _run_qflex(
        pm,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
        ext = extract_params!(pm),
    )
    return result
end

function run_pandamodels_multi_qflex(json_path)
    pm = load_pm_from_json(json_path)
    active_powermodels_silence!(pm)
    pm = check_powermodels_data!(pm)
    model = get_model(pm["pm_model"])
    solver = get_solver(pm)
    mn = set_pq_values_from_timeseries(pm)

    result = _run_multi_qflex(
        mn,
        model,
        solver,
        setting = Dict("output" => Dict("branch_flows" => true)),
        ext = extract_params!(pm),
    )
    return result
end