text:
##############################################################################
#
# Symmetric triangular distribution from Distributions Handbook
#
##############################################################################

immutable Triangular <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Triangular(l::Real, s::Real)
        s > zero(s) || error("scale must be positive")
        new(float64(l), float64(s))
    end
end

Triangular(location::Real) = Triangular(location, 1.0)
Triangular() = Triangular(0.0, 1.0)

function cdf(d::Triangular, x::Real)
    a, b, c = d.location - d.scale, d.location + d.scale, d.location
    if x <= a
        return 0.0
    elseif a <= x <= c
        return (x - a)^2 / ((b - a) * (c - a))
    elseif c < x <= b
        return 1.0 - (b - x)^2 / ((b - a) * (b - c))
    else
        return 1.0
    end
end

entropy(d::Triangular) = 0.5 + log(d.scale)

function insupport(d::Triangular, x::Number)
    return isreal(x) && isfinite(x) &&
           d.location - d.scale <= x <= d.location + d.scale
end

kurtosis(d::Triangular) = -0.6

mean(d::Triangular) = d.location

median(d::Triangular) = d.location

mode(d::Triangular) = d.location
modes(d::Triangular) = [d.location]

function pdf(d::Triangular, x::Real)
    if insupport(d, x)
        return -abs(x - d.location) / (d.scale^2) + 1.0 / d.scale
    else
        return 0.0
    end
end

function quantile(d::Triangular, p::Real)
    a, b, c = d.location - d.scale, d.location + d.scale, d.location
    if p <= 0.0
        return a
    elseif p < 0.5
        return a + sqrt(p * 2.0 * d.scale^2)
    elseif p >= 0.5
        return b -  sqrt((1.0 - p) * 2.0 * d.scale^2)
    else
        return b
    end
end

function rand(d::Triangular)
    ξ1, ξ2 = rand(), rand()
    return d.location + (ξ1 - ξ2) * d.scale
end

function skewness(d::Triangular)
    a = d.location - d.scale
    b = d.location + d.scale
    c = (b - a) / 2 + a
    den = sqrt(2.0) * (a + b - 2.0 * c) *
                      (2.0 * a - b - c) *
                      (a - 2.0 * b + c)
    num = 5.0 * (a^2 + b^2 + c^2 - a * b - a * c - b * c)^1.5
    return den / num
end

var(d::Triangular) = d.scale^2 / 6.0


