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

using Markdown
using InteractiveUtils

# ╔═╡ ae4c61cc-3149-4b10-a61f-125a6fd803a6
begin
	import Pkg;	Pkg.activate(Base.current_project())

	using LaTeXStrings
	using Makie, GLMakie
	using LinearAlgebra, SpecialFunctions
	using Statistics, StatsBase
	using Optim, NLSolversBase
end

# ╔═╡ e09c6169-4c0d-4534-abe9-cf791f185ffb
using NMF

# ╔═╡ f90042ca-4b48-4bd6-aa01-2a2889e3b1b2
import GZip

# ╔═╡ d0b7bd92-8605-11ec-2be3-3be13a4679ff
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 55dbd93d-1d78-4ea0-af78-9410a2a36630
md"""
# scRNAseq Normalization
"""

# ╔═╡ ff88e7f7-e370-41b2-b082-9b250ac0b186
Normalize = ingredients("../src/normalize.jl").Normalize

# ╔═╡ 51c001d2-e85c-4526-b380-044cc26d1e18
scRNA = ingredients("../src/scrna.jl").scRNA

# ╔═╡ b82a3bb4-45b8-4a90-8633-14b2074c9a45
begin

function counts(io::IO)
	genes = String[]
	lines = [
		let
			entry = split(line)
			push!(genes, entry[1])
			parse.(Int, entry[2:end])
		end for line in eachline(io)
	]
	count = Matrix(reduce(hcat,lines)')
	return count, genes
end

function counts(dir::String, subdir::String)
	count = reduce(∪,
		scRNA.load("$dir/$d") for d ∈ readdir(dir) if occursin(subdir,d)
	)

	markers  = (
		yolk = scRNA.locus(count, "sisA", "CG8129", "Corp", "CG8195", "CNT1", "ZnT77C"),
		pole = scRNA.locus(count, "pgc"),
		dvir = scRNA.searchloci(count, "Dvir_")
	)

	count = scRNA.filtercell(count) do cell, _
		sum(cell[markers.dvir]) < .10*sum(cell)
	end

	count = scRNA.filtergene(count) do count, gene
		!occursin("Dvir_", gene) && sum(count) > 20
	end

	markers  = (
		yolk = scRNA.locus(count, "sisA", "CG8129", "Corp", "CG8195", "CNT1", "ZnT77C"),
		pole = scRNA.locus(count, "pgc"),
	)

	count = scRNA.filtercell(count) do cell, _
		(sum(cell[markers.yolk]) < 10
	  && sum(cell[markers.pole]) < 3
	  && sum(cell) > 1e4)
	end

	count = scRNA.filtergene(count) do count, _
		sum(count) > 20
	end

	return count.data, count.gene
end

end

# ╔═╡ deb992f3-72ad-48cb-92f4-b1abd04cfdcb
#raw, gene = counts("/home/nolln/mnt/data/drosophila/raw", "rep"); size(raw)
raw, gene = open(counts,"../data/raw/dge_raw.txt"); size(raw)

# ╔═╡ 8210f158-1d58-44e0-8733-586c584be4e1
md"""
## Negative binomial fit
As seen below, we have huge variability across sequencing depth and gene expression.
We fit a GLM negative-binomial model to each gene, controlling for the confounding efffects of sequencing depth.
"""

# ╔═╡ ba552cce-953f-4f62-bec1-d7c2e9358e6f
let
	fig = Figure(font="Latin Modern Math", fontsize=26)

	ax1 = Axis(
		fig[1,1],
		xscale=Makie.log10,
		xlabel="sequencing depth",
		ylabel="cumulative"
	)

	ax2 = Axis(
		fig[1,2],
		xscale=Makie.log10,
		xlabel="gene counts",
		ylabel="cumulative"
	)

	depth = vec(sum(raw,dims=1))
	level = vec(sum(raw,dims=2))

	lines!(ax1, sort(depth), range(0,1,size(raw,2)),
		linewidth=2,
		color=:black
	)

	lines!(ax2, sort(level), range(0,1,size(raw,1)),
		linewidth=2,
		color=:black
	)

	fig
end

# ╔═╡ 51cde4a1-28ab-413f-8627-d96f2c55be12
model = Normalize.glm(raw);

# ╔═╡ f1ee2764-c1af-4c9d-ad1e-35d629348fe7
md"""
### Filtering

We see a non-monotonic relationship between the uncertainty of our fits and the average expression for each gene.
We filter on the total uncertainty to remove poorly estimated genes
"""

# ╔═╡ 65e87444-f8eb-4596-8af2-4c08995db91d
let
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(raw) ]

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.log10,
		yscale=Makie.log10,
		xlabel=L"\chi",
		ylabel=L"Tr[\delta\Theta]",
	)

	scatter!(ax1, χ, model.δΘ₁ .+ model.δΘ₂ .+ model.δΘ₃,
		color=model.likelihood,
		colormap=:inferno
	)
	hlines!(ax1, [1e1], color=:turquoise, linewidth=4, linestyle=:dashdot)
	fig
end

# ╔═╡ 593f3381-ac7a-4148-8733-99e0791ea587
md"""
#### Componentwise
For completeness, we break it down into components.
Each gene is colored by $\chi$, the scalar measure of sequencing depth
"""

# ╔═╡ 3b3c02da-d5a2-4304-acef-2c1aca63cb03
let
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(raw) ]

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		yscale=Makie.log10,
		xlabel=L"\Theta_1",
		ylabel=L"\delta\Theta_1",
	)
	scatter!(ax1, model.Θ₁, model.δΘ₁, color=log.(χ), colormap=:inferno)

	fig
end

# ╔═╡ b7d31f65-3b5e-4f64-8382-91ecd2d19c38
let
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(raw) ]

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		yscale=Makie.log10,
		xlabel=L"\Theta_2",
		ylabel=L"\delta\Theta_2",
	)
	scatter!(ax1, model.Θ₂, model.δΘ₂, color=log.(χ), colormap=:inferno)
	fig
end

# ╔═╡ fea31b0e-052d-401c-8814-68da136bdfc7
let
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(raw) ]

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.log10,
		yscale=Makie.log10,
		xlabel=L"\Theta_3",
		ylabel=L"\delta\Theta_3",
	)
	scatter!(ax1, model.Θ₃, model.δΘ₃, color=log.(χ), colormap=:inferno)
	fig
end

# ╔═╡ 9d08d1ea-dc65-41ec-991e-e9229ce40cb4
md"""
#### Intuition
We see that the poorly fit genes are basically bernouilli trials.
Below we plot a representative sample of high uncertainty genes.
As seen, these genes have rather uninformative distributions and thus we can "throw" away the associated rows with little loss of information.
"""

# ╔═╡ def1d949-c5ef-4f5a-999c-42db1e6fa305
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="gene expression",
		ylabel="cumulative"
	)

	for g in findall(model.δΘ₃ .> 1e5)
		lines!(ax1,
			sort(vec(raw[g,:])), range(0,1,size(raw,2)),
			color=(:black,.5),
		)
	end
	fig
end

# ╔═╡ fa4c4080-b794-4b69-bdee-b941d83e80c1
md"""
#### Scaling
We see an odd scaling with respect to gene counts and the estimated sequencing depth.
It scales like $n \sim \chi^{1.25}$ instead of proportionally.
"""

# ╔═╡ 4a53093f-702d-4c9d-85b0-e5f06bffada8
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="Θ₂",
		ylabel="Cumulative",
	)

	ι = findall(model.δΘ₂ .< 1e-1)
	lines!(ax1,
		sort(model.Θ₂[ι]), range(0,1,length(ι)),
		color=:black,
		linewidth=2,
	)
	fig
end

# ╔═╡ f9469369-cf4e-4e56-8785-44d63e2af91e
md"""
#### Large variablility in Θ₃
Varies over 4 orders of magnitude
"""

# ╔═╡ 4610b310-f2f8-47aa-8158-ecc453c6e145
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="Θ₃",
		xscale=Makie.log10,
		ylabel="Cumulative",
	)

	ι = findall(model.δΘ₃ .< 1e5)
	lines!(ax1,
		sort(model.Θ₃[ι]), range(0,1,length(ι)),
		color=:black,
		linewidth=2,
	)
	fig
end

# ╔═╡ 1fc1bdf0-a03b-4799-9a9c-ff2f7e95a36f
function filter_uncertain(count, gene, model; cutoff=50)
	uncertainty = model.δΘ₁ .+ model.δΘ₂ .+ model.δΘ₃
	keep = findall(
		(uncertainty .<= cutoff) .&
		(vec(sum(count .== 0, dims=2)) .<= .99*size(count,2))
	)

	return count[keep,:], gene[keep], (
		likelihood = model.likelihood[keep],
		residual   = model.residual[keep, :],

		Θ₁ = model.Θ₁[keep],
		Θ₂ = model.Θ₂[keep],
		Θ₃ = model.Θ₃[keep],

		δΘ₁ = model.δΘ₁[keep],
		δΘ₂ = model.δΘ₂[keep],
		δΘ₃ = model.δΘ₃[keep],
	)
end

# ╔═╡ 77e68c01-447c-4bd9-b0a8-f6b82269c8e4
Counts, Genes, Model = filter_uncertain(raw, gene, model; cutoff=10); size(Counts)

# ╔═╡ d43e0889-6b07-48a5-bbc1-20c37089b50b
md"""
## Verification
To ensure our fits are adequate, we verify the parameter estimates via bootstrap.
"""

# ╔═╡ 5e4f0f97-80aa-4b09-8840-285428bb1faa
bootstrap, sampled_indices = let
	χ = [ cell |> vec |> Normalize.logmean for cell in eachcol(Counts) ]
	ι = 1:10:size(Counts,1)
	data = Array{Tuple}(undef,length(ι))

	Threads.@threads for (n,i) in collect(enumerate(ι))
		data[n] = Normalize.bootstrap(vec(Counts[i,:]), χ)
	end

	data,ι
end;

