text:
using Cat

# "Real" numbers
const R = Float64

# The only objects in Smooth are Nothing, R, and products of these
@category Smooth #where objects = obj -> obj in (Nothing, R, Tuple)

# Basic algebra
@arrow Smooth Plus :: (R, R) --> R
@arrow Smooth Mult :: (R, R) --> R
@arrow Smooth Neg :: R --> R
# Calc 1
@arrow Smooth Exp :: R --> R
@arrow Smooth Sin :: R --> R
@arrow Smooth Cos :: R --> R

# Just an unbound input is equivalent to the identity morphism
Placeholder = Smooth.Identity{R}

@alias Smooth Plus Base.:+ 2
@alias Smooth Neg Base.:- 1
@alias Smooth Mult Base.:* 2

Base.exp(a::Smooth.Arrow) = Exp(a)
Base.sin(a::Smooth.Arrow) = Sin(a)
Base.cos(a::Smooth.Arrow) = Cos(a)

# Evaluate an expression tree
@interpretation Eval (=>) Smooth

@interpret function (e::Eval)(m::Plus, a, b)
    a + b
end

@interpret function (e::Eval)(m::Mult, a, b)
    a * b
end

@interpret function (e::Eval)(m::Neg, a)
    -a
end

@interpret function (e::Eval)(m::Exp, a)
    exp(a)
end

@interpret function (e::Eval)(m::Sin, a)
    sin(a)
end

@interpret function (e::Eval)(m::Cos, a)
    cos(a)
end

# Functor is a map from arrow -> arrow
# @functor defines a function Diff for all of the
# primitive arrows
@functor Diff :: Smooth => Smooth (T -> Tuple{T, T})

function split_inputs(m::Smooth.Arrow)
    A = Diff(source(m))
    p = A.parameters
    tuple([Smooth.Proj{A, p[k]}(k) for k in 1:length(p)]...)
end

function Diff(m::Plus)
    a, b = split_inputs(m)
    Smooth.Product(a[1] + b[1],
            a[2] + b[2])
end

function Diff(m::Mult)
    a, b = split_inputs(m)
    Smooth.Product(a[1]*b[1],
            a[1]*b[2] + a[2]*b[1])
end

function Diff(m::Smooth.Constant)
    Smooth.Product(Smooth.Constant(m.val),
        Smooth.Constant(0.0))
end

function Diff(m::Exp)
    a, da = split_inputs(m)
    Smooth.Product(exp(a), da*exp(a))
end

function Diff(m::Sin)
    a, da = split_inputs(m)
    Smooth.Product(sin(a), da*cos(a))
end

function Diff(m::Cos)
    a, da = split_inputs(m)
    Smooth.Product(cos(a), -da*sin(a))
end

# Here is the actual "user code"
x = Placeholder()
y = sin(x) + 0.5*cos(x)
y = exp(y*y) + 2.0
dy = Diff(y)

println(Eval()(y, 0.0))
println(Eval()(dy, 0.0, 1.0))


text:
module TestConvenience

include("preamble.jl")

@testset "Convenience" begin
    @test_nothrow IPython.envinfo(devnull)
    @test IPython.pyversion("IPython") isa String
    if occursin(".dev", IPython._pyversion("IPython"))
        @test startswith(IPython.pyversion("IPython"),
                         IPython._pyversion("IPython"))
    else
        @test IPython.pyversion("IPython") == IPython._pyversion("IPython")
    end
    @test IPython.pyversion("__NON_EXISTING__") isa Nothing

    println("vvv DRY RUN vvv")
    @test_nothrow IPython.install_dependency("ipython"; dry_run=true)
    @test_nothrow IPython.install_dependency("spam"; dry_run=true)
    println("^^^ DRY RUN ^^^")

    @test IPython.yes_or_no(input=IOBuffer("yes\n"), output=devnull)
    @test IPython.yes_or_no(input=IOBuffer("no\n"), output=devnull) == false
    @test IPython.yes_or_no(input=IOBuffer("spam\n"), output=devnull) == false

    @test IPython.condajl_installation("IPython") isa Tuple
    @test IPython.conda_installation("IPython") isa Tuple
    @test IPython.pip_installation("IPython") isa Tuple
end

end  # module

ext: jl


text:
@testset "$TEST $G" begin

if !@isdefined(readcentrality)
    function readcentrality(f::AbstractString)
        f = open(f,"r")
        c = Vector{Float64}()
        while !eof(f)
            line = chomp(readline(f))
            push!(c, parse(Float64,line))
        end
        return c
    end
end

g = readgraph(joinpath(testdir,"testdata","graph-50-500.gml"),:gml)

c = readcentrality(joinpath(testdir,"testdata","graph-50-500-bc.txt"))
z = betweenness_centrality(g)

@test map(Float32, z) == map(Float32, c)

y = betweenness_centrality(g, endpoints=true, normalize=false)
@test round.(y[1:3], digits=4) ==
    round.([122.10760591498584, 159.0072453120582, 176.39547945994505], digits=4)
x = betweenness_centrality(g; approx=3)
@test length(x) == 50

# self loops
s2 = DG(3)
add_edge!(s2,1,2); add_edge!(s2,2,3); add_edge!(s2,3,3)
s1 = graph(s2)

@test betweenness_centrality(s1) == [0, 1, 0]
@test betweenness_centrality(s2) == [0, 0.5, 0]

g = G(2)
add_edge!(g,1,2)
z = betweenness_centrality(g; normalize=true)
@test z[1] == z[2] == 0.0

g3 = PathGraph(5, G)
z = betweenness_centrality(g3; normalize=false)
@test z[1] == z[5] == 0.0

end # testset


text:
export KLDivergence, kl_divergence
import Distances: KLDivergence, kl_divergence

const TuMvNormal = Union{DistributionsAD.TuringDenseMvNormal,
                         DistributionsAD.TuringDiagMvNormal,
                         DistributionsAD.TuringScalMvNormal}

function _kld_gaussian(μ1::AbstractArray, σ1::AbstractArray, μ2::AbstractArray, σ2::AbstractArray)
    k  = size(μ1, 1)
    m1 = sum(σ1 ./ σ2, dims=1)
    m2 = sum((μ2 .- μ1).^2 ./ σ2, dims=1)
    m3 = sum(log.(σ2 ./ σ1), dims=1)
    (m1 .+ m2 .+ m3 .- k) ./ 2
end

function _kld_gaussian(p, q)
    (μ1, σ1) = mean(p), var(p)
    (μ2, σ2) = mean(q), var(q)
    _kld_gaussian(μ1, σ1, μ2, σ2)
end

(m::KLDivergence)(p::MvNormal, q::MvNormal) = _kld_gaussian(p,q)
(m::KLDivergence)(p::TuMvNormal, q::TuMvNormal) = _kld_gaussian(p,q)
(m::KLDivergence)(p::ConditionalDists.BMN, q::MvNormal) = _kld_gaussian(p,q)
(m::KLDivergence)(p::ConditionalDists.BMN, q::TuMvNormal) = _kld_gaussian(p,q)

function (m::KLDivergence)(p::ConditionalMvNormal, q::MvNormal, z::AbstractArray)
    _kld_gaussian(condition(p,z), q)
end
function (m::KLDivergence)(p::ConditionalMvNormal, q::TuMvNormal, z::AbstractArray)
    _kld_gaussian(condition(p,z), q)
end

kl_divergence(p::ConditionalMvNormal, q::MvNormal, z::AbstractArray) = KLDivergence()(p,q,z)
kl_divergence(p::ConditionalMvNormal, q::TuMvNormal, z::AbstractArray) = KLDivergence()(p,q,z)


text:
using Pkg
Pkg.add("Zygote")
Pkg.add("IterativeSolvers")
Pkg.add("Arpack")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export WWn

"""
    WWn(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the water waves expansion proposed by Dommermuth and Yue (doi:10.1017/s002211208700288x), West et al. (doi:10.1029/jc092ic11p11803), Craig and Sulem (doi:10.1006/jcph.1993.1164)
(see also the account by Choi : Fifth-order nonlinear spectral model for surface gravity waves: From pseudo-spectral to spectral formulations, RIMS Kokyuroku, 2019)
with the "rectification" method proposed by Duchêne and Melinand.

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`

## Optional keyword arguments
- `ν`: shallow/deep water multiplication factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0` (or `μ=Inf`).
- `n :: Int`: the order of the expansion; linear system if `1`, quadratic if `2`, cubic if `3`, quartic if `4` (default and other values yield `2`);
- `δ` and `m`: parameters of the rectifier operator, set as `k->min(1,|δ*k|^m)` or `k->min(1,|δ*k|^m[1]*exp(1-|δ*k|^m[2]))` if `m` is a couple
(by default is `δ=0`, i.e. no regularization and `m=-1`. Notice `m=-Inf` and `δ>0` yields a cut-off filter);
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!` and in particular
1. a function `WWn.f!` to be called in the time-integration solver;
2. a function `WWn.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `WWn.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential.

"""
mutable struct WWn <: AbstractModel

    label   :: String
	f!		:: Function
	f1!		:: Function
	f2!		:: Function
	mapto	:: Function
	mapfro	:: Function
	param	:: NamedTuple
	kwargs	:: NamedTuple

    function WWn(param::NamedTuple;
							ν=nothing,n=2,δ=0,m=-1,ktol=0,dealias=0, verbose=true)

		if n in [1,2,3,4] && verbose
			@info string("Spectral model with nonlinearity of order ", n)
		elseif verbose
			n = 2
			@info string("Spectral model with nonlinearity of order ", n)
		end
		label = string("WW",n)

		μ 	= param.μ
		ϵ 	= param.ϵ
		if ν == nothing
			if μ > 1
				ν = 1/sqrt(μ)
			else
				ν = 1
			end
		end

		n 	= n
		mesh = Mesh(param)

		param = ( ϵ = ϵ, μ = μ, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )
		kwargs = (ν=ν,n=n,δ=δ,m=m,dealias=dealias,ktol=ktol,verbose=verbose)

		x = mesh.x
		k = copy(mesh.k)

		if length(m) == 1
			rectifier = k -> min(1,abs(k)^m)
		else
			rectifier = k -> min(1,abs(k)^(m[1])*exp(1-abs(k)^m[2]))
		end
		Π = rectifier.(δ*k)   # regularizing rectifier
		K = mesh.kmax * (1-dealias/(2+dealias))
		Π⅔ 	= abs.(mesh.k) .<= K # Dealiasing low-pass filter
		if dealias == 0
			if verbose @info "no dealiasing" end
			Π⅔ 	= ones(size(mesh.k))
		elseif verbose
			@info string("dealiasing : spectral scheme for power ", dealias + 1," nonlinearity ")
		end
		if μ == Inf || ν==0
			∂ₓF₀ 	= 1im * sign.(mesh.k)
			G₀ 	= abs.(mesh.k)
			μ = 1
			ν = 1
		else
        	∂ₓF₀ 	= 1im* sign.(mesh.k) .* tanh.(sqrt(μ)*abs.(mesh.k))
			G₀ 	= sqrt(μ)*abs.(mesh.k).*tanh.(sqrt(μ)*abs.(mesh.k))
		end
		∂ₓ	=  1im * sqrt(μ)* mesh.k            # Differentiation
		z = zeros(Complex{Float64}, mesh.N)
		η = copy(z) ; v = copy(z) ; Lphi = copy(z) ; LzLphi = copy(z) ; dxv = copy(z) ;
		fftη = copy(z) ; fftv = copy(z) ; Q = copy(z) ; R = copy(z) ;

		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[ abs.(U).< ktol ].=0
			return U
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real( ifft(U[:,1]) ),real( ifft(U[:,2]) )
		end

		# Evolution equations are ∂t U = f!(U)
		function f!(U)

			fftη .= U[:,1]
		    fftv .= U[:,2]
			Q .= -∂ₓF₀.*fftv
			R .= -fftη*ν

			if n >= 2
				η  .= ifft(Π.*U[:,1])
				v  .= ifft(U[:,2])
				Lphi .= -ifft(Q)
				Q += -ϵ *∂ₓ.*fft(η.*ifft(fftv)) .+ ϵ*G₀.*fft(η.*Lphi)
				R += ϵ/2*Π.*fft(-v.^2 .+ Lphi.^2)
			end
			if n >= 3
				LzLphi .= ifft(-G₀.* fft(η.*Lphi))
				dxv .= ifft(∂ₓ.*fftv)
				Q += ϵ^2*G₀.*fft(η.* LzLphi + 1/2 * η.^2 .* dxv ) .- ϵ^2*∂ₓ.*∂ₓ.*fft( 1/2 * η.^2  .*Lphi)
				R += ϵ^2*Π.*fft(Lphi .* ( LzLphi .+ η.* dxv ) )
			end
			if n >= 4
				Q += ϵ^3 * G₀.*fft(η.*ifft(-G₀.* fft(η.*LzLphi + 1/2 * η.^2 .* dxv ) )
							.+ 1/2 * η.^2 .* ifft(∂ₓ.*∂ₓ.*fft( η  .* Lphi ) )
							.- 1/6 * η.^3 .* ifft(∂ₓ.*∂ₓ.*fft( Lphi ) ) ) .-
					   ϵ^3 * ∂ₓ.*∂ₓ.*fft( 1/2 * η.^2  .*LzLphi .+ 1/3 * η.^3  .* dxv )
				R += ϵ^3 * Π.*fft( Lphi .*  ifft(-G₀.* fft(η.*LzLphi + 1/2 * η.^2 .* dxv ) )
						.+ 1/2* (LzLphi .+ η .* dxv ).^2
						.+ 1/2* η.* Lphi.^2 .* ifft(∂ₓ.*∂ₓ.* fftη)
						.- 1/2* (η.^2).* (ifft(∂ₓ .* fft(Lphi))).^2 ) .+
						1/4*ϵ^3 * ∂ₓ.*∂ₓ.*fft((η .* Lphi).^2)
			end
		   	U[:,1] .= Π⅔.*Q/sqrt(μ)/ν
		   	U[:,2] .= Π⅔.*∂ₓ.*R/sqrt(μ)/ν
			U[abs.(U).< ktol ].=0

		end

		# Evolution equations are ∂t U = (f1!(U) f2!(U))
		function f1!(U1,U2)

			fftη .= U1
		    fftv .= U2
			Q .= -∂ₓF₀.*fftv

			# attention,  G₀=-L dans Choi
			if n >= 2
				η  .= ifft(Π.*U1)
				v  .= ifft(U2)
				Lphi .= -ifft(Q)
				Q += -ϵ *∂ₓ.*fft(η.*ifft(fftv)) .+ ϵ*G₀.*fft(η.*Lphi)
			end
			if n >= 3
				LzLphi .= ifft(-G₀.* fft(η.*Lphi))
				dxv .= ifft(∂ₓ.*fftv)
				Q += ϵ^2*G₀.*fft(η.* LzLphi + 1/2 * η.^2 .* dxv ) .-
						ϵ^2*∂ₓ.*∂ₓ.*fft( 1/2 * η.^2  .*Lphi)
			end
			if n >= 4
				Q += ϵ^3 * G₀.*fft(η.*ifft(-G₀.* fft(η.*LzLphi + 1/2 * η.^2 .* dxv ) )
							.+ 1/2 * η.^2 .* ifft(∂ₓ.*∂ₓ.*fft( η  .* Lphi ) )
							.- 1/6 * η.^3 .* ifft(∂ₓ.*∂ₓ.*fft( Lphi ) ) ) .-
					   ϵ^3 * ∂ₓ.*∂ₓ.*fft( 1/2 * η.^2  .*LzLphi .+ 1/3 * η.^3  .* dxv )
			end
		   	U1 .= Π⅔.*Q/sqrt(μ)/ν
			U1[abs.(U1).< ktol ].=0

		end

		function f2!(U1,U2)

			fftη .= U1
		    fftv .= U2
			Q .= -∂ₓF₀.*fftv
			R .= -fftη*ν

			# attention,  G₀=-L dans Choi
			if n >= 2
				η  .= ifft(Π.*U1)
				v  .= ifft(U2)
				Lphi .= -ifft(Q)
				R += ϵ/2*Π.*fft(-v.^2 .+ Lphi.^2)
			end
			if n >= 3
				LzLphi .= ifft(-G₀.* fft(η.*Lphi))
				dxv .= ifft(∂ₓ.*fftv)
				R += ϵ^2*Π.*fft(Lphi .* ( LzLphi .+ η.* dxv ) )
			end
			if n >= 4
				R += ϵ^3 * Π.*fft( Lphi .*  ifft(-G₀.* fft(η.*LzLphi + 1/2 * η.^2 .* dxv ) )
						.+ 1/2* (LzLphi .+ η .* dxv ).^2
						.+ 1/2* η.* Lphi.^2 .* ifft(∂ₓ.*∂ₓ.* fftη)
						.- 1/2* (η.^2).* (ifft(∂ₓ .* fft(Lphi))).^2 ) .+
						1/4*ϵ^3 * ∂ₓ.*∂ₓ.*fft((η .* Lphi).^2)
			end
		   	U2 .= Π⅔.*∂ₓ.*R/sqrt(μ)/ν
			U2[abs.(U2).< ktol ].=0

		end

        new(label, f!, f1!, f2! , mapto, mapfro, param, kwargs )
    end