text:
function modulation_score(n::NeuronResponse, decoders::Array{ResponseDecoder,1},
                          num_spikes, num_trials, num_surrogates)
    results = Dict(d.name=>Float64[] for d in decoders)
    for _ in 1:num_surrogates
        spikes = [n.get_spikes(num_spikes) for _ in 1:num_trials]
        for decoder in decoders
            activity_matrix = reduce(hcat,fit(Histogram, s, decoder.τ).weights for s in spikes)
            # split into before and after "trials"
            mid = floor(Int,size(activity_matrix,1)/2)
            activity_matrix_split    = [activity_matrix[1:mid,:]';
                                        activity_matrix[mid+1:end,:]']
            before_after_identifiers = [zeros(size(activity_matrix,2));
                                        ones(size(activity_matrix,2))]
            performance = decoder.decode(activity_matrix_split,
                                         before_after_identifiers)
            push!(results[decoder.name],performance)
        end
    end
    return results
end

function modulation_score(responses::Array{NeuronResponse, 1},
                          decoders::Array{ResponseDecoder,1},
                          num_spikes, num_trials, num_surrogates)
    results = Dict(r.name =>
                   modulation_score(r, decoders,  num_spikes, num_trials, num_surrogates)
                   for r in responses)
end

ext: jl


text:
# developed with Julia 1.4.2
#
# arguments parsing for examples

using ArgParse

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table! s begin

        "--save"
            help = "folder to save scores and/or models"
            arg_type = String
            default = joinpath(@__DIR__, "../results")

        "--metadata"
            help = "metadata.csv - site and battery parameters"
            arg_type = String
            default = joinpath(@__DIR__, "../data/metadata.csv")

        "--train"
            help = "train data folder"
            arg_type = String
            default = joinpath(@__DIR__, "../data/train")

        "--test"
            help = "test data folder"
            arg_type = String
            default = joinpath(@__DIR__, "../data/test")

        "--price"
            help = "price folder or .csv file"
            arg_type = String
            default = joinpath(@__DIR__, "../data/prices/edf_prices.csv")

        "--calibrate"
            help = "perform model calibration - SDP only"
            action = :store_true

        "--simulate"
            help = "perform model simulation"
            action = :store_true

    end

    parsed_args = parse_args(s)

    for (key, val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end

    return parse_args(s)

end


text:
# Test sparse Jacobian implementation (even though the Jacobian of the function
# is actually not sparse...)

@testset "sparse" begin

function f_sparse!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function g_sparse!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

df = OnceDifferentiable(f_sparse!, g_sparse!, rand(2), rand(2), spzeros(2, 2))

# Test trust region
r = nlsolve(df, [ -0.5; 1.4], method = :trust_region, autoscale = true)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8

# Test Newton
r = nlsolve(df, [ -0.5; 1.4], method = :newton, linesearch = LineSearches.BackTracking(), ftol = 1e-6)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-6

# Test MCP solver with smooth reformulation
r = mcpsolve(df, [-Inf;-Inf], [Inf; Inf], [-0.5; 1.4], reformulation = :smooth)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8

# Test MCP solver with minmax reformulation
r = mcpsolve(df, [-Inf;-Inf], [Inf; Inf], [-0.5; 1.4], reformulation = :minmax)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8

# Test given sparse
df = OnceDifferentiable(f_sparse!, g_sparse!, rand(2), rand(2), spzeros(2, 2))
r = nlsolve(df, [ -0.5; 1.4], method = :trust_region, autoscale = true)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8

end


text:
export SymmetricEmbedding
export Reflection, Rotation, Symmetry

"""
	Symmetry
Supertype of all symmetries used in [`SymmetricEmbedding`](@ref).
All symmetries are initialized like `Symmetry(n1, n2, ...)` with
`ni` being the indices of the spatial dimensions that have the said symmetry.

Notice that the symmetries are defined with respect to the center point
of the embedding.
"""
abstract type Symmetry end

"""
	Rotation <: Symmetry
Index sets at equal distance from the center point are equivalent,
for the given input dimensions. E.g. for `Rotation(1,3)` and given center
point `u[i,j,k,...]` all indices `m, n` that satisfy
```
|i-m|^2 + |k-n|^2 = r^2
```
for a given `r`, are equivalent.
The same process generalizes to any number of input dimensions.
"""
struct Rotation <: Symmetry
	d::Vector{Int}
	function Rotation(args::Vector{Int})
		@assert length(args) > 1 "Rotation symmetry needs at least 2 dimensions"
		new(args)
	end
end

"""
	Reflection <: Symmetry
Reflection symmetry: x and -x are equivalent (for all given dimensions).
"""
struct Reflection <: Symmetry
	d::Vector{Int}
end

for sym in (:Rotation, :Reflection)
	@eval $(sym)(x::Int, args...) = $(sym)([x, args...])
end

_smallstr(::Rotation) = "rot"
_smallstr(::Reflection) = "refl"
function Base.show(io::IO, sym::Symmetry)
	s = _smallstr(sym)*string(sym.d)
	print(io, s)
end

# internal translation to the nested vectors of integers
# that will remain until we have time to properly re-write
# (if we ever encounter a different symmetry we want to use)
function _nestedvec(syms::Tuple)
	# This part of the code ensures no duplicate symmetries
	allidx = vcat(s.d for s in syms)
	@assert unique(allidx) == allidx
	# Now convert to nestedvec
	v = Vector{Vector{Int}}()
	for s in syms
		push!(v, s.d)
	end
	return v
end

"""
	SymmetricEmbedding(ste::SpatioTemporalEmbedding, sym::Tuple)

A `SymmetricEmbedding` is intended as a means of dimension reduction
for a [`SpatioTemporalEmbedding`](@ref) by exploiting spatial symmetries in
the system, listed as a `Tuple` of `<:Symmetry` (see [`Symmetry`](@ref)) for
all possible symmetries.

All points at a time step equivalent to each other according to the symmetries
passed in `sym` will be **averaged** to a single entry! For example,
the symmetry `Reflection(2)` means that the embedding won't have two entries
`u[i, j+1], u[i, j-1]` but instead a single entry
`(u[i, j+1] + u[i, j-1])/2`, with `i,j` being indices *relative to the
central point of the embedding*. (the same process is
done for any index offset `j+2, j+3`, etc., depending on how large the
spatial radius `r` is)

The resulting structure from `SymmetricEmbedding`
can be used for reconstructing datasets in
the same way as a [`SpatioTemporalEmbedding`](@ref).
"""
struct SymmetricEmbedding{Φ,BC,X} <: AbstractSpatialEmbedding{Φ,BC,X}
    τ::Vector{Int}
    β_groups::Vector{Vector{CartesianIndex{Φ}}}
	inner::Region{Φ}
	whole::Region{Φ}
    boundary::BC
end

function SymmetricEmbedding(ste::SpatioTemporalEmbedding, sym::Tuple)
	nv = _nestedvec(sym)
	_SymmetricEmbedding(ste, nv)
end

# Internal function, uses the nested vectors until we change it /
# clean it up
function _SymmetricEmbedding(ste::SpatioTemporalEmbedding{Φ,BC}, sym) where {Φ,BC}
	check_symmetry(sym, Φ)
    # Mapping dimension for dimension
    # single dim symmetry : x₁ → |x₁|
    # two dim symmetry : (x₁, x₂) → (y₁, y₂) = (|x₁|, |x₂|)
    # (y₁, y₂) → y₁ > y₂ ? id : (y₂, y₁)

    # arbitrary dim symmetry : x = (x₁, ...) → sort(abs.(x), reverse=true)
    τ = Vector{Int}()
    β_groups = []
	#CartesianIndices are nice but difficult to manipulate
	β = map(x-> x.I ,ste.β)
    #Repeat for every timestep
    for t in unique(ste.τ)
		#Filter for β at delay t
        β_t = β[ ste.τ .== t ]
		#Map selected β to their symmetry reduced space
        #and combine into groups according to their mapped value
        groups_t = groupbymap( x->fold_up(x, sym) , β_t)
		#Append results to new lists
        append!(τ, fill(t, length(groups_t)))
        append!(β_groups, groups_t)
    end
	# convert β_groups back to cartesian indices
	β_new = map(β_groups) do b
		map(x-> CartesianIndex(x...), b)
	end
	X = length(τ)
    return SymmetricEmbedding{Φ,BC,X}(τ, β_new, ste.inner, ste.whole, ste.boundary)
end

function fold_up(β, sym)
    b = [β...]
    for s in sym
        b[s] .= sort(abs.(b[s]), rev=true)
    end
    return tuple(b...)
end

function groupbymap(f, v; by=norm)
    d = Dict()
    for e ∈ v
        w = f(e)
        haskey(d, w) ? push!(d[w], e) : d[w] = [e]
    end
    map(k-> d[k], sort!(collect(keys(d)); by=by) )
end

function check_symmetry(sym, Φ)
	s = vcat(sym...)
	unique(s) == s || throw(ArgumentError("Dimension may not be passed more than once."))
	maximum(s) > Φ && throw(ArgumentError("Invalid dimension was passed."))
	minimum(s) < 1 && throw(ArgumentError("Invalid dimension was passed."))
end

function Base.show(io::IO, em::SymmetricEmbedding{Φ,BC, X}) where {Φ,BC,X}
	print(io, "$(Φ)D Spatio-Temporal Delay Embedding with $X Entries")
    if BC == PeriodicBoundary
        println(io, " and PeriodicBoundary condition.")
    else
        println(io, " and ConstantBoundary condition with c = $(em.boundary.c).")
    end
    println(io, "The included neighboring points are (forward embedding):")
    for (τ,β) in zip(em.τ,em.β_groups)
        println(io, "τ = $τ , β = ", getproperty.(β,:I))
    end
end

function (r::SymmetricEmbedding{Φ,ConstantBoundary{T},X})(rvec, s, t, α) where {T,Φ,X}
	if α in r.inner
		@inbounds for n=1:X
			β = r.β_groups[n]
			s_t = s[t + r.τ[n]]
			rvec[n] = zero(T)
			for m in eachindex(β)
				rvec[n] += s_t[ α + β[m] ]
			end
			rvec[n] /= length(β)
		end
	else
		@inbounds for n=1:X
			β = r.β_groups[n]
			s_t = s[t + r.τ[n]]
			rvec[n] = zero(T)
			for m in eachindex(β)
				rvec[n] += α + β[m] in r.whole ? s_t[ α + β[m] ] : r.boundary.b
			end
			rvec[n] /= length(β)
		end
	end
	return nothing
end

function (r::SymmetricEmbedding{Φ,PeriodicBoundary,X})(rvec, s, t, α) where {Φ,X}
	if α in r.inner
		@inbounds for n=1:X
			β = r.β_groups[n]
			s_t = s[t + r.τ[n]]
			rvec[n] = zero(eltype(rvec))
			for m in eachindex(β)
				rvec[n] += s_t[ α + β[m] ]
			end
			rvec[n] /= length(β)
		end
	else
		@inbounds for n=1:X
			β = r.β_groups[n]
			s_t = s[t + r.τ[n]]
			rvec[n] = zero(eltype(rvec))
			for m in eachindex(β)
				rvec[n] = s_t[ project_inside(α + β[m], r.whole) ]
			end
			rvec[n] /= length(β)
		end
	end
	return nothing
end


text:
function methods!(
    KS::KineticSolver,
    psi::Array{Float64,3},
    psiOld::Array{Float64,3},
    k::Int64,
)

    if (KS.ThisQuadratureType == "octa" || KS.ThisQuadratureType == "ico") &&
       KS.rotationmagnitude > 0 # rSN rotation step
        u = randn(3)
        u = u / norm(u)
        theta = KS.rotationmagnitude * KS.dt / KS.norder
        rotateandinterpolate!(u, theta, KS.Q, psiOld, psi, "forward") # writes into psiOld
        @inbounds for i = 1:length(psi)
            psi[i] = psiOld[i]
        end

    elseif KS.ThisQuadratureType == "tensorized" && KS.rotationmagnitude > 0 # Rotation just around the z-axis
        if k % 2 == 0
            theta = KS.rotationmagnitude * KS.dt / KS.norder       # and ensure it lies in Delta psi
            alpha = theta / (pi / KS.norder) % 1.0
            rotateOnly!([0.0; 0.0; 1.0], theta, KS.Q, "forward")
            psiInterp = zeros(Int(KS.norder / 2), 2 * KS.norder)
            for i = 1:size(psi, 2)
                for j = 1:size(psi, 3)
                    psiMat = reshape(psi[:, i, j], Int(KS.norder / 2), 2 * KS.norder)
                    for q = 1:(2*KS.norder-1)
                        psiInterp[:, q] =
                            (1 - alpha) * psiMat[:, q] + alpha * psiMat[:, q+1]
                    end
                    psiInterp[:, 2*KS.norder] =
                        (1 - alpha) * psiMat[:, 2*KS.norder] + alpha * psiMat[:, 1]
                    psi[:, i, j] = psiInterp[:]
                end
            end
        else
            theta = KS.rotationmagnitude * KS.dt / KS.norder       # and ensure it lies in Delta psi
            alpha = theta % (pi / KS.norder)
            rotateOnly!([0.0; 0.0; -1.0], theta, KS.Q, "forward")
            psiInterp = zeros(Int(KS.norder / 2), 2 * KS.norder)
            for i = 1:size(psi, 2)
                for j = 1:size(psi, 3)
                    psiMat = reshape(psi[:, i, j], Int(KS.norder / 2), 2 * KS.norder)
                    for q = 2:2*KS.norder
                        psiInterp[:, q] =
                            (1 - alpha) * psiMat[:, q] + alpha * psiMat[:, q-1]
                    end
                    psiInterp[:, 1] =
                        (1 - alpha) * psiMat[:, 1] + alpha * psiMat[:, 2*KS.norder]
                    psi[:, i, j] = psiInterp[:]
                end
            end
        end
    end

    if KS.convolutionmagnitude > 0 # Idea of Cory Hauck, perform convolution in the angular domain

        @inbounds @views for j = 1:size(psi, 3)
            for i = 1:size(psi, 2)
                if norm(psi[:, i, j]) > 1e-12
                    mul!(
                        KS.tmp,
                        KS.ScatteringKernelConvolution,
                        KS.Q.weights .* psi[:, i, j],
                    )
                    psi[:, i, j] .+=
                        KS.dt .* KS.convolutionmagnitude .*
                        (KS.tmp - KS.OutscatteringVector .* KS.Q.weights .* psi[:, i, j])
                end
            end
        end
    end
end

function solveFluxUpwind!(KS::KineticSolver, psi::Array{Float64,3}, flux::Array{Float64,3})
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))

    # PosPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosPos
        s0 = psi[q, i, j-3]
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]
        s3 = psi[q, i, j]
        northflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        southflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        s0 = psi[q, i-3, j]
        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]
        s3 = psi[q, i, j]
        eastflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        westflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end
    #PosNeg
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosNeg
        s0 = psi[q, i, j+3]
        s1 = psi[q, i, j+2]
        s2 = psi[q, i, j+1]
        s3 = psi[q, i, j]
        northflux = 7 / 4 * s2 - s1 + 1 / 4 * s0
        southflux = 7 / 4 * s3 - s2 + 1 / 4 * s1

        s0 = psi[q, i-3, j]
        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]
        s3 = psi[q, i, j]
        eastflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        westflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end

    # NegPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxNegPos
        s0 = psi[q, i, j-3]
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]
        s3 = psi[q, i, j]
        northflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        southflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        s0 = psi[q, i+3, j]
        s1 = psi[q, i+2, j]
        s2 = psi[q, i+1, j]
        s3 = psi[q, i, j]
        westflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        eastflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end

    # NegNeg
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxNegNeg
        s0 = psi[q, i, j+3]
        s1 = psi[q, i, j+2]
        s2 = psi[q, i, j+1]
        s3 = psi[q, i, j]
        northflux = 7 / 4 * s2 - s1 + 1 / 4 * s0
        southflux = 7 / 4 * s3 - s2 + 1 / 4 * s1

        s0 = psi[q, i+3, j]
        s1 = psi[q, i+2, j]
        s2 = psi[q, i+1, j]
        s3 = psi[q, i, j]
        westflux = 7 / 4 * s3 - s2 + 1 / 4 * s1
        eastflux = 7 / 4 * s2 - s1 + 1 / 4 * s0

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end
end

function Euler!(
    KS::KineticSolver,
    psi::Array{Float64,3},
    flux::Array{Float64,3},
    dt::Float64,
)
    # computes the Euler step
    @inbounds for j in KS.rangex, i in KS.rangey
        integral = 0.0
        @inbounds for q in KS.rangequad
            integral += psi[q, i, j] .* KS.Q.weights[q]
        end
        integral = integral * KS.SigmaS[i, j] / 4.0 / pi
        sigmaT = KS.SigmaT[i, j]
        @inbounds for q in KS.rangequad
            psi[q, i, j] += dt .* (-psi[q, i, j] .* sigmaT - flux[q, i, j] + integral)
        end
    end
