text:
export delta_image
"""
		$(SIGNATURES)

Function to generate a phantom with discrete points. The `distanceOfPoints` argument takes two functions
that take the number of the point to generate and return an integer. This makes the phantoms to generate highly
customizable.

# Arguments
- `size::Tuple{Integer, Integer}`: The size of the phantom
- `numOfPoints::Integer`: The number of points to generate
- `sizeOfPoint::Tuple{Integer, Integer}`: The size of the points in the phantom
- `distanceOfPoints::Tuple{Function, Function}`: The distance to add between each points in x and y direction
- `pivot::Tuple{Integer, Integer}`: The starting point to generate points towards (size, size)
- `circularShape::Bool`: If true, points are generated as circular
# Examples

## Two simple dots
```jldoctest
julia> image = delta_image((8, 8), 2; sizeOfPoint=(3, 2), distanceOfPoints=(x -> 0, x -> 4), pivot=(3, 3))
8×8 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```

## L-shaped arrangement
```jldoctest
julia> image = delta_image((8, 8), 3; sizeOfPoint=(2, 2), distanceOfPoints=(x -> x == 2 ? 3 : 0, x -> x == 3 ? -3 : 3), pivot=(3, 3))
8×8 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  0.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@testimage_gen function delta_image(
		size::Tuple{Integer, Integer},
		numOfPoints::Integer;
		sizeOfPoint::Tuple{Integer, Integer}=(1, 1),
		distanceOfPoints::Tuple{Function, Function}=(x -> 1, x -> 1),
		pivot::Tuple{Integer, Integer}=(1, 1),
		circularShape::Bool=false
	)

	image = zeros(Float64, size)

	(xPos, yPos) = pivot

	# Calculate the pixels that will get filled
	shape = []
	if !circularShape || (sizeOfPoint[1] <= 2 && sizeOfPoint[2] <= 2)
		shape = ones(sizeOfPoint)
	else
		shape = zeros(sizeOfPoint)
		halfX = (sizeOfPoint[1] + 1) / 2
		halfY = (sizeOfPoint[2] + 1) / 2

		for x in 1:sizeOfPoint[1]
			for y in 1:sizeOfPoint[2]
				if ((x - halfX)^2/(halfX - 1.0)^2) + ((y - halfY)^2/(halfY - 1.0)^2) <= 1.25
					shape[x, y] = 1
				end
			end
		end
	end

	for i in 1:numOfPoints
		if yPos + sizeOfPoint[2] > size[2] + 1 || xPos + sizeOfPoint[1] > size[1] + 1
			break
		end

		image[xPos:(xPos+sizeOfPoint[1] - 1), yPos:(yPos+sizeOfPoint[2] - 1)] = shape

		try
			xPos += distanceOfPoints[1](i+1)
			yPos += distanceOfPoints[2](i+1)
		catch
			println("The function to calculate the distances is not valid.")
			return image
		end

	end

	return image
end

export checker_image
"""
		$(SIGNATURES)

Function to generate a phantom with a checker board pattern. This function uses a best effort approach, meaning
that it is tried to cover most of the phantom with the pattern using the specified parameters.

# Arguments
- `size::Tuple{Integer, Integer}`: The size of the phantom
- `checkersCount::Tuple{Integer, Integer}`: How many squares to generate along each axis
- `stripeWidth::Tuple{Integer, Integer}`: By default `(1, 1)`. Sets the width of the lines between the squares