# ╔═╡ 5cf21f9c-0f18-4544-910a-53639859373a
let
	ι = sampled_indices
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(Counts) ]

	function kernel(i)
		Θ̂  = [ median(b[i]) for b in bootstrap ]
		δΘ̂ = [ var(b[i])  for b in bootstrap ]
		return Θ̂, δΘ̂
	end

	Θ₁, δΘ₁ = kernel(1)
	Θ₂, δΘ₂ = kernel(2)
	Θ₃, δΘ₃ = kernel(3)

	fig = Figure(font="Latin Modern Math", fontsize=26)

	ax1 = Axis(
		fig[1,1],
		xlabel="MLE Θ₁",
		ylabel="bootstrap Θ₁",
	)
	scatter!(ax1, Model.Θ₁[ι], Θ₁, color=log.(χ[ι]), colormap=:inferno)

	ax2 = Axis(
		fig[1,2],
		xlabel="MLE δΘ₁",
		xscale=Makie.log10,
		ylabel="bootstrap δΘ₁",
		yscale=Makie.log10,
	)
	scatter!(ax2, Model.δΘ₁[ι], δΘ₁, color=log.(χ[ι]), colormap=:inferno)

	ax3 = Axis(
		fig[2,1],
		xlabel="MLE Θ₂",
		ylabel="bootstrap Θ₁",
	)
	scatter!(ax3, Model.Θ₂[ι], Θ₂, color=log.(χ[ι]), colormap=:inferno)

	ax4 = Axis(
		fig[2,2],
		xlabel="MLE δΘ₂",
		ylabel="bootstrap δΘ₂",
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	scatter!(ax4, Model.δΘ₂[ι], δΘ₂, color=log.(χ[ι]), colormap=:inferno)

	ax5 = Axis(
		fig[3,1],
		xlabel="MLE Θ₃",
		ylabel="bootstrap Θ₃",
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	scatter!(ax5, Model.Θ₃[ι], Θ₃, color=log.(χ[ι]), colormap=:inferno)
	ylims!(ax5, (1e-2, 1e2))

	ax6 = Axis(
		fig[3,2],
		xlabel="MLE δΘ₃",
		xscale=Makie.log10,
		ylabel="bootstrap δΘ₃",
		yscale=Makie.log10,
	)
	scatter!(ax6, Model.δΘ₃[ι], δΘ₃, color=log.(χ[ι]), colormap=:inferno)
	ylims!(ax6, (1e-5, 1e2))

	fig
end

# ╔═╡ ebd96e2d-1146-4f57-9fbd-8b0c7496b553
md"""
## Dimensional reduction
With the mean and variance estimated reliably, we now rescale with sinkhorn and provide a rank estimation.

### Random matrix theory
"""

# ╔═╡ c3820cbf-70e8-49e7-bfe7-032a5b39515a
Utility = ingredients("../src/util.jl").Utility

# ╔═╡ 36ba7c41-bc1f-480d-bd3c-dc64c795aea6
μ, σ², u, v = let
	σ² = Counts.*(Counts.+Model.Θ₃) ./ (1 .+ Model.Θ₃)
	u², v², _ = Utility.sinkhorn(σ²)

	μ  = Diagonal(.√u²)*Counts*Diagonal(.√v²)
	σ² = Diagonal(u²)*σ²*Diagonal(v²)
	μ, σ², .√u², .√v²
end;

# ╔═╡ cf324d59-c127-479f-8878-6bd2b17bce2e
Λ = svd(μ);

# ╔═╡ cd718366-107a-440f-a2bc-9f9a6ec7c96d
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		ylabel="rank",
		xlabel="singular value",
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	λ = sum(sqrt.(size(Counts)))
	k = sum(Λ.S .>= λ)

	lines!(ax1, Λ.S, 1:length(Λ.S),
		color=:black,
		linewidth=2,
		label="spectrum",
	)
	vlines!(ax1, [λ],
		color=:red,
		linestyle=:dashdot,
		linewidth=2,
		label="MP max value",
	)

	axislegend("Rank≈$(k)", position = :rt)

	fig
end

# ╔═╡ 7dee363d-b606-4d74-b73c-80978de0ef58
μ̃ = let
	r = nnmf(μ,40;alg=:multmse,maxiter=200,init=:nndsvda)
	r.W*r.H
	#=
	m = Λ.U[:,1:40] * Diagonal(Λ.S[1:40]) * Λ.Vt[1:40,:]
	m[m .< 0] .= 0
	m
	=#
end;

# ╔═╡ 0326c120-1a3f-4ed6-aad8-dab05a1e1adf
MeanModel = Normalize.glm(μ̃.+1e-8; stochastic=Normalize.gamma, ϵ=1e-8);

# ╔═╡ 7acdbbcf-2e99-4664-943f-89787d7901d8
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		ylabel="empirical quantile",
		xlabel="theoretical quantile",
	)

	Z = erf.( (MeanModel.residual) ./ sqrt(2) )
	for i in 1:50:size(MeanModel.residual,1)
		χ = .5*(1 .+ Z[i,:])
		lines!(ax1, sort(χ), range(0,1,length(χ)), color=(:black,.05))
	end
	fig
end

# ╔═╡ af2f9edf-80b9-4b15-8056-b55482c10500
Z = let
#=
	F = svd(Model.residual)
	F.U[:,1:30]*Diagonal(F.S[1:30])*F.Vt[1:30,:]
=#
	#.5*(1 .+ erf.(MeanModel.residual ./ sqrt(2)))
	MeanModel.residual
end;

# ╔═╡ 79b02fc1-fb56-4814-ac19-882decc2f137
μ̂ = let
	μ = Λ.U[:,1:40]*Diagonal(Λ.S[1:40])*Λ.Vt[1:40,:]

# sinkhorn normalization
	u, v, _ = Utility.sinkhorn(μ)
	Diagonal(u)*μ*Diagonal(v)

#= unscaled pearson normalization
	χ = log.([ cell |> vec |> Normalize.logmean for cell in eachcol(n) ])
	norms = [
		let
			n = (1 ./ u[i]) .* gene .* (1 ./ v)
			m = exp.( Model.Θ₁[i] .+ Model.Θ₂[i].*χ )
			s = m.*(1 .+ m./Model.Θ₃[i])

			v = (n .- m) ./ sqrt.(s)
			v[v.<-5]  .= -5
			v[v.>+5]  .= +5
			v
		end for (i,gene) in enumerate(eachrow(μ))
	]
	Matrix(
		reduce(hcat,norms)'
	)
=#
end;

# ╔═╡ 3f2d6b35-7cfc-4c2a-a15e-8ce69afe0bfb
let
	fig = Figure()
	ax1 = Axis(fig[1,1],
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	ax2 = Axis(fig[1,2],
		xscale=Makie.log10,
		yscale=Makie.log10,
	)

	scatter!(ax1, sum(μ,dims=1) |> vec, -Λ.Vt[1,:]|> vec)
	scatter!(ax2, sum(μ,dims=2) |> vec, -Λ.U[:,1]|> vec)

	fig
end

# ╔═╡ 90c0c19f-99bd-4d42-95a9-3c4b872cfcf0
md"""
### Pearson residuals
Repeat the same analysis, but don't use random matrix theory first.
Use our estimated pearson residuals.
"""

# ╔═╡ dff854f7-97e7-402c-9476-1717159ab518
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		ylabel="rank",
		xlabel="singular value",
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	F = svd(Model.residual)
	Λ = sum(sqrt.(size(N)))
	k = sum(F.S .>= Λ)

	lines!(ax1, F.S, 1:length(F.S),
		color=:black,
		linewidth=2,
		label="spectrum",
	)
	vlines!(ax1, [Λ],
		color=:red,
		linestyle=:dashdot,
		linewidth=2,
		label="MP max value",
	)

	axislegend("Rank≈$(k)", position = :lb)

	fig
end

# ╔═╡ 0dec4b8a-5402-43b6-ab51-4d911c8fd4de
md"""
## Spatial inference
"""

# ╔═╡ 8a21602f-24fe-40cd-820f-271bbbe40935
Inference = ingredients("../src/infer.jl").Inference

# ╔═╡ 6bf4f2a5-04f2-43cc-995a-78f1c7d1abb3
database = Inference.virtualembryo(directory="../data/raw");

# ╔═╡ 98506181-1e09-432e-9e47-68eeda650507
inv_z = Inference.inversion(Z, Genes; refdb=database);

# ╔═╡ 5b74e210-9b3d-422c-8e89-525ccdd17c47
ψz = minimum(size(Counts))*inv_z.invert(80);

# ╔═╡ d2503eec-c26a-41a8-835c-2ddb7bafdec3
size(ψz)

# ╔═╡ 41e74b6b-8835-4013-b4eb-9cdc9318f8f0
S = let
	p = ψz ./ sum(ψz,dims=2)
	-sum( p .* log.(p), dims=2)
end

# ╔═╡ 7425feb4-5c17-4ece-a6d3-80d412481b56
Xz = ψz'*database.position;

# ╔═╡ 97affcd3-e919-48d9-a263-4fa4f81dc0a2
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax3 = Axis3(fig[1,1],
		aspect=:data,
		elevation=0,
		azimuth=π/2,
	)

	scatter!(ax3,
		-database.position[:,1],database.position[:,2],database.position[:,3],
		color=ψz*vec(Z[findfirst(Genes.=="odd"),:]),
		colormap=:inferno,
		markersize=2000,
	)
	hidedecorations!(ax3)
	fig
end

# ╔═╡ a803fc43-75de-44eb-9d23-6661fe44a01d
seq = let
	u,v = Utility.sinkhorn(μ)
	Diagonal(u)*μ*Diagonal(v)
end

# ╔═╡ ef87beae-4f24-46d1-87db-72b8e4c2e5aa
inv_m = Inference.inversion(seq, Genes; refdb=database);

# ╔═╡ 786ef949-248f-4079-a226-e2c6c827c0b4
ψm = minimum(size(Counts))*inv_m.invert(80);

# ╔═╡ 9e856728-34d2-4583-9299-828178e97d33
Xm = ψm'*database.position;

# ╔═╡ c82ad97b-2230-4431-8b89-8d51211cc8b0
let
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax3 = Axis3(fig[1,1],
		aspect=:data,
		elevation=0,
		azimuth=π/2,
	)

	scatter!(ax3,
		-database.position[:,1],database.position[:,2],database.position[:,3],
		color=ψm*vec(μ̂[findfirst(Genes.=="eve"),:]),
		colormap=:inferno,
		markersize=2000,
	)
	hidedecorations!(ax3)
	fig
end

# ╔═╡ 10af84b3-e81c-4e2f-b518-950d1b4aac81
md"""
## Manifold learning
### Isomap
"""

# ╔═╡ 6ef6791f-65be-47ca-8330-7f5f16e8f222
PointCloud = ingredients("../src/pointcloud.jl").PointCloud

# ╔═╡ d93b478b-a094-4edd-8aea-b84cd91cfa82
function scaling(D, N; δ=2)
    ϕ, Rs = PointCloud.scaling(D, N)

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xlabel="ball size",
		ylabel="points enclosed",
		xscale=Makie.log10,
		yscale=Makie.log10,
	)
	for i in 1:δ:size(ϕ,1)
	    lines!(ax1, Rs, vec(ϕ[i,:]),
	           color=(:black,.01),
	           label="",
	    )
	end
    lines!(ax1, Rs, 10*size(ϕ,1)*(Rs ./ maximum(Rs)).^3,
           color=:red,
           linestyle=:dashdot,
           linewidth=2,
           label="3d scaling",
    )
	ylims!(ax1, (10, size(D,1)))
    return fig
end

# ╔═╡ 29cc3bbf-5da3-4a0e-9bcc-dd5bc0d3b53e
Dz = PointCloud.geodesics(Z,6);

# ╔═╡ ea0af752-58d4-4399-8d34-7eff34967a03
scaling(Dz,100)

# ╔═╡ c7a2a425-be2b-4a0d-969c-a8642d30ab58
ξz = PointCloud.isomap(Z, 10); size(ξz)

# ╔═╡ e62d07c8-37c0-440c-b520-4a57fdbcf5ef
let
	fig = Figure()
	ax1 = Axis3(fig[1,1])
	ax2 = Axis3(fig[1,2])

	scatter!(ax1, -ξz[:,1], ξz[:,2], ξz[:,3], markersize=4000,color=Xm[:,1], colormap=:inferno)
	scatter!(ax2, -ξz[:,1], ξz[:,2], ξz[:,3], markersize=4000,color=atan.(Xm[:,2],Xm[:,3]), colormap=:inferno)
	fig
end

# ╔═╡ 7cec641c-2311-4e2b-888c-03c7af82e07d
begin
	AP = Xz[:,1] .+ .5*ξz[:,3]
	DV = atan.(Xz[:,2],Xz[:,3]) .+ .013*ξz[:,1]
end;

# ╔═╡ 8cd10af6-db26-4846-b65d-66c0fc47a53e
let
	fig = Figure()
	ax1 = Axis(fig[1,1])
	ax2 = Axis(fig[1,2])

	scatter!(ax1, ξz[:,1], ξz[:,3], color=AP, colormap=:inferno)
	scatter!(ax2, ξz[:,1], ξz[:,3], color=DV, colormap=:inferno)
	fig
end

# ╔═╡ 474c8b35-4eef-4995-88f7-7b0355560a14
let
	c = [
		let
			Dest = PointCloud.Distances.euclidean(ξz[:,1:i]')
			cor(PointCloud.upper_tri(Dest), PointCloud.upper_tri(Dz))
		end	for i in 1:size(ξz,2)
	]

	lines(c)
end

# ╔═╡ 35cdebf3-5400-47cb-8b3d-759451adcef9
Dm = PointCloud.geodesics(μ̂,12);

# ╔═╡ 03ef65ec-6155-4705-8ab8-d6515610a535
scaling(Dm,100)

# ╔═╡ 975ad2b4-9a79-4b6f-ac4e-ff78b4c5fe1e
ξm = PointCloud.isomap(μ̂, 10);

# ╔═╡ 5dac3ca1-f173-4885-ab39-748da6989e41
let
	fig = Figure()
	ax1 = Axis3(fig[1,1])
	ax2 = Axis3(fig[1,2])

	scatter!(ax1, ξm[:,1], -ξm[:,2], ξm[:,3], markersize=4000, color=Xm[:,1], colormap=:inferno)
	scatter!(ax2, ξm[:,1], -ξm[:,2], ξm[:,3], markersize=4000, color=atan.(Xm[:,2],Xm[:,3]), colormap=:inferno)
	fig
end

# ╔═╡ 06c01255-ecb7-4623-ab5b-8072f0d16553
let
	c = [
		let
			Dest = PointCloud.Distances.euclidean(ξm[:,1:i]')
			cor(PointCloud.upper_tri(Dest), PointCloud.upper_tri(Dm))
		end	for i in 1:size(ξm,2)
	]

	lines(c)
end

# ╔═╡ Cell order:
# ╠═ae4c61cc-3149-4b10-a61f-125a6fd803a6
# ╠═f90042ca-4b48-4bd6-aa01-2a2889e3b1b2
# ╟─d0b7bd92-8605-11ec-2be3-3be13a4679ff
# ╟─55dbd93d-1d78-4ea0-af78-9410a2a36630
# ╠═ff88e7f7-e370-41b2-b082-9b250ac0b186
# ╟─51c001d2-e85c-4526-b380-044cc26d1e18
# ╠═b82a3bb4-45b8-4a90-8633-14b2074c9a45
# ╠═deb992f3-72ad-48cb-92f4-b1abd04cfdcb
# ╟─8210f158-1d58-44e0-8733-586c584be4e1
# ╠═ba552cce-953f-4f62-bec1-d7c2e9358e6f
# ╠═51cde4a1-28ab-413f-8627-d96f2c55be12
# ╟─f1ee2764-c1af-4c9d-ad1e-35d629348fe7
# ╠═65e87444-f8eb-4596-8af2-4c08995db91d
# ╟─593f3381-ac7a-4148-8733-99e0791ea587
# ╠═3b3c02da-d5a2-4304-acef-2c1aca63cb03
# ╠═b7d31f65-3b5e-4f64-8382-91ecd2d19c38
# ╟─fea31b0e-052d-401c-8814-68da136bdfc7
# ╟─9d08d1ea-dc65-41ec-991e-e9229ce40cb4
# ╟─def1d949-c5ef-4f5a-999c-42db1e6fa305
# ╟─fa4c4080-b794-4b69-bdee-b941d83e80c1
# ╠═4a53093f-702d-4c9d-85b0-e5f06bffada8
# ╟─f9469369-cf4e-4e56-8785-44d63e2af91e
# ╠═4610b310-f2f8-47aa-8158-ecc453c6e145
# ╠═1fc1bdf0-a03b-4799-9a9c-ff2f7e95a36f
# ╠═77e68c01-447c-4bd9-b0a8-f6b82269c8e4
# ╟─d43e0889-6b07-48a5-bbc1-20c37089b50b
# ╠═5e4f0f97-80aa-4b09-8840-285428bb1faa
# ╠═5cf21f9c-0f18-4544-910a-53639859373a
# ╟─ebd96e2d-1146-4f57-9fbd-8b0c7496b553
# ╠═e09c6169-4c0d-4534-abe9-cf791f185ffb
# ╠═c3820cbf-70e8-49e7-bfe7-032a5b39515a
# ╠═36ba7c41-bc1f-480d-bd3c-dc64c795aea6
# ╠═cf324d59-c127-479f-8878-6bd2b17bce2e
# ╠═cd718366-107a-440f-a2bc-9f9a6ec7c96d
# ╠═7dee363d-b606-4d74-b73c-80978de0ef58
# ╠═0326c120-1a3f-4ed6-aad8-dab05a1e1adf
# ╠═7acdbbcf-2e99-4664-943f-89787d7901d8
# ╠═af2f9edf-80b9-4b15-8056-b55482c10500
# ╠═79b02fc1-fb56-4814-ac19-882decc2f137
# ╠═3f2d6b35-7cfc-4c2a-a15e-8ce69afe0bfb
# ╟─90c0c19f-99bd-4d42-95a9-3c4b872cfcf0
# ╟─dff854f7-97e7-402c-9476-1717159ab518
# ╟─0dec4b8a-5402-43b6-ab51-4d911c8fd4de
# ╠═8a21602f-24fe-40cd-820f-271bbbe40935
# ╠═6bf4f2a5-04f2-43cc-995a-78f1c7d1abb3
# ╠═98506181-1e09-432e-9e47-68eeda650507
# ╠═5b74e210-9b3d-422c-8e89-525ccdd17c47
# ╠═d2503eec-c26a-41a8-835c-2ddb7bafdec3
# ╠═41e74b6b-8835-4013-b4eb-9cdc9318f8f0
# ╠═7425feb4-5c17-4ece-a6d3-80d412481b56
# ╟─97affcd3-e919-48d9-a263-4fa4f81dc0a2
# ╠═a803fc43-75de-44eb-9d23-6661fe44a01d
# ╠═ef87beae-4f24-46d1-87db-72b8e4c2e5aa
# ╠═786ef949-248f-4079-a226-e2c6c827c0b4
# ╠═9e856728-34d2-4583-9299-828178e97d33
# ╟─c82ad97b-2230-4431-8b89-8d51211cc8b0
# ╟─10af84b3-e81c-4e2f-b518-950d1b4aac81
# ╟─6ef6791f-65be-47ca-8330-7f5f16e8f222
# ╟─d93b478b-a094-4edd-8aea-b84cd91cfa82
# ╠═29cc3bbf-5da3-4a0e-9bcc-dd5bc0d3b53e
# ╠═ea0af752-58d4-4399-8d34-7eff34967a03
# ╠═c7a2a425-be2b-4a0d-969c-a8642d30ab58
# ╠═e62d07c8-37c0-440c-b520-4a57fdbcf5ef
# ╟─7cec641c-2311-4e2b-888c-03c7af82e07d
# ╠═8cd10af6-db26-4846-b65d-66c0fc47a53e
# ╠═474c8b35-4eef-4995-88f7-7b0355560a14
# ╟─35cdebf3-5400-47cb-8b3d-759451adcef9
# ╠═03ef65ec-6155-4705-8ab8-d6515610a535
# ╠═975ad2b4-9a79-4b6f-ac4e-ff78b4c5fe1e
# ╠═5dac3ca1-f173-4885-ab39-748da6989e41
# ╠═06c01255-ecb7-4623-ab5b-8072f0d16553

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
import Base: copy

import LinearAlgebra: norm, normalize
import Statistics: mean

abstract type MosiSystem{T <: MosiVector} end
abstract type Molecule{T} <: MosiSystem{T} end

natoms(s::MosiSystem) = length(positions(s))
positions(::MosiSystem) = error("Unimplemented")
velocities(::MosiSystem{T}) where T = T[]
periods(::MosiSystem{T}) where T = T[]
pbc_box(::MosiSystem{T}) where T = zero(T)
has_pbc(s::MosiSystem{T}) where T = pbc_box(s) ≢ zero(T)
original_positions(s::MosiSystem) = if has_pbc(s)
    original_vectors(
        positions(s), periods(s), pbc_box(s)
    )
else
    positions(s)
end
distance_function(s::MosiSystem) = if has_pbc(s)
    pbc_distance(pbc_box(s))
else
    default_distance
end
update_periods!(s::MosiSystem) = if has_pbc(s)
    update_periods!(positions(s), periods(s), pbc_box(s))
    s
else
    s
end

struct MolecularSystem{T <: MosiVector, AT <: AbstractVector{T}} <: MosiSystem{T}
    positions::AT
    velocities::AT
    periods::AT
    box::T
end

function MolecularSystem(
    rs, vs, box = nothing
)
    N = length(rs)
    @assert N === length(vs)
    T = eltype(rs)
    periods = if isnothing(box)
        box = zero(T)
        T[]
    else
        periods = zeros(T, N)
        update_periods!(rs, periods, box)
        periods
    end
    ms = MolecularSystem(rs, vs, periods, box)
end

positions(s::MolecularSystem) = s.positions
velocities(s::MolecularSystem) = s.velocities
periods(s::MolecularSystem) = s.periods
pbc_box(s::MolecularSystem) = s.box

struct ConfigurationSystem{
    T <: MosiVector, AT1 <: AbstractVector{T}, AT2 <: AbstractVector{T}
} <: MosiSystem{T}
    positions::AT1
    periods::AT2
    box::T
end

function ConfigurationSystem(rs; box = nothing, update_periods = true)
    T = eltype(rs)
    ps, box = if box ≡ nothing || box ≡ zero(T)
        T[], zero(T)
    else
        zeros(T, length(rs)), box
    end
    conf = ConfigurationSystem(rs, ps, box)
    if update_periods
        update_periods!(conf)
    end
    conf
end
function ConfigurationSystem((rs, ps)::Tuple; box, update_periods = true)
    conf = ConfigurationSystem(rs, ps, box)
    if update_periods
        update_periods!(conf)
    end
    conf
end
ConfigurationSystem(s::MosiSystem; box = pbc_box(s)) = ConfigurationSystem(positions(s), periods(s), box)
copy(s::ConfigurationSystem) = ConfigurationSystem(copy(positions(s)), copy(periods(s)), pbc_box(s))

positions(s::ConfigurationSystem) = s.positions
velocities(::ConfigurationSystem{T}) where T = T[]
periods(s::ConfigurationSystem) = s.periods
pbc_box(s::ConfigurationSystem) = s.box

ext: jl


text:
using GraphLasso, DataFrames, Distributions, RDatasets, Base.Test

S_iris = cov(array(dataset("datasets", "iris")[:,1:4]))
α = 1.0
W_R = [1.68569351 0.000000 0.2743155 0.01969989
       0.00000000 1.189979 0.0000000 0.00000000
       0.27431545 0.000000 4.1162779 0.29560940
       0.01969989 0.000000 0.2956094 1.58100626]
Θ_R = [ 0.59973155 0.0000000 -0.03996709  0.00000000
        0.00000000 0.8403507  0.00000000  0.00000000
       -0.03996708 0.0000000  0.24890787 -0.04604166
        0.00000000 0.0000000 -0.04604166  0.64111722]
W_iris, Θ_iris = graphlasso(S_iris, α; tol=1e-6, penalize_diag=true)
@test_approx_eq_eps W_R W_iris 1e-5
@test_approx_eq_eps Θ_R Θ_iris 1e-5


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
# Test summation-by-parts property for our derivative operator
# See Summation by Parts, Projections, and Stability. I Olsson 1995
#--------------------------------------------------------------------

D  = [-1.0  1.0  0.0  0.0;
      -0.5  0.0  0.5  0.0;
       0.0 -0.5  0.0  0.5;
       0.0  0.0 -1.0  1.0]

W  = [ 0.5  0.0  0.0  0.0;
       0.0  1.0  0.0  0.0;
       0.0  0.0  1.0  0.0;
       0.0  0.0  0.0  0.5]

B  = [-1.0  0.0  0.0  0.0;
       0.0  0.0  0.0  0.0;
       0.0  0.0  0.0  0.0;
       0.0  0.0  0.0  1.0]

u = [x^2 for x in linspace(-1, 1, 4)]
v = [x^3 for x in linspace(-1, 1, 4)]
@test W*D + (W*D)' == B
@test u'*W*D*v + u'*(W*D)'*v == u'*B*v

# Now test the spectral discretization operators
@testset "Summation-by-Parts" begin
    @testset "N = $N" for N in 1:6
        D = Float64[chebd(i,j, N) for i in 1:N+1, j in 1:N+1]
        W = diagm(Float64[chebw(i, N) for i in 1:N+1])
        B = zeros(N+1, N+1)
        B[1,1] = 1
        B[end, end] = -1
        u = [x^2 for x in chebgrid(N)]
        v = [x^3 for x in chebgrid(N)]
        @test u'*W*D*v + u'*(W*D)'*v ≈ u'*B*v atol=1e-15
        @test_broken (W*D) + (W*D)' ≈ B atol=1e-15
    end
end;


text:
struct InterfaceTractionOperatorValues
    nn::Any
    np::Any
    pn::Any
    pp::Any
    nnT::Any
    npT::Any
    pnT::Any
    ppT::Any
    function InterfaceTractionOperatorValues(nntop, nptop, pntop, pptop, eta)
        nn = -0.5 * vec(nntop)
        np = -0.5 * vec(nptop)
        pn = -0.5 * vec(pntop)
        pp = -0.5 * vec(pptop)

        nnT = -0.5 * eta * vec(transpose(nntop))
        npT = -0.5 * eta * vec(transpose(nptop))
        pnT = -0.5 * eta * vec(transpose(pntop))
        ppT = -0.5 * eta * vec(transpose(pptop))

        new(nn, np, pn, pp, nnT, npT, pnT, ppT)
    end
end

function coherent_traction_operators(
    basis,
    quad1,
    quad2,
    normals,
    stiffness1,
    stiffness2,
    dim,
    scalearea,
    jac,
    vectosymmconverter,
    eta,
)
    nntop = surface_traction_operator(
        basis,
        quad1,
        quad1,
        normals,
        stiffness1,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    nptop = surface_traction_operator(
        basis,
        quad1,
        quad2,
        normals,
        stiffness2,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    pntop = surface_traction_operator(
        basis,
        quad2,
        quad1,
        normals,
        stiffness1,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    pptop = surface_traction_operator(
        basis,
        quad2,
        quad2,
        normals,
        stiffness2,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )

    return InterfaceTractionOperatorValues(nntop, nptop, pntop, pptop, eta)
end

function component_traction_operators(
    basis,
    quad1,
    quad2,
    components,
    normals,
    stiffness1,
    stiffness2,
    dim,
    scalearea,
    jac,
    vectosymmconverter,
    eta,
)

    nntop = surface_traction_component_operator(
        basis,
        quad1,
        quad1,
        components,
        normals,
        stiffness1,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    nptop = surface_traction_component_operator(
        basis,
        quad1,
        quad2,
        components,
        normals,
        stiffness2,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    pntop = surface_traction_component_operator(
        basis,
        quad2,
        quad1,
        components,
        normals,
        stiffness1,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
    pptop = surface_traction_component_operator(
        basis,
        quad2,
        quad2,
        components,
        normals,
        stiffness2,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )

    return InterfaceTractionOperatorValues(nntop, nptop, pntop, pptop, eta)
end

struct InterfaceMassOperatorValues
    nn::Any
    np::Any
    pn::Any
    pp::Any
    function InterfaceMassOperatorValues(nnmop, npmop, pnmop, ppmop)
        nn = vec(nnmop)
        np = vec(npmop)
        pn = vec(pnmop)
        pp = vec(ppmop)

        new(nn, np, pn, pp)
    end
end

function coherent_mass_operators(basis, quad1, quad2, dim, facescale)
    nn = mass_matrix(basis, quad1, quad1, dim, facescale)
    np = mass_matrix(basis, quad1, quad2, dim, facescale)
    pn = mass_matrix(basis, quad2, quad1, dim, facescale)
    pp = mass_matrix(basis, quad2, quad2, dim, facescale)

    return InterfaceMassOperatorValues(nn, np, pn, pp)
end

function component_mass_operators(
    basis,
    quad1,
    quad2,
    components,
    dim,
    facescale,
)
    nn = component_mass_matrix(basis, quad1, quad1, components, dim, facescale)
    np = component_mass_matrix(basis, quad1, quad2, components, dim, facescale)
    pn = component_mass_matrix(basis, quad2, quad1, components, dim, facescale)
    pp = component_mass_matrix(basis, quad2, quad2, components, dim, facescale)

    return InterfaceMassOperatorValues(nn, np, pn, pp)
end

function assemble_interface_condition!(
    sysmatrix,
    nodeids1,
    nodeids2,
    dofspernode,
    tractionop,
    massop,
)
    assemble_cell_matrix!(sysmatrix, nodeids1, dofspernode, +tractionop.nn)
    assemble_cell_matrix!(sysmatrix, nodeids1, dofspernode, +tractionop.nnT)
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids1,
        nodeids2,
        dofspernode,
        +tractionop.np,
    )
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids1,
        nodeids2,
        dofspernode,
        -tractionop.pnT,
    )

    assemble_cell_matrix!(sysmatrix, nodeids2, dofspernode, -tractionop.pp)
    assemble_cell_matrix!(sysmatrix, nodeids2, dofspernode, -tractionop.ppT)
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids2,
        nodeids1,
        dofspernode,
        -tractionop.pn,
    )
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids2,
        nodeids1,
        dofspernode,
        +tractionop.npT,
    )

    assemble_cell_matrix!(sysmatrix, nodeids1, dofspernode, massop.nn)
    assemble_cell_matrix!(sysmatrix, nodeids2, dofspernode, massop.pp)
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids1,
        nodeids2,
        dofspernode,
        -massop.np,
    )
    assemble_couple_cell_matrix!(
        sysmatrix,
        nodeids2,
        nodeids1,
        dofspernode,
        -massop.pn,
    )
end


text:
using ReadStat
using DataValues
using Test

@testset "ReadStat: $ext files" for (reader, ext) in
    ((read_dta, "dta"),
     (read_sav, "sav"),
     (read_sas7bdat, "sas7bdat"),
     (read_xport, "xpt"))

    dtafile = joinpath(dirname(@__FILE__), "types.$ext")
    rsdf = reader(dtafile)
    data = rsdf.data

    @test length(data) == 6
    @test rsdf.headers == [:vfloat, :vdouble, :vlong, :vint, :vbyte, :vstring]
    @test data[1] == DataValueArray{Float32}([3.14, 7., NA])
    @test data[2] == DataValueArray{Float64}([3.14, 7., NA])
    @test data[3] == DataValueArray{Int32}([2, 7, NA])
    @test data[4] == DataValueArray{Int16}([2, 7, NA])
    @test data[5] == DataValueArray{Int8}([2, 7., NA])
    @test data[6] == DataValueArray{String}(["2", "7", ""])
end


text:
import ForwardDiff
using ForwardDiff.MAX_CHUNK_SIZE
using Base.Test

const XLEN = MAX_CHUNK_SIZE
const YLEN = div(MAX_CHUNK_SIZE, 2) + 1
const X, Y = rand(XLEN), rand(YLEN)
const CHUNK_SIZES = (1, div(MAX_CHUNK_SIZE, 3), div(MAX_CHUNK_SIZE, 2), MAX_CHUNK_SIZE)
const FINITEDIFF_ERROR = 1e-5

# used to test against results calculated via finite difference
test_approx_eps(a::Array, b::Array) = @test_approx_eq_eps a b EPS

##################
# Test Functions #
##################

# f: Number -> Number #
#---------------------#

deriv_test_1(x) = sin(x)^2 / cos(x)^2
deriv_test_2(x) = 2*x + sqrt(x*x*x)
deriv_test_3(x) = 10.31^(x + x) - x
deriv_test_4(x) = x
deriv_test_5(x) = 1

const NUMBER_TO_NUMBER_FUNCS = tuple(deriv_test_1, deriv_test_2,
                                     deriv_test_3, deriv_test_4,
                                     deriv_test_5)

# f: Number -> Array #
#--------------------#

function deriv_test_6(x)
    return reshape([deriv_test_1(x),
                    deriv_test_2(x),
                    deriv_test_3(x),
                    deriv_test_1(x) - deriv_test_2(x),
                    deriv_test_2(x),
                    deriv_test_3(x),
                    deriv_test_2(x),
                    deriv_test_3(x)], 2, 2, 2)
end

const NUMBER_TO_ARRAY_FUNCS = tuple(deriv_test_6)

# f: Vector -> Number #
#---------------------#

function rosenbrock(x::AbstractVector)
    a = one(eltype(x))
    b = 100 * a
    result = zero(eltype(x))
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end

function ackley(x::AbstractVector)
    a, b, c = 20.0, -0.2, 2.0*π
    len_recip = inv(length(x))
    sum_sqrs = zero(eltype(x))
    sum_cos = sum_sqrs
    for i in x
        sum_cos += cos(c*i)
        sum_sqrs += i^2
    end
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + e)
end

self_weighted_logit(x::AbstractVector) = inv(1.0 + exp(-dot(x, x)))

const VECTOR_TO_NUMBER_FUNCS = tuple(rosenbrock, ackley, self_weighted_logit)

# f: Vector -> Vector #
#---------------------#

# Credit for `chebyquad!`, `brown_almost_linear!`, and `trigonometric!` goes to
# Kristoffer Carlsson (@KristofferC) - I (@jrevels) just ported them over from PR #104.
# I've modified these functions to allow for tunable input/output sizes. These changes
# might make these functions incorrect in terms of their original design, but shouldn't
# be too different computationally (which is what we care about for tests/benchmarks).

function chebyquad!(y::Vector, x::Vector)
    tk = 1/length(x)
    for j = 1:length(x)
        temp1 = 1.0
        temp2 = 2x[j]-1
        temp = 2temp2
        for i = 1:length(y)
            y[i] += temp2
            ti = temp*temp2 - temp1
            temp1 = temp2
            temp2 = ti
        end
    end
    iev = -1.0
    for k = 1:length(y)
        y[k] *= tk
        if iev > 0
            y[k] += 1/(k^2-1)
        end
        iev = -iev
    end
    return y
end

chebyquad(x) = (y = zeros(eltype(x), length(Y)); chebyquad!(y, x); return y)

function brown_almost_linear!(y::Vector, x::Vector)
    c = sum(x) - (length(x) + 1)
    for i = 1:(length(x)-1), j = 1:(length(y)-1)
        y[j] += x[i] + c
    end
    y[length(y)] = prod(x) - 1
    return nothing
end

brown_almost_linear(x) = (y = zeros(eltype(x), length(Y)); brown_almost_linear!(y, x); return y)

function trigonometric!(y::Vector, x::Vector)
    for i in x, j in eachindex(y)
        y[j] = cos(i)
    end
    c = sum(y)
    n = length(x)
    for i in x, j in eachindex(y)
        y[j] = sin(i) * y[j] + n - c
    end
    return y
end

trigonometric(x) = (y = zeros(eltype(x), length(Y)); trigonometric!(y, x); return y)

const VECTOR_TO_VECTOR_FUNCS = tuple(tuple(chebyquad!, chebyquad),
                                     tuple(brown_almost_linear!, brown_almost_linear),
                                     tuple(trigonometric!, trigonometric))


text:
using Statistics

function isopening(char)
    if char == '(' || char == '[' || char == '{' || char == '<'
        return true
    else
        return false
    end
end

function isclosing(char)
    if char == ')' || char == ']' || char == '}' || char == '>'
        return true
    else
        return false
    end
end

function iscounterpart(openingchar, char)
    if openingchar == '(' && char == ')'
        return true
    elseif openingchar == '[' && char == ']'
        return true
    elseif openingchar == '{' && char == '}'
        return true
    elseif openingchar == '<' && char == '>'
        return true
    else
        return false
    end
end

function getclosingsequence(opseq)
    seq = reverse(opseq)
    cs = Char[]
    for el in seq
        if el == '('
            cs = push!(cs, ')')
        elseif el == '['
            cs = push!(cs, ']')
        elseif el == '{'
            cs = push!(cs, '}')
        elseif el == '<'
            cs = push!(cs, '>')
        end
    end
    return cs
end

function getPoints(cseq)
    score = 0
    for el in cseq
        score *= 5
        if el == ')'
            score += 1
        elseif el == ']'
            score += 2
        elseif el == '}'
            score += 3
        elseif el == '>'
            score += 4
        end
    end
    return score
end

function solve(filename)
    values = readlines(filename*".txt")
    illegals = []
    remainingopeners = []
    for line in values
        currentchunkopener = 'b'
        openers = Char[]
        firstillegal = nothing
        for el in line
            #println(el)
            #println(currentchunkopener)
            if isnothing(firstillegal) == false
                break
            elseif currentchunkopener == 'b' && isopening(el) == false
                firstillegal = el
                illegals = push!(illegals, el)
                break
            elseif currentchunkopener == 'b'
                openers = push!(openers, el)
                currentchunkopener = el
            else
                if isopening(el)
                    openers = push!(openers, el)
                    currentchunkopener = el
                else
                    if iscounterpart(currentchunkopener, el)
                        pop!(openers)
                        if length(openers) > 0
                            currentchunkopener = last(openers)
                        else
                            currentchunkopener = 'b'
                        end
                    else
                        firstillegal = el
                        break
                    end
                end
            end
        end
        illegals = push!(illegals, firstillegal)
        remainingopeners = push!(remainingopeners, openers)
    end
    points = 0
    for el in illegals
        if isnothing(el) == false
            if el == ')'
                points += 3
            elseif el == ']'
                points += 57
            elseif el == '}'
                points += 1197
            elseif el == '>'
                points += 25137
            end
        end
    end
    println("Points for corrupted: "*string(points))
    deletionindices = Int[]
    for (i,el) in enumerate(illegals)
        if isnothing(el) == false
            deletionindices = push!(deletionindices, i)
        end
    end
    values = deleteat!(values, deletionindices)
    remainingopeners = deleteat!(remainingopeners, deletionindices)
    closingsequences = []
    points = Int[]
    for el in remainingopeners
        cs = getclosingsequence(el)
        closingsequences = push!(closingsequences, cs)
        points = push!(points, getPoints(cs))
    end
    completionpoints = Int(median(sort(points)))
    println("Points for incomplete: "*string(completionpoints))
end

solve("smalltest")
solve("input")


text:
# Connects ASTInterface with CommonRLInterface, allowing compatibility with standard RL-based solvers.

CommonRLInterface.reset!(mdp::ASTMDP) = reset!(mdp.sim)

CommonRLInterface.actions(mdp::ASTMDP) = environment(mdp.sim)

CommonRLInterface.observe(mdp::ASTMDP) = Float32.(observe(mdp.sim))

function CommonRLInterface.act!(mdp::ASTMDP, action::Vector{<:Real})
	env = environment(mdp.sim)
	value = unflatten(mdp, action)
	sample = create_sample(env, value)

	d1 = distance(mdp.sim)
	step!(mdp.sim, value)
	d2 = distance(mdp.sim)
	r = reward(mdp, sample, (d1, d2))
	return r
end

CommonRLInterface.terminated(mdp::ASTMDP) = isterminal(mdp.sim) || isevent(mdp.sim)


text:
using NumPy
using Test

@testset "NumPy.jl" begin
    # Write your tests here.
end


text:
# Test file for 2-d Chebyshev polynomials

using ChebyshevApprox

n1 = 10
n2 = 15

dom_1 = [2.0,-3.0]
dom_2 = [1.5, 0.5]
dom = [dom_1 dom_2]

nodesr_1 = chebyshev_nodes(n1,dom_1)
nodesr_2 = chebyshev_nodes(n2,dom_2)
nodese_1 = chebyshev_extrema(n1,dom_1)
nodese_2 = chebyshev_extrema(n2,dom_2)

yr = zeros(n1,n2)
ye = zeros(n1,n2)

for i = 1:n1
  for j = 1:n2

    yr[i,j] = (nodesr_1[i]+4.0)^0.5+nodesr_1[i]*sqrt(nodesr_2[j])
    ye[i,j] = (nodese_1[i]+4.0)^0.5+nodese_1[i]*sqrt(nodese_2[j])

  end
end

order_tensor = [5, 6]
order_complete = 6

wr_tensor   = chebyshev_weights(yr,[nodesr_1,nodesr_2],order_tensor,dom)
wr_complete = chebyshev_weights(yr,(nodesr_1,nodesr_2),order_complete,dom)
we_tensor   = chebyshev_weights_extrema(ye,(nodese_1,nodese_2),order_tensor,dom)
we_complete = chebyshev_weights_extrema(ye,(nodese_1,nodese_2),order_complete,dom)

point = [1.748, 0.753]

yr_tensor    = chebyshev_evaluate(wr_tensor,point,order_tensor,dom)
yr_complete  = chebyshev_evaluate(wr_complete,point,order_complete,dom)
ye_tensor    = chebyshev_evaluate(we_tensor,point,order_tensor,dom)
ye_complete  = chebyshev_evaluate(we_complete,point,order_complete,dom)

y_actual = (point[1]+4.0)^0.5+point[1]*sqrt(point[2])

println([y_actual yr_tensor yr_complete ye_tensor ye_complete])

yr_grad_tensor      = chebyshev_gradient(wr_tensor,point,order_tensor,dom)
yr_grad_complete    = chebyshev_gradient(wr_complete,point,order_complete,dom)
yr_deriv_tensor_1   = chebyshev_derivative(wr_tensor,point,1,order_tensor,dom)
yr_deriv_complete_1 = chebyshev_derivative(wr_tensor,point,1,order_complete,dom)
ye_grad_tensor      = chebyshev_gradient(we_tensor,point,order_tensor,dom)
ye_grad_complete    = chebyshev_gradient(we_complete,point,order_complete,dom)
ye_deriv_tensor_1   = chebyshev_derivative(we_tensor,point,1,order_tensor,dom)
ye_deriv_complete_1 = chebyshev_derivative(we_complete,point,1,order_complete,dom)

println(yr_grad_tensor, yr_grad_complete, yr_deriv_tensor_1, yr_deriv_complete_1)
println(ye_grad_tensor, ye_grad_complete, ye_deriv_tensor_1, ye_deriv_complete_1)

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
module TestUtils
using Conda
using PyCall
using Pkg

Conda.pip_interop(true)
Conda.pip("install", "SimpleITK")
Conda.pip("install", "pandas")
Conda.pip("install", "pymia")

sitk = pyimport("SimpleITK")
pym = pyimport("pymia")
np= pyimport("numpy")

####### getting example data

## examples taken from https://pymia.readthedocs.io/en/latest/examples.evaluation.basic.html?highlight=dice#Evaluation-of-results

##### downloading data
py"""
import argparse
import io
import os
import urllib.request as request
import zipfile

def main(url: str, data_dir: str):
    print(f'Downloading... ({url})')
    resp = request.urlopen(url)
    zip_ = zipfile.ZipFile(io.BytesIO(resp.read()))
    print(f'Extracting... (to {data_dir})')
    members = zip_.infolist()
    for member in members:
        if member.filename.startswith('Subject_') or member.filename.endswith('.h5'):
            if not os.path.basename(member.filename):
                continue
            zip_.extract(member, data_dir)
"""
py"main"("https://github.com/rundherum/pymia-example-data/releases/download/v0.1.0/example-data.zip", "C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\exampleForTestsData")

py"""
import glob
import os

import numpy as np
import pymia.evaluation.metric as metric
import pymia.evaluation.evaluator as eval_
import pymia.evaluation.writer as writer
import SimpleITK as sitk
data_dir = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\exampleForTestsData'
result_file = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results.csv'
result_summary_file = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results_summary.csv'

def downloadd(url: str, data_dir: str):
    metrics = [metric.DiceCoefficient(), metric.HausdorffDistance(percentile=95, metric='HDRFDST95'), metric.VolumeSimilarity()]
    labels = {1: 'WHITEMATTER',2: 'GREYMATTER'  }
    # get subjects to evaluate
    subject_dirs = [subject for subject in glob.glob(os.path.join(data_dir, '*')) if os.path.isdir(subject) and os.path.basename(subject).startswith('Subject')]
    return subject_dirs
"""

    for subject_dir in subject_dirs:

        subject_id = os.path.basename(subject_dir);ground_truth = sitk.ReadImage(os.path.join(subject_dir, f'{subject_id}_GT.mha'))

        prediction = ground_truth

        for label_val in labels.keys():
            # erode each label we are going to evaluate
            prediction = sitk.BinaryErode(prediction, [1,1,1], sitk.sitkBall, 0, label_val)
"""

 py"main"()

# print('\nSubject-wise results...')
# writer.ConsoleWriter().write(evaluator.results)
"""

data_dir = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\exampleForTestsData'
result_file = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results.csv'
result_summary_file = 'C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results_summary.csv'

end#TestUtils

# Pkg.build("HDF5")
# using HDF5


text:
using AQFED.Math
const SqrtEpsilon = sqrt(eps())
# totalVariance is vol^2*τ
export blackScholesFormula, blackScholesVega

function blackScholesFormula(isCall::Bool, strike::Number, spot::Number, totalVariance::Number, driftDf::Number, discountDf::Number)
    sign = 1
    if !isCall
        sign = -1
    end
    forward = spot / driftDf
    if totalVariance < eps()
        price = discountDf * max(sign * (forward - strike), 0)
        return price
    elseif spot < eps()
        if isCall
            return 0
        else
            return discountDf * strike
        end
    elseif strike < eps()
        if isCall
            return discountDf * forward
        else
            return 0
        end
    else
        sqrtVar = sqrt(totalVariance)
        d1 = log(forward / strike)/sqrtVar + sqrtVar/2
        d2 = d1 - sqrtVar
        nd1 = normcdf(sign * d1)
        nd2 = normcdf(sign * d2)
        price = sign * discountDf * (forward * nd1 - strike * nd2)
        return price
    end
end

function blackScholesVega(strike::Number,
	spot::Number,
	variance::Number,
	driftDf::Number,
	discountDf::Number,
	tte::Number)::Number
	forward = spot / driftDf
	sqrtVar = sqrt(variance)
	d1 = 1.0/sqrtVar*log(forward/strike) + 0.5*sqrtVar
	nd1 = normpdf(d1)
	vega = discountDf * forward * nd1 * sqrt(tte)
	return vega
end


text:
# using RecipesBase
using ImageFiltering
using ImageCore

"""
```
image = pssm(A; exaggeration, resolution)
```
Perceptually Shaded Slope Map by *Pingel, Clarke. 2014* [^pingel2014].

# Output
- `image::Gray{T,2}` Grayscale image

# Arguments
- `A::Array{Real,2}` Input Array
- `exaggeration::Real=2.3` Factor to exaggerate elevation
- `resolution::Real=1.0` Resolution of cell size

[^pingel2014]: Pingel, Thomas, and Clarke, Keith. 2014. ‘Perceptually Shaded Slope Maps for the Visualization of Digital Surface Models’. Cartographica: The International Journal for Geographic Information and Geovisualization 49 (4): 225–40. <https://doi.org/10/ggnthv>.
"""
function pssm(A::AbstractMatrix{<:Real}; exaggeration=2.3, resolution=1.)
    x, y = imgradients(A * exaggeration, ImageFiltering.sobel)
    G = sqrt.(x.^2 .+ y.^2)

    G /= resolution  # account for horizontal resolution

    f = scaleminmax(0, 1)
    clamped = f.(G)
    Gray.(1 .- clamped)
end


text:
using CUDA: CUDNN

# Repeat cudnn calls that fail due to memory issues
macro cudnn_retry(ex)
    quote
        res = CUDNN.@retry_reclaim x->(x ∈ (CUDNN.CUDNN_STATUS_ALLOC_FAILED, CUDNN.CUDNN_STATUS_EXECUTION_FAILED)) begin
            $ex
        end
        if res != CUDNN.CUDNN_STATUS_SUCCESS
            CUDNN.throw_api_error(res)
        end
    end |> esc
end


text:
module StencilTest1

using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DistributedPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    avg = stencil(x->(x[-1]+2*x[0]+x[1])/4.0, t1[:val2])
    return avg
end

using HDF5
using MPI
using Base.Test

function main()

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("test1_1.hdf5", "/userid", [1, 2, 3, 1, 2])
        h5write("test1_1.hdf5", "/val2", [1.2, 2.1, 3.1, 3.2, 1.9])
    end
    a = test1("test1_1.hdf5")
    println(a)

    @test_approx_eq a [1.2, 2.125, 2.875, 2.85, 1.9]

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        rm("test1_1.hdf5")
    end

end

end

StencilTest1.main()


text:
test_LP() = CoreModel(
    zeros(4, 3),
    zeros(4),
    ones(3),
    ones(3),
    ones(3),
    ["r$x" for x = 1:3],
    ["m$x" for x = 1:4],
)

test_simpleLP() = CoreModel(
    [
        1.0 1.0
        -1.0 1.0
    ],
    [3.0, 1.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_simpleLP2() = CoreModel(
    zeros(2, 2),
    [0.0, 0.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_sparseLP() = CoreModel(
    sprand(4000, 3000, 0.5),
    sprand(4000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    ["r$x" for x = 1:3000],
    ["m$x" for x = 1:4000],
)

test_coupledLP() = CoreModelCoupled(
    CoreModel(
        sprand(4000, 3000, 0.5),
        sprand(4000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        ["r$x" for x = 1:3000],
        ["m$x" for x = 1:4000],
    ),
    sprand(2000, 3000, 0.5),
    sprand(2000, 0.5),
    sprand(2000, 0.5),
)

test_toyModel() = CoreModel(
    [
        -1.0 1.0 0.0 0.0 0.0 0.0 0.0
        -2.0 0.0 1.0 0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0 0.0 0.0 -1.0
        0.0 -1.0 0.0 -1.0 0.0 0.0 0.0
        0.0 0.0 -1.0 0.0 -1.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 -1.0 1.0
    ],
    zeros(6),
    [0, 0, 0, 0, 0, 0, 1.0],
    fill(-1000.0, 7),
    fill(1000.0, 7),
    ["r1", "m1t", "m3t", "EX_m1(e)", "EX_m3(e)", "EX_biomass(e)", "biomass1"],
    ["m1[c]", "m3[c]", "m2[c]", "m1[e]", "m3[e]", "biomass[c]"],
)


text:
"""
Struct for Standard American Option

		amOption=AmericanOption(T::num1,K::num2,isCall::Bool=true) where {num1 <: Number,num2 <: Number}

Where:\n
		T	=	Time to maturity of the Option.
		K	=	Strike Price of the Option.
		isCall  = true for CALL, false for PUT.
"""
mutable struct AmericanOption{num1 <: Number, num2 <: Number, numtype <: Number} <: AmericanPayoff{numtype}
    T::num1
    K::num2
    isCall::Bool
    function AmericanOption(T::num1, K::num2, isCall::Bool = true) where {num1 <: Number, num2 <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        elseif K <= 0.0
            error("Strike Price must be positive")
        else
            zero_typed = zero(num1) + zero(num2)
            return new{num1, num2, typeof(zero_typed)}(T, K, isCall)
        end
    end
end

export AmericanOption;

function payout(Sti::numtype_, amPayoff::AmericanOption) where {numtype_ <: Number}
    iscall = amPayoff.isCall ? 1 : -1
    return ((Sti - amPayoff.K) * iscall > 0.0) ? (Sti - amPayoff.K) * iscall : zero(numtype_)
end


text:
using Distributions

"""
    get_energy(pop::binding_sites, emat::Array{T, 2}) where {T<:Real}

Compute binding energies in simple binding energy model. Energy matrix needs to have
dimensions (n, site_length).
"""
function get_energy(pop::simple_sites, emat::Array{T, 2}) where {T<:Real}
    c = length(pop.seqs)
    energy = zeros(T, c)
    for j in 1:c
        for i in eachindex(1:pop.l)
            energy[j] += emat[pop.seqs[j][i], i]
        end
    end
    return energy
end

"""
    get_energy(pop::driver_trailer, emat::Array{T, 2}) where {T<:Real}

Compute binding energies in driver/trailer population.
"""
function get_energy(pop::driver_trailer, emat::Array{T, 2}) where {T<:Real}
    c = length(pop.seqs)
    energy = zeros(Float64, c)
    for j in 1:c
        for i in eachindex(1:pop.l)
            energy[j] += emat[pop.seqs[j][i], pop.driver[i]]
        end
    end
    return energy
end

"""
    get_energy(pop::driver_trailer_l, emat::Array{T, 2}) where {T<:Real}

Compute binding energies in driver/trailer population with different binding lengths.
"""
function get_energy(pop::driver_trailer_l, emat::Array{T, 2}) where {T<:Real}
    c = length(pop.seqs)
    energy = zeros(Float64, c)
    for j in 1:c
        for i in 1:pop.l[j]
            energy[j] += emat[pop.seqs[j][i], pop.driver[i]]
        end
    end
    return energy
end

"""
    get_energy(pop::mono_pop, emat::Array{T, 2}) where {T<:Real}

Compute binding energies in monomorphic population..
"""
function get_energy(pop::mono_pop, emat::Array{T, 2})::T where {T<:Real}

    energy::T = 0
    for i in eachindex(pop.seqs)
        energy += emat[pop.seqs[i], pop.driver[i]]
    end
    return energy
end

function remove_empty!(array_list...)
    main_arr = array_list[1]
    inds = findall(x->x==0, main_arr)
    for array in array_list
        deleteat!(array, inds)
    end
end

"""
    function sample_gen!(pop::populations, fit::fitness_functions, emat::Array{T, 2}; remove=true) where {T<:Real}

Sample a new generation and remove extinct species.
"""
function sample_gen!(pop::populations, fit::fitness_functions, emat::Array{T, 2}; remove=true) where {T<:Real}
    E = get_energy(pop, emat)
    f::Array{Float64, 1} = fitness.(E, fit)
    mean_fitness::Float64 = sum(1 / pop.N * pop.freqs.* f)
    norm::Float64 = sum(pop.freqs .* exp.(f .- mean_fitness))
    probabilities::Array{Float64, 1} = pop.freqs .* exp.(f .- mean_fitness) ./ norm
    pop.freqs = rand(Multinomial(pop.N, probabilities), 1)[:]
    if remove
        remove_empty!(pop.freqs, pop.seqs)
    end
end

"""
    function sample_gen!(pop::driver_trailer, fit::fitness_functions, emat::Array{T, 2}; remove=true) where {T<:Real}

Sample a new generation and remove extinct species.
"""
function sample_gen!(pop::driver_trailer_l, fit::fitness_functions, emat::Array{T, 2}; remove=true) where {T<:Real}
    E = get_energy(pop, emat)
    f::Array{Float64, 1} = fitness.(E, pop.l, fit)
    mean_fitness::Float64 = sum(1 / pop.N * pop.freqs.* f)
    norm::Float64 = sum(pop.freqs .* exp.(f .- mean_fitness))
    probabilities::Array{Float64, 1} = pop.freqs .* exp.(f .- mean_fitness) ./ norm
    pop.freqs = rand(Multinomial(pop.N, probabilities), 1)[:]
    if remove
        remove_empty!(pop.freqs, pop.seqs, pop.l)
    end
end


text:
using MonteCarloX
using Random
using Plots

""" Self-exciting process. Rate jumps up by x_increase on event and subsequently
decays to the baserate with timeconstant tau."""
function hawkes_process(tau::Float64 = 1.0, baserate::Float64 = 0.2, x_increase::Float64 = 0.001;
    plotflag::Bool = false, nSamples::Int = 30, rng::AbstractRNG = MersenneTwister(1000))

    alg = InhomogeneousPoissonPiecewiseDecreasing()

    decay(t, t0, x0) = x0 * exp(-(t - t0) / tau)
    Lambda(t, t0, x0) = baserate + decay(t, t0, x0)

    samples = zeros(nSamples)

    t0 = 0.0
    x0 = 0.0
    for i in 1:nSamples
        t_previous = t0
        t0 += next_time(alg, rng, dt->Lambda(dt, 0.0, x0))
        # The current exponential decays starting value x0 is
        # the decay up to now (t0) plus the increase
        x0 = decay(t0, t_previous, x0) + x_increase
        samples[i] = t0
    end

    if plotflag
        # For plotting we are inefficient; (t > t0) means only causal influence
        ratefun(t) = sum([(t > t0) * decay(t, t0, x_increase) for t0 in samples])
        maxT = samples[end] + 4 * tau
        plotrange = range(0.0, stop = maxT, length = 2000)
        plot(plotrange, map(ratefun, plotrange), grid = false)
        display(vline!(samples)) # plot lines when events happen
    end

    return samples
end


text:
# Copyright 2021, by the California Institute of Technology. ALL RIGHTS
# RESERVED. United States Government Sponsorship acknowledged. Any commercial
# use must be negotiated with the Office of Technology Transfer at the
# California Institute of Technology.
#
# This software may be subject to U.S. export control laws. By accepting this
# software, the user agrees to comply with all applicable U.S. export laws and
# regulations. User has the responsibility to obtain export licenses, or other
# export authority as may be required before exporting such information to
# foreign countries or providing access to foreign persons.

using CryoFaBs
using Healpix
using BlockDiagonals
using Statistics
using Test

@testset "CryoFaBs.jl" begin
    @testset "AngularCryoFaB" begin
        nside = 8
        npix = nside2npix(nside)
        mask = fill(1.0, npix)
        cfb = AngularCryoFaB(mask; nside_hires=npix2nside(length(mask)))
    end

    @testset "Radial Cryofunks" begin
        ℓ = 3.14
        rmin = 500.0
        rmax = 1500.0
        nbins = 50
        Δr = (rmax - rmin) / nbins
        reff = range(rmin + Δr/2, rmax - Δr/2, length=nbins)
        k, g, ginv = CryoFaBs.get_radial_cryofunks(ℓ, reff, Δr)
    end

    @testset "AngRadCryoFaB" begin
        nside = 4
        rmin = 500.0
        rmax = 1500.0
        nbins = 50

        mask = fill(1.0, nside2npix(nside))

        cfb = AngRadCryoFaB(mask, rmin, rmax, nbins)

        for n=1:1:nbins, iL=1:10
            # the mode to excite:
            NL = (iL - 1) * nbins + n
            @show iL,n,NL

            # excite that mode:
            R = blocks(cfb.RadTrans)[iL]
            fofr = R[:,n]
            fofa = cfb.angfab.Trans[:,iL]
            δr = [fofr[N]*fofa[i] for N=1:nbins for i=1:cfb.angfab.npixsurvey]

            # manual transform
            δrlm = cfb.AngTransInv * δr
            δrlm_sort = δrlm[cfb.sortbyl]
            δnlm_new = cfb.RadTransInv * δrlm_sort

            # transform
            δnlm = cfb \ δr

            @test δnlm == δnlm_new

            idx = (1:length(δnlm))
            idx = idx[idx .!= NL]
            prec = √eps(eltype(δnlm))

            @show prec
            @show δnlm[NL]
            @show extrema(δnlm[idx])

            @test δnlm[NL] ≈ 1
            @test all(abs.(extrema(δnlm[idx])) .< prec)
        end

        # compile test
        @time k, ell, W = CryoFaBs.cryo_window(cfb; do_calc=true, lmax=10)
    end

    @testset "estimate_density_contrast" begin
        nside = 8
        rmin = 500.0
        rmax = 1500.0
        nbins = 50
        Ngals = 10^5

        mask = fill(1.0, nside2npix(nside))

        cfb = AngRadCryoFaB(mask, rmin, rmax, nbins)

        rθϕ = fill(NaN32, 3, Ngals)
        @. rθϕ[1,:] = rmin + (rmax - rmin) * rand()
        @. rθϕ[2,:] = π * rand()
        @. rθϕ[3,:] = 2 * π * rand() - π

        println("First:")
        @time n_g = CryoFaBs.calc_numgals_per_cell(rθϕ, cfb)
        println("Second:")
        @time n_g = CryoFaBs.calc_numgals_per_cell(rθϕ, cfb)

        @time n_g1 = CryoFaBs.calc_numgals_per_cell_v1(rθϕ, cfb)
        @test n_g == n_g1
    end

    @testset "cryobin" begin
        println("No binning:")
        kk = collect(0:0.1:1)
        ll = fill(1, length(kk))
        P = collect(1:length(kk))
        @show kk ll P
        kbin, lbin, Pbin = cryobin(P, kk, ll)
        @show kbin lbin Pbin
        @test kbin == kk
        @test lbin == ll
        @test Pbin == P

        println("k binning, ktol=0.15:")
        kk = collect(0:0.1:1)
        ll = fill(1, length(kk))
        P = collect(1:length(kk))
        @show kk ll P
        kbin, lbin, Pbin = cryobin(P, kk, ll; ktol=0.15)
        kbin_th = [middle.(kk[1:2:end-1],kk[2:2:end]); [kk[end]]]
        lbin_th = fill(1, length(kbin))
        Pbin_th = [middle.(P[1:2:end-1],P[2:2:end]); [P[end]]]
        @show kbin lbin Pbin
        @show kbin_th lbin_th Pbin_th
        @test kbin == kbin_th
        @test lbin == lbin_th
        @test Pbin == Pbin_th

        println("l binning, ltol=1.5:")
        kk = fill(0.1, 10)
        ll = 1:length(kk)
        P = 1:length(kk)
        @show kk ll P
        kbin, lbin, Pbin = cryobin(P, kk, ll; ktol=0, ltol=1.5)
        lbin_th = middle.(ll[1:2:end-1],ll[2:2:end])
        kbin_th = fill(0.1, length(lbin))
        Pbin_th = middle.(P[1:2:end-1],P[2:2:end])
        @show kbin lbin Pbin
        @show kbin_th lbin_th Pbin_th
        @test kbin == kbin_th
        @test lbin == lbin_th
        @test Pbin == Pbin_th

        println("k & l binning:")
        kk = collect(0.1:0.1:1)
        ll = collect(1:length(kk))
        P = collect(1:length(kk))
        @show kk ll P
        kbin, lbin, Pbin = cryobin(P, kk, ll; ktol=0.15, ltol=1.5)
        kbin_th = collect(0.15:0.2:1)
        lbin_th = collect(1.5:2:10)
        Pbin_th = collect(1.5:2:10)
        @show kbin lbin Pbin
        @show kbin_th lbin_th Pbin_th
        @test kbin ≈ kbin_th
        @test lbin ≈ lbin_th
        @test Pbin ≈ Pbin_th
    end
end

# vim: set sw=4 et sts=4 :


text:
function sliwak_heuristic(clique::AbstractVector, edge; merge_kwargs)
    clique1 = clique[edge[1]]
    clique2 = clique[edge[2]]
    size_intersect = length(intersect(clique1, clique2))
    size_clique1 = length(clique1)
    size_clique2 = length(clique2)
    size_fusion_clique = size_clique1 + size_clique2 - size_intersect
    nb_lc_new = merge_kwargs["nb_lc"] - size_intersect * (2 * size_intersect + 1)
    q2 = (merge_kwargs["m"] + nb_lc_new)^3
    cost = merge_kwargs["alpha"] * (- size_clique1^3 - size_clique2^3 + size_fusion_clique^3)
           + merge_kwargs["beta"] * q2
    return cost
end

function update_sliwak(;clique, edge, merge_kwargs)
    size_intersect = length(intersect(clique[edge[1]], clique[edge[2]]))
    merge_kwargs["nb_lc"] -= size_intersect * (2 * size_intersect + 1)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using IterativeSolvers
using FactCheck

facts("idfact") do
    srand(1)
    M = randn(4,5)
    k = 3

    F = idfact(M, k, 3)
    @fact vecnorm(F.B*F.P-M) --> less_than_or_equal(2svdvals(M)[k+1])
end


text:
# ---
# title: 1162. As Far from Land as Possible
# id: problem1162
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search, Graph
# link: <https://leetcode.com/problems/as-far-from-land-as-possible/description/>
# hidden: true
# ---
#
# Given an `n x n` `grid` containing only values `0` and `1`, where `0`
# represents water and `1` represents land, find a water cell such that its
# distance to the nearest land cell is maximized, and return the distance. If no
# land or water exists in the grid, return `-1`.
#
# The distance used in this problem is the Manhattan distance: the distance
# between two cells `(x0, y0)` and `(x1, y1)` is `|x0 - x1| + |y0 - y1|`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/03/1336_ex1.JPG)
#
#
#
#     Input: grid = [[1,0,1],[0,0,0],[1,0,1]]
#     Output: 2
#     Explanation: The cell (1, 1) is as far as possible from all the land with distance 2.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/03/1336_ex2.JPG)
#
#
#
#     Input: grid = [[1,0,0],[0,0,0],[0,0,0]]
#     Output: 4
#     Explanation: The cell (2, 2) is as far as possible from all the land with distance 4.
#
#
#
#
# **Constraints:**
#
#   * `n == grid.length`
#   * `n == grid[i].length`
#   * `1 <= n <= 100`
#   * `grid[i][j]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct NetworkProperties
	directed :: Bool
	numConnections :: Union{Int, Nothing}
	meanConnectivity :: Union{Real, Nothing}

	clusteringCoefficients :: Union{Vector{<:Real}, Nothing}
	transitivity :: Union{Real, Nothing}
	shortestPaths :: Vector{Union{Vector{<:Real}, Nothing}}

	σ :: Union{Float64, Nothing}
	ω :: Union{Float64, Nothing}
end

NetworkProperties(N::Integer, directed::Bool) = NetworkProperties(
	directed,
	nothing,
	nothing,
	nothing,
	nothing,
	repeat([nothing], N),
	nothing,
	nothing,
)

NetworkProperties(N::Integer, directed::Bool, numConnections::Integer) = NetworkProperties(
	directed,
	numConnections,
	nothing,
	nothing,
	nothing,
	repeat([nothing], N),
	nothing,
	nothing,
)


text:
export md_hmac, fp_prime_get

# Module initialization function
function __init__()
    # make sure the native lib was configure correctly
    @assert sizeof(Limb) == LIMB_SIZE
    @assert MD_LEN == 32  # SHA256
    @assert BN_SIZE >= 32  # 256-bits
    @assert BN_ST_SIZE == sizeof(BN)
    @assert FP_ST_SIZE == sizeof(FP)
    @assert FP2_ST_SIZE == sizeof(FP2)
    #@assert FP3_ST_SIZE == sizeof(FP3)
    @assert G1_ST_SIZE == sizeof(EP)
    @assert G2_ST_SIZE == sizeof(EP2)
    Threads.@threads for i = 1:Threads.nthreads()
        iszero(ccall((:core_init, LIB), Cint, ())) || error("core_init failed")
        iszero(ccall((:ep_param_set_any_pairf, LIB), Cint, ())) || error("ep_param_set_any_pairf failed")
        @info "$(basename(LIB)) initialised successfully (on Thread $(Threads.threadid()))!"
    end
end

# the data is always last so this should be safe for both BN and BigInt
unsafe_bnptr(bn::BN) = Ptr{Limb}(pointer_from_objref(bn) + sizeof(bn) - sizeof(bn.dp))

function bn_add_dig!(c::BN, a::BN, b::Limb)
    ccall((:bn_add_dig, LIB), Cvoid, (Ref{BN}, Ref{BN}, Limb), c, a, b)
    return c
end
function bn_add!(c::BN, a::BN, b::BN)
    ccall((:bn_add, LIB), Cvoid, (Ref{BN}, Ref{BN}, Ref{BN}), c, a, b)
    return c
end

function bn_mod_basic!(c::BN, a::BN, m::BN)
    ccall((:bn_mod_basic, LIB), Cvoid, (Ref{BN}, Ref{BN}, Ref{BN}), c, a, m)
    return c
end

function bn_mul_dig!(c::BN, a::BN, b::Limb)
    ccall((:bn_mul_dig, LIB), Cvoid, (Ref{BN}, Ref{BN}, Limb), c, a, b)
    return c
end
function bn_mul_comba!(c::BN, a::BN, b::BN)
    ccall((:bn_mul_comba, LIB), Cvoid, (Ref{BN}, Ref{BN}, Ref{BN}), c, a, b)
    return c
end
function bn_neg!(c::BN, a::BN)
    ccall((:bn_neg, LIB), Cvoid, (Ref{BN}, Ref{BN}), c, a)
    return c
end

function bn_rand!(a::BN, bits::Int)
    if bits < 1 || bits > typemax(Cint)
        error("invalid bits specified: $bits")
    end
    ccall((:bn_rand, LIB), Cvoid, (Ref{BN}, Cint, Cint), a, 0, bits)
    return a
end
function bn_read_bin!(a::BN, bin::Vector{UInt8})
    ccall((:bn_read_bin, LIB), Cvoid, (Ref{BN}, Ptr{UInt8}, Cint), a, bin, length(bin))
    return a
end

bn_size_bin(a::BN) = Int(ccall((:bn_size_bin, LIB), Cint, (Ref{BN},), a))

function bn_sub_dig!(c::BN, a::BN, b::Limb)
    ccall((:bn_sub_dig, LIB), Cvoid, (Ref{BN}, Ref{BN}, Limb), c, a, b)
    return c
end
function bn_sub!(c::BN, a::BN, b::BN)
    ccall((:bn_sub, LIB), Cvoid, (Ref{BN}, Ref{BN}, Ref{BN}), c, a, b)
    return c
end

function bn_write_bin!(bin::Vector{UInt8}, a::BN)
    ccall((:bn_write_bin, LIB), Cvoid, (Ptr{UInt8}, Cint, Ref{BN}), bin, length(bin), a)
    return bin
end

function ep_add_basic!(c::EP, a::EP, b::EP)
    ccall((:ep_add_basic, LIB), Cvoid, (Ref{EP}, Ref{EP}, Ref{EP}), c, a, b)
    return c
end

function ep_dbl_basic!(c::EP, a::EP)
    ccall((:ep_dbl_basic, LIB), Cvoid, (Ref{EP}, Ref{EP}), c, a)
    return c
end

function ep_curve_get_gen!(g::EP)
    ccall((:ep_curve_get_gen, LIB), Cvoid, (Ref{EP},), g)
    return g
end

function ep_curve_get_ord!(n::BN)
    ccall((:ep_curve_get_ord, LIB), Cvoid, (Ref{BN},), n)
    return n
end

function ep_map!(p::EP, msg::Vector{UInt8})
    ccall((:ep_map, LIB), Cvoid, (Ref{EP}, Ptr{UInt8}, Cint), p, msg, length(msg))
    return p
end

function ep_mul_lwnaf!(c::EP, a::EP, b::BN)
    ccall((:ep_mul_lwnaf, LIB), Cvoid, (Ref{EP}, Ref{EP}, Ref{BN}), c, a, b)
    return c
end

function ep_neg_basic!(c::EP, a::EP)
    ccall((:ep_neg_basic, LIB), Cvoid, (Ref{EP}, Ref{EP}), c, a)
    return c
end

function ep_read_bin!(a::EP, bin::Vector{UInt8})
    ccall((:ep_read_bin, LIB), Cvoid, (Ref{EP}, Ptr{UInt8}, Cint), a, bin, length(bin))
    return a
end

function ep_rand!(ep::EP)
    ccall((:ep_rand, LIB), Cvoid, (Ref{EP},), ep)
    return ep
end

ep_size_bin(ep::EP, pack::Bool) = Int(ccall((:ep_size_bin, LIB), Cint, (Ref{EP}, Cint), ep, pack ? one(Cint) : zero(Cint)))

function ep_sub_basic!(c::EP, a::EP, b::EP)
    ccall((:ep_sub_basic, LIB), Cvoid, (Ref{EP}, Ref{EP}, Ref{EP}), c, a, b)
    return c
end

function ep_write_bin!(bin::Vector{UInt8}, a::EP, pack::Bool)
    flags = pack ? one(Cint) : zero(Cint)
    ccall((:ep_write_bin, LIB), Cvoid, (Ptr{UInt8}, Cint, Ref{EP}, Cint), bin, length(bin), a, flags)
    return bin
end

function ep2_add_basic!(c::EP2, a::EP2, b::EP2)
    ccall((:ep2_add_basic, LIB), Cvoid, (Ref{EP2}, Ref{EP2}, Ref{EP2}), c, a, b)
    return c
end

function ep2_dbl_basic!(c::EP2, a::EP2)
    ccall((:ep2_dbl_basic, LIB), Cvoid, (Ref{EP2}, Ref{EP2}), c, a)
    return c
end

function ep2_curve_get_gen!(g::EP2)
    ccall((:ep2_curve_get_gen, LIB), Cvoid, (Ref{EP2},), g)
    return g
end

function ep2_curve_get_ord!(n::BN)
    ccall((:ep2_curve_get_ord, LIB), Cvoid, (Ref{BN},), n)
    return n
end

function ep2_map!(p::EP2, msg::Vector{UInt8})
    ccall((:ep2_map, LIB), Cvoid, (Ref{EP2}, Ptr{UInt8}, Cint), p, msg, length(msg))
    return p
end

function ep2_mul_lwnaf!(c::EP2, a::EP2, b::BN)
    ccall((:ep2_mul_lwnaf, LIB), Cvoid, (Ref{EP2}, Ref{EP2}, Ref{BN}), c, a, b)
    return c
end

function ep2_neg_basic!(c::EP2, a::EP2)
    ccall((:ep2_neg_basic, LIB), Cvoid, (Ref{EP2}, Ref{EP2}), c, a)
    return c
end

function ep2_rand!(ep::EP2)
    ccall((:ep2_rand, LIB), Cvoid, (Ref{EP2},), ep)
    return ep
end

function ep2_read_bin!(a::EP2, bin::Vector{UInt8})
    ccall((:ep2_read_bin, LIB), Cvoid, (Ref{EP2}, Ptr{UInt8}, Cint), a, bin, length(bin))
    return a
end

ep2_size_bin(ep::EP2, pack::Bool) = Int(ccall((:ep2_size_bin, LIB), Cint, (Ref{EP2}, Cint), ep, pack ? one(Cint) : zero(Cint)))

function ep2_sub_basic!(c::EP2, a::EP2, b::EP2)
    ccall((:ep2_sub_basic, LIB), Cvoid, (Ref{EP2}, Ref{EP2}, Ref{EP2}), c, a, b)
    return c
end

function ep2_write_bin!(bin::Vector{UInt8}, a::EP2, pack::Bool)
    flags = pack ? one(Cint) : zero(Cint)
    ccall((:ep2_write_bin, LIB), Cvoid, (Ptr{UInt8}, Cint, Ref{EP2}, Cint), bin, length(bin), a, flags)
    return bin
end

function fp_add_basic!(c::FP, a::FP, b::FP)
    ccall((:fp_add_basic, LIB), Cvoid, (Ref{FP}, Ref{FP}, Ref{FP}), c, a, b)
    return c
end

function fp_exp_slide!(c::FP, a::FP, b::BN)
    ccall((:fp_exp_slide, LIB), Cvoid, (Ref{FP}, Ref{FP}, Ref{BN}), c, a, b)
    return c
end

function fp_inv_lower!(c::FP, a::FP)
    ccall((:fp_inv_lower, LIB), Cvoid, (Ref{FP}, Ref{FP}), c, a)
    return c
end

function fp_mul_comba!(c::FP, a::FP, b::FP)
    ccall((:fp_mul_comba, LIB), Cvoid, (Ref{FP}, Ref{FP}, Ref{FP}), c, a, b)
    return c
end

function fp_neg_basic!(c::FP, a::FP)
    ccall((:fp_neg_basic, LIB), Cvoid, (Ref{FP}, Ref{FP}), c, a)
    return c
end

function fp_prime_get()
    fp = ccall((:fp_prime_get, LIB), Ptr{Limb}, ())
    bn = BN(used=FP_SIZE)
    unsafe_copyto!(unsafe_bnptr(bn), fp, bn.used)
    return bn
end

function fp_rand!(fp::FP)
    ccall((:fp_rand, LIB), Cvoid, (Ref{FP},), fp)
    return fp
end

function fp_sqr_comba!(c::FP, a::FP)
    ccall((:fp_sqr_comba, LIB), Cvoid, (Ref{FP}, Ref{FP}), c, a)
    return c
end

function fp_sub_basic!(c::FP, a::FP, b::FP)
    ccall((:fp_sub_basic, LIB), Cvoid, (Ref{FP}, Ref{FP}, Ref{FP}), c, a, b)
    return c
end

function fp_zero(c::FP)
    ccall((:fp_zero, LIB), Cvoid, (Ref{FP},), c)
    return c
end

function fp2_add_integ!(c::FP2, a::FP2, b::FP2)
    ccall((:fp2_add_integ, LIB), Cvoid, (Ref{FP2}, Ref{FP2}, Ref{FP2}), c, a, b)
    return c
end

function fp2_exp!(c::FP2, a::FP2, b::BN)
    ccall((:fp2_exp, LIB), Cvoid, (Ref{FP2}, Ref{FP2}, Ref{BN}), c, a, b)
    return c
end

function fp2_inv!(c::FP2, a::FP2)
    ccall((:fp2_inv, LIB), Cvoid, (Ref{FP2}, Ref{FP2}), c, a)
    return c
end

function fp2_mul_integ!(c::FP2, a::FP2, b::FP2)
    ccall((:fp2_mul_integ, LIB), Cvoid, (Ref{FP2}, Ref{FP2}, Ref{FP2}), c, a, b)
    return c
end

function fp2_neg!(c::FP2, a::FP2)
    ccall((:fp2_neg, LIB), Cvoid, (Ref{FP2}, Ref{FP2}), c, a)
    return c
end

function fp2_rand!(a::FP2)
    ccall((:fp2_rand, LIB), Cvoid, (Ref{FP2},), a)
    return a
end

function fp2_sub_integ!(c::FP2, a::FP2, b::FP2)
    ccall((:fp2_sub_integ, LIB), Cvoid, (Ref{FP2}, Ref{FP2}, Ref{FP2}), c, a, b)
    return c
end

function fp2_sqr_integ!(c::FP2, a::FP2)
    ccall((:fp2_sqr_integ, LIB), Cvoid, (Ref{FP2}, Ref{FP2}), c, a)
    return c
end

function fp2_zero(c::FP2)
    ccall((:fp2_zero, LIB), Cvoid, (Ref{FP2},), c)
    return c
end

function fp6_add!(c::FP6, a::FP6, b::FP6)
    ccall((:fp6_add, LIB), Cvoid, (Ref{FP6}, Ref{FP6}, Ref{FP6}), c, a, b)
    return c
end

function fp6_exp!(c::FP6, a::FP6, b::BN)
    ccall((:fp6_exp, LIB), Cvoid, (Ref{FP6}, Ref{FP6}, Ref{BN}), c, a, b)
    return c
end

function fp6_inv!(c::FP6, a::FP6)
    ccall((:fp6_inv, LIB), Cvoid, (Ref{FP6}, Ref{FP6}), c, a)
    return c
end

function fp6_mul_lazyr!(c::FP6, a::FP6, b::FP6)
    ccall((:fp6_mul_lazyr, LIB), Cvoid, (Ref{FP6}, Ref{FP6}, Ref{FP6}), c, a, b)
    return c
end

function fp6_neg!(c::FP6, a::FP6)
    ccall((:fp6_neg, LIB), Cvoid, (Ref{FP6}, Ref{FP6}), c, a)
    return c
end

function fp6_rand!(a::FP6)
    ccall((:fp6_rand, LIB), Cvoid, (Ref{FP6},), a)
    return a
end

function fp6_sub!(c::FP6, a::FP6, b::FP6)
    ccall((:fp6_sub, LIB), Cvoid, (Ref{FP6}, Ref{FP6}, Ref{FP6}), c, a, b)
    return c
end

function fp6_sqr_lazyr!(c::FP6, a::FP6)
    ccall((:fp6_sqr_lazyr, LIB), Cvoid, (Ref{FP6}, Ref{FP6}), c, a)
    return c
end

function fp6_zero(c::FP6)
    ccall((:fp6_zero, LIB), Cvoid, (Ref{FP6},), c)
    return c
end

function fp12_add!(c::FP12, a::FP12, b::FP12)
    ccall((:fp12_add, LIB), Cvoid, (Ref{FP12}, Ref{FP12}, Ref{FP12}), c, a, b)
    return c
end

function fp12_exp!(c::FP12, a::FP12, b::BN)
    ccall((:fp12_exp, LIB), Cvoid, (Ref{FP12}, Ref{FP12}, Ref{BN}), c, a, b)
    return c
end

function fp12_inv!(c::FP12, a::FP12)
    ccall((:fp12_inv, LIB), Cvoid, (Ref{FP12}, Ref{FP12}), c, a)
    return c
end

function fp12_mul_lazyr!(c::FP12, a::FP12, b::FP12)
    ccall((:fp12_mul_lazyr, LIB), Cvoid, (Ref{FP12}, Ref{FP12}, Ref{FP12}), c, a, b)
    return c
end

function fp12_neg!(c::FP12, a::FP12)
    ccall((:fp12_neg, LIB), Cvoid, (Ref{FP12}, Ref{FP12}), c, a)
    return c
end

function fp12_rand!(a::FP12)
    ccall((:fp12_rand, LIB), Cvoid, (Ref{FP12},), a)
    return a
end

function fp12_sub!(c::FP12, a::FP12, b::FP12)
    ccall((:fp12_sub, LIB), Cvoid, (Ref{FP12}, Ref{FP12}, Ref{FP12}), c, a, b)
    return c
end

function fp12_sqr_lazyr!(c::FP12, a::FP12)
    ccall((:fp12_sqr_lazyr, LIB), Cvoid, (Ref{FP12}, Ref{FP12}), c, a)
    return c
end

function fp12_zero(c::FP12)
    ccall((:fp12_zero, LIB), Cvoid, (Ref{FP12},), c)
    return c
end

function md_hmac(in::Vector{UInt8}, key::Vector{UInt8})
    mac = Vector{UInt8}(undef, MD_LEN)
    ccall((:md_hmac, LIB), Cvoid, (Ref{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint),
        mac, in, length(in), key, length(key))
    return mac
end

function md_map_sh256(msg::Vector{UInt8})
    hash = Vector{UInt8}(undef, 32)
    ccall((:md_map_sh256, LIB), Cvoid, (Ref{UInt8}, Ptr{UInt8}, Cint), hash, msg, length(msg))
    return hash
end

function pp_exp_k12!(c::FP12, a::FP12)
    ccall((:pp_exp_k12, LIB), Cvoid, (Ref{FP12}, Ref{FP12}), c, a)
    return c
end

function pp_map_oatep_k12!(r::FP12, p::EP, q::EP2)
    ccall((:pp_map_oatep_k12, LIB), Cvoid, (Ref{FP12}, Ref{EP}, Ref{EP2}), r, p, q)
    return r
end


text:
function einsum(index,A,B)
  """ Einstein sumation. Limited to produce fourth order tensor from
      two order tensors."""
  # tensor product of two second order tensor A_ij B_kl
  if index == "ijkl"
    C=zeros(Float64,3,3,3,3)
    for i in 1:3
      for j in 1:3
        for k in 1:3
          for l in 1:3
            C[i,j,k,l] += A[i,j]*B[k,l]
          end
        end
      end
    end
  # tensor product of two second order tensor A_ik B_jl with permutation
  elseif index == "ikjl"
    C=zeros(Float64,3,3,3,3)
    for i in 1:3
      for j in 1:3
        for k in 1:3
          for l in 1:3
            C[i,j,k,l] += A[i,k]*B[j,l]
          end
        end
      end
    end
  # tensor product of two second order tensor A_il B_jk with permutation
  elseif index == "iljk"
    C=zeros(Float64,3,3,3,3)
    for i in 1:3
      for j in 1:3
        for k in 1:3
          for l in 1:3
            C[i,j,k,l] += A[i,l]*B[j,k]
          end
        end
      end
    end
  else
    println("Operation not understood.")
    C = Nothing
  end
  return C
end

function Voigt(C)
  """Transform of a fourth-order tensor into a square matrix under Voigt notation.
     Just working for 3dimension."""
  VoigtA = zeros(Float64,6,6)
  VoigtA[1,1] = C[1,1,1,1]
  VoigtA[1,2] = C[1,1,2,2]
  VoigtA[1,3] = C[1,1,3,3]
  VoigtA[1,4] = 0.5*(C[1,1,1,2]+C[1,1,2,1])
  VoigtA[1,5] = 0.5*(C[1,1,1,3]+C[1,1,3,1])
  VoigtA[1,6] = 0.5*(C[1,1,2,3]+C[1,1,3,2])
  VoigtA[2,1] = VoigtA[1,2]
  VoigtA[2,2] = C[2,2,2,2]
  VoigtA[2,3] = C[2,2,3,3]
  VoigtA[2,4] = 0.5*(C[2,2,1,2]+C[2,2,2,1])
  VoigtA[2,5] = 0.5*(C[2,2,1,3]+C[2,2,3,1])
  VoigtA[2,6] = 0.5*(C[2,2,2,3]+C[2,2,3,2])
  VoigtA[3,1] = VoigtA[1,3]
  VoigtA[3,2] = VoigtA[2,3]
  VoigtA[3,3] = C[3,3,3,3]
  VoigtA[3,4] = 0.5*(C[3,3,1,2]+C[3,3,2,1])
  VoigtA[3,5] = 0.5*(C[3,3,1,3]+C[3,3,3,1])
  VoigtA[3,6] = 0.5*(C[3,3,2,3]+C[3,3,3,2])
  VoigtA[4,1] = VoigtA[1,4]
  VoigtA[4,2] = VoigtA[2,4]
  VoigtA[4,3] = VoigtA[3,4]
  VoigtA[4,4] = 0.5*(C[1,2,1,2]+C[1,2,2,1])
  VoigtA[4,5] = 0.5*(C[1,2,1,3]+C[1,2,3,1])
  VoigtA[4,6] = 0.5*(C[1,2,2,3]+C[1,2,3,2])
  VoigtA[5,1] = VoigtA[1,5]
  VoigtA[5,2] = VoigtA[2,5]
  VoigtA[5,3] = VoigtA[3,5]
  VoigtA[5,4] = VoigtA[4,5]
  VoigtA[5,5] = 0.5*(C[1,3,1,3]+C[1,3,3,1])
  VoigtA[5,6] = 0.5*(C[1,3,2,3]+C[1,3,3,2])
  VoigtA[6,1] = VoigtA[1,6]
  VoigtA[6,2] = VoigtA[2,6]
  VoigtA[6,3] = VoigtA[3,6]
  VoigtA[6,4] = VoigtA[4,6]
  VoigtA[6,5] = VoigtA[5,6]
  VoigtA[6,6] = 0.5*(C[2,3,2,3]+C[2,3,3,2])
  return VoigtA
end

function in1d(ar1,ar2)
  """Return a boolean 1D-array with trues where an element of ar2 is presented."""
  ar3 = falses(size(ar1,1))
  for i in 1:size(ar1,1)
    for j in 1:size(ar2,1)
      if ar1[i]==ar2[j]
        ar3[i] = true
      end
    end
  end
  return ar3
end

function delete1d(ar1,ar2)
  """Remove elements in ar1 from de inices in ar2"""
  ar3 = trues(size(ar1,1))
  for i in ar2
    ar3[i] = false
  end
  ar1 = ar1[ar3]
  return ar1
end

max_stars_repo_name: jdlaubrie/Kiltro

max_issues_repo_name: jdlaubrie/Kiltro

max_forks_repo_name: jdlaubrie/Kiltro


text:
using Test
using ICD10Utilities
@testset "ICD-10" begin
  @testset "ICD10" begin
    @test_throws DomainError ICD10("A00.000", true)
    @test_throws DomainError ICD10("A00.A", true)
    @test_throws DomainError ICD10("000.00", true)
    @test_throws DomainError ICD10("A00000", true)
    @test_throws DomainError ICD10("A00A", true)
    @test_throws DomainError ICD10("00000", true)

    @test ICD10("A00") < ICD10("A00.01")
    @test ICD10("A00.01") < ICD10("A0002")

    @test ICD10("A00.01") == ICD10("A0001")

    @test string(ICD10("A00.00")) == "A00.00"
    @test string(ICD10("A0000")) == "A00.00"

    ICDOPTS[:punct] = false
    @test string(ICD10("A00.00")) != "A00.00"
    @test string(ICD10("A0000")) == "A0000"
    ICDOPTS[:punct] = true

    @test isvalidcode(ICD10("A00"), ICD10["A00.01", "A00.02", "A00"])
  end

  @testset "Other ICD-10 versions" begin
    @test isdefined(ICD10Utilities, :ICD10AM)
    @test isdefined(ICD10Utilities, :ICD10CA)
    @test isdefined(ICD10Utilities, :ICD10GM)
  end
end

@testset "ICD-10-CM" begin
  @test isdefined(ICD10Utilities, :ICD10CM)
  @test_throws DomainError ICD10CM("A00.00000", true)
  @test string(ICD10CM("C441092")) == "C44.1092"
  ICDOPTS[:punct] = false
  @test string(ICD10CM("C44.1092")) == "C441092"
  ICDOPTS[:punct] = true
  @test string(ICD10(ICD10CM("C441092"))) == "C44.10"
end

@testset "ICD10AM/ACHI" begin
  @testset "ICD10AM" begin
    @test_throws DomainError ICD10("A00.000", true)
    @test_throws DomainError ICD10("A00.A", true)
    @test_throws DomainError ICD10("000.00", true)
    @test_throws DomainError ICD10("A00000", true)
    @test_throws DomainError ICD10("A00A", true)
    @test_throws DomainError ICD10("00000", true)
  end
  @testset "ACHI" begin
    @test_throws DomainError ACHI("12345-678", true)
    @test_throws DomainError ACHI("12345678", true)
    @test_throws DomainError ACHI("A2345-67", true)
    @test_throws DomainError ACHI("A234567", true)
    @test_throws DomainError ACHI("1234-567", true)

    @test ACHI("12345-67") == ACHI("1234567")

    ICDOPTS[:punct] = true
    @test string(ACHI("12345-67")) == "12345-67"
    ICDOPTS[:punct] = false
    @test string(ACHI("12345-67")) == "1234567"
    ICDOPTS[:punct] = true
  end

  @testset "ICD10AM age coding" begin
    @test_throws DomainError ICD10AMAge(7)
    @test_throws DomainError ICD10AMAge("15")
  end

  @testset "ICD10AM/ACHI functions requiring electronic codes lists" begin
    if isdefined(ICD10Utilities, :_ICD10AMcodes_)
      @test isvalidcode(ICD10AM("A00.1"))
      @test !isvalidcode(ICD10AM("A00"))
    else
      @test_throws ErrorException isvalidcode(ICD10AM("A00.1"))
    end
    if isdefined(ICD10Utilities, :_ACHIcodes_)
      @test isvalidcode(ACHI("10801-00"))
      @test !isvalidcode(ACHI("12345-67"))
    else
      @test_throws ErrorException isvalidcode(ACHI("10801-00"))
    end
  end
end


text:
function loglikelihood(dists, p::Path, consts::Vector{Float64})
    like = 0.0
    xx = reshape(p.x, (length(p.x), 1))
    a=xx[1] * p.θ[5]
    b=xx[2] * p.θ[6]
    for i in eachindex(dists)
        like += pdf(dists[i], [a,b])[1] * consts[i]
    end
    log(like)
end

function getCost!(sys::System, pathEnsemble::Ensemble, i::Int64)
    dists = components(sys.data.dd[sys.times[i]])
    n = ncomponents(sys.data.dd[sys.times[i]])
    consts = (2*pi).^(length.(dists)/2).*sqrt.(det.(diagm.(var.(dists))))./n
    for k in eachindex(pathEnsemble)
        costll!(dists, pathEnsemble[k], consts)
    end
end

function costll!(dists, p::Path, consts::Vector{Float64})
    p.dm -= loglikelihood(dists, p, consts)
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
##==============================================================================
## Abstract Types
##==============================================================================

# TODO consider changing this to AbstractFactor
abstract type AbstractFactor end
abstract type AbstractPackedFactor end

const FunctorInferenceType = AbstractFactor       # will eventually deprecate
const PackedInferenceType = AbstractPackedFactor  # will eventually deprecate

abstract type AbstractPrior <: AbstractFactor end
abstract type AbstractRelative <: AbstractFactor end
abstract type AbstractRelativeRoots <: AbstractRelative end
abstract type AbstractRelativeMinimize <: AbstractRelative end

# NOTE DF, Convolution is IIF idea, but DFG should know about "FactorOperationalMemory"
# DF, IIF.CommonConvWrapper <: FactorOperationalMemory #
# NOTE was `<: Function` as unnecessary
abstract type FactorOperationalMemory end   # <: Function
# TODO to be removed from DFG,
# we can add to IIF or have IIF.CommonConvWrapper <: FactorOperationalMemory directly
# abstract type ConvolutionObject <: FactorOperationalMemory end

##==============================================================================
## GenericFunctionNodeData
##==============================================================================

"""
$(TYPEDEF)

Notes
- S::Symbol

Designing (WIP)
- T <: Union{FactorOperationalMemory, PackedInferenceType}
- in IIF.CCW{T <: DFG.FunctorInferenceType}
- in DFG.AbstractRelativeMinimize <: FunctorInferenceType
- in Main.SomeFactor <: AbstractRelativeMinimize
"""
mutable struct GenericFunctionNodeData{T<:Union{<:PackedInferenceType, <:AbstractFactor, <:FactorOperationalMemory}}
    eliminated::Bool
    potentialused::Bool
    edgeIDs::Vector{Int}
    fnc::T
    multihypo::Vector{Float64} # TODO re-evaluate after refactoring w #477
    certainhypo::Vector{Int}
    nullhypo::Float64
    solveInProgress::Int
    inflation::Float64
    # inner constructor needed for dual use
    # GenericFunctionNodeData{T}(x1::Bool,x2::Bool,x3,x4::T,args...) where T = new{T}(x1,x2,x3,x4,args...)
    # TODO deprecate all these inner constructors at end of DFG v0.9.x (was added for GFND.nullhypo::Float64 breaking change)
    # GenericFunctionNodeData{T}(el,po,ed,fn,mu::Vector{<:Real},ce::Vector{Int},so::Int) where T =
    #                            new{T}(el,po,ed,fn,mu,ce,0.0,so)
    # GenericFunctionNodeData{T}(el,po,ed,fn,mu::Vector{<:Real},ce::Vector{Int},nu::Real,so::Int,infl::Real) where T = new{T}(el,po,ed,fn,mu,ce,nu,so,infl)
end

## Constructors

GenericFunctionNodeData{T}() where T =
    GenericFunctionNodeData(false, false, Int[], T())

function GenericFunctionNodeData(   eliminated::Bool,
                                    potentialused::Bool,
                                    edgeIDs::Vector{Int},
                                    fnc::T,
                                    multihypo::Vector{<:Real}=Float64[],
                                    certainhypo::Vector{Int}=Int[],
                                    nullhypo::Real=0,
                                    solveInP::Int=0,
                                    inflation::Real=0.0 ) where T
    return GenericFunctionNodeData{T}(eliminated, potentialused, edgeIDs, fnc, multihypo, certainhypo, nullhypo, solveInP, inflation)
end

##------------------------------------------------------------------------------
## PackedFunctionNodeData and FunctionNodeData

const PackedFunctionNodeData{T} = GenericFunctionNodeData{T} where T <: AbstractPackedFactor
PackedFunctionNodeData(args...) = PackedFunctionNodeData{typeof(args[4])}(args...)

const FunctionNodeData{T} = GenericFunctionNodeData{T} where T <: Union{<:AbstractFactor, <:FactorOperationalMemory}
FunctionNodeData(args...) = FunctionNodeData{typeof(args[4])}(args...)

# PackedFunctionNodeData(x2, x3, x4, x6::T, multihypo::Vector{Float64}=[], certainhypo::Vector{Int}=Int[], x9::Int=0) where T <: PackedInferenceType =
#     GenericFunctionNodeData{T}(x2, x3, x4, x6, multihypo, certainhypo, x9)
# FunctionNodeData(x2, x3, x4, x6::T, multihypo::Vector{Float64}=[], certainhypo::Vector{Int}=Int[], x9::Int=0) where T <: Union{FunctorInferenceType, FactorOperationalMemory} =
#     GenericFunctionNodeData{T}(x2, x3, x4, x6, multihypo, certainhypo, x9)

##==============================================================================
## Factors
##==============================================================================
#
# |                   | label | tags | timestamp | solvable | solverData |
# |-------------------|:-----:|:----:|:---------:|:--------:|:----------:|
# | SkeletonDFGFactor |   X   |   x  |           |          |            |
# | DFGFactorSummary  |   X   |   X  |     X     |          |            |
# | DFGFactor         |   X   |   X  |     X     |     X    |      X     |

## DFGFactor lv2

"""
$(TYPEDEF)
Complete factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGFactor{T, N} <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::ZonedDateTime
    """Nano second time, for more resolution on timestamp (only subsecond information)"""
    nstime::Nanosecond
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Solver data.
    Accessors: [`getSolverData`](@ref), [`setSolverData!`](@ref)"""
    solverData::Base.RefValue{GenericFunctionNodeData{T}}
    """Solvable flag for the factor.
    Accessors: [`getSolvable`](@ref), [`setSolvable!`](@ref)"""
    solvable::Base.RefValue{Int}
    """Internal cache of the ordering of the neighbor variables. Rather use getVariableOrder to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::NTuple{N,Symbol}
    # Inner constructor
    function DFGFactor{T}(label::Symbol,
                 timestamp::Union{DateTime,ZonedDateTime},
                 nstime::Nanosecond,
                 tags::Set{Symbol},
                 solverData::GenericFunctionNodeData{T},
                 solvable::Int,
                 _variableOrderSymbols::NTuple{N,Symbol}) where {T,N}

        #TODO Deprecate remove in v0.10.
        if timestamp isa DateTime
            Base.depwarn("DFGFactor timestamp field is now a ZonedTimestamp", :DFGFactor)
            return new{T,N}(label, ZonedDateTime(timestamp, localzone()), nstime, tags, Ref(solverData), Ref(solvable), _variableOrderSymbols)
        end
        return new{T,N}(label, timestamp, nstime, tags, Ref(solverData), Ref(solvable), _variableOrderSymbols)
    end

end

##------------------------------------------------------------------------------
## Constructors

"""
$(SIGNATURES)

Construct a DFG factor given a label.
"""
DFGFactor(label::Symbol,
          timestamp::Union{DateTime,ZonedDateTime},
          nstime::Nanosecond,
          tags::Set{Symbol},
          solverData::GenericFunctionNodeData{T},
          solvable::Int,
          _variableOrderSymbols::Tuple) where {T} =
          DFGFactor{T}(label, timestamp, nstime, tags, solverData, solvable, _variableOrderSymbols)

DFGFactor{T}(label::Symbol, variableOrderSymbols::Vector{Symbol}, timestamp::Union{DateTime,ZonedDateTime}=now(localzone()), data::GenericFunctionNodeData{T} = GenericFunctionNodeData{T}()) where {T} =
                DFGFactor(label, timestamp, Nanosecond(0), Set{Symbol}(), data, 1, Tuple(variableOrderSymbols))

DFGFactor(label::Symbol,
          variableOrderSymbols::Vector{Symbol},
          data::GenericFunctionNodeData{T};
          tags::Set{Symbol}=Set{Symbol}(),
          timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
          solvable::Int=1,
          nstime::Nanosecond = Nanosecond(0)) where {T} =
                DFGFactor{T}(label, timestamp, nstime, tags, data, solvable, Tuple(variableOrderSymbols))

Base.getproperty(x::DFGFactor,f::Symbol) = begin
    if f == :solvable || f == :solverData
        getfield(x, f)[]
    elseif f == :_variableOrderSymbols
        [getfield(x,f)...]
    else
        getfield(x,f)
    end
end

Base.setproperty!(x::DFGFactor,f::Symbol, val) = begin
    if f == :solvable || f == :solverData
        getfield(x,f)[] = val
    else
        setfield!(x,f,val)
    end
end
##------------------------------------------------------------------------------
## DFGFactorSummary lv1
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Read-only summary factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGFactorSummary <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref)"""
    timestamp::ZonedDateTime
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Internal cache of the ordering of the neighbor variables. Rather use getNeighbors to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::Vector{Symbol}
end

##------------------------------------------------------------------------------
## SkeletonDFGFactor lv0
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Skeleton factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct SkeletonDFGFactor <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Internal cache of the ordering of the neighbor variables. Rather use getNeighbors to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::Vector{Symbol}
end

##------------------------------------------------------------------------------
## Constructors

#NOTE I feel like a want to force a variableOrderSymbols
SkeletonDFGFactor(label::Symbol, variableOrderSymbols::Vector{Symbol} = Symbol[]) = SkeletonDFGFactor(label, Set{Symbol}(), variableOrderSymbols)

##==============================================================================
## Define factor levels
##==============================================================================
const FactorDataLevel0 = Union{DFGFactor, DFGFactorSummary, SkeletonDFGFactor}
const FactorDataLevel1 = Union{DFGFactor, DFGFactorSummary}
const FactorDataLevel2 = Union{DFGFactor}

##==============================================================================
## Convertion constructors
##==============================================================================

DFGFactorSummary(f::DFGFactor) =
    DFGFactorSummary(f.label, f.timestamp, deepcopy(f.tags), deepcopy(f._variableOrderSymbols))

SkeletonDFGFactor(f::FactorDataLevel1) =
    SkeletonDFGFactor(f.label, deepcopy(f.tags), deepcopy(f._variableOrderSymbols))


text:
using WSMP
using Base.Test

const ia = Int32[1,5,9,13,17,21,25,27,29,30];
const ja = Int32[1,3,7,8,2,3,8,9,3,7,8,9,4,6,7,8,5,6,8,9,6,7,8,9,7,8,8,9,9];
const avals = [14.,-1.,-1.,-3.,14.,-1.,-3.,-1.,16.,-2.,-4.,-2.,14.,-1.,-1.,
               -3.,14.,-1.,-3.,-1.,16.,-2.,-4.,-2.,16.,-4.,71.,-4.,16.];
w = Wssmp(ia,ja,avals)
x = w\ones(9)
@test_approx_eq x [0.10094824180165943,0.10094824180165943,0.11951797708415644,
                   0.10094824180165943,0.10094824180165943,0.11951797708415644,
                   0.11951797708415646,0.058079810351639664,0.11951797708415648]
V,R,M = wsversion()
@test V == 14


text:
"""

Create an array of the difference between the event times.

# Notes

* The event times must be unique.

"""
function event_difference_list(eventTimes::Array{<:Number, 1})
    map(x -> event_diff(x, eventTimes), eventTimes)
end

function event_diff(event::Number, eventTimes::Array{<:Number, 1})
    subEvents = eventTimes[eventTimes .< event]
    event .- subEvents
end


text:
"""
    get_data(dir::String; mouse_tag = nothing)

Read directory `dir` and make a list of all the csv files present.
Additionally if a `mouse_tag` is provided filter out all the files where the tag doesn't occur in the filename
"""
function get_data(dir; mouse_tag = nothing)
    f = readdir(dir)
    filter!(x -> occursin(".csv", x), f)
    if mouse_tag != nothing
        filter!(x -> occursin(mouse_tag,x), f)
    end
    return f
end

"""
    get_DataIndex(dir::String; mouse_tag = nothing)

Read and filters files in directory `dir` using `get_data` and extract informations from the filename in a DataFrame.
It returns a DataFrame with File path, Session, MouseID, Day, Turn.
"""
function get_DataIndex(dir; mouse_tag = nothing)
    f = get_data(dir;mouse_tag = mouse_tag)
    DataIndex = DataFrame(Path = joinpath.(dir,f), Session = replace.(f, ".csv"=>""))
    DataIndex.MouseID = [match(r"[a-zA-Z]{1,2}\d{1,2}",file).match for file in f]
    DataIndex.Day = Date.(["20" * match(r"\d{6}",file).match for file in f], "yyyymmdd")
    DataIndex.Turn = [file[end] for file in DataIndex.Session]
    return DataIndex
end


text:
using FixedSizeArrays
include(Pkg.dir("FixedSizeArrays/test/runtests.jl"))


text:
"""
$(TYPEDEF)

Defines a fluctuator ensemble controller

# Fields

$(FIELDS)
"""
mutable struct FluctuatorLiouvillian <: AbstractLiouvillian
    "system-bath coupling operator"
    coupling::Any
    "waitting time distribution for every fluctuators"
    dist::Any
    "cache for each fluctuator value"
    b0::Any
    "index of the fluctuator to be flipped next"
    next_idx::Any
    "time interval for next flip event"
    next_τ::Any
    "noise value"
    n::Any
end

function (F::FluctuatorLiouvillian)(du, u, p, t)
    s = p(t)
    H = sum(F.n .* F.coupling(s))
    du .= -1.0im * (H*u - u*H)
    # gemm! does not work for static matrix
    #gemm!('N', 'N', -1.0im, H, u, 1.0 + 0.0im, du)
    #gemm!('N', 'N', 1.0im, u, H, 1.0 + 0.0im, du)
end

function update_cache!(cache, F::FluctuatorLiouvillian, p, t)
    s = p(t)
    cache .+= -1.0im * sum(F.n .* F.coupling(s))
end

function update_vectorized_cache!(cache, F::FluctuatorLiouvillian, p, t)
    s = p(t)
    hmat = sum(F.n .* F.coupling(s))
    iden = one(hmat)
    cache .+= 1.0im * (transpose(hmat) ⊗ iden - iden ⊗ hmat)
end

function next_state!(F::FluctuatorLiouvillian)
    next_τ, next_idx = findmin(rand(F.dist, size(F.b0, 2)))
    F.next_τ = next_τ
    F.next_idx = next_idx
    F.b0[next_idx] *= -1
    nothing
end

function reset!(F::FluctuatorLiouvillian, initializer)
    F.b0 = abs.(F.b0) .* initializer(length(F.dist), size(F.b0, 2))
    F.n = sum(F.b0, dims = 1)[:]
    next_state!(F)
end