end

# the first minmod code is the fast version of the second minmod  below that is commented
@inline minmod(x::Float64, y::Float64) = ifelse(x < 0, clamp(y, x, 0.0), clamp(y, 0.0, x))
#@inline function minmod(x::Float64, y::Float64)
#    return sign(x) * max(0.0, min(abs(x),y*sign(x) ) );
#end

@inline function slopefit(left::Float64, center::Float64, right::Float64)
    tmp = minmod(0.5 * (right - left), 2.0 * (center - left))
    return minmod(2.0 * (right - center), tmp)
end

function computeMoments(KS::KineticSolver, psi::Array{Float64,3})
    rho = zeros(KS.ny, KS.nx)
    for q in KS.rangequad, j = 1:KS.nx, i = 1:KS.ny#KS.rangex, i=KS.rangey
        rho[i, j] += psi[q, i, j] * KS.Q.weights[q]
    end
    return rho
end

function forcePeriodicity!(KS::KineticSolver, psi::Array{Float64,3})
    if KS.periodicXflag
        for q = 1:KS.nquadpoints
            for i = 1:KS.ny
                psi[q, i, 1] = psi[q, i, end-3]
                psi[q, i, 2] = psi[q, i, end-2]
                psi[q, i, end-1] = psi[q, i, 3]
                psi[q, i, end] = psi[q, i, 4]
            end
        end
    end
    if KS.periodicYflag
        for q = 1:KS.nquadpoints
            for j = 1:KS.nx
                psi[q, 1, j] = psi[q, end-3, j]
                psi[q, 2, j] = psi[q, end-2, j]
                psi[q, end-1, j] = psi[q, 3, j]
                psi[q, end, j] = psi[q, 4, j]
            end
        end
    end
end

function copyFromToNoGhosts!(
    psiFrom::Array{Float64,3},
    psiTo::Array{Float64,2},
    KS::KineticSolver,
)
    for q in KS.rangequad
        k = 0
        for i = 3:KS.ny-2
            for j = 3:KS.nx-3
                k = k + 1
                psiTo[q, k] = psiFrom[q, i, j]
            end
        end
    end
end

function copyFromToWithGhosts!(
    psiFrom::Array{Float64,2},
    psiTo::Array{Float64,3},
    KS::KineticSolver,
)
    for q in KS.rangequad
        k = 0
        for i = 1:KS.ny-4
            for j = 1:KS.nx-4
                k = k + 1
                psiTo[q, i+2, j+2] = psiFrom[q, k]
            end
        end
    end
end

function solveFlux!(KS::KineticSolver, phi::Array{Float64,3}, flux::Array{Float64,3})
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))

    # PosPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosPos
        s1 = phi[q, i, j-2]
        s2 = phi[q, i, j-1]
        s3 = phi[q, i, j]
        s4 = phi[q, i, j+1]
        northflux = s3 + 0.5 .* slopefit(s2, s3, s4)
        southflux = s2 + 0.5 .* slopefit(s1, s2, s3)

        s1 = phi[q, i-2, j]
        s2 = phi[q, i-1, j]
        s3 = phi[q, i, j]
        s4 = phi[q, i+1, j]
        eastflux = s3 + 0.5 .* slopefit(s2, s3, s4)
        westflux = s2 + 0.5 .* slopefit(s1, s2, s3)

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end
    #PosNeg
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosNeg
        s1 = phi[q, i, j-1]
        s2 = phi[q, i, j]
        s3 = phi[q, i, j+1]
        s4 = phi[q, i, j+2]
        northflux = s3 - 0.5 .* slopefit(s2, s3, s4)
        southflux = s2 - 0.5 .* slopefit(s1, s2, s3)

        s1 = phi[q, i-2, j]
        s2 = phi[q, i-1, j]
        s3 = phi[q, i, j]
        s4 = phi[q, i+1, j]
        eastflux = s3 + 0.5 .* slopefit(s2, s3, s4)
        westflux = s2 + 0.5 .* slopefit(s1, s2, s3)

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end

    # NegPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxNegPos
        s1 = phi[q, i, j-2]
        s2 = phi[q, i, j-1]
        s3 = phi[q, i, j]
        s4 = phi[q, i, j+1]
        northflux = s3 + 0.5 .* slopefit(s2, s3, s4)
        southflux = s2 + 0.5 .* slopefit(s1, s2, s3)

        s1 = phi[q, i-1, j]
        s2 = phi[q, i, j]
        s3 = phi[q, i+1, j]
        s4 = phi[q, i+2, j]
        eastflux = s3 - 0.5 .* slopefit(s2, s3, s4)
        westflux = s2 - 0.5 .* slopefit(s1, s2, s3)

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end

    # NegNeg
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxNegNeg
        s1 = phi[q, i, j-1]
        s2 = phi[q, i, j]
        s3 = phi[q, i, j+1]
        s4 = phi[q, i, j+2]
        northflux = s3 - 0.5 .* slopefit(s2, s3, s4)
        southflux = s2 - 0.5 .* slopefit(s1, s2, s3)

        s1 = phi[q, i-1, j]
        s2 = phi[q, i, j]
        s3 = phi[q, i+1, j]
        s4 = phi[q, i+2, j]
        eastflux = s3 - 0.5 .* slopefit(s2, s3, s4)
        westflux = s2 - 0.5 .* slopefit(s1, s2, s3)

        flux[q, i, j] =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux - westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux - southflux)
    end
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
include("cpu_stencils.jl")

## Stencil Definition
radius = 4
l = 2*radius + 1
coefs = zeros(Float64, (l,l,l))

for i = 1:radius
    global radius
    posi = radius + i + 1
    negi = radius -i + 1
    mid = radius + 1
    coefs[posi,mid,mid] = 1
    coefs[negi,mid,mid] = 1
    coefs[mid, posi, mid] = 1
    coefs[mid, negi, mid] = 1
    coefs[mid, mid, posi] = 1
    coefs[mid, mid, negi] = 1

end
mid = radius + 1
coefs[mid,mid,mid] = 1
st_def = CreateStencilDefinition(coefs, uses_vsq=true)
st_inst = NewStencilInstance(st_def)
## Input Data size Definition
radius = 4
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 1
## GPU Implementation
original_data = CreateData(dx,dy,dz)
data = PadData(radius, original_data)
vsq = ones(Float32, dx, dy, dz)
vsq[10:15,:,:] .= 0.1
dout = ApplyStencil(st_inst, data, t_steps,vsq=vsq)
gpu_out = Array(dout)

## CPU Implementation
coefs = Float64.([1 1 1 1 1])
cpu_out =  cpu_star_stencil_v(original_data, vsq, radius, coefs, t_steps)

## Check Result
mean_er = sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out)))
@assert(sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out))) < 0.001)
## Visual Check
# using Plots
# shift = -10
# Int(round(dz/2))+shift
# heatmap(cout[:,:,Int(round(dz/2))+shift])
# heatmap(o[:,:,Int(round(dz/2))+shift])


text:
module SymbolicMDPs

export SymbolicMDP, SymbolicRLEnv

using Random

import PDDL: PDDL, Domain, Problem, State, Term, Compound, Const, @pddl_str
import POMDPs: POMDPs, MDP
import POMDPModelTools: Deterministic, Uniform, SparseCat
import POMDPModelTools: CommonRLInterface, MDPCommonRLEnv

## Utility functions ##

typerange(T) = nothing
typerange(T::Type{<:Integer}) = typemin(T):typemax(T)
typerange(T::Type{<:Enum}) = instances(T)

compare_terms(x::Term, y::Term) = x.name < y.name
compare_terms(x::Compound, y::Compound) = string(x) < string(y)

## SymbolicStateSpace ##

"Product representation of a symbolic state space for ground PDDL domains."
struct SymbolicStateSpace{S <: State, I, FT, FR}
    state::S # Reference state with static values of fluents
    ftypes::FT # Map from (non-static) fluent names to their Julia types
    franges::FR # Map from (non-static) fluent names to their value ranges
    fluents::Vector{Term} # List of all (non-static) ground fluents
    iter::I # Iterator over states
end

function SymbolicStateSpace(domain::Domain, state::State)
    # Make copy of state for reference
    state = deepcopy(state)
    # Infer static fluents so that we can avoid enumerating over their values
    statics = PDDL.infer_static_fluents(domain)
    # Extract and order fluent types and value ranges
    sigs = sort(collect(pairs(PDDL.get_fluents(domain))), by=first)
    ftypes = (; (f => PDDL.GLOBAL_DATATYPES[sig.type]
                 for (f, sig) in sigs if !(f in statics))...)
    franges = (; (f => typerange(ty) for (f, ty) in pairs(ftypes))...)
    # Construct list of all (non-static) ground fluents
    fluents = Term[]
    for (fname, sig) in sigs
        if fname in statics continue end
        if PDDL.arity(sig) == 0
            push!(fluents, Const(fname))
        else
            fs = [Compound(fname, collect(Term, args))
                  for args in PDDL.groundargs(domain, state, fname)] |> vec
            append!(fluents, sort!(fs, lt=compare_terms))
        end
    end
    # Construct iterator over states
    if (any(!(ftypes[f.name] <: Union{Integer,Enum}) for f in fluents) ||
        length(fluents) > Sys.WORD_SIZE)
        # Not iterable if some fluents are continuous, or there are too many
        state_iter = nothing
    else
        # Construct product iterator over ground fluent values
        val_iters = (franges[f.name] for f in fluents)
        state_iter = Base.Generator(Iterators.product(val_iters...)) do vals
            s = copy(state)
            for (f, val) in zip(fluents, vals)
                s[f] = val
            end
            return s
        end
    end
    return SymbolicStateSpace(state, ftypes, franges, fluents, state_iter)