# Examples
```jldoctest
julia> image = checker_image((8, 8), (2, 3), (2, 1))
8×8 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@testimage_gen function checker_image(size::Tuple{Integer, Integer}=(8, 8), checkersCount::Tuple{Integer, Integer}=(2, 2), stripeWidth::Tuple{Integer, Integer}=(1, 1))
	image = zeros(Float64, size)

	# Calculate the space of each square
	xSpace = convert(Integer, floor((size[1] - (checkersCount[1]+1)*stripeWidth[1]) / checkersCount[1]))
	ySpace = convert(Integer, floor((size[2] - (checkersCount[2]+1)*stripeWidth[2]) / checkersCount[2]))

	# Calculate the rest of the division that will be added to the border of the checker board
	xBorder = convert(Integer, round((size[1] - (checkersCount[1]+1)*stripeWidth[1] - xSpace * checkersCount[1]) / 2))
	yBorder = convert(Integer, round((size[2] - (checkersCount[2]+1)*stripeWidth[2] - ySpace * checkersCount[2]) / 2))

	# Generate the stripes
	for x in 0:checkersCount[1]-1
		for y in 0:checkersCount[2]-1
			image[xBorder+stripeWidth[1]+1+x*(xSpace+stripeWidth[1]):xBorder+x*(xSpace+stripeWidth[1])+xSpace+stripeWidth[1],
				  yBorder+stripeWidth[2]+1+y*(ySpace+stripeWidth[2]):yBorder+y*(ySpace+stripeWidth[2])+ySpace+stripeWidth[2]] .= 1
		end
	end

	return image
end


text:
module TestBroadcast

using BangBang
using BangBang.Extras
using Test

@testset "broadcast!!" begin
    @testset "mutate" begin
        xs = fill(1234, 3)
        @test broadcast!!(+, xs, zeros(Int, 3), 1) === xs == ones(Int, 3)
    end

    @testset "widen" begin
        xs = fill(123, 3)
        @test broadcast!!(+, xs, [1, 2, 3], 0.5) == [1.5, 2.5, 3.5]
        @test all(==(123), xs)
    end
end

@testset "broadcast_inplace!!" begin
    @testset "mutate" begin
        xs = zeros(Int, 3)
        @test broadcast_inplace!!(+, xs, 1) === xs == ones(Int, 3)
    end

    @testset "widen" begin
        xs = ones(Int, 3)
        @test broadcast_inplace!!(+, xs, [1, 2, 3], 0.5) == [2.5, 3.5, 4.5]
        @test all(==(1), xs)
    end
end

end  # module

ext: jl


text:
mutable struct ADM{U, O} <: AbstractSubspaceOptimizer
    λ::U
    R::O
end

get_threshold(opt::ADM) = opt.λ
set_threshold!(opt::ADM, λ) = (opt.λ = λ; return)

"""
    ADM()
    ADM(λ = 0.1)

Optimizer for finding a sparse basis vector in a subspace based on [this paper](https://arxiv.org/pdf/1412.4659.pdf).
`λ` is the weight for the soft-thresholding operation.
"""
function ADM(λ::T = 0.1) where T <: Real
    f = NormL1(λ)
    return ADM(λ, f)
end

# ADM algorithm
function fit!(q::AbstractArray{T, 1}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10, tol::T = eps(eltype(q))) where T <: Real

    x = Y*q
    q_ = deepcopy(q)
    iters_ = 0

    while iters_ <= maxiter
        iters_ += 1
        prox!(x, opt.R, Y*q)
        mul!(q, Y', x)
        normalize!(q, 2)

        if norm(q - q_) < tol
            break
        else
            q_ .= q
        end
    end

    return iters_
end

# ADM initvary
function fit!(q::AbstractArray{T, 2}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10, tol::T = eps(eltype(q))) where T <: Real
    iters_ = Inf
    i_ = Inf

    @inbounds for i in 1:size(q, 2)
        i_ = fit!(q[:, i], Y, opt, maxiter = maxiter, tol = tol)
        if iters_ > i_
            iters_ = i_
        end
    end

    return iters_
end

# ADM pareto
function fit!(X::AbstractArray, A::AbstractArray, Y::AbstractArray, opt::ADM; rtol = 0.99, maxiter::Int64 = 1, convergence_error::T = eps(), f::Function = (xi, theta)->[norm(xi, 0); norm(theta'*xi, 2)], g::Function = x->norm(x),) where T <: Real
    # Return just the best candidate for the subspace optimization

    θ = zeros(eltype(A), size(X, 1), size(A, 2))
    θ[size(A, 1)+1:end, :] .= A

    fg(xi, theta) = (g∘f)(xi, theta)

    iters = Inf

    @inbounds for i in 1:size(Y, 1)
        for j in 1:size(A, 2)
            @views θ[1:size(A, 1), j] .= Y[i, j].*A[:, j]
        end

        N = nullspace(θ', rtol = rtol)
        Q = deepcopy(N) # Deepcopy for inplace
        # Find sparse vectors in nullspace
        # Calls effectively the ADM algorithm with varying initial conditions
        fit!(Q, N', opt, maxiter = maxiter)

        # Find sparse vectors in nullspace
        # Calls effectively the ADM algorithm with varying initial conditions
        iters_ = fit!(Q, N', opt, maxiter = maxiter, tol = convergence_error)
        iters_ < iters ? iters = iters_ : nothing

        # Compute pareto front
        for (j, ξ) in enumerate(eachcol(Q))
            if j == 1
                X[:, i] .= ξ
            else
                evaluate_pareto!(view(X, :, i), view(ξ, :), fg, view(θ, :, :))
            end
        end

        X[abs.(X[:, i]) .< get_threshold(opt), i] .= zero(eltype(X))
        @views X[:, i] .= X[:, i] ./ maximum(abs, X[:, i])
    end

    return iters
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ConstValueShape{T} <: AbstractValueShape

A `ConstValueShape` describes the shape of constant values of type `T`.

Constructor:

    ConstValueShape(value)

`value` may be of arbitrary type, e.g. a constant scalar value or array:

    ConstValueShape(4.2)
    ConstValueShape([11 21; 12 22])

Shapes of constant values have zero degrees of freedom (see
[`totalndof`](@ref)).

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ConstValueShape{T, strict} <: AbstractValueShape
    value::T
end

export ConstValueShape

ConstValueShape{T}(x::T) where T = ConstValueShape{T,true}(x)
ConstValueShape(x::T) where T = ConstValueShape{T,true}(x)

Base.:(==)(a::ConstValueShape, b::ConstValueShape) = a.value == b.value
Base.isequal(a::ConstValueShape, b::ConstValueShape) = isequal(a.value, b.value)
Base.isapprox(a::ConstValueShape, b::ConstValueShape; kwargs...) = isapprox(a.value, b.value; kwargs...)
Base.hash(x::ConstValueShape, h::UInt) = hash(x.value, hash(:ConstValueShape, hash(:ValueShapes, h)))

@inline Base.size(shape::ConstValueShape) = size(shape.value)
@inline Base.length(shape::ConstValueShape) = length(shape.value)

@inline Base.:(<=)(a::ConstValueShape{T}, b::ConstValueShape{U}) where {T,U} = T<:U && a.value ≈ b.value

@inline default_unshaped_eltype(shape::ConstValueShape) = Int32

@inline shaped_type(shape::ConstValueShape, ::Type{T}) where {T<:Real} = typeof(shape.value)

@inline totalndof(::ConstValueShape) = 0

# ToDo/Decision: Return a copy instead?
(shape::ConstValueShape)(::UndefInitializer) = shape.value

function unshaped(x::Any, shape::ConstValueShape)
    x == shape.value || throw(ArgumentError("Given value does not match value of ConstValueShape"))
    Float32[]
end

@inline _valshapeoftype(::Type{Nothing}) = ConstValueShape(nothing)

"""
    const_zero_shape(shape::ConstValueShape)

Get the equivalent of a constant zero shape for shape `shape` that will
only allow zero values to be set via an accessor.
"""
const_zero_shape(shape::ConstValueShape) = ConstValueShape(const_zero(shape.value))

"""
    nonstrict_const_zero_shape(shape::ConstValueShape)

Get the equivalent of a constant zero shape for shape `shape` that will
ignore any attempt to set a value via an accessor.

Useful as a gradient/tangent varshape of constants, as they can ignore
attempts to set non-zero values.
"""
function nonstrict_const_zero_shape(shape::ConstValueShape)
    x = const_zero(shape.value)
    ConstValueShape{typeof(x),false}(x)
end

replace_const_shapes(f::Function, shape::ConstValueShape) = f(shape)

const ConstAccessor{T,strict} = ValueAccessor{ConstValueShape{T,strict}}

@inline vs_getindex(data::AbstractVector{<:Real}, va::ConstAccessor) = va.shape.value

@inline vs_unsafe_view(::AbstractVector{<:Real}, va::ConstAccessor) = va.shape.value

# Zygote has a generic `@adjoint getindex(x::AbstractArray, inds...)` and same for view that
# will result in overwriting va.shape.value with dy without these custom adjoints:
ZygoteRules.@adjoint function getindex(x::AbstractVector{<:Real}, va::ConstAccessor)
    getindex(x, va), dy -> nothing, nothing
end
ZygoteRules.@adjoint function view(x::AbstractVector{<:Real}, va::ConstAccessor)
    view(x, va), dy -> nothing, nothing
end

function vs_setindex!(data::AbstractVector{<:Real}, v, va::ConstAccessor{T,true}) where T
    v == va.shape.value || throw(ArgumentError("Cannot set constant value to a different value"))
    data
end

function vs_setindex!(data::AbstractVector{<:Real}, v, va::ConstAccessor{T,false}) where T
    data
end

@inline _bcasted_view(data::AbstractArrayOfSimilarVectors{<:Real,N}, va::ConstAccessor) where N =
    Fill(va.shape.value, size(data)...)

Base.Broadcast.broadcasted(::typeof(getindex), A::AbstractArrayOfSimilarVectors{<:Real,N}, acc::Ref{<:ConstAccessor}) where N =
    _bcasted_view(A, acc[])

Base.Broadcast.broadcasted(::typeof(view), A::AbstractArrayOfSimilarVectors{<:Real,N}, acc::Ref{<:ConstAccessor}) where N =
    _bcasted_view(A, acc[])


text:
type TPTSystem{Tr <: ReferenceSystem, Tp <: Perturbation}
  ncomp::Int # number of components
  density::Float64 # total number density
  composition::Vector{Float64}
  temp::Float64 # temperature
  mass::Vector{Float64} # mass of components
  ref::Tr # reference system
  pert::Tp # perturbing system
end

function TPTSystem(ref::ReferenceSystem, pert::Perturbation = NullPerturbation(ncomp(ref)); T = InvalTemp, m = InvalMass)
  N::Int = ncomp(ref)
  ρ::Float64 = totalnumberdensity(ref)
  c::Vector{Float64} = composition(ref)

  temp::Float64 = T

  if temp == InvalTemp
    temp = temperature(ref) # may be also InvalTemp
  end

  mass::Vector{Float64} = []

  if typeof(m) <: Vector
    mass = m
  else m > 0
    mass = [m]
  end

  TPTSystem(N, ρ, c, temp, mass, ref, pert)
end

paircorrelation(sys::TPTSystem) = paircorrelation(sys.ref)
structurefactor(sys::TPTSystem) = structurefactor(sys.ref)
pairpotential(sys::TPTSystem) = pairpotential(sys.pert)
ncomp(sys::TPTSystem) = sys.ncomp
totalnumberdensity(sys::TPTSystem) = sys.density
composition(sys::TPTSystem) = sys.composition
temperature(sys::TPTSystem) = sys.temp
mass(sys::TPTSystem) = sys.mass

function kinetic(sys::TPTSystem)::Float64
  T = temperature(sys)
  K = 3/2 * kB*T
end

function entropy(sys::TPTSystem)::Float64
  T::Float64 = temperature(sys)

  S_gas = entropy_gas(sys)
  S_conf = entropy_conf(sys)
  S_ref = entropy(sys.ref) # reference system entropy
  S_pert = entropy(sys.pert, sys.ref, T) # perturbation entropy

  S = S_gas + S_conf + S_ref + S_pert
end

function entropy_gas(sys::TPTSystem)::Float64
  N::Int = ncomp(sys)
  ρ::Float64 = totalnumberdensity(sys)
  c::Vector{Float64} = composition(sys)
  T::Float64 = temperature(sys)
  m::Vector{Float64} = mass(sys)

  # Ideal gas entropy
  if m == InvalMass || T == InvalTemp
    return S = 0
  else
    P::Float64 = 1
    for i in 1:N
      if c[i] > 0
        P *= m[i]^c[i]
      end
    end
    return S = 5/2 + log((P * kB*T / 2π)^(3/2) / ρ)
  end
end

function entropy_conf(sys::TPTSystem)::Float64
  N::Int = ncomp(sys)
  c::Vector{Float64} = composition(sys)
  T::Float64 = temperature(sys)

  if T == InvalTemp
    return S = 0
  else
    S::Float64 = 0
    for i in 1:N
      if c[i] > 0
        S -= c[i]*log(c[i])
      end
    end
    return S
  end
end

function internal(sys::TPTSystem)::Float64
  U = internal(sys.pert, sys.ref)
end

function helmholtz(sys::TPTSystem)::Float64
  T = temperature(sys)
  K = kinetic(sys)
  U = internal(sys)
  S = entropy(sys)

  F = K + U - kB*T*S
end


text:
# space_promotions.jl

# Find the first space in a list of arguments that is not AnySpace
lcd(::Type{GSpace{T}}) where {T} = GSpace{T}
lcd(::Type{AnySpace}, ::Type{AnySpace}) = AnySpace
lcd(::Type{AnySpace}, ::Type{GSpace{T}}) where {T} = GSpace{T}
lcd(::Type{GSpace{T}}, ::Type{AnySpace}) where {T} = GSpace{T}
lcd(::Type{GSpace{T}}, ::Type{GSpace{S}}) where {T,S} = GSpace{T}
lcd(A::Type{GSpace{T}}, B::Type{GSpace{S}}, C::Type{GSpace{U}}, D...) where {T,S,U} = lcd(A, lcd(B, C, D...))

###############
# Isomorphisms
###############

"""
Return `True` if `A` and `B` are isomorphic. This implies that any point in `A`
can be converted to a point in `B`, and vice-versa, with an exact inverse.

The functionality is implemented recursively in terms of the
`isomorphism_reduction` function.
"""
isomorphism(::Type{GSpace{T}}, ::Type{GSpace{T}}) where {T} = True

isomorphism(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} =
    isomorphism_reduction_result(A, B, isomorphism_reduction(A,B)..., isomorphism_reduction(B,A)...)

"""
Indicate a reduction rule for isomorphism. The function either returns `()` or
a 2-tuple. If `isomorphism_reduction(A,B)` returns `(C,D)`, this means that `A`
and `B` are isomorphic if `C` and `D` are. Here, `C` and `D` are the subeltypes
of `A` and `B` respectively.

For example, `VectorSpace{2,T}` and `ComplexSpace{S}` are isomorphic if
`GeometricSpace{T}` and `GeometricSpace{S}` are isomorphic.

The isomorphism only has to be declared in one direction, and is automatically
valid in both directions. For space promotions, in case of ambiguity the first
space `A` here is chosen as the preferred promotion space.
"""
isomorphism_reduction(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} = ()

isomorphism_reduction_result(A, B) = False
isomorphism_reduction_result(A, B, C, D) = isomorphism(C, D)
isomorphism_reduction_result(A, B, C, D, E, F) = _duplicate_reduction_result(A, B, C, D, E, F)
# We have a situation here: isomorphism_reduction leads to a result in both
# directions.This is fine, as long as the outcome is the same. This actually
# happens for composite types, because they may differ only in their subtypes.
# Return an error if the results are not consistent, because that may never be
# the case.
_duplicate_reduction_result(A, B, C::Type{T}, D::Type{T}, E::Type{T}, F::Type{T}) where {T} = isomorphism(C, D)
_duplicate_reduction_result(A, B, C::Type{T}, D::Type{S}, E::Type{T}, F::Type{S}) where {T,S} = isomorphism(C, D)
_duplicate_reduction_result(A, B, C::Type{T}, D::Type{S}, E::Type{S}, F::Type{T}) where {T,S} = isomorphism(C, D)
_duplicate_reduction_result(A, B, C, D, E, F) = error("Duplicate isomorphism_reduction defined for ", A, " and ", B, ", but with inconsistent results.")

"""
Returns true if `A` is isomorphic to `B`, and false otherwise.
"""
isomorphic(a::GSpace, b::GSpace) = isomorphic(typeof(a), typeof(b))
isomorphic(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} = result(isomorphism(A, B))

"The symbol ≅ (\\cong) is a synonym for `isomorphic`."
≅ = isomorphic

#############
# Embeddings
#############

"""
The function `embedding` describes in the type domain whether a space with
type `T` is embedded into a space with type `B`. If so, it returns True,
otherwise it returns False.

Embeddings are the result of the following rules:
1) A{T} is embedded in A{S} if T promotes to S in the Julia type system
2) If a rule embedding_reduction(A,B) = (C,D) has been defined, then A is
   embedded in B if C is embedded in D.
3) Say A and B are isomorphic if C and D are. This automatically results in
   embedding reduction rules from A to B and from B to A.
4) A is embedded in B if the superspace of A is embedded in B.
"""
embedding(::Type{GeometricSpace{T}}, ::Type{GeometricSpace{T}}) where {T} = True

# Rules for AnySpace, necessary for termination of the recursive procedure using
# superspace in what follows:
# All spaces are embedded in AnySpace, but AnySpace is only embedded in itself
embedding(::Type{GeometricSpace{T}}, ::Type{AnySpace}) where {T} = True
embedding(::Type{AnySpace}, ::Type{AnySpace}) = True
embedding(::Type{AnySpace}, ::Type{GeometricSpace{T}}) where {T} = False

# We try all different cases. For isomorphisms, we have to check isomorphism_reduction
# with both possible orders of the arguments (A,B) and (B,A).
embedding(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} =
    one_of(
        embedding_via_promotion(A, B),
        embedding_via_reduction(A, B),
        embedding_via_superspace(A, B))

embedding_via_promotion(A, B) =
    _embedding_via_promotion(A, B, promote_type(eltype(A),eltype(B)))
_embedding_via_promotion(::Type{GSpace{T}}, ::Type{GSpace{S}}, ::Type{S}) where {T,S} = True
_embedding_via_promotion(::Type{GSpace{T}}, ::Type{GSpace{S}}, ::Type{U}) where {T,S,U} = False

# By default, two spaces have no embedding reduction rules.
# However, if there is an isomorphic_reduction rule, then there is automatically
# an embedding_reduction possible in both directions.
embedding_reduction(::Type{GSpace{T}}, ::Type{GSpace{T}}) where {T} = ()
embedding_reduction(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} =
   _embedding_reduction(A, B, isomorphism_reduction(A,B)..., 1, isomorphism_reduction(B,A)...)
_embedding_reduction(A, B, ::Int) = ()
_embedding_reduction(A, B, C, D, ::Int) = (C,D)
_embedding_reduction(A, B, ::Int, E, F) = (F,E)
_embedding_reduction(A, B, C, D, ::Int, E, F) = (C,D)

embedding_via_reduction(A, B) = _embedding_via_reduction(A, B, embedding_reduction(A, B)...)
# There was no specific rule for A and B
_embedding_via_reduction(A, B) = False
# A is embedded in B if C is embedded in D
_embedding_via_reduction(A, B, C, D) = embedding(C, D)

embedding_via_superspace(A, B) = embedding(superspace(A), B)

"""
Returns true if `A` is embedded in `B`.
"""
embedded(a::GSpace, b::GSpace) = embedded(typeof(a), typeof(b))
embedded(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} = result(embedding(A, B))

"The symbol ↪ (\\hookrightarrow) is a synonym for `embedded`."
↪ = embedded

#############
# Conversion
#############

"""
Convert the variable `x` to an element of the space `B`. This is possible if
the space of `x` is embedded in `B`.
"""
# We don't need to do anything if the space of x is B
convert_space(B::Type{GSpace{T}}, x::T) where {T} = x

# If it isn't, then we dispatch on the type of embedding
convert_space(B::Type{GSpace{T}}, x) where {T} = convert_spaces(x, spaceof(x), B)

convert_spaces(x, A, B) = convert_spaces(x, A, B,
    embedding_via_promotion(A, B),
    embedding_via_reduction(A, B),
    embedding_via_superspace(A, B))

convert_spaces(x, A, B, ::Type{True}, d2, d3) =
    convert_space_via_promotion(x, A, B)
convert_spaces(x, A, B, ::Type{False}, ::Type{True}, d3) =
    convert_space_via_reduction(x, A, B)
convert_spaces(x, A, B, ::Type{False}, ::Type{False}, ::Type{True}) =
    convert_space_via_superspace(x, A, B)
convert_spaces(x, A, B, ::Type{False}, ::Type{False}, ::Type{False}) =
    throw(InexactError(:convert_spaces, B, x))

# Embedding via promotion: promote the type of x using convert
convert_space_via_promotion(x, A, B) = convert(eltype(B), x)

# Embedding via reduction: if A is embedded in B if C is embedded in D,
# then we have to upgrade the space A{C} to A{D}
convert_space_via_reduction(x, A, B) = _convert_space_via_reduction(x, A, B, embedding_reduction(A, B)...)
_convert_space_via_reduction(x, A, B, C, D) = convert_space(B, convert_space(similar_space(A, eltype(D)), x))

# Embedding via superspace: we convert x first to the superspace of A and then continue
convert_space_via_superspace(x, A, B) = convert_space(B, convert_space(superspace(A), x))

## And now for the converse: restrict_space

"""
Restrict the variable `x` to an element of the space `B`. This is possible if
`B` is embedded in the space of `x`.

Mathematically, `restrict_space(A, y)` for `y` in space `B` is a left inverse
of `convert_space(B, x)` for `x` in space `A`. This means that
`restrict_space(A, convert_space(B, x)) == x` for any `x` in space `A`. However,
`restrict_space(A, y)` for a `y` not in the range of `convert_space(B, x)` could
take any value.
"""# We don't need to do anything if the space of x is B
restrict_space(B::Type{GSpace{T}}, x::T) where {T} = x

# If it isn't, then we dispatch on the type of embedding
restrict_space(B::Type{GSpace{T}}, x) where {T} = restrict_spaces1(x, spaceof(x), B)

# If A and B are isomorphic, then we can use convert_space rather than restrict_space.
restrict_spaces1(x, A, B) = _restrict_spaces1(x, A, B, isomorphism(A, B))
_restrict_spaces1(x, A, B, ::Type{True}) = convert_space(B, x)
_restrict_spaces1(x, A, B, ::Type{False}) = restrict_spaces2(x, A, B)

# We have to follow the inverse logic compared to the case of conversions.
# Here, B is embedded in A and we figure out by which path.
restrict_spaces2(x, A, B) = restrict_spaces2(x, A, B,
    embedding_via_promotion(B, A),
    embedding_via_reduction(B, A),
    embedding_via_superspace(B, A))

restrict_spaces2(x, A, B, ::Type{True}, d2, d3) =
    restrict_space_via_promotion(x, A, B)
restrict_spaces2(x, A, B, ::Type{False}, ::Type{True}, d3) =
    restrict_space_via_reduction(x, A, B)
restrict_spaces2(x, A, B, ::Type{False}, ::Type{False}, ::Type{True}) =
    restrict_space_via_superspace(x, A, B)
restrict_spaces2(x, A, B, ::Type{False}, ::Type{False}, ::Type{False}) =
    throw(InexactError(:restrict_spaces, B, x))

# Embedding via promotion
restrict_space_via_promotion(x, A, B) = demote(eltype(B), x)

"""
The function `demote(S, y::T)` is a left inverse of `convert(T, x::S)`, for the
case where type `S` promotes to type `T` (i.e. `promote_type(S,T) == T`).

This means that `demote(S, convert(T, x::S)) == x`, while
`demote(S, y)` may be anything for `y::T` not in the range of
`convert(T, x::S)`.

Note that Julia's `convert` function is its own left inverse, in the sense that
`convert(S, convert(T,x)) == x` usually holds if `x::S`. However, the `convert`
function generally throws an `InexactError()` for elements of type `T` that are
not in (or close to) the range of `convert(T, x::S)`. In those cases, `demote`
does not throw an error, but there may be an arbitrarily large difference between
`convert(T, demote(S, y))` and `y` itself. In other words, `demote` may differ
wildly from a right inverse of `convert`.
"""
demote(::Type{T}, x) where {T} = convert(T, x)

# We need this rule, because Julia only converts a complex number to a real
# number if the imaginary part is less than some threshold.
demote(::Type{T}, x::Complex{T}) where {T} = real(x)

# Embedding via reduction: if B is embedded in A if C is embedded in D,
# then we have to downgrade the space A{D} to A{C}
restrict_space_via_reduction(x, A, B) = _restrict_space_via_reduction(x, A, B, embedding_reduction(B, A)...)
_restrict_space_via_reduction(x, A, B, C, D) = restrict_space(B, restrict_space(similar_space(A, eltype(C)), x))

# Embedding via superspace: we convert x first to the superspace of A and then continue
restrict_space_via_superspace(x, A, B) = restrict_space(B, restrict_space(superspace(B), x))

#############
# Promotions
#############

# We mimick Julia's promotion system, with the key difference that the types
# a space can be promoted to are not actual supertypes, but superspaces.

# Implementation:
# - the user calls promote_space(a,b)
# - the space to promote to is decided by promote_space_type(spaceof(a),spaceof(b))
# - promote_space_type is implemented in terms of promote_space_rule. This function
#   checks whether one of the spaces is embedded into the other.
# - If that is not the case, then both spaces (one at a time) are upgraded to
#   their superspaces and promote_space_type is called again.
# - It ends as soon as a joined promotable type is found. If that is not the case,
#   all branches end in AnySpace.

"""
Promote the geometric spaces of the arguments to a joined space, to which all
arguments can be converted using embeddings. If no such concrete space exists,
this leads to `AnySpace`.
"""
promote_space() = ()
promote_space(x) = (x,)
promote_space(x::T, y::T) where {T} = (x,y)

promote_space(x, y) = _promote_space(x, y, promote_space_type(spaceof(x),spaceof(y)))
_promote_space(x, y, G) = (convert_space(G, x), convert_space(G, y))

promote_space_type(::Type{AnySpace}, ::Type{AnySpace}) = AnySpace
promote_space_type(::Type{AnySpace}, ::Type{GSpace{T}}) where {T} = AnySpace
promote_space_type(::Type{GSpace{T}}, ::Type{AnySpace}) where {T} = AnySpace
promote_space_type(::Type{GSpace{T}}, ::Type{GSpace{T}}) where {T} = GSpace{T}

promote_space_type(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} =
    lcd(
        promote_via_promotion(A,B),
        promote_via_embedding_reduction(A,B),
        promote_via_isomorphism_reduction(A,B),
        promote_via_superspace(A,B))

promote_via_promotion(A::Type{GSpace{T}}, B::Type{GSpace{S}}) where {T,S} =
    _promote_via_promotion(A, B, T, S, promote_type(T,S))

# S promotes to T: return A
_promote_via_promotion(A, B, ::Type{T}, ::Type{S}, ::Type{T}) where {T,S} = A
# T promotes to S: return B
_promote_via_promotion(A, B, ::Type{T}, ::Type{S}, ::Type{S}) where {T,S} = B
# T and S are equal: this should have been caught before!
_promote_via_promotion(A, B, ::Type{T}, ::Type{T}, ::Type{T}) where {T} = A
# neither is the case: return AnySpace
_promote_via_promotion(A, B, ::Type{T}, ::Type{S}, ::Type{U}) where {T,S,U} = AnySpace

promote_via_superspace(A, B) = lcd(promote_space_type(superspace(A),B), promote_space_type(A,superspace(B)))

promote_via_embedding_reduction(A, B) = _promote_via_embedding_reduction(A, B, embedding_reduction(A,B)..., 1, embedding_reduction(B,A)...)
# We passed a dummy Int to distinguish between the cases below
_promote_via_embedding_reduction(A, B, ::Int) = AnySpace
# This line: A is embedded in B via reduction if C is in D: we try to promote C and D
_promote_via_embedding_reduction(A, B, C, D, ::Int) = _promote_via_embedding_reduction1(A, B, C, D, promote_space_type(C,D))
# - Not succesful: we bail out
_promote_via_embedding_reduction1(A, B, C, D, E::Type{AnySpace}) = AnySpace
# - Successful: we promote B{D} to B{E}
_promote_via_embedding_reduction1(A, B, C, D, E::Type{GSpace{T}}) where {T} = similar_space(B, T)
# This case is similar, but with B embedded in A
_promote_via_embedding_reduction(A, B, ::Int, C, D) = _promote_via_embedding_reduction2(A, B, C, D, promote_space_type(C,D))
_promote_via_embedding_reduction2(A, B, C, D, E::Type{AnySpace}) = AnySpace
# - in case promotion of C and D is succesful, we promote A{C} to A{E}
_promote_via_embedding_reduction2(A, B, C, D, E::Type{GSpace{T}}) where {T} = similar_space(A, T)
# Third case: A is embedded in B and B is in A: we rely on isomorphism_reduction
# in order to get the right preferential promotion type.
# The only risk we take here is that A is embedded in B and B is embedded in A, but
# for some reason the isomorphism between A and B is not defined.
_promote_via_embedding_reduction(A, B, C, D, ::Int, E, F) = AnySpace

promote_via_isomorphism_reduction(A, B) = _promote_via_isomorphism_reduction(A, B, isomorphism_reduction(A,B)..., 1, isomorphism_reduction(B,A)...)
_promote_via_isomorphism_reduction(A, B, ::Int) = AnySpace
# Here, we can choose either to promote A or B. We choose A because it was listed first in the reduction rule
_promote_via_isomorphism_reduction(A, B, C, D, ::Int, args...) = similar_space(A, eltype(promote_space_type(C,D)))
# In this second case B was listed first, hence we continue promotion with B.
_promote_via_isomorphism_reduction(A, B, ::Int, C, D) = similar_space(B, eltype(promote_space_type(C,D)))


text:
const lowercaseletters = collect('a':'z')
const uppercaseletters = collect('A':'Z')
const letters = append!(collect('a':'z'),collect('A':'Z'))

seed(seed::Int=83)::Random.MersenneTwister = Random.seed!(seed)
random_int(;min::Int=0, max::Int=9999)::String = string(rand(min : max))
random_digit()::String = string(rand(0:9))
random_digit_not_null()::String= string(rand(1:9))
random_digit_or_empty()::String = string(rand(append!([' '],collect('0':'9'))))
random_digit_not_null_or_empty()::String = string(rand(append!([' '],collect('1':'9'))))
random_number(;digits::Int= 1 )::String = string(rand(0 : ((10^digits) - 1)))
random_letter()::String = string(rand(letters))

random_lowercase_letter()::String = string(rand(lowercaseletters))
random_uppercase_letter()::String = string(rand(uppercaseletters))
random_special_caracteres()::String = string(rand('!':'+'))

random_element(elements::Tuple=("a", "b", "b")) = rand(elements)

function randomize_nb_elements(;number::Int=10, le::Bool=false, ge::Bool=false)::String
    if le==true && ge==true
        return string(number)
    end
    min::Int=60
    max::Int=140

    if ge==true
        min = 100
    end

    if le==true
        max = 100
    end
    string(((rand(min : max))*number / 100) + 1 )
end

numerify( text::String = "####")::String = map( x -> x=='#' ? rand('0':'9') : x  , text )
lexify( text::String="????")::String = map( x -> x=='?' ? rand(letters) : x  , text )
bothify( text::String="## ??")::String = lexify(numerify(text))
unicode()::String = string(Char(rand(0x2200:0x22FF)))

generator(;items::Tuple=(user_name(),first_name(),city(),postcode(), user_name(),"@",domain_name()), separator=stringWithComa) = reduce(separator, items )


text:
using Test, Random, Flux

@testset "Weights" begin
    @testset "uniform" begin
        @test size(OperatorLearning.cglorot_uniform(128,64,10)) == (128, 64, 10)
        @test eltype(OperatorLearning.cglorot_uniform(128,64,10)) == ComplexF32
    end

    @testset "normal" begin
        @test size(OperatorLearning.cglorot_normal(128,64,10)) == (128, 64, 10)
        @test eltype(OperatorLearning.cglorot_normal(128,64,10)) == ComplexF32
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# Example of 2d linear elasticity + volume load + surface load

# Dictionary containing node coordinates
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])

# Create new problem of type `Elasticity`
block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = false
block.properties.geometric_stiffness = false

# Add volume element
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
update!(element, "displacement load 2", 576.0)
add_element!(block, element)

# Add boundary element for tractoin force
traction_element = Element(Seg2, (3, 4))
update!(traction_element, "geometry", X)
update!(traction_element, "displacement traction force 2", 288.0)
add_element!(block, traction_element)

# Define boundary conditions
bc = Problem(Dirichlet, "symmetry boundary conditions", 2, "displacement")
bc_elements = [Element(Seg2, (1, 2)), Element(Seg2, (4, 1))]
update!(bc_elements, "geometry", X)
update!(bc_elements[1], "displacement 2", 0.0)
update!(bc_elements[2], "displacement 1", 0.0)
add_elements!(bc, bc_elements)

# Create analysis, add problems to it and run analysis
analysis = Analysis(Linear)
add_problems!(analysis, block, bc)
run!(analysis)

# Analytical solution is known:
f = 288.0
g = 576.0
E = 288.0
nu = 1/3
u3 = block("displacement", 0.0)[3]
u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
@test isapprox(u3, u3_expected)


text:
using Symbolics, Test
@variables t, x(t), y(t), z(t)
@test Symbolics.islinear(x + y,[x,y])
@test Symbolics.islinear(x,[x,y])
@test Symbolics.islinear(y,[x,y])
@test !Symbolics.islinear(z,[x,y])

@test Symbolics.isaffine(x + y,[x,y])
@test Symbolics.isaffine(x,[x,y])
@test Symbolics.isaffine(y,[x,y])
@test Symbolics.isaffine(z,[x,y])
@test Symbolics.isaffine(x + y + z,[x,y])

@test Symbolics.isaffine(x + z * y,[x,y])
@test Symbolics.islinear(x + z * y,[x,y])
@test Symbolics.islinear(z * x + z * y,[x,y])
@test Symbolics.islinear(z * (x + y),[x,y])
@test Symbolics.isaffine(z * (x + y),[x,y])


text:
"""
    Update the clusters using a Rao-Blackwellized Gibbs algorithm,
    see Sudderth (2006, p. 94, algorithm 2.2)
"""
function update_z!(c::ChainState, s::SuffStats{A, B, C}) where {A, B, C}
    @unpack refresh_rate, α, z, τ, O, U, rng = c
    @unpack n = s
    randperm!(rng, τ)

    counter = 0
    @inbounds @fastmath for i ∈ τ
        # Recompute from scratch each 30 iterations
        (counter += 1) % refresh_rate == 0 && suffstats0!(s, c)

        # Initialize the search
        z̄ = z[i]
        p̄ = -Inf
        k̄ = first(U)

        # Update zi
        for k ∈ O
            p = log_pl(s, c, i, k) + log(n[k] - (z̄ == k))
            p = p - log(-log(rand(rng)))
            p < p̄ && continue
            z[i] = k
            p̄ = p
        end
        p = log_pl(s, c, i, k̄) + log(α[1])
        p = p - log(-log(rand(rng)))
        p > p̄ && (z[i] = k̄)

        # Update the sufficient statistics
        z̄ != z[i] && suffstats!(s, c, z̄, z[i], i)
    end
end

function update_θ!(c::ChainState, s::SuffStats{A, B, C}) where {A, B, C}
    @unpack J, u, ν, r, S, Σ, μ = s
    @unpack O, rng = c
    for k in O, j ∈ 1:J
        Σ[j, k] = rand(rng, InverseWishart(ν[j, k], S[j, k]))
        Σ[j, k] = (Σ[j, k] + Σ[j, k]') / 2
        μ[j, k] = rand(rng, MvNormal(u[j, k], Σ[j, k] / r[j, k]))
    end
end

function update_γ!(
    c::ChainState,
    s::SuffStats{A, B, C},
    pγ0::OffsetArray{Float64, 1, Vector{Float64}}
) where {A, B, C}
    @unpack n, J = s
    @unpack γ, O, rng = c

    # Resample γ[g], given the other γ's
    for g = 2:J
        # log-odds (numerator)
        γ[g] = 1
        suffstats0!(s, c)
        log_num = log(pγ0[sum(γ)])
        for k ∈ O, j ∈ (1, g)
            log_num += log_ml(s, j, k)
        end

        # log-odds (denominator)
        γ[g] = 0
        suffstats0!(s, c)
        log_den = log(pγ0[sum(γ)])
        for k ∈ O, j ∈ (1)
            log_den += log_ml(s, j, k)
        end

        # log-odds and new γ[g]
        log_odds = log_num - log_den
        γ[g] = rand(rng) <= 1 / (1 + exp(-log_odds))
    end
end

function update_α!(c::ChainState, s::SuffStats{A, B, C}) where {A, B, C}
    @unpack m, N = s
    @unpack a0, b0 = m
    @unpack α, K, rng = c

    ϕ = rand(rng, Beta(α[1] + 1.0, N))
    ψ = (a0 + K[1] - 1.0) / (N * (b0 - log(ϕ))); ψ = ψ / (1 + ψ)
    α[1] = rand(rng, Gamma(a0 + K[1] - (rand(rng) > ψ), 1.0 / (b0 - log(ϕ))))
end

function fit(
    m::Model{A, B},
    y::Vector{C},
    x::Vector{Int};
    K::Int = 5,
    iter::Int = 4000,
    warmup::Int = iter ÷ 2,
    refresh_rate::Int = length(x),
    rng::MersenneTwister = MersenneTwister()
) where {A, B, C}
    # Initialization
    s = SuffStats(m = m, y = y, x = x)
    c = ChainState(N = s.N, J = s.J, rng = rng, K = [K])
    pγ1 = zeros(2^(s.J - 1))
    pγ0 = ph0(s.J - 1, 1.0)

    # Updating
    @fastmath for t ∈ 1:iter
        suffstats0!(s, c)
        update_z!(c, s)
        update_α!(c, s)
        update_γ!(c, s, pγ0)
        t > warmup || continue
        pγ1[γcode(c.γ)] += 1
    end
    return pγ1 / (iter - warmup)
end

function fit(
    m::Model{A, B},
    y::Vector{C},
    x::Vector{Int},
    grid::Vector{Float64};
    K::Int = 5,
    iter::Int = 4000,
    warmup::Int = iter ÷ 2,
    rng::MersenneTwister = MersenneTwister()
) where {A, B, C}
    # Initialization
    D = length(y[1])
    J = length(unique(x))
    s = SuffStats(m = m, y = y, x = x)
    c = ChainState(N = s.N, J = s.J, rng = rng, K = [K])
    pγ1 = zeros(2^(s.J - 1))
    pγ0 = ph0(s.J - 1, 1.0)
    ygrid = Iterators.product(fill(grid, 2)...)
    ygrid = collect.(ygrid)[:] |> x -> hcat(x...)
    fgrid = [zeros(size(ygrid, 2)) for j = 1:J, z1 in 1:D, z2 in 1:D]
    M = size(ygrid, 2)

    # Updating
    @unpack γ, α, O = c
    @unpack n, N, μ, Σ = s
    @unpack r0, ν0 = m
    r1 = r0 + 1
    ν1 = ν0 + 1
    u0 = zeros(2)
    u1 = zeros(2)
    S0 = Cholesky(m.S0.factors[1:2, 1:2], :U, 0)
    S1 = deepcopy(S0)
    ggrid = zeros(M)
    yi = zeros(2)
    @fastmath for t ∈ 1:iter
        suffstats0!(s, c)
        update_z!(c, s)
        update_θ!(c, s)
        update_α!(c, s)
        update_γ!(c, s, pγ0)
        t > warmup || continue
        pγ1[γcode(c.γ)] += 1

        # Compute a base contribution to f(y)
        for i in 1:M
            S1.factors .= S0.factors
            yi .= ygrid[:, i]
            u1 .= (r0 * u0 + yi) / r1
            u1 .= (yi - u1) * √(r1 / r0)
            lowrankupdate!(S1, u1)
            ggrid[i] =
                0.5ν0 * logdet(S0) -
                0.5ν1 * logdet(S1) +
                logmvgamma(2, ν1 / 2) -
                logmvgamma(2, ν0 / 2) +
                0.5 * 2 * log(r0 / r1) -
                0.5 * 2 * log(π) * (r1 - r0)
            ggrid[i] = exp(ggrid[i])
        end

        # Accumulate f(y)
        μsub = [0.0; 0.0]
        Σsub = [1.0 0.0; 0.0 1.0]
        for j = 1:J, z1 = 1:D, z2 = z1+1:D
            for k in O
                μsub .= μ[j^γ[j], k][[z1; z2]]
                Σsub .= Σ[j^γ[j], k][[z1; z2], [z1; z2]]
                fgrid[j, z1, z2] += pdf(MvNormal(μsub, Σsub), ygrid) * n[k] / (N + α[1])
            end
            fgrid[j, z1, z2] += ggrid * α[1] / (N + α[1])
        end
    end

    # Convert fgrid into a dataframe
    dfs = [
        DataFrame(
            j    = j,
            var1 = z1,
            var2 = z2,
            y1   = ygrid[1, :],
            y2   = ygrid[2, :],
            f    = fgrid[j, z1, z2] / (iter - warmup)
        ) for j = 1:J, z1 = 1:D, z2 = 1:D
    ]
    df = reduce(vcat, dfs)
    filter!([:var1, :var2] => (x, y) -> x < y, df)
    return pγ1 / (iter - warmup), df
end

function train(
    y::Matrix{Float64},
    x::Vector{Int};
    r0::Int = 1,
    v0::Int = size(y, 2) + 2,
    u0::Vector{Float64} = zeros(size(y, 2)),
    S0::Matrix{Float64} = Matrix{Float64}(I(size(y, 2))),
    a0::Float64 = 1.0,
    b0::Float64 = 1.0,
    z0::Float64 = 1.0,
    iter::Int   = 4000,
    warmup::Int = iter ÷ 2,
    rng::MersenneTwister = MersenneTwister()
)
    N, D = size(y)
    J = length(unique(x))
    m = MANOVABNPTest.Model(
        D  = D,
        r0 = r0,
        ν0 = v0,
        u0 = u0,
        S0 = cholesky(S0),
        a0 = a0,
        b0 = b0,
        ζ0 = z0
    )
    y = standardize(ZScoreTransform, y, dims = 1)
    y = [SVector{D}(y[i, :]) for i ∈ 1:N]
    ps = MANOVABNPTest.fit(m, y, x; iter = iter, warmup = warmup, rng = rng)
    γs = [γvector(J, u)[2:end] for u in 1:length(ps)]
    DataFrame(hypothesis = γs, prob = ps)
end

function fit_γ(
    m::Model{A, B},
    y::Vector{C},
    x::Vector{Int};
    K::Int = 5,
    iter::Int = 4000,
    warmup::Int = iter ÷ 2,
    refresh_rate::Int = length(x),
    rng::MersenneTwister = MersenneTwister()
) where {A, B, C}
    # Initialization
    s = SuffStats(m = m, y = y, x = x)
    c = ChainState(N = s.N, J = s.J, rng = rng, K = [K])
    pγ1 = zeros(2^(s.J - 1))
    pγ0 = ph0(s.J - 1, 1.0)
    γchain = zeros(iter - warmup, s.J)

    # Updating
    @fastmath for t ∈ 1:iter
        suffstats0!(s, c)
        update_z!(c, s)
        update_α!(c, s)
        update_γ!(c, s, pγ0)
        t > warmup || continue
        γchain[t - warmup, :] = c.γ
    end
    return γchain
end

function train_γ(
    y::Matrix{Float64},
    x::Vector{Int};
    r0::Int = 1,
    v0::Int = size(y, 2) + 2,
    u0::Vector{Float64} = zeros(size(y, 2)),
    S0::Matrix{Float64} = Matrix{Float64}(I(size(y, 2))),
    a0::Float64 = 1.0,
    b0::Float64 = 1.0,
    z0::Float64 = 1.0,
    iter::Int   = 4000,
    warmup::Int = iter ÷ 2,
    rng::MersenneTwister = MersenneTwister()
)
    N, D = size(y)
    J = length(unique(x))
    m = MANOVABNPTest.Model(
        D  = D,
        r0 = r0,
        ν0 = v0,
        u0 = u0,
        S0 = cholesky(S0),
        a0 = a0,
        b0 = b0,
        ζ0 = z0
    )
    y = standardize(ZScoreTransform, y, dims = 1)
    y = [SVector{D}(y[i, :]) for i ∈ 1:N]
    MANOVABNPTest.fit_γ(m, y, x; iter = iter, warmup = warmup, rng = rng)
end


text:
function read_ace_header(filename)
    open(filename) do f
        idx = 1
        line = readline(f)
        if length(split(line)) == 3
            readline(f)
            line = readline(f)
        end
        HZ, AW, TZ, HD = split(line, " ", keepempty = false)
        suffix = HZ[end-2 : end]
        if suffix == "00c"
            T = 293.6
        elseif suffix == "01c"
            T = 600.0
        elseif suffix == "02c"
            T = 900.0
        elseif suffix == "03c"
            T = 1200.0
        elseif suffix == "04c"
            T = 2500.0
        elseif suffix == "05c"
            T = 0.1
        elseif suffix == "06c"
            T = 250.0
        end
        return HZ, T, parse(Float64, AW)
    end
end

xs_dir = "C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\Lib80x\\Lib80x"
#=
xs_paths = []
for (root, dirs, files) in walkdir(xs_dir)
    for filename in joinpath.(root, files)
        HZ, T, AW = read_ace_header(filename)
        push!(xs_paths, [HZ, AW, T, filename])
        #println(HZ)
    end
end
=#
ZAs = Int64[]
Ts = Float64[]
xs_paths = String[]
for (root, dirs, files) in walkdir(xs_dir)
    for filename in joinpath.(root, files)
        HZ, T, AW = read_ace_header(filename)
        ZA = parse(Int, HZ[1:end-4])
        push!(ZAs, ZA)
        push!(Ts, T)
        push!(xs_paths, filename)
        #println(HZ)
    end
end

order = sortperm(ZAs)
xs_paths = [ZAs[order] Ts[order] xs_paths[order]]

xs_dir_filename = "C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\Lib80x\\xsdir.txt"
using DelimitedFiles
writedlm(xs_dir_filename, xs_paths, ',')

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
# A geometric graph is one in which vertices are mapped to points and
# two vertices are adjacent if the distance between them is at most 1
# (or a set value).

export GeometricGraph, RandomGeometricGraph

"""
`GeometricGraph(A)` where `A` is an `m` by `n` matrix creates a
geometric graph in which the columns of `A` are the points
representing the vertices `1:n`. Two vertices are adjacent iff the
distance between the points is at most `1` (or an optional second
parameter, `d`).
"""
function GeometricGraph(A::Array{T,2}, d::Real = 1) where {T<:Real}
    r, n = size(A)
    G = IntGraph(n)
    dd::T = d * d

    for i = 1:n-1
        vi = A[:, i]
        for j = i+1:n
            vj = A[:, j]
            dv = vi - vj
            if dot(dv, dv) <= dd
                add!(G, i, j)
            end
        end
    end
    return G
end

"""
`GeometricGraph(f)` where `f` is a `Dict` mapping vertex names to
vectors creates a geometric graph in which two vertices are adjacent
iff distance between their points is at most `1` (or `d` if given as a
second argument).
"""
function GeometricGraph(f::Dict{S,Vector{T}}, d::Real = 1) where {S,T<:Real}
    G = SimpleGraph{S}()
    vtcs = collect(keys(f))
    for v in vtcs
        add!(G, v)
    end

    dd = d * d
    n = length(vtcs)
    for i = 1:n-1
        vi = f[vlist[i]]
        for j = i+1:n
            vj = f[vlist[j]]
            dv = vi - vj
            if dot(dv, dv) <= dd
                add!(G, vlist[i], vlist[j])
            end
        end
    end
    return G
end

"""
`RandomGeometricGraph(n::Int, dim::Int=2, d::Real=1)` creates a random
geometric graph by generating `n` points at random in the unit
`dim`-cube. Vertices are adjacent if their corresponding points are at
distance at most `d`.
"""
function RandomGeometricGraph(n::Int, dim::Int = 2, d::Real = 1)
    return GeometricGraph(rand(dim, n), d)
end


text:
export TrivialMeasure

struct TrivialMeasure <: AbstractMeasure end

sampletype(::TrivialMeasure) = Nothing


text:
# m5.7s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.29

df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
df[!, :neocortex_perc] = parse.(Float64, df[:, :neocortex_perc])
df[!, :lmass] = log.(df[:, :mass])
scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])

stan5_7 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] K; // Outcome
    vector[N] NC; // Predictor
    vector[N] M; // Predictor
}

parameters {
    real a; // Intercept
    real bM; // Slope (regression coefficients)
    real bN; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}

model {
    vector[N] mu;               // mu is a vector
    a ~ normal(0, 0.2);           //Priors
    bN ~ normal(0, 0.5);
    bM ~ normal(0, 0.5);
    sigma ~ exponential(1);
    mu = a + bM * M + bN * NC;
    K ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5_7s = SampleModel("m5.7", stan5_7);

# Input data for cmdstan

m5_7_data = Dict("N" => size(df, 1), "M" => df[!, :lmass_s],
    "K" => df[!, :kcal_per_g_s], "NC" => df[!, :neocortex_perc_s]);

# Sample using StanSample

rc5_7s = stan_sample(m5_7s, data=m5_7_data);

if success(rc5_7s)

    part5_7s = read_samples(m5_7s, :particles)
    part5_7s |> display

    rethinking = "
             mean   sd  5.5% 94.5%
    a      0.07 0.13 -0.15  0.28
    bN     0.68 0.25  0.28  1.07
    bM    -0.70 0.22 -1.06 -0.35
    sigma  0.74 0.13  0.53  0.95
    "

end


text:
myfunc(x) = log(sum(x))
test_compare(myfunc; x=rand(2))

quote
    tmp685 = @get_or_create(mem, :tmp685, zero(Float64))
    tmp686 = @get_or_create(mem, :tmp686, zero(Float64))
    dtmp686!dx = @get_or_create(mem, :dtmp686!dx, zeros(Float64, (2,)))
    dtmp686!dtmp686 = @get_or_create(mem, :dtmp686!dtmp686, zero(Float64))
    tmp694 = @get_or_create(mem, :tmp694, zero(Float64))
    dtmp686!dtmp686 = 1.0
    tmp685 = sum(x)
    dtmp686!dx .= dtmp686!dtmp686 .* (dtmp686!dtmp686 ./ tmp685)
    tmp686 = log(tmp685)
    tmp696 = (tmp686, dtmp686!dx)
end


text:
# Autogenerated wrapper script for glslang_jll for i686-linux-gnu-cxx03
export glslangValidator, slibGenericCodeGen, slibHLSL, slibMachineIndependent, slibOGLCompiler, slibOSDependent, slibglslang

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `glslangValidator`
const glslangValidator_splitpath = ["bin", "glslangValidator"]

# This will be filled out by __init__() for all products, as it must be done at runtime
glslangValidator_path = ""

# glslangValidator-specific global declaration
function glslangValidator(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(glslangValidator_path)
    end
end

# Relative path to `slibGenericCodeGen`
const slibGenericCodeGen_splitpath = ["lib", "libGenericCodeGen.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibGenericCodeGen_path = ""

# slibGenericCodeGen-specific global declaration
# This will be filled out by __init__()
slibGenericCodeGen = ""

# Relative path to `slibHLSL`
const slibHLSL_splitpath = ["lib", "libHLSL.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibHLSL_path = ""

# slibHLSL-specific global declaration
# This will be filled out by __init__()
slibHLSL = ""

# Relative path to `slibMachineIndependent`
const slibMachineIndependent_splitpath = ["lib", "libMachineIndependent.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibMachineIndependent_path = ""

# slibMachineIndependent-specific global declaration
# This will be filled out by __init__()
slibMachineIndependent = ""

# Relative path to `slibOGLCompiler`
const slibOGLCompiler_splitpath = ["lib", "libOGLCompiler.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibOGLCompiler_path = ""

# slibOGLCompiler-specific global declaration
# This will be filled out by __init__()
slibOGLCompiler = ""

# Relative path to `slibOSDependent`
const slibOSDependent_splitpath = ["lib", "libOSDependent.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibOSDependent_path = ""

# slibOSDependent-specific global declaration
# This will be filled out by __init__()
slibOSDependent = ""

# Relative path to `slibglslang`
const slibglslang_splitpath = ["lib", "libglslang.a"]

# This will be filled out by __init__() for all products, as it must be done at runtime
slibglslang_path = ""

# slibglslang-specific global declaration
# This will be filled out by __init__()
slibglslang = ""

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"glslang")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global glslangValidator_path = normpath(joinpath(artifact_dir, glslangValidator_splitpath...))

    push!(PATH_list, dirname(glslangValidator_path))
    global slibGenericCodeGen_path = normpath(joinpath(artifact_dir, slibGenericCodeGen_splitpath...))

    global slibGenericCodeGen = slibGenericCodeGen_path
    global slibHLSL_path = normpath(joinpath(artifact_dir, slibHLSL_splitpath...))

    global slibHLSL = slibHLSL_path
    global slibMachineIndependent_path = normpath(joinpath(artifact_dir, slibMachineIndependent_splitpath...))

    global slibMachineIndependent = slibMachineIndependent_path
    global slibOGLCompiler_path = normpath(joinpath(artifact_dir, slibOGLCompiler_splitpath...))

    global slibOGLCompiler = slibOGLCompiler_path
    global slibOSDependent_path = normpath(joinpath(artifact_dir, slibOSDependent_splitpath...))

    global slibOSDependent = slibOSDependent_path
    global slibglslang_path = normpath(joinpath(artifact_dir, slibglslang_splitpath...))

    global slibglslang = slibglslang_path
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')
end  # __init__()


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

# import ..WaveOptics

# -------------------------------------------------------------------------
#               WIDE FIELD MICROSCOPE APSF (DEFOCUS MODEL)
# -----------------------------------------------------------------------------

mutable struct DefocusModel <: LateralShiftInvariantModel
    recipe::OpticalRecipe
    sim_size_m::Float64
    sim_size_px::Int64

    padding::Int32
    num_observations::Int32

    # Private model state information
    aperture_plane_coordinates::OpticalCoordinates
    microscope_atf::Array{ComplexF32,2}
    user_intensity_scaling_factor::Float32
end

"""
This code implements compute a microscope APSF by simulating the refocusing of
a microscope with a defocus mask in the back aperture plane [1].  We then
Fourier transforming the back aperture ATF to obtain the image plane APSF.  This
is more efficient than carrying out numerical integration, as in the Debye model.
However, this model does utilize the Debye model to produce a baseline ATF, on top
of which defocus is added.

This approach is faster than the Debye models below, although it is not quite as
accurate and it can be prone to sampling artifacts.  Generally this model produces
results that are within 1% (and often better) of the Debye model.  However you should
always compare the two models to be sure.

References:

[1] Hanser, B. M., Gustafsson, M., & Agard, D. A. (2004). Phase-retrieved pupil
functions in wide-field fluorescence microscopy. Journal of Microscopy. Vol.
216. pp. 32-48. 2004.

**Parameters**

      recipe - an optical recipe containing important simulation parameters
      sim_size_m - linear dimension (used for width & height) in meters (or other units
               of your choice) to be used for the PSF simulation.
      sim_size_px - the number of pixel in each linear dimension.  All simulations are square.
"""
function DefocusModel(recipe::OpticalRecipe, sim_size_m::Float64, sim_size_px::Int64)

    # TODO: Add a check to see if sim_size_m and sim_size_px are sufficient to produce
    # an accurate simulation.

    num_observations = 1

    # We will zero pad the FFT to increase the sampling rate in k-space (i.e. the back aperture)
    # and ensure that there are periodic edge wrapping artifacts in the final PSF.  (Multiplication
    # in the Fourier Domain corresponds to convolution in the primal domain, and in the case of the
    # FFT that means the circular convolution due to periodicity of the discrete fourier transform.)
    padding = 2
    f_objective = recipe.f_tubelens / recipe.objective_mag

    aperture_plane_coordinates = fourier_plane_meshgrid(
        sim_size_px * padding,
        sim_size_m / sim_size_px,
        recipe.objective_na,
        recipe.wavelength,
        f_objective,
        NormalizedUnits(),
    )

    # Use the DebyeModel to accurately compute the ATF for the microscope.  We use this in
    # conjunction with the FFT-based defocus calculation below, making the overall computation
    # of the PSf more accurate.
    debye_model = DebyeModel(recipe, sim_size_m, sim_size_px)
    apsf_z0 = apsf(debye_model, (0.0, 0.0, 0.0))
    if recipe.normalize_intensity
        apsf_z0 *= sqrt(1.0 / maximum(abs2(apsf_z0)))
    end
    microscope_atf = im_to_ap_xform(WaveOptics.fftpad(apsf_z0))

    # If a static aperture mask has been supplied, include it in the model here.
    if !isa(recipe.static_aperture_mask, NullCodedAperture)
        static_ap_atf = mask(recipe.static_aperture_mask, aperture_plane_coordinates)
        microscope_atf = microscope_atf .* static_ap_atf
    end

    # Intensities can be scaled further by the user by setting this value explicitly.
    user_intensity_scaling_factor = 1.0

    return DefocusModel(
        recipe,
        sim_size_m,
        sim_size_px,
        padding,
        num_observations,
        aperture_plane_coordinates,
        microscope_atf,
        user_intensity_scaling_factor,
    )
end

function na_limit!(
    model::DefocusModel,
    pupil::AbstractArray{Complex{T}},
) where {T<:AbstractFloat}
    return (model.aperture_plane_coordinates.rho .<= 1.0) .* pupil
end

function na_limit!(
    rho::AbstractArray{T,2},
    pupil::AbstractArray{Complex{T}},
) where {T<:AbstractFloat}
    return (rho .<= 1.0) .* pupil
end

function num_parameters(model::DefocusModel)
    return num_parameters(model.recipe.dynamic_aperture_mask)
end

"""
Set the total number of photons for an in-focus PSF of an equivalent wide field microsocpe
model to be equal to the specified number of total_photons.  This helps to remove the effects
of simulation pixel size and can be used to achive similar intensity scaling across
very different PSF models.  This is especially useful for generating consistent Fisher
information results.
"""
function find_intensity_scaling_factor(
    model::PointSpreadFunctionModel,
    total_photons::T,
) where {T<:Real}

    # Generate an equivalent widefield model
    recipe = model.recipe
    wf_recipe = WaveOptics.WideFieldOpticalRecipe(
        recipe.objective_mag,
        recipe.objective_na,
        recipe.wavelength,
        recipe.medium_index,
        recipe.f_tubelens,
    )
    wf_model = WaveOptics.DefocusModel(wf_recipe, model.sim_size_m, model.sim_size_px)

    # Find the areaof a simulation pixel in um^2
    #sim_pixel_area_um_squared = (model.sim_size_m / model.sim_size_px / 1000)^2
    #area_scaling_factor = 1. / sim_pixel_area_um_squared

    # Find the max intensity without scaling
    WaveOptics.set_intensity_scaling_factor!(wf_model, 1.0)
    psf0 = WaveOptics.psf(wf_model, (0, 0, 0))

    #    alpha = photons_per_um_squared / (maximum(psf0) * area_scaling_factor)
    alpha = total_photons / (sum(psf0))# * area_scaling_factor)

    # Check the result
    WaveOptics.set_intensity_scaling_factor!(wf_model, alpha)
    psf1 = WaveOptics.psf(wf_model, (0, 0, 0))

    println("Adjusting intensity scaling factor...")
    println("   Widefield PSF sum before: $(sum(psf0))")
    println("   After: $(sum(psf1))")
    println("   Intensity scaling factor: $alpha.")
    return alpha
end

"""
Set a scalar multiplier that is used to adjust the final intensity pattern generate
by the PSF model.  This multiplier can be found using find_intensity_scaling_factor()
above.
"""
function set_intensity_scaling_factor!(model::DefocusModel, f::Float64)
    model.user_intensity_scaling_factor = f
end

function atf(
    model::DefocusModel,
    p::NTuple{3,T},
    theta = Float64[];
    na_limit::Bool = true,
    apply_tilt::Bool = true,
) where {T<:Real}
    x, y, z = p  # Unpack p

    # Generate a defocus ATF at the back aperture plane. This is a more
    # direct and accurate method of computing the defocus ATF than the
    # spherical wavefront and FFT method above.
    defocus_atf = defocus_mask(
        z,
        model.aperture_plane_coordinates.rho,
        model.aperture_plane_coordinates.theta,
        model.recipe.objective_na,
        model.recipe.wavelength,
        model.recipe.medium_index,
    )

    # Compute back aperture simulation size directly
    fourier_size_m =
        model.sim_size_px / model.sim_size_m *
        model.recipe.f_tubelens *
        model.recipe.wavelength

    # Compute the apodization appropriate to account for an objective with the abbe
    # sine correction. This would be appropriate for a plan corrected objective.
    #        microscope_atf = abbe_sine_apodization(model.aperture_plane_coordinates.rho, model.aperture_plane_coordinates.theta, model.recipe.objective_na, model.recipe.medium_index)
    intensity = model.user_intensity_scaling_factor
    microscope_atf = sqrt(intensity) .* model.microscope_atf  # Use Debye model for microscope ATF

    # If the user has supplied an aperture code, we add it to the PSF here.
    if apply_tilt && ((x != 0.0) || (y != 0.0))
        translation_atf = translation_mask(
            x,
            y,
            model.aperture_plane_coordinates.fx,
            model.aperture_plane_coordinates.fy,
            model.recipe.objective_na,
            model.recipe.wavelength,
        )
        full_atf = microscope_atf .* defocus_atf .* translation_atf
    else
        full_atf = microscope_atf .* defocus_atf
    end

    if !isa(model.recipe.dynamic_aperture_mask, NullCodedAperture)
        ap_atf = mask(
            model.recipe.dynamic_aperture_mask,
            model.aperture_plane_coordinates,
            theta,
        )
        full_atf .= full_atf .* ap_atf
    end

    if na_limit
        return na_limit!(model, full_atf)
    else
        return full_atf
    end
end

function otf_2d(model::DefocusModel, p, theta = Float64[]; return_freqs = false)
    full_atf = atf(model, p, theta, apply_tilt = true, na_limit = true)

    # Compute the normalized autocorrelation using an FFT method.
    A = fftshift(fft(fftshift(fftpad(full_atf))))
    B = fftunpad(fftshift(ifft(fftshift(A .* conj(A)))))

    normalization = sum(full_atf .* conj(full_atf))
    norm_autocorr = B ./ normalization

    if return_freqs
        sample_period_x = model.sim_size_m / (model.sim_size_px * model.padding)
        freqs = fftshift(fftfreq(model.sim_size_px * model.padding, sample_period_x))
        return norm_autocorr, freqs
    else
        return norm_autocorr
    end
end

"""
Compute the amplitude point spread function.

**Parameters**

    p - position of the point source

    theta - the parameters to pass to the aperture_mask.  If theta = 'None', the aperture's default
            parameters will be used.

**Returns**

    apsf - 2D image of the APSF on the image sensor
"""
function apsf(model::DefocusModel, p, theta = Float64[])
    x, y, z = p  # Unpack p

    # Don't apply tilt in the back aperture, as we will opt to instead apply
    # this as a lateral translation in the image plane.
    full_atf = atf(model, p, theta, na_limit = false, apply_tilt = false)

    # Transform back to get the field at the image plane. Restrict to
    # circular aperture using the na_limit() function.
    fourier_size_m =
        model.sim_size_px / model.sim_size_m *
        model.recipe.f_tubelens *
        model.recipe.wavelength
    camera_size_m, apsf = ap_to_im(
        na_limit!(model, full_atf),
        fourier_size_m,
        model.recipe.f_tubelens,
        model.recipe.wavelength,
    )
    apsf = fftshift(apsf)

    # Unpad the FFT, and return the full result. The translate_apsf() function
    # here translates the point in x and y, returning an interpolated result.
    if x != 0.0 || y != 0.0
        return translate_apsf(
            model.sim_size_m,
            model.sim_size_px,
            WaveOptics.fftunpad(apsf),
            x,
            y,
        )
    else
        return WaveOptics.fftunpad(apsf)
    end
end

function psf_3d(model::DefocusModel, zrange::AbstractArray{T,1}) where {T<:Real}
    # Transform back to get the field at the image plane. Restrict to
    # circular aperture using the na_limit() function.
    fourier_size_m =
        model.sim_size_px / model.sim_size_m *
        model.recipe.f_tubelens *
        model.recipe.wavelength

    c_rho = model.aperture_plane_coordinates.rho
    c_theta = model.aperture_plane_coordinates.theta
    microscope_atf = model.microscope_atf
    c_rho = ThreeDeconv.to_gpu_or_not_to_gpu(c_rho)
    c_theta = ThreeDeconv.to_gpu_or_not_to_gpu(c_theta)
    microscope_atf = ThreeDeconv.to_gpu_or_not_to_gpu(microscope_atf)

    M, N = size(c_rho)

    # Compute back aperture simulation size directly
    fourier_size_m =
        model.sim_size_px / model.sim_size_m *
        model.recipe.f_tubelens *
        model.recipe.wavelength

    # Compute the apodization appropriate to account for an objective with the abbe
    # sine correction. This would be appropriate for a plan corrected objective.
    #        microscope_atf = abbe_sine_apodization(model.aperture_plane_coordinates.rho, model.aperture_plane_coordinates.theta, model.recipe.objective_na, model.recipe.medium_index)
    intensity = model.user_intensity_scaling_factor
    microscope_atf = sqrt(intensity) .* microscope_atf  # Use Debye model for microscope ATF

    psf3d = zeros(Float32, (M, N, length(zrange)))
    for (idx, z) in enumerate(zrange)
        # Don't apply tilt in the back aperture, as we will opt to instead apply
        # this as a lateral translation in the image plane.

        # Generate a defocus ATF at the back aperture plane. This is a more
        # direct and accurate method of computing the defocus ATF than the
        # spherical wavefront and FFT method above.
        defocus_atf = defocus_mask(
            z,
            c_rho,
            c_theta,
            model.recipe.objective_na,
            model.recipe.wavelength,
            model.recipe.medium_index,
        )

        # If the user has supplied an aperture code, we add it to the PSF here.
        full_atf = microscope_atf .* defocus_atf
        camera_size_m, apsf = ap_to_im(
            na_limit!(c_rho, full_atf),
            fourier_size_m,
            model.recipe.f_tubelens,
            model.recipe.wavelength,
        )
        psf3d[:, :, idx] .= collect(abs2.(apsf))
    end

    psf3d = fftshift(psf3d, [1, 2])

    # Unpad the FFT, and return the full result. The translate_apsf() function
    # here translates the point in x and y, returning an interpolated result.
    return WaveOptics.fftunpad2(psf3d)
end


text:
#-------------julia-fractal.jl-------------------------------------------------#
#
# Purpose: This file tests the creation of the julia fractal in julia for fun
#
#------------------------------------------------------------------------------#

# We will need a few libraries for visualization
using PyPlot

# function to normalize fractal
function normalize(carr::Array{Float64,2})
    normalization_factor = sum(carr) / length(carr)
    carr /= normalization_factor
    return carr
end

# Function to write a fractal to a file (assuming square data)
function write_fractal(filename::String, res::Int64, fractal::Array{Float64})
    file = open(filename, "w")
    for i = 1:res
        line = ""
        for j = 1:res
            if j != res
                line *= string(fractal[i,j]) * '\t'
            else
                line *= string(fractal[i,j]) * '\n'
            end
        end
        write(file, line)
    end
    close(file)
end

# Function to calculate and output julia fractal to image
function create_fractal(cutoff::Float64, step::Float64, c_val::Float64,
                        res::Int64, range::Float64)
    # Setting up the array to be used
    carr = zeros(res,res)

    # For now, let's just assume that the imaginary and real ranges are from
    #     -range < x < range
    for i = 1:res
        for j = 1:res
            z = complex((i/res)*range - range/2.0,(j/res)*range - range/2.0)
            c = complex(0,c_val)
            n = 255
            while (abs(z) < cutoff && n > 5)
                z = z*z + c
                n -= 5
            end
            carr[i,j] = n/255.0;
        end
    end

    return normalize(carr)

end

# Fuction to scan through C values
function c_scan(cutoff::Float64, step::Float64, res::Int64,
                range::Float64, max_c::Float64, c_step::Float64)

    id = 0
    for i = 0:c_step:max_c
        println(i)
        carr = create_fractal(cutoff, step, i, res, range)
        filename = string("cdata/c_scan",
                          lpad(string(id), 5, string(0)), ".dat")
        write_fractal(filename, res, carr)
        id += 1
    end
end

# Function to zoom in on fractal
function fractal_zoom(cutoff::Float64, step::Float64, c_val::Float64,
                      res::Int64, max_range::Float64,
                      min_range::Float64)
    id = 0
    range = max_range
    while range > min_range
        println(range)
        carr = create_fractal(cutoff, step, c_val, res, range)
        filename = string("data/fractal_zoom",
                          lpad(string(id), 5, string(0)), ".dat")
        write_fractal(filename, res, carr)
        id += 1
        range -= range*0.05
    end
end

#fractal_zoom(10.0, 0.1, 1.0, 512, 1.0, 0.0001)

c_scan(10.,5., 512, 3.0, 1.0, 0.025)


text:
struct Closereason
    code::Int
    description::String
    valid::Bool
    function Closereason(code::Int, description::String = "")
        if length(description) === 0 && haskey(CLOSE_DESCRIPTIONS, code)
            description = CLOSE_DESCRIPTIONS[code]
        end

        parsed = parsedescription(description)
        parsed[2] !== nothing && (code = parsed[2])
        description = parsed[1]
        new(
            code,
            description,
            validateCloseReason(code),
        )
    end
end
function parsedescription(description::String)::Array{Union{String, Int, Nothing}, 1}
    reg = r"^\[(\d+)*?\]"
    code = nothing
    description = replace(description, reg => newcode -> (
        begin
            result = match(reg, newcode).captures[1]
            result !== nothing && (code = parse(Int, result))
            s""
        end
    ))
    [description, code]
end
function validateCloseReason(code::Int)
    if code < 1000
        #Status codes in the range 0-999 are not used
        return false;
    end
    if code >= 1000 && code <= 2999
        #Codes from 1000 - 2999 are reserved for use by the protocol.  Only
        #a few codes are defined, all others are currently illegal.
        return code in [1000, 1001, 1002, 1003, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015]
    end
    if code >= 3000 && code <= 3999
        #Reserved for use by libraries, frameworks, and applications.
        #Should be registered with IANA.  Interpretation of these codes is
        #undefined by the WebSocket protocol.
        return true;
    end
    if code >= 4000 && code <= 4999
        #Reserved for private use.  Interpretation of these codes is
        #undefined by the WebSocket protocol.
        return true;
    end
    if code >= 5000
        return false;
    end
end
"""
    abstract type WebsocketError <: Exception
Websocket child error types have the following fields:
- msg::String
- log::Function

The `log()` function will internally call: @error msg  exception = (err, trace)

Where `trace` is the backtrace of the exception origin.
"""
abstract type WebsocketError <: Exception end
function msg(err::Exception)
    try
        return err.msg
    catch
        return string(typeof(err))
    end
    # hasfield(typeof(err), :msg) ? err.msg : string(typeof(err)) #hasfield needs julia >= 1.2
end
logError(self::Exception, err::Exception, trace::Array) = @error string(typeof(self))*"(\"$(self.msg)\")" exception = (err, trace)
"""
    struct ConnectError <: WebsocketError
An exception originated while trying to start a server or connect to a server
"""
struct ConnectError <: WebsocketError
    msg::String
    log::Function
    function ConnectError(err::Exception, trace::Array = [])
        self = new(
            msg(err),
            () -> logError(self, err, trace)
        )
    end
end
"""
    struct CallbackError <: WebsocketError
An exception originated in a user provided callback function
"""
struct CallbackError <: WebsocketError
    msg::String
    log::Function
    function CallbackError(err::Exception, trace::Array = [])
        self = new(
            msg(err),
            () -> logError(self, err, trace)
        )
    end
end
"""
    struct FrameError <: WebsocketError
An exception originated while parsing a websocket data frame
"""
struct FrameError <: WebsocketError
    msg::String
    log::Function
    function FrameError(err::Exception, trace::Array = [])
        self = new(
            msg(err),
            () -> logError(self, err, trace)
        )
    end
end

requestHash() = base64encode(rand(UInt8, 16))
function acceptHash(key::String)
    hashkey = "$(key)258EAFA5-E914-47DA-95CA-C5AB0DC85B11"
    base64encode(digest(MD_SHA1, hashkey))
end

textbuffer(data::String)::Array{UInt8,1} = Array{UInt8,1}(data)
textbuffer(data::Number)::Array{UInt8,1} = Array{UInt8,1}(string(data))

modindex(i::Int, m::Int) = ((i-1) % m) + 1
function mask!(mask::IOBuffer, data::Array{UInt8, 1})
    for (i, value) in enumerate(data)
       data[i] = value ⊻ mask.data[modindex(i, mask.size)]
    end
end

function makeHeaders(extend::Dict{String, String})
    headers = Dict{String, String}(
        "Sec-WebSocket-Version" => defaultHeaders["Sec-WebSocket-Version"],
    )
    for (key, value) in extend
        headers[key] = value
    end
    headers["Upgrade"] = defaultHeaders["Upgrade"]
    headers["Connection"] = defaultHeaders["Connection"]
    headers["Sec-WebSocket-Key"] = requestHash()

    headers
end

function getmaxwithheaders(config::NamedTuple)
    !config.fragmentOutgoingMessages && return (Int(config.maxReceivedMessageSize) + 10)
    Int(config.maxReceivedMessageSize + (ceil(config.maxReceivedMessageSize / config.fragmentationThreshold)*10))
end

function explain(object::Any)
    type = typeof(object)
    @show type
    println("----------------------------------")
    for field in fieldnames(type)
        value = getfield(object, field)
        println(field, " | ", typeof(value), " | ", value)
    end
    println("----------------------------------")
end

max_stars_repo_name: citkane/Websocket

max_issues_repo_name: citkane/Websocket

max_forks_repo_name: citkane/Websocket


text:
""" Extends the matrices u,v,η,sst with a halo of ghost points for boundary conditions."""
function add_halo(  u::Array{T,2},
                    v::Array{T,2},
                    η::Array{T,2},
                    sst::Array{T,2},
                    S::ModelSetup) where {T<:AbstractFloat}

    @unpack nx,ny,nux,nuy,nvx,nvy = S.grid
    @unpack halo,haloη,halosstx,halossty = S.grid

    # Add zeros to satisfy kinematic boundary conditions
    u = cat(zeros(T,halo,nuy),u,zeros(T,halo,nuy),dims=1)
    u = cat(zeros(T,nux+2*halo,halo),u,zeros(T,nux+2*halo,halo),dims=2)

    v = cat(zeros(T,halo,nvy),v,zeros(T,halo,nvy),dims=1)
    v = cat(zeros(T,nvx+2*halo,halo),v,zeros(T,nvx+2*halo,halo),dims=2)

    η = cat(zeros(T,haloη,ny),η,zeros(T,haloη,ny),dims=1)
    η = cat(zeros(T,nx+2*haloη,haloη),η,zeros(T,nx+2*haloη,haloη),dims=2)

    sst = cat(zeros(T,halosstx,ny),sst,zeros(T,halosstx,ny),dims=1)
    sst = cat(zeros(T,nx+2*halosstx,halossty),sst,zeros(T,nx+2*halosstx,halossty),dims=2)

    ghost_points!(u,v,η,S)
    ghost_points_sst!(sst,S)
    return u,v,η,sst
end

"""Cut off the halo from the prognostic variables."""
function remove_halo(   u::Array{T,2},
                        v::Array{T,2},
                        η::Array{T,2},
                        sst::Array{T,2},
                        S::ModelSetup) where {T<:AbstractFloat}

    @unpack halo,haloη,halosstx,halossty = S.grid

    uview = u[halo+1:end-halo,halo+1:end-halo]
    vview = v[halo+1:end-halo,halo+1:end-halo]
    ηview = η[haloη+1:end-haloη,haloη+1:end-haloη]
    sstview = sst[halosstx+1:end-halosstx,halossty+1:end-halossty]

    return uview,vview,ηview,sstview
end

""" Copy ghost points for u from inside to the halo in the nonperiodic case. """
function ghost_points_u_nonperiodic!(C::Constants,u::AbstractMatrix)

    @unpack one_minus_α = C

    # tangential boundary condition
    @views @inbounds u[:,1] .= one_minus_α*u[:,4]
    @views @inbounds u[:,2] .= one_minus_α*u[:,3]
    @views @inbounds u[:,end-1] .= one_minus_α*u[:,end-2]
    @views @inbounds u[:,end] .= one_minus_α*u[:,end-3]
end

""" Copy ghost points for u from inside to the halo in the periodic case. """
function ghost_points_u_periodic!(::Type{Tcomm},u::Array{T,2},C::Constants) where {Tcomm,T}

    @unpack one_minus_α = C

    # periodic bc   - mimic MPI communication with (reduced precision) type Tcomm
    @views @inbounds u[1,:] .= T.(Tcomm.(u[end-3,:]))
    @views @inbounds u[2,:] .= T.(Tcomm.(u[end-2,:]))
    @views @inbounds u[end-1,:] .= T.(Tcomm.(u[3,:]))
    @views @inbounds u[end,:] .= T.(Tcomm.(u[4,:]))

    # tangential bc
    @views @inbounds u[:,1] .= one_minus_α*u[:,4]
    @views @inbounds u[:,2] .= one_minus_α*u[:,3]
    @views @inbounds u[:,end-1] .= one_minus_α*u[:,end-2]
    @views @inbounds u[:,end] .= one_minus_α*u[:,end-3]

end

""" Copy ghost points for v from inside to the halo in the nonperiodic case. """
function ghost_points_v_nonperiodic!(C::Constants,v::AbstractMatrix)

    @unpack one_minus_α = C

    # tangential boundary condition
    @views @inbounds v[1,:] .= one_minus_α*v[4,:]
    @views @inbounds v[2,:] .= one_minus_α*v[3,:]
    @views @inbounds v[end-1,:] .= one_minus_α*v[end-2,:]
    @views @inbounds v[end,:] .= one_minus_α*v[end-3,:]
end

""" Copy ghost points for v from inside to the halo in the periodic case. """
function ghost_points_v_periodic!(::Type{Tcomm},v::Array{T,2}) where {Tcomm,T}

    # tangential boundary condition - mimic MPI communication with (reduced precision) type Tcomm
    @views @inbounds v[1,:] .= T.(Tcomm.(v[end-3,:]))
    @views @inbounds v[2,:] .= T.(Tcomm.(v[end-2,:]))
    @views @inbounds v[end-1,:] .= T.(Tcomm.(v[3,:]))
    @views @inbounds v[end,:] .= T.(Tcomm.(v[4,:]))
end

""" Copy ghost points for η from inside to the halo in the nonperiodic case. """
function ghost_points_η_nonperiodic!(η::AbstractMatrix)

    # assume no gradients of η across solid boundaries
    # the 4 corner points are copied twice, but it's faster!
    @views @inbounds η[1,:] .= η[2,:]
    @views @inbounds η[end,:] .= η[end-1,:]

    @views @inbounds η[:,1] .= η[:,2]
    @views @inbounds η[:,end] .= η[:,end-1]
end

""" Copy ghost points for η from inside to the halo in the periodic case. """
function ghost_points_η_periodic!(::Type{Tcomm},η::Array{T,2}) where {Tcomm,T}

    # corner points are copied twice, but it's faster!
    # mimic MPI communication with (reduced precision) type Tcomm
    @views @inbounds η[1,:] .= T.(Tcomm.(η[end-1,:]))
    @views @inbounds η[end,:] .= T.(Tcomm.(η[2,:]))

    @views @inbounds η[:,1] .= η[:,2]
    @views @inbounds η[:,end] .= η[:,end-1]
end

""" Copy ghost points for η from inside to the halo in the nonperiodic case. """
function ghost_points_sst_nonperiodic!(G::Grid,sst::AbstractMatrix)

    @unpack halosstx,halossty = G

    # assume no gradients of η across solid boundaries
    # the 4 corner points are copied twice, but it's faster!
    for i ∈ 1:halosstx
        @views @inbounds sst[i,:] .= sst[halosstx+1,:]
        @views @inbounds sst[end-i+1,:] .= sst[end-halosstx,:]
    end

    for j ∈ 1:halossty
        @views @inbounds sst[:,j] .= sst[:,halossty+1]
        @views @inbounds sst[:,end-j+1] .= sst[:,end-halossty]
    end
end

""" Copy ghost points for η from inside to the halo in the periodic case. """
function ghost_points_sst_periodic!(::Type{Tcomm},sst::Array{T,2},G::Grid) where {Tcomm,T}

    @unpack halosstx,halossty = G

    # corner points are copied twice, but it's faster!
    for i ∈ 1:halosstx
        @views @inbounds sst[i,:] .= T.(Tcomm.(sst[end-2*halosstx+i,:]))
        @views @inbounds sst[end-halosstx+i,:] .= T.(Tcomm.(sst[halosstx+i,:]))
    end

    for j ∈ 1:halossty
        @views @inbounds sst[:,j] .= sst[:,halossty+1]
        @views @inbounds sst[:,end-j+1] .= sst[:,end-halossty]
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points!( u::AbstractMatrix,
                        v::AbstractMatrix,
                        η::AbstractMatrix,
                        S::ModelSetup)

    @unpack bc = S.parameters
    C = S.constants

    if bc == "periodic"
        @unpack Tcomm = S.parameters
        ghost_points_u_periodic!(Tcomm,u,C)
        ghost_points_v_periodic!(Tcomm,v)
        ghost_points_η_periodic!(Tcomm,η)
    else
        ghost_points_u_nonperiodic!(C,u)
        ghost_points_v_nonperiodic!(C,v)
        ghost_points_η_nonperiodic!(η)
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points_sst!(sst::AbstractMatrix,S::ModelSetup)

    @unpack bc = S.parameters
    G = S.grid

    if bc == "periodic"
        @unpack Tcomm = S.parameters
        ghost_points_sst_periodic!(Tcomm,sst,G)
    else
        ghost_points_sst_nonperiodic!(G,sst)
    end
end


text:
# Autogenerated wrapper script for MPICH_jll for i686-linux-musl-libgfortran3
export generate_compile_time_mpi_constants, libload_time_mpi_constants, libmpi, libmpicxx, libmpifort, mpiexec

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MPICH")
JLLWrappers.@declare_library_product(libload_time_mpi_constants, "libload_time_mpi_constants.so")
JLLWrappers.@declare_library_product(libmpi, "libmpi.so.12")
JLLWrappers.@declare_library_product(libmpicxx, "libmpicxx.so.12")
JLLWrappers.@declare_library_product(libmpifort, "libmpifort.so.12")
JLLWrappers.@declare_executable_product(generate_compile_time_mpi_constants)
JLLWrappers.@declare_executable_product(mpiexec)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libload_time_mpi_constants,
        "lib/libload_time_mpi_constants.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmpi,
        "lib/libmpi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmpicxx,
        "lib/libmpicxx.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmpifort,
        "lib/libmpifort.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        generate_compile_time_mpi_constants,
        "bin/generate_compile_time_mpi_constants",
    )

    JLLWrappers.@init_executable_product(
        mpiexec,
        "bin/mpiexec",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export rate_session, ses_mean, ses_std, rate_event, rate_window, rate_window_pop, zscore, zscore_pop, psth_per_con,plot_raster, plot_psth, plot_psth_raster, plot_raster_psth, plot_pop_psth, plot_pop_psth_cb, psth_pop

#=
Histogram Methods

=#

function rate_session(myrate::rate_bin,n::Int64)
    #hist(myrate.spikes[n].ts,collect(myrate.spikes[n].trials[1].time:myrate.binsize:myrate.spikes[n].trials[end].time))[2]./myrate.binsize
h=fit(Histogram,myrate.spikes[n].ts,collect(myrate.spikes[n].trials[1].time:myrate.binsize:myrate.spikes[n].trials[end].time))
h.weights ./ myrate.binsize
end

function rate_window(myrate::rate_bin,ind::Int64,ts::FloatRange{Float64},n::Int64)
    spikehist(myrate,ind,ts,n)./myrate.binsize
end

function rate_window(myrate::rate_bin,ind::Int64,ts::Array{Float64,1},n::Int64)
    spikehist(myrate,ind,ts,n)./myrate.binsize
end

function rate_event(myrate::rate_bin,inds::Array{Int64,1},ts::FloatRange{Float64},n::Int64)
    collect_ts(myrate,inds,ts,n)./(ts.step/ts.divisor)
end

function rate_event(myrate::rate_bin,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64)
    collect_ts(myrate,inds,ts,n)./(ts[2]-ts[1])
end

function rate_event(myrate::rate_bin,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64,output::Array{Float64,1})
    collect_ts(myrate,inds,ts,n,output)
    @inbounds for i=1:length(output)
        output[i]=output[i]/(ts[i+1]-ts[i])
    end
    nothing
end

function rate_event(myrate::rate_bin,inds::Array{Int64,1},ts::FloatRange{Float64},n::Int64,output::Array{Float64,1})
    collect_ts(myrate,inds,collect(ts),n,output)
end

#=
TODO
Shimazaki H. and Shinomoto S., A method for selecting the bin size of a time histogram
=#

#=
Kernel Density Methods

=#
function rate_session(myrate::rate_KD,n::Int64)
    binned=hist(myrate.spikes[n].ts,myrate.spikes[n].trials[1].time:myrate.binsize:myrate.spikes[n].trials[end].time)[2]
    imfilter_gaussian(binned,[myrate.kern/myrate.binsize])./myrate.binsize
end

function rate_window(myrate::rate_KD,ind::Int64,time::FloatRange{Float64},n::Int64)
    binned=spikehist(myrate,ind,time,n)
    imfilter_gaussian(binned,[myrate.kern/myrate.binsize])./myrate.binsize
end

function rate_event(myrate::rate_KD,inds::Array{Int64,1},time::FloatRange{Float64},n::Int64)
    binned=collect_ts(myrate,inds,time,n)
    imfilter_gaussian(binned,[myrate.kern/myrate.binsize])./myrate.binsize
end

#=
Shimazaki H. and Shinomoto S., Kernel Bandwidth Optimization in Spike Rate Estimation
=#
function rate_event_kdopt(spikes::SpikeTrain,inds::Array{Int64,1},time::Array{Float64,1})
    spikes_temp=sort(vcat([spikes.ts[spikes.trials[i].inds]-spikes.center[i,1] for i in inds]...))

    isi=diff(spikes_temp)
    dt=minimum(isi)

    #t=linspace(time[1],time[2],round(Int64,min((time[2]-time[1])/dt+.5,1000)))
    t=linspace(time[1],time[2],1000)

    dt=minimum(diff(t))

    binned=hist(spikes_temp, [t-dt/2; t[end]+dt/2])[2]

    L=length(binned)
    N=length(spikes_temp)

    wmin=2*dt
    wmax=time[2]-time[1]
    imax=20

    w=zeros(Float64,imax)
    c=zeros(Float64,imax)
    optw=0.0
    tol=1e-5

    phi=.5*(sqrt(5)+1)

    a=ilogexp(wmin)
    b=ilogexp(wmax)

    c1= (phi-1) * a + (2-phi) * b
    c2= (2 - phi)*a*(phi-1) *b

    f1=cost_function(binned, N, logexp(c1), dt)
    f2=cost_function(binned, N, logexp(c1), dt)

    k=1

    while ((abs(b-a)) > (tol * (abs(c1) + abs(c2)))) & (k<imax)
        if f1<f2
            b=c2
            c2=c1
            c1 = (phi-1)*a + (2-phi) * b
            f2=f1
            f1=cost_function(binned,N,logexp(c1),dt)
            w[k]=logexp(c1)
            c[k]=f1
            optw=logexp(c1)

        else
            a=c1
            c1=c2
            c2=(2-phi)*a+(phi-1)*b
            f1=f2
            f2=cost_function(binned,N,logexp(c1),dt)
            w[k]=logexp(c2)
            c[k]=f2
            optw=logexp(c2)
        end
        k=k+1
    end

    yh=abs(imfilter_gaussian(binned, [optw/dt]))./(length(inds)*dt)

    (yh, optw/dt)
end

function ilogexp(x::Float64)
   if x < 1e2
        y = log(exp(x)-1)
    else
        y = x
    end

    y
end

function logexp(x::Float64)
    if x < 1e2
        y = log(1 + exp(x))
    else
        y = x
    end
    y
end

function cost_function(binned::Array{Int64,1},N::Int64,w::Float64,dt::Float64)
    #The cost function
    #Cn(w) = sum_{i,j} int k(x - x_i) k(x - x_j) dx - 2 sum_{i~=j} k(x_i - x_j)

    yh=abs(imfilter_gaussian(binned, [w/dt]))

    C=sum(yh.*2)*dt - 2 * sum(yh.*binned) * dt + 2 / sqrt(2*pi) / w / N
    C = C * N * N

    C
end

#=
Helper Methods
=#

#Mean firing rate for entire session
function ses_mean(myrate::rate,n::Int64)
    mean(rate_session(myrate,n))
end

#STD of firing rate for whole session
function ses_std(myrate::rate,n::Int64)
    std(rate_session(myrate,n))
end

#
function rate_trials(myrate::rate,inds::Array{Int64,1},time::FloatRange{Float64},n::Int64)

    spikes_temp=zeros(Float64,length(inds),length(time[1]:myrate.binsize:time[end])-1)
    for i=1:length(inds)
        spikes_temp[i,:]=rate_window(myrate,inds[i],time,n)
    end
    spikes_temp
end

function spikehist(myrate::rate,ind::Int64,ts::FloatRange{Float64},n::Int64)
	 h=fit(Histogram,myrate.spikes[n].ts[myrate.spikes[n].trials[ind].inds]-myrate.spikes[n].center[ind,1],ts[1]:myrate.binsize:ts[end])
    #hist(myrate.spikes[n].ts[myrate.spikes[n].trials[ind].inds]-myrate.spikes[n].center[ind,1],ts[1]:myrate.binsize:ts[end])[2]
h.weights
end

function spikehist(myrate::rate,ind::Int64,ts::Array{Float64,1},n::Int64)
    hist(myrate.spikes[n].ts[myrate.spikes[n].trials[ind].inds]-myrate.spikes[n].center[ind,1],ts)[2]
end

function collect_ts(myrate::rate,inds::Array{Int64,1},ts::FloatRange{Float64},n::Int64,output=zeros(Float64,length(ts)-1))
    collect_ts(myrate,inds,collect(ts),n,output)
end

function collect_ts(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64,output=zeros(Float64,length(ts)-1))
    @inbounds for i=1:length(output)
        output[i]=0.0
    end
   @inbounds for i in inds
        count=1
        mycent=myrate.spikes[n].center[i,1]
        for j=1:length(myrate.spikes[n].trials[i].inds)
            myt=myrate.spikes[n].ts[myrate.spikes[n].trials[i].inds[j]]-mycent
            if myt>ts[count]
                k=count+1
                while k < length(ts)+1
                    if myt<ts[k]
                        output[k-1]+=1.0
                        count=k-1
                        break
                    end
                    k+=1
                end
            end
        end
    end
    @inbounds for i=1:length(output)
        output[i]=output[i]/length(inds)
    end
    output
end

function zscore(myrate::rate,inds::Array{Int64,1},ts::FloatRange{Float64},n::Int64)
    mypsth=rate_event(myrate,inds,ts,n)
    wholetrain=rate_session(myrate,n)
    (mypsth-mean(wholetrain))./std(wholetrain)
end

function zscore(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64)
    mypsth=rate_event(myrate,inds,ts,n)
    wholetrain=rate_session(myrate,n)
    (mypsth-mean(wholetrain))./std(wholetrain)
end

function rate_window_pop(myrate::rate,ind::Int64,ts::FloatRange{Float64})
     rate_window_pop(myrate,ind,collect(ts))
end

function rate_window_pop(myrate::rate,ind::Int64,ts::Array{Float64,1})
    raster=zeros(Float64,length(time[1]:myrate.binsize:time[end])-1,length(myrate.spikes))

    for i=1:size(raster,2)
        raster[:,i]=rate_window(myrate,ind,time,i)
    end
    raster
end

zscore_pop(myrate::rate,inds::Array{Int64,1},ts::FloatRange{Float64})=zscore_pop(myrate,inds,collect(ts))

function zscore_pop(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1})
    raster=zeros(Float64,length(myrate.spikes),length(ts)-1)

    for k=1:size(raster,1)
        raster[k,:]=zscore(myrate,inds,ts,k)
    end
    raster
end

psth_pop(myrate::rate,inds::Array{Int64,1},ts::FloatRange{Float64},raster=zeros(Float64,length(myrate.spikes),length(ts)-1))=psth_pop(myrate,inds,collect(ts),raster)

function psth_pop(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},raster=zeros(Float64,length(myrate.spikes),length(ts)-1))
    tempmat=zeros(Float64,length(ts)-1)
    for k=1:size(raster,1)
        rate_event(myrate,inds,ts,k,tempmat)
        @inbounds for j=1:length(tempmat)
            raster[k,j]=tempmat[j]
        end
    end
    raster
end

function total_spikes(myrate::rate,inds::Array{Int64,1},n::Int64)
    mysize=0
    for i in inds
        @inbounds mysize+=length(myrate.spikes[n].trials[i].inds)
    end
    mysize
end

function total_spikes(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64)

    mysize=0
    @inbounds for i in inds
        mycent=myrate.spikes[n].center[i,1]
        for j in myrate.spikes[n].trials[i].inds
            if myrate.spikes[n].ts[j]-mycent > ts[end]
                break
            elseif myrate.spikes[n].ts[j]-mycent > ts[1]
                mysize+=1
            end
        end
    end
    mysize
end

function psth_per_con(myrate::rate,ts::Array{Float64,1},inds::Array{Int64,1},n::Int64,cons::Array{Int64,1})

    psths=zeros(Float64,length(ts)-1,length(cons))

    for i=1:length(cons)
        con_inds=find(inds.==cons[i])
        psths[:,i]=rate_event(myrate,con_inds,ts,n)
    end

    psths
end

#=
Plotting
=#

#=
Raster
=#

function plot_raster(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64,ax; event_ts=zeros(Float64,0))

    plot_raster(myrate,inds,ts,n,ax)

    if length(event_ts)>0

        xy=[zeros(Float64,2,2) for i=1:(size(event_ts,2)*length(inds))]

        mycount=1.0
        myinds=1
        for i in inds
            for j=1:size(event_ts,2)
                xy[myinds][1,1]=event_ts[i]-.015
                xy[myinds][1,2]=mycount+.5
                xy[myinds][2,1]=event_ts[i]+.015
                xy[myinds][2,2]=mycount+.5
                myinds+=1
            end
            mycount+=1.0
        end
        c2=C.LineCollection(xy,color="black",linewidth=2.0)
        ax[:add_collection](c2)
    end
    nothing
end

function plot_raster(myrate::rate,inds::Array{Int64,1},ts::Array{Float64,1},n::Int64,ax)

    mycount=1.0
    myspikes=1
    xy=[zeros(Float64,2,2) for i=1:total_spikes(myrate,inds,ts,n)]
    @inbounds for i in inds
        mycenter=myrate.spikes[n].center[i,1]
        for j in myrate.spikes[n].trials[i].inds
            if myrate.spikes[n].ts[j]-mycenter > ts[end]
                break
            elseif myrate.spikes[n].ts[j] - mycenter > ts[1]
                xy[myspikes][1,1]=myrate.spikes[n].ts[j]-mycenter
                xy[myspikes][1,2]=mycount
                xy[myspikes][2,1]=myrate.spikes[n].ts[j]-mycenter
                xy[myspikes][2,2]=mycount+1.0
                myspikes+=1
            end
        end
        mycount+=1.0
    end
    c1=C.LineCollection(xy,color="black",linewidth=1.0)
    ax[:add_collection](c1)
    ax[:set_yticks]([])
    ax[:set_ylabel]("")
    ax[:set_ylim]([1,mycount])
    ax[:set_xticks]([])
    ax[:set_xticklabels]([])
    ax[:set_xlim]([ts[1],ts[end]])
    nothing
end

#=
PSTH
=#

plot_psth(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},n::Int64)=plot_psth(myrate,collect(ts),inds,n)

plot_psth(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},n::Int64,ax)=plot_psth(myrate,collect(ts),inds,n,ax)

function plot_psth(myrate::rate,ts::Array{Float64,1},inds::Array{Int64,1},n::Int64)

    (fig,ax)=subplots(1,1)

    plot_psth(myrate,ts,inds,n,ax)

    (fig,ax)
end

function plot_psth(myrate::rate,ts::Array{Float64,1},inds::Array{Int64,1},n::Int64,ax)

    psth=rate_event(myrate,inds,ts,n)

    ylimit=round(Int64,maximum(psth)/10)*10
    ax[:set_yticks]([0, ylimit])
    ax[:set_yticklabels]([0,ylimit],size=6)
    ax[:set_ylabel]("Rate (spikes/s)", size=6)

    ax[:set_xticks]([])
    ax[:set_xticklabels]([])
    ax[:bar](2:length(ts),psth,width = 1.0, color="black")
    ax[:set_ylim]([0, ylimit])
    nothing
end

plot_psth_raster(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},n::Int64)=plot_psth_raster(myrate,collect(ts),inds,n)

function plot_psth_raster(myrate::rate,ts::Array{Float64,1},inds::Array{Int64,1},n::Int64)

    (fig,ax)=subplots(2,1)
    subplots_adjust(hspace=0.0)

    plot_psth(myrate,ts,inds,n,ax[1,1])
    plot_raster(myrate,inds,ts,n,ax[2,1])
    ax[2,1][:set_xlabel]("Time (s)", size=8)
    (fig,ax)
end

function plot_psth_raster(myrate::rate,ts::Array{Float64,1},inds::Array{Int64,1},n::Int64,ax_x::Int64,ax_y::Int64,ax;event_ts=zeros(Float64,0))
    plot_psth(myrate,ts,inds,n,ax[ax_x,ax_y])
    plot_raster(myrate,inds,ts,n,ax[ax_x+1,ax_y],event_ts=event_ts)
    nothing
end

function plot_raster_psth(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},n::Int64)

    (fig,ax)=subplots(2,1)
    subplots_adjust(hspace=0.0)

    plot_raster(myrate,inds,ts,n,ax[1,1])
    plot_psth(myrate,ts,inds,n,ax[2,1])
    ax[2,1][:set_xlabel]("Time (s)", size=8)
    (fig,ax)
end

function plot_pop_psth(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},lims::NTuple{2,Float64})

    (fig,ax)=subplots(1,1)
    myplot=plot_pop_psth(myrate,ts,inds,lims,ax)
    (fig,ax)
end

function plot_pop_psth(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},lims::NTuple{2,Float64},ax)

    popth=zscore_pop(myrate,inds,ts)

    ax[:set_yticks]([])
    ax[:set_yticklabels]([])
    ax[:set_ylabel]("Neuron", size=8)

    ax[:set_xlabel]("Time (s)", size=8)

    myplot=ax[:pcolor](popth,cmap="jet",vmin=lims[1],vmax=lims[2])
end

function plot_pop_psth_cb(myrate::rate,ts::FloatRange{Float64},inds::Array{Int64,1},lims::NTuple{2,Float64})
    (fig,ax)=subplots(1,1)
    myplot=plot_pop_psth(myrate,ts,inds,lims,ax)
    cbar_ax = fig[:add_axes]([.85,.15,.02,.7])
    fig[:colorbar](myplot,cax=cbar_ax,label="Z score",ticks=[lims[1],lims[2]])
    (fig,ax)
end