end


text:
using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval, SpecifiedTimes
using Oceananigans.TimeSteppers: Clock

@testset "Schedules" begin
    @info "Testing schedules..."

    # Some fake models
    fake_model_at_iter_3 = (; clock=Clock(time=0.0, iteration=3))
    fake_model_at_iter_5 = (; clock=Clock(time=0.0, iteration=5))

    fake_model_at_time_2 = (; clock=Clock(time=2.0, iteration=0))
    fake_model_at_time_3 = (; clock=Clock(time=3.0, iteration=0))
    fake_model_at_time_4 = (; clock=Clock(time=4.0, iteration=0))
    fake_model_at_time_5 = (; clock=Clock(time=5.0, iteration=0))

    # TimeInterval
    ti = TimeInterval(2)
    @test ti.interval == 2.0
    @test ti(fake_model_at_time_2)
    @test !(ti(fake_model_at_time_3))

    # IterationInterval
    ii = IterationInterval(3)

    @test !(ii(fake_model_at_iter_5))
    @test ii(fake_model_at_iter_3)

    # WallTimeInterval
    wti = WallTimeInterval(1e-9)

    @test wti.interval == 1e-9
    @test wti(nothing)

    # SpecifiedTimes
    st = st_list = SpecifiedTimes(2, 5, 6)
    st_vector = SpecifiedTimes([2, 5, 6])
    @test st_list.times == st_vector.times
    @test st.times == [2.0, 5.0, 6.0]

    # Times are sorted
    st = SpecifiedTimes(5, 2, 6)
    @test st.times == [2.0, 5.0, 6.0]

    @test st(fake_model_at_time_2)

    @test !(st(fake_model_at_time_4))
    @test st(fake_model_at_time_5)
end


text:
boundary_conditions(::LandModel) = (1, 2)