end

Base.eltype(space::SymbolicStateSpace{S}) where {S} = S

Base.iterate(space::SymbolicStateSpace{S,Nothing}) where {S} =
    error("Non-iterable space.")
Base.iterate(space::SymbolicStateSpace{S,Nothing}, i) where {S} =
    error("Non-iterable space.")
Base.length(space::SymbolicStateSpace{S,Nothing}) where {S} =
    error("Non-iterable space.")

Base.iterate(space::SymbolicStateSpace{S}) where {S} =
    iterate(space.iter)
Base.iterate(space::SymbolicStateSpace{S}, i) where {S} =
    iterate(space.iter, i)
Base.length(space::SymbolicStateSpace{S}) where {S} =
    prod(length(space.franges[f.name]) for f in space.fluents)

function Base.rand(rng::AbstractRNG, space::SymbolicStateSpace)
    state = copy(space.state)
    val_iters = for f in space.fluents
        ty = space.ftypes[f.name]
        state[f] = Base.rand(rng, ty)
    end
    return state
end

function POMDPs.stateindex(space::SymbolicStateSpace{S}, state::S) where {S}
    if space.iter === nothing
        error("State space is not discrete.")
    end
    idx = 0
    for f in Iterators.reverse(space.fluents)
        ty = space.ftypes[f.name]
        vrng = space.franges[f.name]
        val = state[f]
        offset = (ty <: Enum ? Int(val) : findfirst(==(val), vrng)) - 1
        idx = idx * length(vrng) + offset
    end
    return idx + 1
end

function vectorize(space::SymbolicStateSpace, state::State)
    T = typejoin(space.ftypes...)
    return broadcast(f -> state[f]::T, space.fluents)
end

## SymbolicMDP ##

"""
    SymbolicMDP

MDP wrapper for a ground PDDL domain.
"""
struct SymbolicMDP{D<:Domain,S<:State,SS<:SymbolicStateSpace,C} <: MDP{S,Term}
    domain::D
    init::S
    goal::Term
    metric::Union{Term,Nothing}
    discount::Float64
    states::SS
    actions::Vector{Compound}
    a_cache::C
end

"""
    SymbolicMDP(domain, state, [goal, metric, discount]; cache_actions=false)

Construct a symbolic MDP from a PDDL `domain` and initial `state`.
"""
function SymbolicMDP(domain::Domain, state::State,
                     goal=pddl"(true)", metric=nothing, discount=1.0;
                     cache_actions::Bool=true)
    s_space = SymbolicStateSpace(domain, state)
    actions = Compound[act.term for act in PDDL.groundactions(domain, state)]
    actions = sort!(actions, lt=compare_terms)
    a_cache = cache_actions ? Dict{UInt64,Vector{Compound}}() : nothing
    return SymbolicMDP(domain, state, goal, metric, discount,
                       s_space, actions, a_cache)
end

"""
    SymbolicMDP(domain, problem; cache_actions=false)

Construct a symbolic MDP from a PDDL `domain` and `problem`.
"""
function SymbolicMDP(domain::Domain, problem::Problem; options...)
    state = PDDL.initstate(domain, problem)
    goal = PDDL.get_goal(problem)
    metric = PDDL.get_metric(problem)
    if metric !== nothing # Extract metric formula to minimize
        metric = metric.name == :minimize ?
            metric.args[1] : Compound(:-, metric.args)
    end
    return SymbolicMDP(domain, state, goal, metric; options...)
end

POMDPs.states(m::SymbolicMDP) = m.states
POMDPs.actions(m::SymbolicMDP) = m.actions
POMDPs.actions(m::SymbolicMDP, s) = m.a_cache isa Nothing ?
    collect(PDDL.available(m.domain, s)) :
    get!(() -> collect(PDDL.available(m.domain, s)), m.a_cache, hash(s))

POMDPs.initialstate(m::SymbolicMDP) =
    Deterministic(m.init)
POMDPs.transition(m::SymbolicMDP, s, a) =
    Deterministic(PDDL.transition(m.domain, s, a))
POMDPs.reward(m::SymbolicMDP, s, a, sp) = m.metric === nothing ?
    -1 : m.domain[s => m.metric] - m.domain[sp => m.metric]
POMDPs.discount(m::SymbolicMDP) =
    m.discount
POMDPs.isterminal(m::SymbolicMDP, s) =
    PDDL.satisfy(m.domain, s, m.goal)

POMDPs.stateindex(m::SymbolicMDP, s) =
    POMDPs.stateindex(m.states, s)
POMDPs.actionindex(m::SymbolicMDP, a) =
    searchsortedfirst(m.actions, a; lt=compare_terms)

POMDPs.convert_s(::Type{Any}, s::State, m::SymbolicMDP) =
    s
POMDPs.convert_s(::Type{<:AbstractArray}, s::State, m::SymbolicMDP) =
    vectorize(m.states, s)
POMDPs.convert_s(::Type{Int}, s::State, m::SymbolicMDP) =
    POMDPs.stateindex(m, s)

POMDPs.convert_a(::Type{Any}, a::Term, m::SymbolicMDP) =
    a
POMDPs.convert_a(::Type{<:AbstractArray}, a::Term, m::SymbolicMDP) =
    [POMDPs.actionindex(m, a)]
POMDPs.convert_a(::Type{Int}, a::Term, m::SymbolicMDP) =
    POMDPs.actionindex(m, a)

## Implement additional common RL environment interface functions ##

const CRL = CommonRLInterface
const SymbolicMDPCommonRLEnv =
    MDPCommonRLEnv{RLO, M} where {RLO, M <: SymbolicMDP}

"""
    SymbolicRLEnv([S::Type], mdp::SymbolicMDP)
    SymbolicRLEnv([S::Type], domain, problem)
    SymbolicRLEnv([S::Type], domain, state, [goal, metric, discount])

Construct a reinforcement learning environment that adheres to the
[`CommonRLInterface`](https://git.io/J9OzR). Environments are constructed
in a similar manner to [`SymbolicMDP`](@ref), with the additional option
of specifiying a type `S` which the underlying PDDL state will be converted to.

By default, `S` is `Any`, implying no conversion. Other supported types are:
- `Int`: states are converted using [`POMDPs.stateindex`](@ref).
- `AbstractArray`: states are converted via [`SymbolicMDPs.vectorize`](@ref).

Conversion is implemented via [`POMDPs.convert_s`](@ref). To support
custom state converters, define a new type `T`, and then a new
method `POMDPs.convert_s(::Type{T}, s::State, ::SymbolicMDP)` that
performs the relevant conversion. Note that `T` is allowed to be abstract,
and need not correspond to the type actually returned.
"""
SymbolicRLEnv(S::Type, mdp::SymbolicMDP) =
    MDPCommonRLEnv{S}(mdp)
SymbolicRLEnv(S::Type, domain::Domain, args...; kwargs...) =
    SymbolicRLEnv(S, SymbolicMDP(domain, args...; kwargs...))
SymbolicRLEnv(domain::Domain, args...; kwargs...) =
    SymbolicRLEnv(Any, domain, args...; kwargs...)

CRL.@provide CRL.valid_action_mask(env::SymbolicMDPCommonRLEnv) =
    broadcast(in(CRL.valid_actions(env)), env.m.actions)

end


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

"""
	Representation of a 3D vortex particle in CVortex

You do not need this to use the CVortex API.

coord is a particle's position.
vorticity is the particle's vorticity
volume is the volume of the particle. This is only important for
viscous vortex particle strength exchange methods (not included in this
wrapper)
"""
struct VortexParticle3D
    coord :: Vec3f
    vorticity :: Vec3f
    volume :: Float32
end

function VortexParticle3D(coord::Vec3f, vort::Vec3f)
    return VortexParticle3D(coord, vort, 0.0)
end

function VortexParticle3D(coord::Vector{<:Real}, vort::Vector{<:Real}, vol::Real)
	@assert(length(coord)==3)
	@assert(length(vort)==3)
    return VortexParticle3D(Vec3f(coord), Vec3f(vort), vol)
end

function VortexParticle3D(coord::Vector{<:Real}, vort::Vector{<:Real})
	@assert(length(coord)==3)
	@assert(length(vort)==3)
    return VortexParticle3D(coord, vort, 0.0)
end

"""
# 3D vortex particle functions
Single particle -> single velocity measurement

	particle_induced_velocity(
		inducing_particle_position :: Vector{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		measurement_point :: Vector{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

Multiple particles -> single velocity measurement

	particle_induced_velocity(
		inducing_particle_positions :: Matrix{<:Real},
		inducing_particle_vorticities :: Matrix{<:Real},
		measurement_point :: Vector{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

Multiple particles -> multiple velocity measurements

	particle_induced_velocity(
		inducing_particle_positions :: Matrix{<:Real},
		inducing_particle_vorticities :: Matrix{<:Real},
		measurement_points :: Matrix{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

# 2D vortex particle functions
Single particle -> single velocity measurement

	particle_induced_velocity(
		inducing_particle_position :: Vector{<:Real},
		inducing_particle_vorticity :: Real,
		measurement_point :: Vector{<:Real},
		kernel :: RegularisationFunction,
		regularisation_radius :: Real)

Multiple particles -> single velocity measurement

	particle_induced_velocity(
		inducing_particle_positions :: Matrix{<:Real},
		inducing_particle_vorticities :: Vector{<:Real},
		measurement_point :: Vector{<:Real},
		kernel :: RegularisationFunction,
		regularisation_radius :: Real)

Multiple particles -> multiple velocity measurements

	particle_induced_velocity(
		inducing_particle_positions :: Matrix{<:Real},
		inducing_particle_vorticities :: Vector{<:Real},
		measurement_points :: Matrix{<:Real},
		kernel :: RegularisationFunction,
		regularisation_radius :: Real)

Compute the velocity induced in the flow field by vortex particles. The
third (multiple-multiple) method may be GPU accelerated.

# Arguments
- `inducing_particle_position` : Position of inducing particles
- `inducing_particle_vorticity` : Vorticity of inducing particles
- `mesurement_points` : Measurement points
- `kernel :: RegularisationFunction` : Regularisation function (winckelmans_regularisation() for example)
- `regularisation_radius :: Real` : Regularisation distance

In 3D, matrix arguments are expected to have size N by 3.
In 2D, matrix arguments are expected to have size N by 2.

The method will return velocity as:
- 3D - S2S: Vector length 3
- 3D - M2S: Vector length 3
- 3D - M2M: Matrix size N by 3
- 2D - S2S: Vector length 2
- 2D - M2S: Vector length 2
- 2D - M2M: Matrix size N by 2
"""
function particle_induced_velocity(
    inducing_particle_position :: Vector{<:Real},
    inducing_particle_vorticity :: Vector{<:Real},
	measurement_point :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec3f_vect(measurement_point, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

    inducing_particle = VortexParticle3D(
        inducing_particle_position,
        inducing_particle_vorticity, 0.0)
    mes_pnt = Vec3f(measurement_point)
	ret = Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_vel(
		const cvtx_P3D *self,
		const bsv_V3f mes_point,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ret = ccall(
			("cvtx_P3D_S2S_vel", libcvortex),
			Vec3f,
			(Ref{VortexParticle3D}, Vec3f, Ref{RegularisationFunction}, Cfloat),
			inducing_particle, mes_pnt, kernel, regularisation_radius
			)
	return [ret.x, ret.y, ret.z]
end

function particle_induced_velocity(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Matrix{<:Real},
	measurement_point :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec3f_vect(measurement_point, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(inducing_particle_position)[1]
    inducing_particles = map(
        i->VortexParticle3D(
            inducing_particle_position[i, :],
            inducing_particle_vorticity[i, :], 0.0),
        1:np)
    mes_pnt = Vec3f(measurement_point)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret =Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_M2S_vel(
		const cvtx_P3D **array_start,
		const int num_particles,
		const bsv_V3f mes_point,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ret = ccall(
			("cvtx_P3D_M2S_vel", libcvortex),
			Vec3f,
			(Ref{Ptr{VortexParticle3D}}, Cint, Vec3f,
				Ref{RegularisationFunction}, Cfloat),
			pargarr, np, mes_pnt, kernel,	regularisation_radius)
	return [ret.x, ret.y, ret.z]
end

function particle_induced_velocity(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Matrix{<:Real},
	measurement_points :: Matrix{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec3f_vect(measurement_points, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(inducing_particle_position)[1]
	ni = size(measurement_points)[1]
    inducing_particles = map(
        i->VortexParticle3D(
            inducing_particle_position[i, :],
            inducing_particle_vorticity[i, :], 0.0),
        1:np)
    mes_pnt = map(i->Vec3f(measurement_points[i,:]), 1:ni)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret = Vector{Vec3f}(undef, ni)
	#=
	CVTX_EXPORT void cvtx_P3D_M2M_vel(
		const cvtx_P3D **array_start,
		const int num_particles,
		const bsv_V3f *mes_start,
		const int num_mes,
		bsv_V3f *result_array,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ccall(
		("cvtx_P3D_M2M_vel", libcvortex),
		Cvoid,
		(Ptr{Ptr{VortexParticle3D}}, Cint, Ptr{Vec3f},
			Cint, Ref{Vec3f}, Ref{RegularisationFunction}, Cfloat),
		pargarr, np, mes_pnt, ni, ret, kernel, regularisation_radius)
	return Matrix{Float32}(ret)
end

"""

	particle_induced_dvort(
		inducing_particle_position :: Vector{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Vector{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

	particle_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Matrix{<:Real},
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Vector{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

	particle_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Matrix{<:Real},
		induced_particle_position :: Matrix{<:Real},
		induced_particle_vorticity :: Matrix{<:Real},
		kernel :: RegularisationFunction, regularisation_radius :: Real)

Rate of change of vorticity induced on vortex particles by vortex stretching.
The third multiple-multiple variant may be GPU accelerated.
In 2D, vortex stretching does not occur, so there is no 2D counterpart.

# Arguments
- `inducing_particle_position` : Position of inducing particles
- `inducing_particle_vorticity` : Vorticity of inducing particles
- `induced_particle_position` : Position of induced particles
- `induced_particle_vorticity` : Vorticity of induced particles
- `kernel :: RegularisationFunction` : Regularisation function
(VortFunc_winckelmans for example)
- `regularisation_radius :: Real` : Regularisation distance

Vector arguments are expected to have length 3. Matrix arguments are
expected to have size N by 3.
"""
function particle_induced_dvort(
    inducing_particle_position :: Vector{<:Real},
    inducing_particle_vorticity :: Vector{<:Real},
    induced_particle_position :: Vector{<:Real},
    induced_particle_vorticity :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity)
	convertable_to_F32(regularisation_radius, "regularisation_radius")

    inducing_particle = VortexParticle3D(
        inducing_particle_position,
        inducing_particle_vorticity, 0.0)
	induced_particle = VortexParticle3D(
		induced_particle_position,
		induced_particle_vorticity, 0.0)
	ret = Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_dvort(
		const cvtx_P3D *self,
		const cvtx_P3D *induced_particle,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ret = ccall(
			("cvtx_P3D_S2S_dvort", libcvortex),
			Vec3f,
			(Ref{VortexParticle3D}, Ref{VortexParticle3D},
				Ref{RegularisationFunction}, Cfloat),
			inducing_particle, induced_particle, kernel, regularisation_radius
			)
	return ret
end

function particle_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Matrix{<:Real},
    induced_particle_position :: Vector{<:Real},
    induced_particle_vorticity :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity)
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(induced_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], 0.0),
		1:np)
	induced_particle = VortexParticle3D(
		induced_particle_position,
		induced_particle_vorticity, 0.0)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret = Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_M2S_dvort(
		const cvtx_P3D **array_start,
		const int num_particles,
		const cvtx_P3D *induced_particle,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ret = ccall(
			("cvtx_P3D_M2S_dvort", libcvortex),
			Vec3f,
			(Ref{Ptr{VortexParticle3D}}, Cint, Ref{VortexParticle3D},
				Ref{RegularisationFunction}, Cfloat),
			pargarr, np, induced_particle, kernel, regularisation_radius
			)
	return Vector{Float32}(ret)
end

function particle_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Matrix{<:Real},
    induced_particle_position :: Matrix{<:Real},
    induced_particle_vorticity :: Matrix{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity)
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(inducing_particle_position)[1]
	ni = size(induced_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], 0.0),
		1:np)
	induced_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], 0.0),
		1:ni)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	indarg = Vector{Ptr{VortexParticle3D}}(undef, ni)
	for i = 1 : length(indarg)
		indarg[i] = Base.pointer(induced_particles, i)
	end
	ret = Vector{Vec3f}(undef, ni)
	#=
	CVTX_EXPORT void cvtx_P3D_M2M_dvort(
		const cvtx_P3D **array_start,
		const int num_particles,
		const cvtx_P3D **induced_start,
		const int num_induced,
		bsv_V3f *result_array,
		const cvtx_VortFunc *kernel,
		float regularisation_radius);
	=#
	ccall(
		("cvtx_P3D_M2M_dvort", libcvortex),
		Cvoid,
		(Ptr{Ptr{VortexParticle3D}}, Cint, Ptr{Ptr{VortexParticle3D}}, Cint,
			Ptr{Vec3f}, Ref{RegularisationFunction}, Cfloat),
		pargarr, length(inducing_particles), indarg, length(induced_particles),
			ret, kernel, regularisation_radius
		)
	return Matrix{Float32}(ret)
end

"""
# 3D functions

	particle_visc_induced_dvort(
		inducing_particle_position :: Vector{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		inducing_particle_volume :: Real,
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Vector{<:Real},
		induced_particle_volume :: Real,
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

	particle_visc_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Matrix{<:Real},
		inducing_particle_volume :: Vector{<:Real},
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Vector{<:Real},
		induced_particle_volume :: Real,
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

	particle_visc_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Matrix{<:Real},
		inducing_particle_volume :: Vector{<:Real},
		induced_particle_position :: Matrix{<:Real},
		induced_particle_vorticity :: Matrix{<:Real},
		induced_particle_volume :: Vector{<:Real},
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

# 2D functions

	particle_visc_induced_dvort(
		inducing_particle_position :: Vector{<:Real},
		inducing_particle_vorticity :: Real,
		inducing_particle_area :: Real,
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Real,
		induced_particle_area :: Real,
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

	particle_visc_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		inducing_particle_area :: Vector{<:Real},
		induced_particle_position :: Vector{<:Real},
		induced_particle_vorticity :: Real,
		induced_particle_area :: Real,
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

	particle_visc_induced_dvort(
		inducing_particle_position :: Matrix{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		inducing_particle_area :: Vector{<:Real},
		induced_particle_position :: Matrix{<:Real},
		induced_particle_vorticity :: Vector{<:Real},
		induced_particle_area :: Vector{<:Real},
		kernel :: RegularisationFunction,
		regularisation_radius :: Real,
		kinematic_visc :: Real)

The rate of change of vorticity induced on vortex particles by elements in the
flowfield due to viscosity.
The third multiple-multiple variant may be GPU accelerated.

# Arguments
- `inducing_particle_position` : Position of inducing particles
- `inducing_particle_vorticity` : Vorticity of inducing particles
- `inducing_particle_volume` : Volume of inducing particles
- `induced_particle_position` : Position of induced particles
- `induced_particle_vorticity` : Vorticity of induced particles
- `induced_particle_volume` : Volume of induced particles
- `kernel :: RegularisationFunction` : Regularisation function (winckelmans_regularisation() for example)
- `regularisation_radius :: Real` : Regularisation distance
- `kinematic_visc :: Real` : Kinematic viscosity
"""
function particle_visc_induced_dvort(
    inducing_particle_position :: Vector{<:Real},
	inducing_particle_vorticity :: Vector{<:Real},
	inducing_particle_volume :: Real,
    induced_particle_position :: Vector{<:Real},
	induced_particle_vorticity :: Vector{<:Real},
	induced_particle_volume :: Real,
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_volume)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity, induced_particle_volume)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_fn!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

    inducing_particle = VortexParticle3D(
        inducing_particle_position,
        inducing_particle_vorticity, induced_particle_volume)
	induced_particle = VortexParticle3D(
		induced_particle_position,
		induced_particle_vorticity, induced_particle_volume)
	ret = Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_visc_dvort(
		const cvtx_P3D *self,
		const cvtx_P3D *induced_particle,
		const cvtx_VortFunc *kernel,
		float regularisation_radius,
		float kinematic_visc);
	=#
	ret = ccall(
			("cvtx_P3D_S2S_visc_dvort", libcvortex),
			Vec3f,
			(Ref{VortexParticle3D}, Ref{VortexParticle3D},
				Ref{RegularisationFunction}, Cfloat, Cfloat),
			inducing_particle, induced_particle, kernel,
				regularisation_radius, kinematic_visc
			)
	return ret
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
	inducing_particle_vorticity :: Matrix{<:Real},
	inducing_particle_volume :: Vector{<:Real},
    induced_particle_position :: Vector{<:Real},
	induced_particle_vorticity :: Vector{<:Real},
	induced_particle_volume :: Real,
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_volume)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity, induced_particle_volume)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_fn!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	np = size(induced_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], inducing_particle_volume[i]),
		1:np)
	induced_particle = VortexParticle3D(
		induced_particle_position,
		induced_particle_vorticity, induced_particle_volume)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret = Vec3f(0., 0., 0.)
	#=
	CVTX_EXPORT bsv_V3f cvtx_P3D_M2S_visc_dvort(
		const cvtx_P3D **array_start,
		const int num_particles,
		const cvtx_P3D *induced_particle,
		const cvtx_VortFunc *kernel,
		float regularisation_radius,
		float kinematic_visc);
	=#
	ret = ccall(
			("cvtx_P3D_M2S_visc_dvort", libcvortex),
			Vec3f,
			(Ref{Ptr{VortexParticle3D}}, Cint, Ref{VortexParticle3D},
				Ref{RegularisationFunction}, Cfloat, Cfloat),
			pargarr, np, induced_particle, kernel,
				regularisation_radius, kinematic_visc
			)
	return Vector{Float32}(ret)
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
	inducing_particle_vorticity :: Matrix{<:Real},
	inducing_particle_volume :: Vector{<:Real},
    induced_particle_position :: Matrix{<:Real},
	induced_particle_vorticity :: Matrix{<:Real},
	induced_particle_volume :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real)

	check_particle_definition_3D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_volume)
	check_particle_definition_3D(induced_particle_position,
		induced_particle_vorticity, induced_particle_volume)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_3D!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	np = size(inducing_particle_position)[1]
	ni = size(induced_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], inducing_particle_volume[i]),
		1:np)
	induced_particles = map(
		i->VortexParticle3D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i, :], induced_particle_volume[i]),
		1:ni)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	indarg = Vector{Ptr{VortexParticle3D}}(undef, ni)
	for i = 1 : length(indarg)
		indarg[i] = Base.pointer(induced_particles, i)
	end
	ret = Vector{Vec3f}(undef, ni)
	#=
	CVTX_EXPORT void cvtx_P3D_M2M_visc_dvort(
		const cvtx_P3D **array_start,
		const int num_particles,
		const cvtx_P3D **induced_start,
		const int num_induced,
		bsv_V3f *result_array,
		const cvtx_VortFunc *kernel,
		float regularisation_radius,
		float kinematic_visc);
	=#
	ccall(
		("cvtx_P3D_M2M_visc_dvort", libcvortex),
		Cvoid,
		(Ptr{Ptr{VortexParticle3D}}, Cint, Ptr{Ptr{VortexParticle3D}}, Cint,
			Ptr{Vec3f}, Ref{RegularisationFunction}, Cfloat, Cfloat),
		pargarr, length(inducing_particles), indarg, length(induced_particles),
			ret, kernel, regularisation_radius, kinematic_visc
		)
	return Matrix{Float32}(ret)