function boundary_state!(
    nf,
    bctype,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    args = (state⁺, diff⁺, aux⁺, n̂, state⁻, diff⁻, aux⁻, t)
    soil_boundary_state!(nf, bctype, land, land.soil, land.soil.water, args...)
    soil_boundary_state!(nf, bctype, land, land.soil, land.soil.heat, args...)
end

function boundary_state!(
    nf,
    bctype,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    args = (state⁺, aux⁺, n̂, state⁻, aux⁻, t)
    soil_boundary_state!(nf, bctype, land, land.soil, land.soil.water, args...)
    soil_boundary_state!(nf, bctype, land, land.soil, land.soil.heat, args...)
end

max_forks_repo_name: jm-c/CLIMA


text:
# An MTR is an MTRBasis together with a coefficient vector
#   note this is an MTR for one treatment arm (d = 0 or d = 1)
#   so in practice we are carrying around 2 element tuples of MTR objects
#
# θ[j,k] corresponds to a_{j}(z)b_{k}(u)
@with_kw struct MTR
    basis::MTRBasis
    θ::Matrix{<:Real} # coefficient vector, arranged in a matrix

    function MTR(basis, θ)
        @assert length(basis.a) == size(θ, 1)
        @assert length(basis.b) == size(θ, 2)
        new(basis, θ)
    end
end
export MTR

function evaluate_mtr(mtr::MTR, ev::DataFrame)
    result = fill(NaN, nrow(ev))
    for (i,r) in enumerate(eachrow(ev))
        result[i] = sum([aj(r[:z]) * bk(r[:u]) for
                            aj in mtr.basis.a, bk in mtr.basis.b]
                        .* mtr.θ)
    end
    return result
end
export evaluate_mtr

function evaluate_mtr(mtrs::Tuple{MTR, MTR}, ev::DataFrame)
    result = fill(NaN, nrow(ev))
    for d in 0:1
        result[(ev.d == d)] = evaluate_mtr(mtrs[d+1], ev[(ev.d == d)])
    end
    return result
end


text:
module AbstractAlgorithms

export
    AbstractCTAlgorithm,
    AbstractProjectionAlgorithm,
    AbstractReconstructionAlgorithm,
    AbstractIRadonAlgorithm
export AbstractParams
export AlgorithmInfo
export project_image, reconstruct_image
export radon, iradon
export alg_name, alg_params

using SimpleTraits
using ..Monads
using ..Geometry:
    AbstractGeometry, AbstractParallelBeamGeometry, AbstractFanBeamGeometry
using ..FanBeam: fan2para
using ..Coordinates
using IntervalSets: Interval

@inline function _alg_progress(
    ::Type{T},
    desc::AbstractString,
    n::Integer,
    enabled::Bool,
    dt::Real=0.2
) where T
    T(n; dt, desc, enabled)
end

abstract type AbstractCTAlgorithm end
abstract type AbstractProjectionAlgorithm <: AbstractCTAlgorithm end
abstract type AbstractReconstructionAlgorithm <: AbstractCTAlgorithm end
abstract type AbstractIRadonAlgorithm <: AbstractReconstructionAlgorithm end
abstract type AlgorithmInfo{A<:AbstractCTAlgorithm} end
abstract type AbstractParams end

alg_name(::Type{A}) where {A <: AbstractCTAlgorithm} = nameof(A)
alg_name(a::AbstractCTAlgorithm) = alg_name(typeof(a))
alg_name(ainfo::AlgorithmInfo) = alg_name(ainfo.algorithm)
@inline function alg_params(ainfo::AInfo) where {AInfo <: AlgorithmInfo}
    :params ∈ fieldnames(AInfo) && return Some(ainfo.params)
    return nothing
end

function radon end
function iradon end
function project_image end
function reconstruct_image end

const _alg_fns =
    (:radon, :iradon, :project_image, :reconstruct_image)
const _algfn_map = (;
    :AbstractProjectionAlgorithm => (:radon, :project_image),
    :AbstractReconstructionAlgorithm => (:iradon, :reconstruct_image),
)

for (T, fns) ∈ pairs(_algfn_map)
    for nm ∈ fns
        @eval begin
            @inline $nm(; kwargs...) = x -> $nm(x; kwargs...)
            @inline $nm(alg::$T; kwargs...) = x -> $nm(x, alg; kwargs...)
            @traitfn @inline function $nm(
                x::M, alg::$T; kwargs...
            ) where {M; Monad{M}}
                x ↣ $nm(alg; kwargs...)
            end
            @traitfn @inline function $nm(x::M; kwargs...) where {M; Monad{M}}
                x ↣ $nm(; kwargs...)
            end
            @inline function $nm(
                ainfo::AlgorithmInfo{A}; kwargs...
            ) where {A<:$T}
                x -> $nm(x, ainfo; kwargs...)
            end
            @inline function $nm(
                x::AbstractMatrix,
                ainfo::AlgorithmInfo{A};
                kwargs...
            ) where {A<:$T}
                g = ainfo.geometry
                p = alg_params(ainfo)
                a = ainfo.algorithm
                isnothing(p) && return $nm(x, g, a; kwargs...)
                $nm(x, g, something(p), a; kwargs...)
            end
            @traitfn @inline function $nm(
                x::M,
                ainfo::AlgorithmInfo{A};
                kwargs...
            ) where {M,A<:$T; Monad{M}}
                x ↣ $nm(ainfo; kwargs...)
            end
        end
    end
end

macro _defapi(f::Symbol, g::Symbol, A::Symbol)
    esc(quote
        @inline function $f(
            image::AbstractMatrix,
            alg::$A;
            kwargs...
        )
            $g(image, alg; kwargs...)
        end
        @inline function $f(
            image::AbstractMatrix,
            geometry::AbstractGeometry,
            alg::$A;
            kwargs...
        )
            $g(image, geometry, alg; kwargs...)
        end
        @inline function $f(
            image::AbstractMatrix,
            geometry::AbstractGeometry,
            params::AbstractParams,
            alg::$A;
            kwargs...
        )
            $g(image, geometry, params, alg; kwargs...)
        end
    end)
end

@_defapi project_image radon AbstractProjectionAlgorithm
@_defapi reconstruct_image iradon AbstractReconstructionAlgorithm

"""
    iradon(sinog::AbstractMatrix[, algorithm::AbstractProjectionAlgorithm]; <keyword arguments>)

Compute the inverse Radon transform of `sinog` with parameters given as keyword
arguments. The parameters depend on the algorithm, please see the relative
documentation. The default algorithm is [`FBP`](@ref).

See Also: [`reconstruct_image`](@ref)
"""
@inline function iradon(
    sinog::AbstractMatrix,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    alg(sinog; kwargs...)
end

"""
    iradon(sinog::AbstractMatrix, xs[, ys[, algorithm::AbstractProjectionAlgorithm[, coo::AbstractCoordinates]]]; <keyword arguments>)

Compute the inverse Radon transform of `sinog` on the points given by the
vectors `xs` and `ys`. If `ys` is omitted, the reconstruction is performed on
the square with `xs == ys`. The default reconstruction algorithm is
[`FBP`](@ref). Please refer to the respective documentation for additional
parameters.

See Also: [`reconstruct_image`](@ref)
"""
@inline function iradon(
    sinog::AbstractMatrix,
    xs::Union{AbstractVector,Interval},
    ys::Union{AbstractVector,Interval},
    alg::AbstractReconstructionAlgorithm,
    coo::AbstractCoordinates = Cartesian();
    kwargs...
)
    alg(sinog, xs, ys, coo; kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    xs::AbstractVector,
    alg::AbstractReconstructionAlgorithm,
    coo::AbstractCoordinates = Cartesian();
    kwargs...
)
    alg(sinog, xs, coo; kwargs...)
end

"""
    iradon(sinog::AbstractMatrix[, geometry::AbstractGeometry[, params::AbstractParams[, algorithm::AbstractProjectionAlgorithm]]]; <keyword arguments>)

Compute the inverse Radon transform of `sinog` with explicit geometry. If the
algorithm needs specific parameters, these can be passed with `params`.
Additional parameters can be passed to the algorithm through keyword arguments.
Please see the respective documentation for more details. The default algorithm
is [`FBP`](@ref).

See Also: [`reconstruct_image`](@ref)
"""
@inline function iradon(
    sinog::AbstractMatrix,
    geometry::AbstractParallelBeamGeometry,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    rows = geometry.rows
    cols = geometry.cols
    α = geometry.α
    α₀ = geometry.α₀
    iradon(sinog, alg; rows, cols, α, α₀, kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    geometry::AbstractFanBeamGeometry,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    g′, para_sinog = fan2para(sinog, geometry)
    iradon(para_sinog, g′, alg; kwargs...)
end

@inline function iradon(
    g::AbstractGeometry,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    x -> iradon(x, g, alg; kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    geometry::AbstractParallelBeamGeometry,
    params::AbstractParams,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    rows = geometry.rows
    cols = geometry.cols
    α = geometry.α
    α₀ = geometry.α₀
    alg(sinog, params; rows, cols, α, α₀, kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    geometry::AbstractFanBeamGeometry,
    params::AbstractParams,
    alg::AbstractReconstructionAlgorithm;
    kwargs...
)
    g′, para_sinog = fan2para(sinog, geometry)
    iradon(sinog, g′, params, alg; kwargs...)
end

end # module


text:
"""
    get_interaction_matrix(::LinearOptics)
"""
function get_interaction_matrix(problem::LinearOptics)
    @debug "start: get_interaction_matrix"

    G = get_empty_matrix(problem.physic, problem.atoms)
    problem.kernelFunction(problem.atoms, problem.laser, G)

    @debug "end  : get_interaction_matrix"
    return G
end

function get_empty_matrix(physic::Scalar, atoms::Atom{<:ThreeD})
    Array{ComplexF64}(undef, atoms.N, atoms.N)
end
function get_empty_matrix(physic::Vectorial, atoms::Atom{<:ThreeD})
    Array{ComplexF64}(undef, 3atoms.N, 3atoms.N)
end

function get_interaction_matrix(problem::NonLinearOptics{MeanField})
    @debug "start: get_interaction_matrix"

    if haskey(problem.data, :G)
        return problem.data[:G]
    end

    temp_scalar_problem = LinearOptics(Scalar(), problem.atoms, problem.laser)
    G = get_empty_matrix(temp_scalar_problem.physic, temp_scalar_problem.atoms)
    temp_scalar_problem.kernelFunction(temp_scalar_problem.atoms, problem.laser, G)

    temp_scalar_problem = 1
    GC.gc()
    problem.data[:G] = G
    @debug "end  : get_interaction_matrix"
    return G
end

"""
    green_scalar!(atoms, laser, G)

Computes:
    @. G = -(Γ/2)*exp(1im*k₀ * R_jk) / (1im * k₀ * R_jk)
    G[diagind(G)] .= 1im * laser.Δ - Γ/2
"""
function green_scalar!(atoms, laser, G)
    @debug "start: green_scalar!"

    G[:] = get_pairwise_matrix(atoms.r) # R_jk

    Threads.@threads for j in eachindex(G) #
        @inbounds G[j] = -(Γ / 2) * cis(k₀ * G[j]) / (1im * k₀ * G[j])
    end
    G[diagind(G)] .= 1im * laser.Δ - Γ / 2

    @debug "end  : green_scalar!"
    return nothing
end

"""
    green_vectorial!(atoms, laser, G)

check bencmark file for details
"""
function green_vectorial!(atoms, laser, G)
    @debug "start: green_vectorial!"

    δ(x, y) = Int(==(x, y))
    N = atoms.N
    Δ = laser.Δ

    Xt, Yt, Zt = atoms.r[1, :], atoms.r[2, :], atoms.r[3, :]

    Xjn = Xt * ones(1, N) - ones(N, 1) * Xt'
    Yjn = Yt * ones(1, N) - ones(N, 1) * Yt'
    Zjn = Zt * ones(1, N) - ones(N, 1) * Zt'

    array_XYZ_jn = [Xjn, Yjn, Zjn]
    Rjn = sqrt.(Xjn .^ 2 + Zjn .^ 2 + Yjn .^ 2)

    P(x) = 1 - 1 / x + 1 / x^2
    Q(x) = -1 + 3 / x - 3 / x^2

    α_range = β_range = [-1, +1, 0]

    term2 = (3 / 2) * exp.(im * k₀ * Rjn) ./ (k₀ * Rjn)
    term2[findall(isnan.(term2))] .= 0
    term2[findall(isinf.(term2))] .= 0

    P_Rjn = P.(im * k₀ * Rjn)
    P_Rjn[findall(isnan.(P_Rjn))] .= 0

    Q_Rjn_over_Rjn_squared = Q.(im * k₀ * Rjn) ./ (k₀ * Rjn) .^ 2
    Q_Rjn_over_Rjn_squared[findall(isnan.(Q_Rjn_over_Rjn_squared))] .= 0

    A = []
    for (α_idx, α) ∈ enumerate(α_range)
        B = []
        for (β_idx, β) ∈ enumerate(β_range)
            term1 = im * I(N) .* δ(α, β)
            # term2 = (3/2)*exp.(im*Rjn)./Rjn  ## Defined outside
            term3 = P_Rjn .* δ(α, β)
            term4 = Q_Rjn_over_Rjn_squared .* (array_XYZ_jn[α_idx] .* array_XYZ_jn[β_idx])

            K = term1 + term2 .* (term3 .+ term4)
            push!(B, K)
        end
        push!(A, vcat(B[1:length(β_range)]...))
    end

    G[:] = (-Γ / 2) * hcat(A[1:length(α_range)]...) #join all matrices
    G .= -im * G # For my code to work, I need this imaginary number. Credits to Sheila
    G[diagind(G)] .= (-Γ / 2 + 1im * Δ)

    @debug "end  : green_vectorial!"
    return nothing
end


text:
__precompile__(false)
module CodeGeneration

import JSON, Calculus
using StaticArrays

using PyCall: pyimport
using JuLIP: save_dict
using ACEtb.SlaterKoster: max_symbol_idx
using ACEtb.SlaterKoster: SKBond, get_l, get_bidx, StandardSigns, sksign, signmod

# load the python part of the code generation
# const _pysys = pyimport("sys")
# push!(_pysys."path", @__DIR__())
# const _codegen = pyimport("codegen")

# TODO: switch this to SymPy.jl ???

sympy = pyimport("sympy")
sp_spin = pyimport("sympy.physics.quantum.spin")

const spsin = sympy.sin
const spcos = sympy.cos
const Rotation = sp_spin.Rotation

_lookupkey(l1, l2, m1, m2, sym) = "$l1,$l2,$m1,$m2,$sym"
_fname_sktable() = joinpath(@__DIR__(), "sktable.json")

"""
`small_d(l, m, M, θ)` : Numerical calculation for Wigner small d-function
"""
function small_d(l, m, mp, θ)
    fc1 = factorial(l+m)
    fc2 = factorial(l-m)
    fc3 = factorial(l+mp)
    fc4 = factorial(l-mp)
    fcm1 = sqrt(fc1 * fc2 * fc3 * fc4)

    cosb = cos(θ / 2.0)
    sinb = sin(θ / 2.0)

    p = m - mp
    lo = max(0,p)
    hi = min(0,l+m,l-mp)

    rtn_sum = 0.0
    for k = lo:hi
       fc5 = factorial(k)
       fc6 = factorial(l+m-k)
       fc7 = factorial(l-mp-k)
       fc8 = factorial(k-p)
       fcm2 = fc5 * fc6 * fc7 * fc8
       pow1 = 2 * l - 2 * k + p
       pow2 = 2 * k - p
       rtn_sum += (-1)^(k+p) * cosb^pow1 * sinb^pow2 / fcm2
    end
    rtn_sum *= fcm1

    return rtn_sum
end

"""
`sk_table(L::Integer)` : read or create a table of SK matrix element expressions
"""
function sk_table(L::Integer; WriteAllow::Bool=true)
   filepath = _fname_sktable()
   try
      tbl = JSON.parsefile(filepath)
      if tbl["L"] >= L
         # return the existing table
         return tbl
      end
      if WriteAllow
         @info("""The existing lookup table does not contain a high enough degree.
                  I'm now going to create a new one which can take a little while.""")
      end
   catch
      if WriteAllow
         @info("""Reading SK lookup table was unsuccesful, I'm now going to create
                  a new one. This could take a while (O(minutes))""")
      end
   end

   # create a new lookup table
   tbl = Dict{String, Any}("L" => L)
   for l1 = 0:L, l2=0:L, m1=-l1:l1, m2=-l2:l2
      for sym = 0:max_symbol_idx(l1,l2)
         tbl[_lookupkey(l1,l2,m1,m2,sym)] = Gsym(l1,l2,m1,m2,sym)
      end
   end
   if WriteAllow
       save_dict(filepath, tbl)
   end

   # return the new table
   return tbl
end

"""
INPUTS: l1: Angular quantum number of atom 1.
        l2: Angular quantum number of atom 2.
        m1: Magnetic quantum number of atom 1.
        m2: Magnetic quantum number of atom 2.
        sym: A float which represents the type of symmetry bond we are considering, e.g. M = 0
           for sigma, M = 1 for pi, M = 2 for delta, etc.
        alpha: Azithmuthal coordinate.
        beta: Polar coordinate.
RETURNS: This function returns the relevant coefficient that should be are used in writing the
         Slater-Koster transformations.
"""
function Gsym(l1, l2, m1, m2, mm)

   φ, θ = sympy.symbols("phi, theta")

   τ(m) = (m >= 0) ? 1 : 0

   A(m, φ) = ( (m == 0) ? sympy.sqrt(2)/2 :
               (-1)^m * (τ(m) * spcos(abs(m)*φ) + τ(-m) * spsin(abs(m)*φ)) )

   B(m, φ) = ( (m == 0) ? 0 :
               (-1)^m * (τ(-m) * spcos(abs(m)*φ) - τ(m) * spsin(abs(m)*φ) ) )

   rot(l, m, sym, θ) = Rotation.d(l, m, sym, θ).doit()

   S(l, m, sym, φ, θ) = A(m, φ) * ( (-1)^sym * rot(l, abs(m),  sym, θ)
                                             + rot(l, abs(m), -sym, θ) )

   T(l, m, sym, φ, θ) = B(m, φ) * ( (-1)^sym * rot(l, abs(m),  sym, θ)
                                             - rot(l, abs(m), -sym, θ) )

   if mm == 0
      expr = (2 * A(m1, φ) * rot(l1, abs(m1), 0, θ)
                * A(m2, φ) * rot(l2, abs(m2), 0, θ) )
   else
      expr = (   S(l1, m1, mm, φ, θ) * S(l2, m2, mm, φ, θ)
               + T(l1, m1, mm, φ, θ) * T(l2, m2, mm, φ, θ) )
   end

   #expr = expr.simplify()

   # x, y, z = sympy.symbols("x y z")
   # expr = expr.simplify()
   # expr = expr.subs(spsin(φ) * spsin(θ), y)
   # expr = expr.subs(spcos(φ) * spsin(θ), x)
   # expr = expr.subs(spcos(θ), z)
   # # expr = expr.subs(spsin(φ)^2, 1.0 - spcos(φ)^2)
   # expr = expr.simplify()

   # convert the sympy expression into a string
   return py2jlcode(expr.__str__())
end

function py2jlcode(str)
   # replace some symbols to turn python into julia code
   str = replace(str, "**" => "^")
   # φ : Azithmuthal coordinate
   str = replace(str, "phi" => "φ")   # the python code uses alpha here
   # θ : Polar angle
   str = replace(str, "theta" => "θ")   # the python code uses beta here
end

@generated function sk_gen(::SKBond{O1, O2, SYM}, φ, θ; W::TF=Val{true}) where {O1, O2, SYM, TF}
   tf = TF == Val{true} ? true : false
   # get the SK expressions table
   l1, l2 = get_l(Val{O1}()), get_l(Val{O2}())
   tbl = sk_table(max(l1,l2), WriteAllow=tf)
   expr_str = "SMatrix{$(2*l1+1),$(2*l2+1)}("
   for m2 = -l2:l2, m1 = -l1:l1
      ex = tbl[_lookupkey(l1, l2, m1, m2, get_bidx(Val{SYM}()))]
      expr_str *= ex * ", "
   end
   code = Meta.parse(expr_str[1:end-2] * ")")
   quote
      $code
   end
end

tau(m) = (m >= 0) ? 1 : 0

AA(m, phi) = ( (m == 0) ? sqrt(2)/2 :
             (-1)^m * (tau(m) * cos(abs(m)*phi) + tau(-m) * sin(abs(m)*phi)) )

BB(m, phi) = ( (m == 0) ? 0 :
             (-1)^m * (tau(-m) * cos(abs(m)*phi) - tau(m) * sin(abs(m)*phi) ) )

function small_d_rot_sym(ll, mm, pm)
   θ = sympy.symbols("theta")
   rot(l, m, sym) = Rotation.d(l, m, sym, θ).doit()
   rot_ex = rot(ll, mm, pm)
   return py2jlcode(rot_ex.__str__())
end

#@generated function small_d_rot(l, m, sym, θ)
function small_d_rot(l, m, sym, θ)
   expr_str = small_d_rot_sym(l, m, sym)
   expr_str = replace(expr_str, "θ" => θ)
   code = Meta.parse(expr_str)
   #eval(quote
   #   $code
   #end)
   eval(code)
end

function d_rot(l, m, sym, θ)
  sd = small_d_rot(l, m, sym, θ)
  return sd
end

SS(l, m, sym, phi, theta) = AA(m, phi) * ( (-1)^sym * d_rot(l, abs(m),  sym, theta)
                                                    + d_rot(l, abs(m), -sym, theta) )

TT(l, m, sym, phi, theta) = BB(m, phi) * ( (-1)^sym * d_rot(l, abs(m),  sym, theta)
                                                    - d_rot(l, abs(m), -sym, theta) )

function Gnum(l1, l2, m1, m2, sym, phi, theta)

   rtn = 0.0

   if sym == 0
      rtn += (2 * AA(m1, phi) * d_rot(l1, abs(m1), 0, theta)
                * AA(m2, phi) * d_rot(l2, abs(m2), 0, theta) )
   else
      rtn += (   SS(l1, m1, sym, phi, theta) * SS(l2, m2, sym, phi, theta)
               + TT(l1, m1, sym, phi, theta) * TT(l2, m2, sym, phi, theta) )
   end

   return rtn
end

function sk_num(::SKBond{O1, O2, SYM}, φ, θ) where {O1, O2, SYM}
   l1, l2 = get_l(Val{O1}()), get_l(Val{O2}())
   E = zeros(2*l1+1,2*l2+1)
   b_l = get_bidx(Val{SYM}())
   for (mj, m2) = enumerate(-l2:l2), (mi, m1) = enumerate(-l1:l1)
      E[mi, mj] = Gnum(l1, l2, m1, m2, b_l, φ, θ)
   end
   return E
end

end

# TODO: create in-place type-stable versions of sk_gen
#       taking l, m as input
# _parse_key(k) = eval(Meta.parse("[" * k * "]"))
#
# # generate code for sk blocks
# # require at least spd
# tbl = sk_table(2)
# lmax = 2
# for (key, expr) in tbl
#    lmax = max(lmax, _parse_key(key)[1])
# end
# _sk_f_tbl =


text:
module LLVMCGUtils

using LLVM

include("misc.jl")

# AddressSpace
const GENERIC = 0
const TRACKED = 10
const DERIVED = 11
const CALLEEROOTED = 12
const LOADED = 13

# jlvalue types
@inline get_jlvalue_ty(ctx::Context) = LLVM.StructType(LLVMType[]; ctx)
@inline get_pjlvalue_ty(ctx::Context) = LLVM.PointerType(get_jlvalue_ty(ctx), GENERIC)
@inline get_prjlvalue_ty(ctx::Context) = LLVM.PointerType(get_jlvalue_ty(ctx), TRACKED)
@inline get_ppjlvalue_ty(ctx::Context) = LLVM.PointerType(get_pjlvalue_ty(ctx), GENERIC)
@inline get_pprjlvalue_ty(ctx::Context) = LLVM.PointerType(get_prjlvalue_ty(ctx), GENERIC)

@inline function get_jlfunc_ty(ctx::Context)
    T_prjlvalue = get_prjlvalue_ty(ctx)  # function
    T_pprjlvalue = LLVM.PointerType(T_prjlvalue, GENERIC) # args[]
    T_nargs = LLVM.Int32Type(ctx) # nargs
    return LLVM.FunctionType(T_prjlvalue, [T_prjlvalue, T_pprjlvalue, T_nargs])
end

@inline function get_jlfuncparams_ty(ctx::Context)
    T_prjlvalue = get_prjlvalue_ty(ctx)  # function
    T_pprjlvalue = LLVM.PointerType(T_prjlvalue, GENERIC) # args[]
    T_nargs = LLVM.Int32Type(ctx) # nargs
    T_sparam_vals = T_pprjlvalue # linfo->sparam_vals
    return LLVM.FunctionType(T_prjlvalue, [T_prjlvalue, T_pprjlvalue, T_nargs, T_sparam_vals])
end

@inline get_voidfunc_ty(ctx::Context) = LLVM.FunctionType(LLVM.VoidType(ctx))
@inline get_pvoidfunc_ty(ctx::Context) = LLVM.PointerType(get_voidfunc_ty(ctx), GENERIC)

# important functions
# TODO...

# placeholder functions
function create_gcroot_flush(mod::LLVM.Module)
    ctx = context(mod)
    ft = LLVM.FunctionType(get_voidfunc_ty(ctx))
    fn = LLVM.Function(mod, "julia.gcroot_flush", ft)
    return fn
end

function create_gc_preserve_begin(mod::LLVM.Module)
    ctx = context(mod)
    ft = LLVM.FunctionType(LLVM.TokenType(ctx); vararg=true)
    fn = LLVM.Function(mod, "llvm.julia.gc_preserve_begin", ft)
    return fn
end

function create_gc_preserve_end(mod::LLVM.Module)
    ctx = context(mod)
    ft = LLVM.FunctionType(get_voidfunc_ty(ctx), [LLVM.TokenType(ctx)]; vararg=true)
    fn = LLVM.Function(mod, "llvm.julia.gc_preserve_end", ft)
    return fn
end

function create_except_enter(mod::LLVM.Module)
    ctx = context(mod)
    ft = LLVM.FunctionType(LLVM.Int32Type(ctx); vararg=true)
    fn = LLVM.Function(mod, "julia.except_enter", ft)
    fn_attrs = function_attributes(fn)
    push!(fn_attrs, StringAttribute("returns_twice", ""; ctx))
    return fn
end

function create_pointer_from_objref(mod::LLVM.Module)
    ctx = context(mod)
    ft = LLVM.FunctionType(get_pjlvalue_ty(ctx), [LLVM.PointerType(get_jlvalue_ty(ctx), DERIVED)])
    fn = LLVM.Function(mod, "julia.pointer_from_objref", ft)
    fn_attrs = function_attributes(fn)
    push!(fn_attrs, StringAttribute("readnone", ""; ctx))
    push!(fn_attrs, StringAttribute("nounwind", ""; ctx))
    ret_attrs = return_attributes(fn)
    push!(ret_attrs, StringAttribute("nonnull", ""; ctx))
    return fn
end

# utils
function track_pjlvalue(builder::Builder, v::Value)
    ctx = context(v)
    T_jlvalue = get_jlvalue_ty(ctx)
    T_pjlvalue = LLVM.PointerType(T_jlvalue, GENERIC)
    @assert llvmtype(v) == T_pjlvalue
    T_prjlvalue = LLVM.PointerType(T_jlvalue, TRACKED)
    return addrspacecast!(builder, v, T_prjlvalue)
end

function maybe_decay_untracked(builder::Builder, v::Value)
    ctx = context(v)
    T_jlvalue = get_jlvalue_ty(ctx)
    T_pjlvalue = LLVM.PointerType(T_jlvalue, GENERIC)
    T_prjlvalue = LLVM.PointerType(T_jlvalue, TRACKED)
    if llvmtype(v) == T_pjlvalue
        return addrspacecast!(builder, v, T_prjlvalue)
    end
    @assert llvmtype(v) == T_prjlvalue
    return v
end

function decay_derived(builder::Builder, v::Value)
    ty = llvmtype(v)
    addrspace(ty) == DERIVED && return v
    new_ty = LLVM.PointerType(eltype(ty), DERIVED) # TODO:[#44310][LLVM 13] use getWithSamePointeeType instead
    return addrspacecast!(builder, v, new_ty)
end

function maybe_decay_tracked(builder::Builder, v::Value)
    ty = llvmtype(v)
    addrspace(ty) != TRACKED && return v
    new_ty = LLVM.PointerType(eltype(ty), DERIVED) # TODO:[#44310][LLVM 13] use getWithSamePointeeType instead
    return addrspacecast!(builder, v, new_ty)
end

function mark_callee_rooted(builder::Builder, v::Value)
    ctx = context(v)
    T_jlvalue = get_jlvalue_ty(ctx)
    T_pjlvalue = LLVM.PointerType(T_jlvalue, GENERIC)
    T_prjlvalue = LLVM.PointerType(T_jlvalue, TRACKED)
    ty = llvmtype(v)
    @assert ty == T_pjlvalue || ty == T_prjlvalue
    T_callee_rooted = LLVM.PointerType(T_jlvalue, CALLEEROOTED)
    return addrspacecast!(builder, v, T_callee_rooted)
end

function emit_pointer_from_objref(mod::LLVM.Module, builder::Builder, v::Value)
    address_space = addrspace(llvmtype(v))
    if address_space != TRACKED && address_space != DERIVED
        return v
    end
    v = decay_derived(builder, v)
    ctx = context(v)
    T_jlvalue = get_jlvalue_ty(ctx)
    ty = LLVM.PointerType(T_jlvalue, DERIVED)
    if llvmtype(v) != ty
        v = bitcast!(builder, v, ty)
    end
    funcs = functions(mod)
    if haskey(funcs, "julia.pointer_from_objref")
        fn = funcs["julia.pointer_from_objref"]
    else
        fn = create_pointer_from_objref(mod)
    end
    callinst = call!(builder, fn, [v])
    # TODO: sync attributes
    return callinst
end

function emit_bitcast(builder::Builder, v::Value, jl_value::Type)
    v_as = addrspace(llvmtype(v))
    if jl_value isa LLVM.PointerType && v_as != addrspace(jl_value)
        jl_value_addr = LLVM.PointerType(jl_value, v_as) # TODO:[#44310][LLVM 13] use getWithSamePointeeType instead
        return bitcast!(builder, v, jl_value_addr)
    else
        return bitcast!(builder, v, jl_value)
    end
end

maybe_bitcast(builder::Builder, v::Value, to::Type) = to != llvmtype(v) ? emit_bitcast(builder, v, to) : v

end


text: my_f(x,y) = 2x + 1


text:
using Test
using ChainPartitioners
using SparseArrays
using Random
using Base.Iterators
include("matrices.jl")

macro ignore(ex)
    nothing
end

@testset "ChainPartitioners" begin
    Random.seed!(0xDEADBEEF)
    include("test_util.jl")
    include("test_SparsePrefixMatrices.jl")
    include("test_EnvelopeMatrices.jl")
    include("test_Partitions.jl")
    include("test_SparseColorArrays.jl")
    include("test_Costs.jl")
    include("test_Partitioners.jl")
    include("test_CuthillMcKee.jl")
end


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, MLDatasets, Plots
using Printf, Debugger, IterativeSolvers, CuArrays
train_x, yy_train_all = MLDatasets.MNIST.traindata()
test_x,  yy_test_all  = MLDatasets.MNIST.testdata()
xx_train=MLDatasets.MNIST.convert2features(MLDatasets.MNIST.traintensor())'
xx_test=MLDatasets.MNIST.convert2features(test_x)'
xx_train=convert(Array{Float64,2},xx_train)
xx_test=convert(Array{Float64,2},xx_test)
yy_train_all=convert(Array{Float64,1},yy_train_all)
yy_test_all=convert(Array{Float64,1},yy_test_all)
i_train=size(xx_train,1)
i_test=size(xx_test,1)
vars=size(xx_train,2)

path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

maetrs=Array{Float64}(undef,0)
maetes=Array{Float64}(undef,0)
dts=Array{Float64}(undef,0)
predl1_train_all=Array{Float64}(undef,i_train,0)
predl1_test_all=Array{Float64}(undef,i_test,0)
@time for num_to_test=0:9
    global yy_train_all,yy_test_all,predl1_train_all,predl1_test_all
    yy_train=copy(yy_train_all)
    yy_test=copy(yy_test_all)
    tol1=0.02
    yy_train=copy(yy_train_all)
    i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
    mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
    yy_test=copy(yy_test_all)
    i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
    yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1
    inds_all=1:i_train;neurons=1000;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
    t1=now()
    a_all,a_layer1,layer1,mat1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
    yy_train=cat(yy_train,dims=[1,2]);yy_test=cat(yy_test,dims=[1,2])
    nodes=neurons.*ones(Int64,1,10)
    nodes=cat(nodes,dims=[1,2])
    n_weigths = vars*nodes[1]
    for i=2:length(nodes) n_weigths+=nodes[i]*nodes[i-1] end; n_weigths+=sum(nodes)
    for i=1:size(yy_train,2) n_weigths+=1+nodes[end] end; opti_ww=rand(n_weigths)
    opti_ww,layer_prev=ANNBN.initialize_deep_weights_per_leyer(neurons,vars,a_all,nodes,a_layer1,n_weigths,xx_train,yy_train,i_train,layer1,mat1)
    t2=now();dt=convert(Float64,Dates.value(t2-t1))/1000;push!(dts,dt)
    predl1=ANNBN.deep_nnm(opti_ww,xx_train,nodes,i_train,vars,size(yy_train,2))
    predl1_test=ANNBN.deep_nnm(opti_ww,xx_test,nodes,i_test,vars,size(yy_test,2))
    maetr=100mean(abs.(round.(predl1).-round.(yy_train)).<0.1)
    maete=100mean(abs.(round.(predl1_test).-round.(yy_test)).<0.1)

    predl1_train_all=[predl1_train_all predl1]
    predl1_test_all=[predl1_test_all predl1_test]

    push!(maetrs,maetr)
    push!(maetes,maete)
    println(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>",num_to_test," ",maete)
end
println(dts)
println(maetes)

yy_train_pred=zeros(i_train)
for i=1:i_train
    yy_train_pred[i]=sortperm(predl1_train_all[i,:],rev=true)[1]-1
end
CC=100-100mean(yy_train_all-yy_train_pred.>0.01)

yy_test_pred=zeros(i_test)
for i=1:i_test
    yy_test_pred[i]=sortperm(predl1_test_all[i,:],rev=true)[1]-1
end
CC=100-100mean(yy_test_all-yy_test_pred.>0.01)


text:
# """
#    simulate(x0::Vector{Float64}, p::Parameters; t_start=0.0, t_end=1.0)

# Function to do simple simulation given starting abundances `x0` and a time range `t_start` to `t_end`.
# """
# function simulate(x0::Vector{Float64}, p::Parameters; t_start=0.0, t_end=1.0)
#     prob = OrdinaryDiffEq.ODEProblem(dx!,x0,(t_start,t_end),deepcopy(p))
#     sol = DiffEqBase.solve(prob, OrdinaryDiffEq.AutoTsit5(OrdinaryDiffEq.Rosenbrock23()) )
#     return(sol)
# end

# """
# simulate(x0::Vector{Float64}, p::Parameters, assembly::Any; t_start=0.0, t_end=1.0)

# Function to do a simulation with assembly rules as dictated by the callback function `assembly`. The assembly parameter must be generated by either the `add_at_equilibrium` or `add_at_t` functions.
# """
# function simulate(x0::Vector{Float64}, p::Parameters, assembly_rules::Any; t_start=0.0, t_end=1.0)
#     #for equilibrium assembly
#     if assembly_rules[1] == :equilibrium
#         prob = OrdinaryDiffEq.ODEProblem(dx!,x0,(t_start,t_end),deepcopy(p))
#         sol = DiffEqBase.solve(prob,
#             OrdinaryDiffEq.AutoTsit5(OrdinaryDiffEq.Rosenbrock23()),
#             callback = assembly_rules[2])

#     #for given t_stops assembly
#     elseif assembly_rules[1] == :at_time_stop
#         prob = DiffEqCallbacks.ODEProblem(dx!,x0,(t_start,t_end),deepcopy(p))
#         sol = DiffEqBase.solve(prob,
#             OrdinaryDiffEq.AutoTsit5(OrdinaryDiffEq.Rosenbrock23()),
#             callback = assembly_rules[2],
#             tstops = assembly_rules[3])
#     else
#         error("Incorrect assembly process specfified. Please see ?simulate")
#     end

#     return(sol)
# end


text:
const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# abstract type AsbtractBilevelOptimizer end
# struct SOS1Optimizer{O} <: AsbtractBilevelOptimizer
#     solver::O
#     # options
# end
# function SOS1Optimizer(solver::O) where O
#     return SOS1Optimizer{O}(solver)
# end

struct Complement#{M1 <: MOI.ModelLike, M2 <: MOI.ModelLike, F, S}
    # primal::M1
    func_w_cte#::F
    set_w_zero#::S
    # dual::M2
    variable#::VI
end

abstract type BilevelSolverMode{T} end

mutable struct SOS1Mode{T} <: BilevelSolverMode{T}
    epsilon::T
    function SOS1Mode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ComplementMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ComplementMode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ComplementWithSlackMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ComplementWithSlackMode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ProductMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ProductMode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ProductWithSlackMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ProductWithSlackMode()
        return new{Float64}(zero(Float64))
    end
end

function get_canonical_complements(primal_model, primal_dual_map)
    map = primal_dual_map.primal_con_dual_var
    out = Complement[]
    T = Float64
    for ci in keys(map)
        func = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
        set = MOI.get(primal_model, MOI.ConstraintSet(), ci)
        i = 1
        func.constant = Dualization.set_dot(i, set, T)*Dualization.get_scalar_term(primal_model, i, ci)
        # todo - set dot on function
        con = Complement(func, set_with_zero(set), map[ci][1])
        push!(out, con)
    end
    return out
end

function set_with_zero(set::Union{MOI.LessThan{T},
    MOI.GreaterThan{T}, MOI.EqualTo{T}}) where T
    return typeof(set)(0.0)
end
function set_with_zero(set)
    return copy(set)
end

function build_bilevel(
    upper::MOI.ModelLike, lower::MOI.ModelLike,
    link::Dict{VI,VI}, upper_variables::Vector{VI},
    mode)

    # Start with an empty problem
    moi_mode = MOIU.AUTOMATIC
    m = MOIU.CachingOptimizer(MOIU.Model{Float64}(), moi_mode)

    # add the first level
    copy_names = false
    # key are from src, value are from dest
    upper_idxmap = MOIU.default_copy_to(m, upper, copy_names)
    pass_names(m, upper, upper_idxmap)

    # append the second level primal
    lower_idxmap = MOIU.IndexMap() #

    for (upper_key, lower_val) in link
        lower_idxmap[lower_val] = upper_idxmap[upper_key]
    end

    append_to(m, lower, lower_idxmap, copy_names)
    pass_names(m, lower, lower_idxmap)

    # dualize the second level
    dual_problem = dualize(lower,
        dual_names = DualNames("dual_","dual_"),
        variable_parameters = upper_variables,
        ignore_objective = true)
    lower_dual = dual_problem.dual_model
    lower_primal_dual_map = dual_problem.primal_dual_map

    # append the second level dual
    lower_dual_idxmap = MOIU.IndexMap()

    append_to(m, lower_dual, lower_dual_idxmap, copy_names)
    pass_names(m, lower_dual, lower_dual_idxmap)

    # complete KKT
    # 1 - primal dual equality (quadratic equality constraint)
    #     use quadratic or log expansion (requires bounds on variables get from outside or from prob)
    # 1a - no slacks
    # 1b - use slack
    # 2 - complementarity
    # 2a - actual complementarity constraints w and w/o slack (ok)
    # 2b - SOS1 (slack * duals) (ok)
    # 2c - Big-M formulation (Fortuny-Amat and McCarl, 1981)
    # 2d - Product of variables  w and w/o slack (ok)
    # 3 - NLP (y = argmin_y problem)

    # futurely add slacks to primal model
    comps = get_canonical_complements(lower, lower_primal_dual_map)
    for comp in comps
        add_complement(mode, m, comp, lower_idxmap, lower_dual_idxmap)
    end

    return m, upper_idxmap, lower_idxmap#, lower_primal_dual_map, lower_dual_idxmap
end

function add_complement(mode::ComplementMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)

    dual = idxmap_dual[v]

    new_f = MOIU.operate(vcat, T, f_dest, MOI.SingleVariable(dual))

    c1 = MOI.add_constraint(m,
        new_f,
        MOI.Complements(1))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.ConstraintName(), c1, "compl_compl_($(nm))")

    return c1
end

function add_complement(mode::ComplementWithSlackMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    dual = idxmap_dual[v]
    c1 = MOI.add_constraint(m,
        MOI.VectorOfVariables([slack, dual]),
        MOI.Complements(1))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_complWslk_($(nm))")

    return slack, slack_in_set, equality, c1
end

function add_complement(mode::SOS1Mode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    dual = idxmap_dual[v]
    c1 = MOI.add_constraint(m,
        MOI.VectorOfVariables([slack, dual]),
        MOI.SOS1([1.0, 2.0]))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_sos1_($(nm))")

    return slack, slack_in_set, equality, c1
end

function add_complement(mode::ProductMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)

    dual = idxmap_dual[v]

    new_f = MOIU.operate(*, T, f_dest, MOI.SingleVariable(dual))

    c1 = MOI.add_constraint(m,
        new_f,
        MOI.EqualTo(zero(T)))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.ConstraintName(), c1, "compl_prod_($(nm))")

    return c1
end

function add_complement(mode::ProductWithSlackMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    dual = idxmap_dual[v]

    prod_f = MOIU.operate(*, T, MOI.SingleVariable(slack), MOI.SingleVariable(dual))

    c1 = MOI.add_constraint(m,
        prod_f,
        MOI.EqualTo(zero(T)))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_prodWslk_($(nm))")

    return slack, slack_in_set, equality, c1
end

function pass_names(dest, src, map)
    for vi in MOI.get(src, MOI.ListOfVariableIndices())
        name = MOI.get(src, MOI.VariableName(), vi)
        MOI.set(dest, MOI.VariableName(), map[vi], name)
    end
    for (F,S) in MOI.get(src, MOI.ListOfConstraints())
        for con in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
            name = MOI.get(src, MOI.ConstraintName(), con)
            MOI.set(dest, MOI.ConstraintName(), map[con], name)
        end
    end
end

function append_to(dest::MOI.ModelLike, src::MOI.ModelLike, idxmap, copy_names::Bool)
    # MOI.empty!(dest)

    # idxmap = MOIU.IndexMap()

    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    constraint_types = MOI.get(src, MOI.ListOfConstraints())
    single_variable_types = [S for (F, S) in constraint_types
                             if F == MOI.SingleVariable]
    vector_of_variables_types = [S for (F, S) in constraint_types
                                 if F == MOI.VectorOfVariables]

    # The `NLPBlock` assumes that the order of variables does not change (#849)
    if MOI.NLPBlock() in MOI.get(src, MOI.ListOfModelAttributesSet())
        error("not nlp for now")
        vector_of_variables_not_added = [
            MOI.get(src, MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}())
            for S in vector_of_variables_types
        ]
        single_variable_not_added = [
            MOI.get(src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
            for S in single_variable_types
        ]
    else
        vector_of_variables_not_added = [
            MOIU.copy_vector_of_variables(dest, src, idxmap, S)
            for S in vector_of_variables_types
        ]
        single_variable_not_added = [
            MOIU.copy_single_variable(dest, src, idxmap, S)
            for S in single_variable_types
        ]
    end

    MOIU.copy_free_variables(dest, idxmap, vis_src, MOI.add_variables)

    # Copy variable attributes
    # MOIU.pass_attributes(dest, src, copy_names, idxmap, vis_src)

    # Copy model attributes
    # pass_attributes(dest, src, copy_names, idxmap)

    # Copy constraints
    MOIU.pass_constraints(dest, src, copy_names, idxmap,
                     single_variable_types, single_variable_not_added,
                     vector_of_variables_types, vector_of_variables_not_added)

    return idxmap
end


text:
using PyPlot

phi = linspace(-pi,pi,1000)
dphi = collect(linspace(-2pi,2pi,100))

spot = cos(phi)
spot = 0.5.*(spot+abs(spot))

amp = zeros(100)
for i=1:100
  sspot = cos(phi+dphi[i])
  sspot = 0.5.*(sspot+abs(sspot))
  plot(phi,spot)
  plot(phi,sspot)
  plot(phi,spot.*sspot)
  amp[i] = sum(spot.*sspot)
end

clf()
plot(dphi,amp)
gamma = 1.7
plot(dphi,(exp(-gamma.*sin(.5.*dphi).^2)-exp(-gamma))./(1.-exp(-gamma)).*maximum(amp))

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
export liftFactor

function liftFactorVars(A,B,D)
    Base.depwarn("liftFactor and liftFactorVars are no longer necessary and will be removed in a future release. Use `O\\x` instead.", :liftFactorVars)
    A  = sparse(A)
    B  = sparse(B)
    Di = sparse(inv(D))
    n  = size(A,1)
    k  = size(B,2)
    M = [A    B   ;
         B'  -Di ];
    M = lu(M) # ldltfact once it's available.
    return x -> (M\[x; zeros(k,1)])[1:n,:]
end

function liftFactorVars(A,B,D::SparseMatrixCSC)
  liftFactorVars(A,B,Matrix(D))
end

# This could be optimized to avoid the extra allocation of a single element matrix
function liftFactorVars(A,B,D::Real)
  liftFactorVars(A,B,ones(1,1)*D)
end

"""
    liftFactor(A)

More stable version of inv(A).  Returns a function which computes the inverse
on evaluation, i.e. `liftFactor(A)(x)` is the same as `inv(A)*x`.
"""
liftFactor(O::SymWoodbury) = liftFactorVars(O.A,O.B,O.D)


text:
n = rand(2:10)

kappa = rand(2:6)
mode = rand(1:5)

A = matrixdepot("randsvd", n, kappa, mode)
B = matrixdepot("randsvd", n, kappa)

if mode == 5
    @test cond(A) <= kappa
else
    @test_approx_eq cond(A) kappa
end

A5 = matrixdepot("randsvd", n, kappa, 5)
A4 = matrixdepot("randsvd", n, kappa, 4)
A3 = matrixdepot("randsvd", n, kappa, 3)
A2 = matrixdepot("randsvd", n, kappa, 2)
A1 = matrixdepot("randsvd", n, kappa, 1)
θ = matrixdepot("randsvd", 1)
try
    matridepot("randsvd", 5, kappa, 6)
catch ArgumentError
    println("randsvd: invalid mode value")
end

println("'randsvd' passed test...")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailabilitySetProperties <: SwaggerModel
    platformUpdateDomainCount::Any # spec type: Union{ Nothing, Int32 } # spec name: platformUpdateDomainCount
    platformFaultDomainCount::Any # spec type: Union{ Nothing, Int32 } # spec name: platformFaultDomainCount
    virtualMachines::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: virtualMachines
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function AvailabilitySetProperties(;platformUpdateDomainCount=nothing, platformFaultDomainCount=nothing, virtualMachines=nothing, statuses=nothing)
        o = new()
        validate_property(AvailabilitySetProperties, Symbol("platformUpdateDomainCount"), platformUpdateDomainCount)
        setfield!(o, Symbol("platformUpdateDomainCount"), platformUpdateDomainCount)
        validate_property(AvailabilitySetProperties, Symbol("platformFaultDomainCount"), platformFaultDomainCount)
        setfield!(o, Symbol("platformFaultDomainCount"), platformFaultDomainCount)
        validate_property(AvailabilitySetProperties, Symbol("virtualMachines"), virtualMachines)
        setfield!(o, Symbol("virtualMachines"), virtualMachines)
        validate_property(AvailabilitySetProperties, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type AvailabilitySetProperties

const _property_map_AvailabilitySetProperties = Dict{Symbol,Symbol}(Symbol("platformUpdateDomainCount")=>Symbol("platformUpdateDomainCount"), Symbol("platformFaultDomainCount")=>Symbol("platformFaultDomainCount"), Symbol("virtualMachines")=>Symbol("virtualMachines"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_AvailabilitySetProperties = Dict{Symbol,String}(Symbol("platformUpdateDomainCount")=>"Int32", Symbol("platformFaultDomainCount")=>"Int32", Symbol("virtualMachines")=>"Vector{SubResource}", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ AvailabilitySetProperties }) = collect(keys(_property_map_AvailabilitySetProperties))
Swagger.property_type(::Type{ AvailabilitySetProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AvailabilitySetProperties[name]))}
Swagger.field_name(::Type{ AvailabilitySetProperties }, property_name::Symbol) =  _property_map_AvailabilitySetProperties[property_name]

function check_required(o::AvailabilitySetProperties)
    true
end

function validate_property(::Type{ AvailabilitySetProperties }, name::Symbol, val)
end


text:
using Base.Test
using SparseRF

tests = ["classification_rand.jl", "regression_rand.jl",
         "classification_scikitlearn.jl", "regression_scikitlearn.jl"]

println("Running tests...")
for curtest in tests
    println("Test: $curtest")
    include(curtest)
    println("=" ^ 50)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "Frei0r"
version = v"1.6.1"
# Collection of sources required to build imagemagick

sources = [
    "https://files.dyne.org/frei0r/frei0r-plugins-1.6.1.tar.gz" =>
    "e0c24630961195d9bd65aa8d43732469e8248e8918faa942cfb881769d11515e",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd frei0r-plugins-1.6.1/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make -j${ncore}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # Windows
    Windows(:i686),
    Windows(:x86_64),

    # linux
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),

    # musl
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl),

    # The BSD's
    FreeBSD(:x86_64),
    MacOS(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "addition", :addition),
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
abstract type Runner; end

function target_nbits(target::AbstractString)
    if startswith(target, "i686-") || startswith(target, "arm-")
        return "32"
    else
        return "64"
    end
end

function target_proc_family(target::AbstractString)
    if startswith(target, "arm") || startswith(target, "aarch")
        return "arm"
    elseif startswith(target, "power")
        return "power"
    else
        return "intel"
    end
end

function target_dlext(target::AbstractString)
    if endswith(target, "-mingw32")
        return "dll"
    elseif contains(target, "-apple-")
        return "dylib"
    else
        return "so"
    end
end

"""
    target_envs(target::String)

Given a `target` (this term is used interchangeably with `triplet`), generate a
`Dict` mapping representing all the environment variables to be set within the
build environment to force compiles toward the defined target architecture.
Examples of things set are `PATH`, `CC`, `RANLIB`, as well as nonstandard
things like `target`.
"""
function target_envs(target::AbstractString)
    global use_ccache

    # Helper function to generate paths such as /opt/x86_64-apple-darwin14/bin/llvm-ar
    tool = x -> "/opt/$(target)/bin/$(x)"
    # Helper function to generate paths such as /opt/x86_64-linux-gnu/bin/x86_64-linux-gnu-gcc
    target_tool = x -> tool("$(target)-$(x)")

    # Start with the default musl ld path:
    lib_path = "/usr/local/lib64:/usr/local/lib:/lib:/usr/local/lib:/usr/lib"

    # Then add on our target-specific locations
    lib_path *= ":/opt/$(target)/lib64:/opt/$(target)/lib"
    lib_path *= ":/opt/$(target)/$(target)/lib64:/opt/$(target)/$(target)/lib"

    # Finally add on our destination location
    lib_path *= ":/workspace/destdir/lib64:/workspace/destdir/lib"

    # Start with the standard PATH:
    path = "/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin"

    # Slip our tools into the front
    path = "/opt/super_binutils/bin:/opt/$(target)/bin:" * path

    mapping = Dict(
        # Activate the given target via `PATH` and `LD_LIBRARY_PATH`
        "PATH" => path,
        "LD_LIBRARY_PATH" => lib_path,

        # binutils/toolchain envvars
        "STRIP" => target_tool("strip"),
        "LIBTOOL" => target_tool("libtool"),
        "LIPO" => target_tool("lipo"),

        # Useful tools for our buildscripts
        "target" => target,
        "nproc" => "$(Sys.CPU_CORES)",
        "nbits" => target_nbits(target),
        "proc_family" => target_proc_family(target),
        "dlext" => target_dlext(target),
        "TERM" => "screen",

        # We should always be looking for packages already in the prefix
        "PKG_CONFIG_PATH" => "/workspace/destdir/lib/pkgconfig",
        "PKG_CONFIG_SYSROOT_DIR" => "/workspace/destdir",

        # We like to be able to get at our .bash_history afterwards. :)
        "HISTFILE"=>"/meta/.bash_history",

        # There is no broad agreement of what this should be called,
        # so we set all the environment variables that we've seen it called as:
        "CC_FOR_BUILD" => "/opt/x86_64-linux-gnu/bin/gcc",
        "BUILD_CC"     => "/opt/x86_64-linux-gnu/bin/gcc",
        "HOSTCC"       => "/opt/x86_64-linux-gnu/bin/gcc",
    )

    # If we're on MacOS or FreeBSD, we default to LLVM tools instead of GCC.
    # On all of our clangy platforms we actually have GCC tools available as well,
    # they're just not used by default.  Override these environment variables in
    # your scripts if you want to use them.
    if contains(target, "-apple-") || contains(target, "-freebsd")
        mapping["AR"] = tool("llvm-ar")
        mapping["AS"] = tool("llvm-as")
        mapping["CC"] = tool("clang")
        mapping["CXX"] = tool("clang++")
        # flang isn't a realistic option yet, so we still use gfortran here
        mapping["FC"] = target_tool("gfortran")
        mapping["LD"] = tool("llvm-ld")
        mapping["NM"] = tool("llvm-nm")
        mapping["OTOOL"] = target_tool("otool")
        mapping["INSTALL_NAME_TOOL"] = target_tool("install_name_tool")
        mapping["LLVM_TARGET"] = target
        mapping["RANLIB"] = tool("llvm-ranlib")
    else
        mapping["AR"] = target_tool("ar")
        mapping["AS"] = target_tool("as")
        mapping["CC"] = target_tool("gcc")
        mapping["CXX"] = target_tool("g++")
        mapping["FC"] = target_tool("gfortran")
        mapping["LD"] = target_tool("ld")
        mapping["NM"] = target_tool("nm")
        mapping["RANLIB"] = target_tool("ranlib")
    end

    # On OSX, we need to do a little more work.
    if contains(target, "-apple-")
        # First, tell CMake what our deployment target is, so that it tries to
        # set -mmacosx-version-min and such appropriately
        mapping["MACOSX_DEPLOYMENT_TARGET"] = "10.8"

        # Also put this into LDFLAGS because some packages are hard of hearing
        mapping["LDFLAGS"] = "-mmacosx-version-min=10.8"
    end

    # If we're using `ccache`, prepend it to `CC`, `CXX`, `FC`, etc....
    if use_ccache
        for k in ("CC", "CXX", "FC")
            mapping[k] = string("ccache ", mapping[k])
        end
    end

    return mapping
end

runner_override = ""
function preferred_runner()
    global runner_override
    if runner_override != ""
        if runner_override in ["userns", "privileged"]
            return UserNSRunner
        elseif runner_override in ["qemu"]
            return QemuRunner
        end
    end

    @static if Compat.Sys.islinux()
        # If runner_override is not yet set, let's probe to see if we can use
        # unprivileged containers, and if we can't, switch over to privileged.
        if !probe_unprivileged_containers()
            msg = strip("""
            Unable to run unprivileged containers on this system!
            This may be because your kernel does not support mounting overlay
            filesystems within user namespaces. To work around this, we will
            switch to using privileged containers. This requires the use of
            sudo. To choose this automatically, set the BINARYBUILDER_RUNNER
            environment variable to "privileged" before starting Julia.
            """)
            Compat.@warn(replace(msg, "\n" => " "))
            runner_override = "privileged"
        end

        return UserNSRunner
    else
        return QemuRunner
    end
end

"""
    runshell(platform::Platform = platform_key())

Launch an interactive shell session within the user namespace, with environment
setup to target the given `platform`.
"""
function runshell(platform::Platform = platform_key(); verbose::Bool = false)
    runshell(preferred_runner(), platform; verbose=verbose)
end


text:
using Optim

# Find the posterior mode
function reginv(y_data::Vector{Float64}, mat::Matrix{Float64}, sigma::Float64=1.0)
	x = ones(size(mat)[1])

	# logprob shorthand that only depends on x
	function badness(x)
		return -logprob(x, y_data, mat, sigma)
	end
	return optimize(badness, x, method=:cg, show_trace=true).minimum
end

function logprob(x::Vector{Float64}, y_data::Vector{Float64}, mat::Matrix{Float64}, sigma::Float64)
	y_model = mat*x
	return -3*sum(x.^2) - 0.5*sum((y_data - y_model).^2/sigma^2)
end

max_stars_repo_name: eggplantbren/MeanImage

max_issues_repo_name: eggplantbren/MeanImage

max_forks_repo_name: eggplantbren/MeanImage


text:
abstract Bootstrap <: StreamStat

## Double-or-nothing online bootstrap
type BernoulliBootstrap{S <: ContinuousUnivariateStreamStat} <: Bootstrap
    replicates::Vector{S}            # replicates of base stat
    cached_state::Vector{Float64}    # cache of replicate states
    n::Int                           # number of observations
    cache_is_dirty::Bool
end

## Poisson weighted online bootstrap
type PoissonBootstrap{S <: ContinuousUnivariateStreamStat} <: Bootstrap
    replicates::Vector{S}           # replicates of base stat
    cached_state::Vector{Float64}  # cache of replicate states
    n::Int                          # number of observations
    cache_is_dirty::Bool
end

# Frozen bootstrap object are generated when two bootstrap distributions
# are combined, e.g., if they are differenced.
immutable FrozenBootstrap <: Bootstrap
    cached_state::Vector{Float64}  # cache of replicate states
    n::Int                          # number of observations
end

# Double or nothing bootstrap
function BernoulliBootstrap{S <: ContinuousUnivariateStreamStat}(
    stat::S,
    R::Int = 1_000,
    α::Real = 0.05,
)
    replicates = S[copy(stat) for i in 1:R]
    cached_state = Array(Float64, R)
    return BernoulliBootstrap(replicates, cached_state, 0, true)
end

function PoissonBootstrap{S <: ContinuousUnivariateStreamStat}(
    stat::S,
    R::Int = 1_000,
    α::Float64 = 0.05,
)
    replicates = S[copy(stat) for i in 1:R]
    cached_state = Array(Float64, R)
    return PoissonBootstrap(replicates, cached_state, 0, true)
end

function update!(stat::BernoulliBootstrap, args::Any...)
    stat.n += 1

    for replicate in stat.replicates
        if rand() > 0.5
            update!(replicate, args...)
            update!(replicate, args...)
        end
    end
    stat.cache_is_dirty = true

    return
end

function update!(stat::PoissonBootstrap, args::Any...)
    stat.n += 1

    for replicate in stat.replicates
        repetitions = rand(Distributions.Poisson(1))
        for repetition in 1:repetitions
            update!(replicate, args...)
        end
    end
    stat.cache_is_dirty = true

    return
end

function update!(stat::FrozenBootstrap, args::Any...)
    error("Cannot update a FrozenBootstrap object")
    return
end

# TODO: Make this work with non-univariate statistics by taking marginal
#       quantiles
function replicates(stat::Bootstrap)
    return stat.replicates
end

function replicates(stat::FrozenBootstrap)
    return stat.cached_state
end

function Base.mean(stat::Bootstrap)
  return mean(cached_state(stat))
end

function Base.std(stat::Bootstrap)
  return std(cached_state(stat))
end

state(stat::Bootstrap) = stat.replicates

# update cached_state' states if necessary and return their values
function cached_state(stat::Bootstrap)
    if stat.cache_is_dirty
        for (i, replicate) in enumerate(stat.replicates)
            stat.cached_state[i] = state(replicate)
        end
        stat.cache_is_dirty = false
    end
    return stat.cached_state
end

function cached_state(stat::FrozenBootstrap)
    return stat.cached_state
end

function ci(stat::Bootstrap, α=0.05, method=:quantile)
    states = cached_state(stat)

    # If any NaN, return NaN, NaN
    if any(isnan, states)
        return (NaN, NaN)
    else
        if method == :quantile
            return (
                quantile(states, α / 2),
                quantile(states, 1 - α / 2),
            )
        elseif method == :normal
            norm_approx = Distributions.Normal(
                mean(states),
                std(states)
            )
            return (
                quantile(norm_approx, α / 2),
                quantile(norm_approx, 1 - α / 2)
            )
        else
            error("Unrecognized confidence interval type: ", ci)
        end
    end
end

nobs(stat::Bootstrap) = stat.n

# Assumes a and b are independent.
function Base.(:-)(a::Bootstrap, b::Bootstrap)
    return FrozenBootstrap(
        cached_state(a) - cached_state(b),
        nobs(a) + nobs(b)
    )
end

# Base.copy(stat::Bootstrap)

# function Base.merge(stat1::Bootstrap, stat2::Bootstrap)

# function Base.empty!(stat::Bootstrap)

function Base.rand(stat::Bootstrap)
    cs = cached_state(stat)
    cs[rand(1:length(cs))]
end

function Base.show(io::IO, stat::Bootstrap)
    @printf("%s:\n", typeof(stat))
    @printf(" * Replicates: %d\n", length(stat.cached_state))
    @printf(" * Observations: %d\n", stat.n)
    lower, upper = ci(stat, 0.05)
    @printf(" * Confidence Interval: [%f, %f]", lower, upper)
end


text:
using FastBEAST
using Pkg
Pkg.add("Documenter")
using Documenter

makedocs(
         sitename = "FastBEAST.jl",
         modules  = [FastBEAST],
         pages=[
                "Home" => "index.md",
                "Getting Started" => "gstarted.md",
                "Manual" => Any[
                    "man/clustering.md",
                    "man/hmatrix.md",
                    "man/aca.md"
                    ],
                "Types and Functions" => "functions.md"
               ])
deploydocs(
    repo="github.com/sbadrian/FastBEAST",
    devbranch = "main"
)

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
const DNASeq = BioSequence{DNAAlphabet{2}}
const DNASeqUnion = Union{DNASequence, DNASeq}
const Phred = Int8
const Prob = Float64
const LogProb = Float64
const Score = Float64


text:
# Alternative build process where the package provides a `deps/Project.toml`
# project that Pkg activates and instantiate before running the build script.
#
# The code in this file would be internal to Pkg, and packages would provide
# what's in `make.jl`.

append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg

mktempdir() do tmp
    build_project = joinpath(@__DIR__, "Project.toml")
    build_manifest = joinpath(@__DIR__, "Manifest.toml")
    makefile = joinpath(@__DIR__, "make.jl")

    tmp_build_project = joinpath(tmp, "Project.toml")
    tmp_build_manifest = joinpath(tmp, "Manifest.toml")

    # Copy build project to tmp directory
    cp(build_project, tmp_build_project)
    if isfile(build_manifest)
        cp(build_manifest, tmp_build_manifest)
    end
    # Activate and instantiate
    Pkg.activate(tmp_build_project)
    Pkg.instantiate()
    # Run build script
    withenv("JULIA_LOAD_PATH" => tmp_build_project) do
        run(`$(Base.julia_cmd()) --startup-file=no $(makefile)`)
    end
end


text:
@testset "Test creating `Lattice` with units" begin
    a = 4u"nm"
    b = 180u"bohr"
    c = 3u"angstrom"
    @test Lattice(a, b, c, 90, 90, 90) == Lattice(
        [
            4u"nm" 0u"m" 0.0u"cm"
            0u"cm" 180.0u"bohr" 0u"m"
            0u"bohr" 0u"nm" (3//1)*u"angstrom"
        ],
    )
end


text:
"""
requires git

returns location of all the models (every version), and the

for a particular version just do

    `filter(x->occursin("l1v2", x), fns)`
"""
function sbml_test_suite(repo_path=joinpath(datadir, "sbml-test-suite"))
    isdir(repo_path) && return nothing
    run(`git clone "https://github.com/anandijain/sbml-test-suite" $(repo_path)`)
end

function get_sbml_suite_fns(repo_path=joinpath(datadir, "sbml-test-suite"); sfx="l3v2.xml")
    ds = filter(isdir, readdir(joinpath(repo_path, "semantic"); join=true))
    fns = reduce(vcat, glob.("*.xml", ds))
    # fs = map(x -> splitdir(x)[end], fns)
    sfx !==nothing && filter!(endswith(sfx), fns) # latest level and version
    fns
end

function setup_settings_txt(fn)
    ls = readlines(fn)
    spls = split.(ls, ": ")
    filter!(x->length(x) == 2, spls)
    Dict(map(x -> x[1] => Meta.parse(x[2]), spls))
end

"""
dir = "data/sbml-test-suite/semantic/00001/"
"""
function verify_case(dir;saveplot=false)
    res = false
    atol = 0
    err = ""
    try
        fns = readdir(dir;join=true)
        model_fn = filter(endswith("l2v3.xml"), fns)[1]
        case_no = basename(dirname(model_fn))
        settings = setup_settings_txt(filter(endswith("settings.txt"), fns)[1])
        results = CSV.read(filter(endswith("results.csv"), fns)[1], DataFrame)

        ml = SBML.readSBML(model_fn, doc -> begin
            set_level_and_version(3, 2)(doc)
            convert_simplify_math(doc)
            end)
        rs = ReactionSystem(ml)
        sys = ODESystem(ml)
        statenames = [string(s.f.name) for s in sys.states]

        ts = LinRange(settings["start"], settings["duration"], settings["steps"]+1)
        prob = ODEProblem(sys, Pair[], (settings["start"], Float64(settings["duration"])); saveat=ts)
        sol = solve(prob, CVODE_BDF(); abstol=settings["absolute"], reltol=settings["relative"])
        solm = Array(sol)'
        m = Matrix(results[1:end, 2:end])[:, sortperm(sortperm(statenames))]
        res = isapprox(solm, m; atol=1e-9, rtol=3e-2)
        diff = m .- solm
        atol = maximum(diff)
        saveplot && !res && verify_plot(case_no, rs, solm, m)
        return [dir, res, atol, err]
    catch e
        err = string(e)
        return [dir, res, atol, err]
    end
end

function verify_all(ds;verbose=true, saveplot=false)
    df = DataFrame(dir=String[], retcode=Bool[], atol=Float64[], error=String[])
    for dir in ds
        ret = verify_case(dir; saveplot=saveplot)
        verbose && @info ret
        push!(df, ret)
    end
    verbose && print(df)
    df
end

"plots the difference between the suites' reported solution and DiffEq's sol"
function verify_plot(case_no, rs, solm, m)
    sys = convert(ODESystem, rs)
    open(joinpath(logdir, case_no*".txt"), "w") do file
        write(file, "Reactions:\n")
        write(file, repr(equations(rs))*"\n")
        write(file, "ODEs:\n")
        write(file, repr(equations(sys))*"\n")
    end
    plt = plot(solm)
    plt = plot!(m, linestyle=:dot)
    savefig(joinpath(logdir, case_no*".png"))
end


text:
using GAP_pkg_cvec
using Test

import GAP_pkg_cvec
println("TODO")


text:
# This file was generated, do not modify it. # hide
rms(y[test], predict(tree, rows=test))


text:
abstract type AbstractChargeDensityModel{T <: SSDFloat} end

@inline function ChargeDensityModel(T::DataType, dict::Union{Dict{String, Any}, Dict{Any, Any}})
    return ChargeDensityModel(T, Val{Symbol(dict["name"])}(), dict)
end

"""
    struct LinearChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}

Simple charge density model which assumes a linear gradient in charge density in each dimension.
`offsets::NTuple{3, T}` are the charge densities at 0 and `gradients::NTuple{3, T}` are the linear
slopes in each dimension.
"""
struct LinearChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}
    offsets::NTuple{3, T}
    gradients::NTuple{3, T}
end

function get_charge_density(lcdm::LinearChargeDensityModel{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    ρ::T = 0
    for i in eachindex(lcdm.offsets)
        ρ += (lcdm.offsets[i] + pt[i] * lcdm.gradients[i]) * T(1e16) # * T(1e10) * T(1e6) -> 1/cm^3 -> 1/m^3
    end
    return ρ
end

function ChargeDensityModel(T::DataType, t::Val{:linear}, dict::Union{Dict{String, Any}, Dict{Any, Any}})
    return LinearChargeDensityModel{T}( dict )
end

function LinearChargeDensityModel{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}})::LinearChargeDensityModel{T} where {T <: SSDFloat}
    offsets, gradients = zeros(T,3), zeros(T,3)
    if haskey(dict, "r")     offsets[1] = dict["r"]["init"];     gradients[1] = dict["r"]["gradient"]    end
    if haskey(dict, "phi")   offsets[2] = dict["phi"]["init"];   gradients[2] = dict["phi"]["gradient"]  end
    if haskey(dict, "z")     offsets[3] = dict["z"]["init"];     gradients[3] = dict["z"]["gradient"]    end
    if haskey(dict, "x")     offsets[1] = dict["x"]["init"];     gradients[1] = dict["x"]["gradient"]    end
    if haskey(dict, "y")     offsets[2] = dict["y"]["init"];     gradients[2] = dict["y"]["gradient"]    end
    LinearChargeDensityModel{T}( NTuple{3,T}(offsets), NTuple{3,T}(gradients) )

end


text:
using TableLogger
using Test

@testset "TableLogger.jl" begin
    include("unit/main.jl")
end


text:
"""
    testpattern() -> Matrix{RGBA{N0f8}}

Load and return the provided 300x400 test image.

The returned image was specifically designed to be informative
about the effects of the applied augmentation operations. It is
thus well suited to prototype an augmentation pipeline, because it
makes it easy to see what kind of effects one can achieve with it.
"""
testpattern() = load(joinpath(@__DIR__, "..", "resources", "testpattern.png"))

function use_testpattern()
    info("No custom image specifed. Using \"testpattern()\" for demonstration.")
    testpattern()
end

# --------------------------------------------------------------------
# rand() is not threadsafe (https://discourse.julialang.org/t/4683)

# Because we only require random numbers to sample parameters
# and not the actual expensive computation, this seems like a better
# approach than using separate RNG per thread.
const rand_mutex = Ref{Threads.Mutex}()

# constant overhead of about 80 ns compared to unsafe rand
function safe_rand(args...)
    lock(rand_mutex[])
    result = rand(args...)
    unlock(rand_mutex[])
    result
end

# --------------------------------------------------------------------

@inline maybe_copy(A::OffsetArray) = A
@inline maybe_copy(A::Array) = A
@inline maybe_copy(A::SArray) = A
@inline maybe_copy(A::MArray) = A
@inline maybe_copy(A::AbstractArray) = match_idx(collect(A), indices(A))
@inline maybe_copy(A::Tuple) = map(maybe_copy, A)

# --------------------------------------------------------------------

@inline _plain_array(A::OffsetArray) = parent(A)
@inline _plain_array(A::Array) = A
@inline _plain_array(A::SArray) = A
@inline _plain_array(A::MArray) = A
@inline _plain_array(A::Tuple) = map(_plain_array, A)
# avoid recursion
@inline plain_array(A) = _plain_array(maybe_copy(A))

# --------------------------------------------------------------------

@inline plain_indices(A::Array) = A
@inline plain_indices(A::OffsetArray) = parent(A)
@inline plain_indices(A::AbstractArray) = _plain_indices(A, indices(A))
@inline plain_indices(A::SubArray) = _plain_indices(A, A.indexes)

@inline function _plain_indices(A::AbstractArray{T,N}, ids::NTuple{N,Base.OneTo}) where {T, N}
    A
end

@inline function _plain_indices(A::AbstractArray, ids::Tuple{Vararg{Any}})
    view(A, indices(A)...)
end

@inline function _plain_indices(A::SubArray{T,N}, ids::NTuple{N,IdentityRange}) where {T, N}
    view(parent(A), indices(A)...)
end

# --------------------------------------------------------------------

@inline match_idx(buffer::AbstractArray, inds::Tuple) = buffer
@inline match_idx(buffer::Union{Array,SubArray}, inds::NTuple{N,UnitRange}) where {N} =
    OffsetArray(buffer, inds)

# --------------------------------------------------------------------

function indirect_indices(::Tuple{}, ::Tuple{})
    throw(MethodError(indirect_indices, ((),())))
end

@inline function indirect_indices(O::NTuple{N,Base.OneTo}, I::NTuple{N,AbstractUnitRange}) where N
    map(IdentityRange, I)
end

@inline function indirect_indices(O::NTuple{N,Base.OneTo}, I::NTuple{N,StepRange}) where N
    I
end

function indirect_indices(O::NTuple{N,AbstractUnitRange}, I::NTuple{N,AbstractUnitRange}) where N
    map((i1,i2) -> IdentityRange(UnitRange(i1)[i2]), O, I)
end

function indirect_indices(O::NTuple{N,AbstractUnitRange}, I::NTuple{N,StepRange}) where N
    map((i1,i2) -> UnitRange(i1)[i2], O, I)
end

function indirect_indices(O::NTuple{N,StepRange}, I::NTuple{N,Range}) where N
    map((i1,i2) -> i1[i2], O, I)
end

# --------------------------------------------------------------------

function indirect_view(A::AbstractArray, I::Tuple)
    view(A, indirect_indices(indices(A), I)...)
end

function indirect_view(A::SubArray{T,N,TA,<:NTuple{N,Range}}, I::Tuple) where {T,N,TA}
    view(parent(A), indirect_indices(A.indexes, I)...)
end

# --------------------------------------------------------------------

function direct_indices(::Tuple{}, ::Tuple{})
    throw(MethodError(direct_indices, ((),())))
end

# TODO: Figure out why this method exists
function direct_indices(O::NTuple{N,IdentityRange}, I::NTuple{N,StepRange}) where N
    throw(MethodError(direct_indices, (O, I)))
end

@inline function direct_indices(O::NTuple{N,Range}, I::NTuple{N,AbstractUnitRange}) where N
    map(IdentityRange, I)
end

@inline function direct_indices(O::NTuple{N,Range}, I::NTuple{N,StepRange}) where N
    I
end

# --------------------------------------------------------------------

function direct_view(A::AbstractArray{T,N}, I::NTuple{N,Range}) where {T,N}
    view(A, direct_indices(indices(A), I)...)
end

function direct_view(A::SubArray{T,N,TA,<:NTuple{N,Range}}, I::NTuple{N,Range}) where {T,N,TA}
    view(A, direct_indices(A.indexes, I)...)
end

# --------------------------------------------------------------------

@inline vectorize(A::AbstractVector) = A
@inline vectorize(A::Real) = A:A

@inline round_if_float(num::Integer, d) = num
round_if_float(num::AbstractFloat, d) = round(num,d)
round_if_float(nums::Tuple, d) = map(num->round_if_float(num,d), nums)

function unionrange(i1::AbstractUnitRange, i2::AbstractUnitRange)
    map(min, first(i1), first(i2)):map(max, last(i1), last(i2))
end

@inline _showcolor(io::IO, T::Type{<:Number}) = print(io, T)
@inline _showcolor(io::IO, T) = ColorTypes.colorant_string_with_eltype(io, T)

# --------------------------------------------------------------------

function _2dborder!(A::AbstractArray{T,3}, val::T) where T
    ndims, h, w = size(A)
    @inbounds for j = (1,w), i = 1:h
        for d = 1:ndims
            A[d,i,j] = val
        end
    end
    @inbounds for j = 1:w, i = (1,h)
        for d = 1:ndims
            A[d,i,j] = val
        end
    end
    A
end


text:
export Model

"""
    abstract type Model end

Abstract type for the representation of a physical model for which calculations can be performed.
Needs to be implemented system-specifically, in particular the following functions need to be defined: (for now)

- the one-particle matrix element, `energy(::Model, ::Orbital)`
- the two-particle matrix element, `w(::Model, ::Orbital, ::Orbital, ::Orbital, ::Orbital)`

"""
abstract type Model end

function energy(m::Model, o::Orbital)
    error("missing implementation of energy(::$(typeof(m)), ::$(typeof(o)))")
end

function w(m::Model, o1::Orbital, o2::Orbital, o3::Orbital, o4::Orbital)
    error("missing implementation of w(::$(typeof(m)), ::$(typeof(o1)), ::$(typeof(o2)))")
end

@doc raw"""
    wminus(::Model, i,j,k,l)

return the difference of the two-particle matrix element with the same but the last two indices transposed
antisymmetric difference of the two-particle matrix elements:

$w^-_{ijkl} = w_{ijkl} - w_{ijlk}$

This is called the antisymmetrized two-particle matrix element.
This is an abbreviation for these terms arising in the
Slater-Condon rules for the calculation of the many-body matrix elements
via the one- and two-particle matrix elements of the underlying single-particle basis.
"""
wminus(m::Model, i,j,k,l) = w(m,i,j,k,l) - w(m,i,j,l,k)

@doc raw"""
    Woffdiag_element(::Model, ::Ensemble, ::Orbital, ::Orbital, ::Orbital, ::Orbital)

Return the offdiagonal many body matrix element of the interaction operator

$\frac{1}{4} ( w_{ijkl} - w_{ijlk} ) (\pm \langle\{\tilde{n}\}|\{n\}^{ij}_{kl}\rangle)$

for an excitation given by creating orbitals i,j and annihilating orbitals k, l
as given by the Slater-Condon rules.
"""
function Woffdiag_element(m::Model, e::Ensemble, i::Orbital, j::Orbital, k::Orbital, l::Orbital)
    # We sample with the weight of antisymmetrized matrix element but we do not restrict
    # the order of indices of our possible kinks. We therefor need an extra factor 1/4 in the weight-function
    1/4 * e.λ * wminus(m, i,j,k,l)
end

Woffdiag_element(m::Model, e::Ensemble, kink::T4) = Woffdiag_element(m, e, kink.i, kink.j, kink.k, kink.l)

"""
    ΔWoffdiag_element(::Model, e::Ensemble, itr, itr)

return the change in the offdiagonal many body matrix element of the interaction operator
given by adding the kinks in the first iterable and removing the kinks in the second iterable

both arguments `itr` are required to be iterables containing kinks
"""
function ΔWoffdiag_element(m::Model, e::Ensemble, add_kinks, drop_kinks)
    if isempty(drop_kinks)
        prod(Woffdiag_element(m, e, k) for k in add_kinks)
    elseif isempty(add_kinks)
        1.0 / prod(Woffdiag_element(m, e, k) for k in drop_kinks)
    else
        prod(Woffdiag_element(m, e, k) for k in add_kinks) / prod(Woffdiag_element(m, e, k) for k in drop_kinks)
    end
end

ΔWoffdiag_element(m::Model, e::Ensemble, add_kinks::Kinks, drop_kinks::Kinks) = ΔWoffdiag_element(m, e, excitations(add_kinks), excitations(drop_kinks))

"""
    ΔW_diag(m::Model, i, j, k, l, occ)

change in the diagonal interaction matrix element due to a change in the occupation `occ`
in a periodic interval (τ1,τ2) where no kinks occur
the change in the occupation is assumed to consist in a creation of two orbitals i, j
and in the annihlation of two orbitals k, l
"""
function ΔW_diag(m::Model, i, j, k, l, occ)
    @assert (i ∉ occ) & (j ∉ occ) "Calculation of the change in the many-body diagonal interaction matrix element: This function assumes that the first two orbitals\n\t $(i)\n and\n\t $(j) given are the creator orbitals and thus that they are not occupied in the given occupation\n\t $(occ). "
    # contributions due to mean field interactions of the annihilated orbitals
    Δ = sum( wminus(m, ν,k,k,ν) + wminus(m, ν,l,l,ν) for ν in drop_orbs(occ, [k,l]) )# interactions of mean field with k and l
    Δ += wminus(m, k,l,l,k)# interaction between k and l
    # contributions due to mean field interactions of the created orbitals
    # note: the annihilator orbitals k, l are not in the new occupation
    Δ -= sum( wminus(m, ν,i,i,ν) + wminus(m, ν,j,j,ν) for ν in drop_orbs(occ, [k,l]) )# interactions of mean field with i and j
    Δ -= wminus(m, i,j,j,i)# interaction between i and j
    return Δ
end

"""
    ΔW_diag(m::Model, i, j, occ)

change in the diagonal interaction matrix element due to a change in the occupation `occ`
in a periodic interval (τ1,τ2) where no kinks occur
the change in the occupation is assumed to consist in a creation of one orbitals i
and in the annihlation of one orbitals j
"""
function ΔW_diag(m::Model, i, j, occ)
    @assert (i ∉ occ) "Calculation of the change in the many-body diagonal interaction matrix element: This function assumes that the first two orbitals\n\t $(i)\n and\n\t $(j) given are the creator orbitals and thus that they are not occupied in the given occupation\n\t $(occ). "
    # contributions due to mean field interactions of the annihilated orbitals
    Δ = sum( wminus(m, ν,j,j,ν) for ν in drop_orbs(occ, (j,)) )# interactions of mean field with k
    # contributions due to mean field interactions of the created orbitals
    Δ -= sum( wminus(m, ν,i,i,ν) for ν in drop_orbs(occ, (j,)) )
end

### convention: all ΔX_element represent only the difference in the matrix elements will be used as exp(-Δ) for the weight change
###             thus, contributions from the new (proposed) configuration (creators) appear positive (+)
###             and the contributions from the old configuration (annihilators) appear negative (-)

"""
    ΔT_element(::Model, i,j,k,l)

return the change in the kinetic many body matrix element due
to creating orbitals i, j and annihilating orbitals k, l
"""
ΔT_element(m::Model, i,j,k,l) = energy(m, i) + energy(m, j) - energy(m, k) - energy(m, l)

"""
    ΔWdiag_element(::Model, ::Ensemble, ::Configuration, i, j, k, l, τ1, τ2)

Calculate the change in the diagonal interaction many-body matrix element
due to a change in the occupations given by creating two orbitals i and j
and annihilating two orbitals k, l in the interval (τ1, τ2).
This interval may be periodically extended over the bounds (0,1) if τ1 > τ2,
i.e. the change in the occupation is considered for (τ1,1] ∪ [0,τ2) in that case.
We do not need to evaluate the diagonal interaction
between all orbitals in all time-intervalls, but it is sufficient to evaluate
the full diagonal interaction with the occupations at the start of the intervall
and then consider only contributions of orbitals that are changed by kinks in the intervall.
"""
function ΔWdiag_element(m::Model, e::Ensemble, c::Configuration, i, j, k, l, τ1, τ2)# TODO: assuming that i, j are creators and k, l are annihilators. Use Step instead ?
    @assert τ1 != τ2 " The diagonal interaction matrix element changes when kinks are added at different times and thus the occupations between the kinks are altered. It has no meaning to calculate this matrix element (or to add kinks) at equal times τ1=$(τ1), τ2=$(τ2). "

    # get all kinks between τ1 and τ2
    Ks = kinks_from_periodic_interval(c.kinks, τ1, τ2)

    # calculate Wdiag with the occupation at the start of the interval
    ΔWdiag = ΔW_diag(m, i,j,k,l, occupations_at(c,τ1)) * Δ(τ1,τ2)

    if !isempty(Ks)
        # calculate contrubutions to ΔWdiag from the orbitals changed by kinks in the intervall:
        # add a contribution if an orbital is created and
        # remove a contribution if an orbital is annilated
        # this is more efficient than calculating the occupation for each consecutive time-interval
        # via `occupation(occ,t)` since this function applies all kinks up to t::ImgTime
        τs = times_from_periodic_interval(c.kinks, τ1, τ2)# get a time-ordered list of the times of the kinks between τ1 and τ2
        ΔWdiag += sum( (ΔW_diag(m, i,j,k,l, creators(Ks[t1])) - ΔW_diag(m, i,j,k,l, annihilators(Ks[t1]))) * Δ(t1,τ2) for t1 in τs)
    end
    e.λ * ΔWdiag
end

"""
    ΔWdiag_element(::Model, ::Ensemble, ::Configuration, i, j, τ1, τ2)

Calculate the change in the diagonal interaction many-body matrix element
due to a change in the occupations given by creating an orbital i
and annihilating an orbital j in the interval (τ1, τ2).
This interval may be periodically extended over the bounds (0,1) if τ1 > τ2,
i.e. the change in the occupation is considered for (τ1,1] ∪ [0,τ2) in that case.
We do not need to evaluate the diagonal interaction
between all orbitals in all time-intervalls, but it is sufficient to evaluate
the full diagonal interaction with the occupations at the start of the intervall
and then consider only contributions of orbitals that are changed by kinks in the intervall.
"""
function ΔWdiag_element(m::Model, e::Ensemble, c::Configuration, i, j, τ1, τ2)# TODO: assuming that i is creator and j is annihilator. Use Step instead ?
    @assert τ1 != τ2 " The diagonal interaction matrix element changes when kinks are added at different times and thus the occupations between the kinks are altered. It has no meaning to calculate this matrix element (or to add kinks) at equal times τ1=$(τ1), τ2=$(τ2). "

    # get all kinks between τ1 and τ2
    Ks = kinks_from_periodic_interval(c.kinks, τ1, τ2)

    # calculate Wdiag with the occupation at the start of the interval
    ΔWdiag = ΔW_diag(m, i,j, occupations_at(c,τ1)) * Δ(τ1,τ2)
    if !isempty(Ks)
        # calculate contrubutions to ΔWdiag from the orbitals changed by kinks in the intervall:
        # add a contribution if an orbital is created and
        # remove a contribution if an orbital is annilated
        # this is more efficient than calculating the occupation for each consecutive time-interval
        # via `occupation(occ,t)` since this function applies all kinks up to t::ImgTime
        τs = times_from_periodic_interval(c.kinks, τ1, τ2)# get a time-ordered list of the times of the kinks between τ1 and τ2
        ΔWdiag += sum( (ΔW_diag(m, i,j, creators(Ks[t1])) - ΔW_diag(m, i,j, annihilators(Ks[t1]))) * Δ(t1,τ2) for t1 in τs)
    end
    e.λ * ΔWdiag
end

@doc raw"""
    sign_offdiagonal_product(::Model, ::Configuration)

Return the sign of the product of two-particle terms in the offdiagonal many-body matrix elements.
These are given by function `wminus`, i.e.

$\langle \{\tilde{n}\} | a^{\dagger}_i a^{\dagger}_j a_k a_l | \{n\} \rangle
    = \pm ( w_{ijkl} - w_{ijlk} ) \text{ for } \{\tilde{n}\} = \{n\}_{kl}^{ij}$

the term in the braces may be negative and this function returns the product of the
sign of all these contributions $w_{ijkl} - w_{ijlk}$ from all kinks.
The sign $\pm$ is determined from the permutation factor of the orbitals i,j,k,l
and is not calculated here.

used for the calculation the sign of the weight function
"""
function sign_offdiagonal_product(m::Model, c::Configuration)
    s = 1
    for κ in excitations(c.kinks)
        s *= sign(wminus(m, κ.i, κ.j, κ.k, κ.l))
    end
    return s
end

"""
    signum(m::Model, c::Configuration)

Calculate the sign of the `Configuration`'s weight.
"""
signum(m::Model, c::Configuration) = ladder_operator_order_factor(c.kinks)*sign_offdiagonal_product(m, c)