end

"""
	redistribute_particles_on_grid(
		particle_positions :: Matrix{<:Real},
		particle_vorticities :: Matrix{<:Real},
		redistribution_function :: RedistributionFunction,
		grid_density :: Real;
		negligible_vort::Real=1e-4,
		max_new_particles::Integer=-1)

	redistribute_particles_on_grid(
		particle_positions :: Matrix{<:Real},
		inducing_particle_vorticity :: Vector{<:Real},
		redistribution_function :: RedistributionFunction,
		grid_density :: Real;
		negligible_vort::Real=1e-4,
		max_new_particles::Integer=-1)

Redistribute vorticity from vortex particles onto a regular grid.

# Arguments
- `particle_positions` : Position of particles to be redistributed.
- `particle_vorticities` : Vorticity of particles to be redistributed
- `redistribution_function :: RedistributionFunction` : Redistribution function (lambda3_redistribution() for example).
- `grid density :: Real` : The distance between newly placed particles on the grid.
- `negligible_vort :: Real` : The vorticity threshold at which to discard vortex particles. Given as abs(vorticity) as a proportion of the average particle's abs(vorticity)
- `max_new_particles::Integer` : The maximum number of particles that the grid can contain. -1 indicates that any number of particles can be created.

returns

	particle_positions, particle_vorts, particle_areas = redistribute_particles_on_grid(...)

where the
number of created particles may be less that the max_new_particles variable, even
when the number of particles is limited by the max_new_particles variable.
"""
function redistribute_particles_on_grid(
    particle_positions :: Matrix{<:Real},
	particle_vorticities :: Matrix{<:Real},
	redistribution_function :: RedistributionFunction,
	grid_density :: Real;
	negligible_vort::Real=1e-4,
	max_new_particles::Integer=-1)

	@assert(0 <= negligible_vort < 1, "The negligible_vort vort must be"*
		" in [0, 1). Given "*string(negligible_vort)*"." )
	@assert((max_new_particles==-1)||(0<=max_new_particles),
		"max_new_particles must be -1 (indicating any number of new particles"*
		") or a positive integer.")
	@assert(0<grid_density, "Grid density must be positive. Was "*
		string(grid_density)*".")
	check_particle_definition_3D(particle_positions,
		particle_vorticities)
	convertable_to_F32(grid_density, "grid_density")
	convertable_to_F32(negligible_vort, "negligible_vort")

	np = size(particle_positions)[1]
	inducing_particles = map(
		i->VortexParticle3D(
			particle_positions[i, :],
			particle_vorticities[i, :], 0.0),
		1:np)

	pargarr = Vector{Ptr{VortexParticle3D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	#=
	CVTX_EXPORT int cvtx_P3D_redistribute_on_grid(
		const cvtx_P3D **input_array_start,
		const int n_input_particles,
		cvtx_P3D *output_particles,		/* input is &(*cvtx_P3D) to write to */
		int max_output_particles,		/* Set to resultant num particles.   */
		const cvtx_RedistFunc *redistributor,
		float grid_density,
		float negligible_vort);
	=#
	println(Base.pointer(pargarr, 1))
	println(Base.pointer(inducing_particles, 1))
	if max_new_particles == -1
		max_new_particles = ccall(
			("cvtx_P3D_redistribute_on_grid", libcvortex),
			Cint,
			(Ptr{Ptr{VortexParticle3D}}, Cint, Ptr{VortexParticle3D},
				Cint, Ref{RedistributionFunction}, Cfloat, Cfloat),
			pargarr, np, C_NULL, 1, redistribution_function,
			grid_density, negligible_vort)
	end

	ret = Vector{VortexParticle3D}(undef, max_new_particles)
	nnp = ccall(
		("cvtx_P3D_redistribute_on_grid", libcvortex),
		Cint,
		(Ptr{Ptr{VortexParticle3D}}, Cint, Ptr{VortexParticle3D},
			Cint, Ref{RedistributionFunction}, Cfloat, Cfloat),
		pargarr, np, ret, max_new_particles, redistribution_function,
		grid_density, negligible_vort)

	# nnp is the number of new particles.
	nvorts = zeros(Float32, nnp, 3)
	nposns = zeros(Float32, nnp, 3)
	nareas = zeros(Float32, nnp)
	for i = 1 : nnp
		nvorts[i, :] = Vector{Float32}(ret[i].vorticity)
		nposns[i, :] = Vector{Float32}(ret[i].coord)
		nareas[i] = ret[i].volume
	end
	return nposns, nvorts, nareas
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# Path constructor
isym = ["a","b","c","d"];
osym = ["α","β","γ","δ"];
W = ProbabilityWeight{Float32};
p = Path(isym,osym,["a","b","c"] => ["α","β","γ"], one(W))
@test get_weight(p) == one(W)
FiniteStateTransducers.update_path!(p, 4, 4, 0.5)
@test get_weight(p) == W(0.5)
println(p)

p = Path(isym,["a","b","c"] => ["c","b","a"])
println(p)
@test typeofweight(p) == TropicalWeight{Float32}
@test get_weight(p) == one(TropicalWeight{Float32})

isym = ["a","b","c","d","e"];
osym = ["α","β","γ","δ","ε"];
W = ProbabilityWeight{Float32}
fst = WFST(isym,osym; W=W)
w = randn(6)
add_arc!(fst, 1, 2, "a", "α", w[1])
add_arc!(fst, 1, 3, "b", "β", w[2])
add_arc!(fst, 1, 4, "c", "γ", w[3])
add_arc!(fst, 2, 4, "d", "δ", w[4])
add_arc!(fst, 3, 4, "e", "ε", w[5])
initial!(fst, 1)
wf = rand()
final!(fst, 4, wf)
paths = collectpaths(fst)
for p in get_paths(fst)
  @test p in paths
end
@test length(paths) == 3
@test all( sort(get.(get_weight.(paths))) .≈ sort([w[1]*w[4]*wf, w[2]*w[5]*wf, w[3]*wf]) )

@test get_isequence(paths[1]) == ["c"]
@test get_isequence(paths[2]) == ["a","d"]
@test get_isequence(paths[3]) == ["b","e"]

@test get_osequence(paths[1]) == ["γ"]
@test get_osequence(paths[2]) == ["α","δ"]
@test get_osequence(paths[3]) == ["β","ε"]

# empty fst
fst = WFST(Dict(1=>1); W = W)
@test_throws ErrorException collectpaths(fst)
@test_throws ErrorException get_paths(fst)

# with eps
I = O = Int32
isym = Dict(i=>i     for i=1:5)
osym = Dict(100*i=>i for i=1:5)
fst = WFST(isym,osym; W = ProbabilityWeight{Float32})
w = randn(6)
ϵ = get_eps(I)
add_arc!(fst, 1, 2, 1, 100*1, w[1])
add_arc!(fst, 1, 3, 2,     ϵ, w[2])
add_arc!(fst, 1, 4, ϵ,     ϵ, w[3])
add_arc!(fst, 2, 4, 4, 100*4, w[4])
add_arc!(fst, 3, 4, 5, 100*5, w[5])
initial!(fst, 1)
wf = rand()
final!(fst, 4, wf)
paths = collectpaths(fst)
for p in get_paths(fst)
  @test p in paths
end
@test all( sort(get.(get_weight.(paths))) .≈ sort([w[1]*w[4]*wf, w[2]*w[5]*wf, w[3]*wf]) )


text:
"""
    LocationScale(μ,σ,ρ)

A location-scale transformed distribution with location parameter `μ`,
scale parameter `σ`, and given univariate distribution `ρ`.

If ``Z`` is a random variable with distribution `ρ`, then the distribution of the random
variable
```math
X = μ + σ Z
```
is the location-scale transformed distribution with location parameter `μ` and scale
parameter `σ`.

If `ρ` is a discrete distribution, the probability mass function of
the transformed distribution is given by
```math
P(X = x) = P\\left(Z = \\frac{x-μ}{σ} \\right).
```
If `ρ` is a continuous distribution, the probability density function of
the transformed distribution is given by
```math
f(x) = \\frac{1}{σ} ρ \\! \\left( \\frac{x-μ}{σ} \\right).
```

```julia
LocationScale(μ,σ,ρ) # location-scale transformed distribution
params(d)            # Get the parameters, i.e. (μ, σ, and the base distribution ρ)
location(d)          # Get the location parameter
scale(d)             # Get the scale parameter
```

External links
[Location-Scale family on Wikipedia](https://en.wikipedia.org/wiki/Location%E2%80%93scale_family)
"""
struct LocationScale{T<:Real, S<:ValueSupport, D<:UnivariateDistribution{S}} <: UnivariateDistribution{S}
    μ::T
    σ::T
    ρ::D
    function LocationScale{T,S,D}(μ::T, σ::T, ρ::D; check_args=true) where {T<:Real, S<:ValueSupport, D<:UnivariateDistribution{S}}
        check_args && @check_args(LocationScale, σ > zero(σ))
        new{T, S, D}(μ, σ, ρ)
    end
end

function LocationScale(μ::T, σ::T, ρ::UnivariateDistribution; check_args=true) where {T<:Real}
    _T = promote_type(eltype(ρ), T)
    D = typeof(ρ)
    S = value_support(D)
    return LocationScale{_T,S,D}(_T(μ), _T(σ), ρ; check_args=check_args)
end

LocationScale(μ::Real, σ::Real, ρ::UnivariateDistribution) = LocationScale(promote(μ, σ)..., ρ)

# aliases
const ContinuousLocationScale{T<:Real,D<:ContinuousUnivariateDistribution} = LocationScale{T,Continuous,D}
const DiscreteLocationScale{T<:Real,D<:DiscreteUnivariateDistribution} = LocationScale{T,Discrete,D}

Base.eltype(::Type{<:LocationScale{T}}) where T = T

minimum(d::LocationScale) = d.μ + d.σ * minimum(d.ρ)
maximum(d::LocationScale) = d.μ + d.σ * maximum(d.ρ)
support(d::LocationScale) = locationscale_support(d.μ, d.σ, support(d.ρ))
function locationscale_support(μ::Real, σ::Real, support::RealInterval)
    return RealInterval(μ + σ * support.lb, μ + σ * support.ub)
end
locationscale_support(μ::Real, σ::Real, support) = μ .+ σ .* support

LocationScale(μ::Real, σ::Real, d::LocationScale) = LocationScale(μ + d.μ * σ, σ * d.σ, d.ρ)

#### Conversions

convert(::Type{LocationScale{T}}, μ::Real, σ::Real, ρ::D) where {T<:Real, D<:UnivariateDistribution} = LocationScale(T(μ),T(σ),ρ)
convert(::Type{LocationScale{T}}, d::LocationScale{S}) where {T<:Real, S<:Real} = LocationScale(T(d.μ),T(d.σ),d.ρ, check_args=false)

#### Parameters

location(d::LocationScale) = d.μ
scale(d::LocationScale) = d.σ
params(d::LocationScale) = (d.μ,d.σ,d.ρ)
partype(::LocationScale{T}) where {T} = T

#### Statistics

mean(d::LocationScale) = d.μ + d.σ * mean(d.ρ)
median(d::LocationScale) = d.μ + d.σ * median(d.ρ)
mode(d::LocationScale) = d.μ + d.σ * mode(d.ρ)
modes(d::LocationScale) = d.μ .+ d.σ .* modes(d.ρ)

var(d::LocationScale) = d.σ^2 * var(d.ρ)
std(d::LocationScale) = d.σ * std(d.ρ)
skewness(d::LocationScale) = skewness(d.ρ)
kurtosis(d::LocationScale) = kurtosis(d.ρ)

isplatykurtic(d::LocationScale) = isplatykurtic(d.ρ)
isleptokurtic(d::LocationScale) = isleptokurtic(d.ρ)
ismesokurtic(d::LocationScale) = ismesokurtic(d.ρ)

entropy(d::ContinuousLocationScale) = entropy(d.ρ) + log(d.σ)
entropy(d::DiscreteLocationScale) = entropy(d.ρ)

mgf(d::LocationScale,t::Real) = exp(d.μ*t) * mgf(d.ρ,d.σ*t)

#### Evaluation & Sampling

pdf(d::ContinuousLocationScale,x::Real) = pdf(d.ρ,(x-d.μ)/d.σ) / d.σ
pdf(d::DiscreteLocationScale, x::Real) = pdf(d.ρ,(x-d.μ)/d.σ)

logpdf(d::ContinuousLocationScale,x::Real) = logpdf(d.ρ,(x-d.μ)/d.σ) - log(d.σ)
logpdf(d::DiscreteLocationScale, x::Real) = logpdf(d.ρ,(x-d.μ)/d.σ)

# additional definitions are required to fix ambiguity errors and incorrect defaults
for f in (:cdf, :ccdf, :logcdf, :logccdf)
    _f = Symbol(:_, f)
    @eval begin
        $f(d::LocationScale, x::Real) = $_f(d, x)
        $f(d::DiscreteLocationScale, x::Real) = $_f(d, x)
        $f(d::DiscreteLocationScale, x::Integer) = $_f(d, x)
        $_f(d::LocationScale, x::Real) = $f(d.ρ, (x - d.μ) / d.σ)
    end
end

quantile(d::LocationScale,q::Real) = d.μ + d.σ * quantile(d.ρ,q)

rand(rng::AbstractRNG, d::LocationScale) = d.μ + d.σ * rand(rng, d.ρ)
cf(d::LocationScale, t::Real) = cf(d.ρ,t*d.σ) * exp(1im*t*d.μ)
gradlogpdf(d::ContinuousLocationScale, x::Real) = gradlogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ

#### Syntactic sugar for simple transforms of distributions, e.g., d + x, d - x, and so on

Base.:+(d::UnivariateDistribution, x::Real) = LocationScale(x, one(x), d)
Base.:+(x::Real, d::UnivariateDistribution) = d + x
Base.:*(x::Real, d::UnivariateDistribution) = LocationScale(zero(x), x, d)
Base.:*(d::UnivariateDistribution, x::Real) = x * d
Base.:-(d::UnivariateDistribution, x::Real) = d + -x
Base.:/(d::UnivariateDistribution, x::Real) = inv(x) * d


text:
#####
##### Calculate the areas, volumes and spacings of different kind of `AbstractGrid`
##### Currently for RectilinearGrid and LatLonGrid
#####

#####
##### RectilinearGrid, halo points included
#####
@inline ΔxC(i, j, k, g::RectilinearGrid) = g.Δx
@inline ΔyC(i, j, k, g::RectilinearGrid) = g.Δy
@inline ΔzC(i, j, k, g::RectilinearGrid) = @inbounds g.dzC[k]

@inline ΔxF(i, j, k, g::RectilinearGrid) = g.Δx
@inline ΔyF(i, j, k, g::RectilinearGrid) = g.Δy
@inline ΔzF(i, j, k, g::RectilinearGrid) = @inbounds g.dzF[k]

@inline Ax(i, j, k, g::RectilinearGrid) = @inbounds g.Δy * g.dzF[k]
@inline Ay(i, j, k, g::RectilinearGrid) = @inbounds g.Δx * g.dzF[k]
@inline Az(i, j, k, g::RectilinearGrid) = g.Δx * g.Δy

@inline volume(i, j, k, g::RectilinearGrid) = @inbounds g.Δx * g.Δy * g.dzF[k]

#####
##### LatLonGrid (degree to meter), halo points included
#####
@inline ΔxC(i, j, k, g::LatLonGrid) = @inbounds g.dxC[i,j]
@inline ΔyC(i, j, k, g::LatLonGrid) = @inbounds g.dyC[i,j]
@inline ΔzC(i, j, k, g::LatLonGrid) = @inbounds g.dzC[k]

@inline ΔxF(i, j, k, g::LatLonGrid) = @inbounds g.dxF[i,j]
@inline ΔyF(i, j, k, g::LatLonGrid) = @inbounds g.dyF[i,j]
@inline ΔzF(i, j, k, g::LatLonGrid) = @inbounds g.dzF[k]

@inline Ax(i, j, k, g::LatLonGrid) = @inbounds g.Ax[i,j,k]
@inline Ay(i, j, k, g::LatLonGrid) = @inbounds g.Ay[i,j,k]
@inline Az(i, j, k, g::LatLonGrid) = @inbounds g.Az[i,j]

@inline volume(i, j, k, g::LatLonGrid) = @inbounds g.Vol[i,j,k]


text:
using NearestNeighbors
using StaticArrays
"""
    Fkmnn(xy, r, window, N)

Return empty space function based on Kaplan-Meier estimator.
"""
function Fkmnn(xy, window, r, N)
    if length(xy) == 0
        return zeros(length(r))
    end
    xy = @. SVector(getx(xy), gety(xy))
    (x1,x2), (y1, y2) = window.x, window.y
    tree = KDTree(xy)
    sx = (x2-x1)/N
    Ny = round(Int, N*(y2-y1)/(x2-x1))
    sy = (y2-y1)/Ny
    rx = range(x1 + 0.5sx, x2 - 0.5sx, length=N)
    ry = range(y1 + 0.5sy, y2 - 0.5sy, length=Ny)
    grid = vec([SVector(x, y) for x in rx, y in ry])
    _, D = knn(tree, grid, 1)
    D = first.(D)
    B = [bdist(window, x) for x in grid]
    ind = D .<= B
    T = min.(D, B)
    km(ind, T, r)
end


text:
module TestSingletonContainers

using BangBang.NoBang: SingletonDict, SingletonVector
using Tables: Tables
using Test

function test_shows(x)
    # smoke tests
    @test sprint(show, x) isa String
    @test sprint(show, "text/plain", x) isa String
    @test string(x) isa String
end

@testset "SingletonVector" begin
    v = SingletonVector(1)
    @test v[1] == 1
    @test_throws BoundsError v[2]
    @test size(v) == (1,)
    @test collect(v)::Vector{Int} == [1]
    test_shows(v)
end

@testset "SingletonVector as a table" begin
    @testset for T in [
        SingletonVector,
        SingletonVector{Int},
        SingletonVector{<:NamedTuple},
        SingletonVector{NamedTuple},
    ]
        @test !Tables.istable(T)
    end

    @testset for T in [
        SingletonVector{<:NamedTuple{(:a,)}},
        SingletonVector{<:NamedTuple{(:a,),Tuple{Int}}},
    ]
        @test Tables.istable(T)
        @test Tables.rowaccess(T)
        @test Tables.columnaccess(T)
    end

    @testset begin
        @test Tables.rows(SingletonVector((a = 1,))) == [(a = 1,)]
        @test Tables.columns(SingletonVector((a = 1,))) == (a = [1],)
    end
end

@testset "SingletonDict" begin
    d = SingletonDict(:a, 1)
    @test d[:a] == 1
    @test_throws BoundsError d[:non_existing_key]
    @test collect(d)::Vector{Pair{Symbol,Int}} == [:a => 1]
    test_shows(d)
end

end  # module


text:
using JLD
using NetCDF

function Array2NetCDF(varname::String, outputname::String, x::AbstractArray{tp, 3}; overwrite::Bool = false) where {tp}
  if overwrite
    if isfile(outputname)
    rm(outputname)
    end
  end
  nccreate(outputname, varname,"dim1_size$(size(x,1))",size(x,1),"dim2_size$(size(x,2))",size(x, 2),"dim3_size$(size(x,3))",size(x,3),t=eltype(x))
  ncwrite(x,outputname, varname)
  ncclose(outputname)
end

function Array2NetCDF(varname::String, outputname::String, x::AbstractArray; overwrite::Bool = false)
  if overwrite
    if isfile(outputname)
      rm(outputname)
    end
  end
  nccreate(outputname,varname,string(varname,"_array__length$(length(x))"),length(x),t=eltype(x))
  ncwrite(reshape(x, length(x)),outputname,varname)
  ncclose(outputname)
end

function Array2NetCDF(varname::String, outputname::String, x::Array{Bool, 3}; overwrite::Bool = false)
  if overwrite
    if isfile(outputname)
      rm(outputname)
    end
  end
  nccreate(outputname, varname,"dim1_size$(size(x,1))",size(x,1),"dim2_size$(size(x,2))",size(x, 2),"dim3_size$(size(x,3))",size(x,3),t=Int32)
  ncwrite(broadcast(Int32, x),outputname,varname)
  ncclose(outputname)
end

function Array2NetCDF(varname::String, outputname::String, x::StepRangeLen; overwrite::Bool = false)
  if overwrite
    if isfile(outputname)
      rm(outputname)
    end
  end
  nccreate(outputname,varname,string(varname,"_array_length$(length(x))"),length(x),t=eltype(x))
  ncwrite(collect(x),outputname,varname)
  ncclose(outputname)
end

function Array2NetCDF(varname::String, outputname::String, x::tp; overwrite::Bool = false) where {tp}
  if overwrite
    if isfile(outputname)
      rm(outputname)
    end
  end
  nccreate(outputname,varname,string(varname,"_array_length$(length(x))"),length(x),t=eltype(x))
  ncwrite([x],outputname,varname)
  ncclose(outputname)
end

function CCobj2NetCDF(outputname, CCobj; overwrite = true, endskip::Int = 3)
  if overwrite
    if isfile(outputname)
      rm(outputname)
    end
  end
  for i = 1:(length(propertynames(CCobj))-endskip)
    varname = string(propertynames(CCobj)[i])
    print("\n write $(varname) ...")
    x = getfield(CCobj, propertynames(CCobj)[i])
    Array2NetCDF(varname, outputname, x, overwrite = false)
  end
  return "-- done -- \n"
end

export CCobj2NetCDF

function NetCDF2CCobj(inputname; endskip::Int = 3)
  CCobj = init_CCobj()
  NetCDF2CCobj!(CCobj, inputname, endskip = endskip)
  return CCobj
end

export NetCDF2CCobj

function NetCDF2CCobjVarStats(inputname, CCobj::CCstats, variable_name; endskip::Int = 1, startskip::Int = 1)
  CCobj.VarStats[variable_name] = init_VariableStats()
  NetCDF2CCobj!(CCobj.VarStats[variable_name], inputname, endskip = endskip, startskip = startskip)
  return CCobj
end

export NetCDF2CCobjVarStats

function NetCDF2CCobj!(CCobj, inputname; endskip::Int = 3, startskip::Int = 0 )
  for i = (1+startskip):(length(propertynames(CCobj))-endskip)
    varname = string(propertynames(CCobj)[i])
    print("\n read $(varname) ...")
    x = ncread(inputname, varname)
    if length(x) == 1 || varname == "Latitudes" || varname == "Longitudes"
        x = broadcast(eltype(getfield(CCobj, propertynames(CCobj)[i])), x)
    else
        x = convert(typeof(getfield(CCobj, propertynames(CCobj)[i])), x)
    end
    if length(x) == 1
      setfield!(CCobj, propertynames(CCobj)[i], x[1])
    elseif varname == "Latitudes"
      setfield!(CCobj, propertynames(CCobj)[i], maximum(x):round(x[2]-x[1], digits=5):minimum(x))
    elseif varname == "Longitudes"
      setfield!(CCobj, propertynames(CCobj)[i], minimum(x):round(x[2]-x[1], digits=5):maximum(x))
    else
      setfield!(CCobj, propertynames(CCobj)[i], x)
    end
  end
  ncclose()
  return CCobj
end

export NetCDF2CCobj!