text:
"""
    mask = flood(f, src, idx, nbrhood_function=diamond_iterator((3,3,...)))

Return an array `mask` with the same axes as `src`, marked `true` for all elements of `src` that:

- satisfy `f(src[i]) == true` and
- are connected by such elements to the starting point `idx` (an integer index or `CartesianIndex`).

This throws an error if `f` evaluates as `false` for the starting value `src[idx]`.
The sense of connectivity is defined by `nbrhood_function`, with two choices being
[`ImageSegmentation.diamond_iterator`](@ref) and [`ImageSegmentation.box_iterator`](@ref.)

# Examples

```jldoctest; setup=:(using ImageSegmentation, ImageCore)
julia> mostly_red(c) = red(c) > green(c) && red(c) > blue(c)
mostly_red (generic function with 1 method)

julia> img = repeat(LinRange(colorant"red", colorant"blue", 4), 1, 2) # red-to-blue
4×2 Array{RGB{Float32},2} with eltype RGB{Float32}:
 RGB{Float32}(1.0,0.0,0.0)            RGB{Float32}(1.0,0.0,0.0)
 RGB{Float32}(0.666667,0.0,0.333333)  RGB{Float32}(0.666667,0.0,0.333333)
 RGB{Float32}(0.333333,0.0,0.666667)  RGB{Float32}(0.333333,0.0,0.666667)
 RGB{Float32}(0.0,0.0,1.0)            RGB{Float32}(0.0,0.0,1.0)

julia> flood(mostly_red, [img; img], 1)   # only first copy of `img` is connected
8×2 BitMatrix:
 1  1
 1  1
 0  0
 0  0
 0  0
 0  0
 0  0
 0  0
```

See also [`flood_fill!`](@ref).
"""
flood(f, src::AbstractArray, idx::Union{Integer,CartesianIndex}, nbrhood_function=diamond_iterator(window_neighbors(src))) =
    flood_fill!(f, falses(axes(src)) #=fill!(similar(src, Bool), false)=#, src, idx, nbrhood_function)

"""
    mask = flood(src, idx, nbrhood_function=diamond_iterator((3,3,...)); thresh)

Return an array `mask` with the same axes as `src`, marked `true` for all connected elements
of `src` for which `src[i]` has a Euclidean distance less than `thresh` from `src[idx]`.
"""
function flood(src::AbstractArray, idx::Union{Integer,CartesianIndex},
                    nbrhood_function=diamond_iterator(window_neighbors(src)); thresh)
    validx = src[idx]
    validx = accum_type(typeof(validx))(validx)
    return let validx=validx
        flood(val -> default_diff_fn(val, validx) < thresh, src, idx, nbrhood_function)
    end
end

"""
    flood_fill!(f, dest, src, idx, nbrhood_function=diamond_iterator((3,3,...)); fillvalue=true, isfilled = isequal(fillvalue))

Set entries of `dest` to `fillvalue` for all elements of `src` that:

- satisfy `f(src[i]) == true` and
- are connected by such elements to the starting point `idx` (an integer index or `CartesianIndex`).

This throws an error if `f` evaluates as `false` for the starting value `src[idx]`.
The sense of connectivity is defined by `nbrhood_function`, with two choices being
[`ImageSegmentation.diamond_iterator`](@ref) and [`ImageSegmentation.box_iterator`](@ref.)

You can optionally omit `dest`, in which case entries in `src` will be set to `fillvalue`.
This is recommended only in cases where `fillvalue` is distinct from the values in `src`,
as otherwise there is a risk that some voxels will be incorrectly identified as having been
visited and their neighbors excluded from further search.
You may also supply `isfilled`, which should return `true` for any value in `dest`
which does not need to be set or visited; one requirement is that `isfilled(fillvalue) == true`.

# Examples

```jldoctest; setup=:(using ImageSegmentation)
julia> a = repeat([1:4; 1:4], 1, 3)
8×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3
 4  4  4
 1  1  1
 2  2  2
 3  3  3
 4  4  4

julia> flood_fill!(>=(3), a, CartesianIndex(3, 2); fillvalue = -1, isfilled = <(0))
8×3 Matrix{Int64}:
  1   1   1
  2   2   2
 -1  -1  -1
 -1  -1  -1
  1   1   1
  2   2   2
  3   3   3
  4   4   4
```

See also [`flood`](@ref).
"""
function flood_fill!(f,
                     dest,
                     src::AbstractArray,
                     idx::Union{Int,CartesianIndex},
                     nbrhood_function = diamond_iterator(window_neighbors(src));
                     fillvalue = true,
                     isfilled = (fillvalue === true) ? identity : isequal(fillvalue))
    R = CartesianIndices(src)
    idx = R[idx]  # ensure cartesian indexing
    f(src[idx]) || throw(ArgumentError("starting point fails to meet criterion"))
    q = [idx]
    fillvalue = convert(eltype(dest), fillvalue)
    axes(dest) == R.indices || throw(DimensionMismatch("$(axes(dest)) do not match $(Tuple(R))"))
    _flood_fill!(f, dest, src, R, q, nbrhood_function, fillvalue, isfilled)
    return dest
end
flood_fill!(f, dest, src::AbstractArray, idx::Integer, args...; kwargs...) =
    flood_fill!(f, dest, src, Int(idx)::Int, args...; kwargs...)
flood_fill!(f, src::AbstractArray, idx::Union{Integer,CartesianIndex}, args...; kwargs...) =
    flood_fill!(f, src, src, idx, args...; kwargs...)

# This is a trivial implementation (just to get something working), better would be a raster implementation
function _flood_fill!(f::F, dest, src, R::CartesianIndices{N}, q, nbrhood_function::FN, fillvalue, isfilled::C) where {F,N,FN,C}
    isfilled(fillvalue) == true || throw(ArgumentError("`isfilled(fillvalue)` must return `true`"))
    if Base.mightalias(dest, src)
        try
            f(fillvalue) && @warn """
                in-place `fillvalue` may not be distinct from array values, and filling may be incomplete.
                To suppress this warning, avoid filling in-place or ensure `f(fillvalue)` is `false`."""
        catch
        end
    end
    while !isempty(q)
        idx = pop!(q)
        dest[idx] = fillvalue
        @inbounds for j in nbrhood_function(idx)
            j ∈ R || continue
            if f(src[j]) && !isfilled(dest[j])
                push!(q, j)
            end
        end
    end
end


text:
struct DiffEqNLSolveTag end

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
    du::T
    dual_du::S
end

Base.@pure function DiffCache(T, size, ::Type{Val{chunk_size}}) where chunk_size
    DiffCache(fill(zero(T), size...), fill(zero(Dual{typeof(ForwardDiff.Tag(DiffEqNLSolveTag(),T)),T,chunk_size}), size...))
end

Base.@pure DiffCache(u::AbstractArray) = DiffCache(eltype(u),size(u),Val{ForwardDiff.pickchunksize(length(u))})
Base.@pure DiffCache(u::AbstractArray,nlsolve) = DiffCache(eltype(u),size(u),Val{get_chunksize(nlsolve)})
Base.@pure DiffCache(u::AbstractArray,T::Type{Val{CS}}) where {CS} = DiffCache(eltype(u),size(u),T)

get_du(dc::DiffCache, ::Type{T}) where {T<:Dual} = dc.dual_du
get_du(dc::DiffCache, T) = dc.du

# Default nlsolve behavior, should move to DiffEqDiffTools.jl

Base.@pure determine_chunksize(u,alg::DEAlgorithm) = determine_chunksize(u,get_chunksize(alg))
Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

struct NLSOLVEJL_SETUP{CS,AD} end
Base.@pure NLSOLVEJL_SETUP(;chunk_size=0,autodiff=true) = NLSOLVEJL_SETUP{chunk_size,autodiff}()
(::NLSOLVEJL_SETUP)(f,u0; kwargs...) = (res=NLsolve.nlsolve(f,u0; kwargs...); res.zero)
function (p::NLSOLVEJL_SETUP{CS,AD})(::Type{Val{:init}},f,u0_prototype) where {CS,AD}
  AD ? autodiff = :forward : autodiff = :central
  OnceDifferentiable(f, u0_prototype, u0_prototype, autodiff,
                     ForwardDiff.Chunk(determine_chunksize(u0_prototype,CS)))
end

get_chunksize(x) = 0
get_chunksize(x::NLSOLVEJL_SETUP{CS,AD}) where {CS,AD} = CS

"""
    calculate_residuals(ũ, u₀, u₁, α, ρ, scalarnorm)

Return element-wise residuals
```math
\\frac{ũ}{α+\\max{scalarnorm(u₀),scalarnorm(u₁)}*ρ}.
```
"""
@inline @muladd function calculate_residuals(ũ::Number, u₀::Number, u₁::Number, α::Real,
                                             ρ::Real, scalarnorm)
    ũ / (α + max(scalarnorm(u₀), scalarnorm(u₁)) * ρ)
end

"""
    calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)

Return element-wise residuals
```math
\\frac{δ E₁ + E₂}{α+\\max{scalarnorm(u₀),scalarnorm(u₁)}*ρ}.
```
"""
@inline @muladd function calculate_residuals(E₁::Number, E₂::Number, u₀::Number, u₁::Number,
                                             α::Real, ρ::Real, δ::Number, scalarnorm)
    (δ * E₁ + E₂) / (α + max(scalarnorm(u₀), scalarnorm(u₁)) * ρ)
end

"""
    calculate_residuals!(out, ũ, u₀, u₁, α, ρ, scalarnorm)

Same as [`calculate_residuals`](@ref) but save result in `out`.
"""
@inline function calculate_residuals!(out, ũ, u₀, u₁, α, ρ, scalarnorm)
  @. out = calculate_residuals(ũ, u₀, u₁, α, ρ, scalarnorm)
  out
end

@inline function calculate_residuals!(out::Array{<:Number}, ũ::Array{<:Number},
                                      u₀::Array{<:Number}, u₁::Array{<:Number}, α::Real,
                                      ρ::Real, scalarnorm)
  @tight_loop_macros @inbounds for i in eachindex(out)
    out[i] = calculate_residuals(ũ[i], u₀[i], u₁[i], α, ρ, scalarnorm)
  end
  out
end

@inline function calculate_residuals(ũ, u₀, u₁, α, ρ, scalarnorm)
  @. calculate_residuals(ũ, u₀, u₁, α, ρ, scalarnorm)
end

@inline function calculate_residuals(ũ::Array{<:Number}, u₀::Array{<:Number},
                                     u₁::Array{<:Number}, α::Real, ρ::Real, scalarnorm)
  out = similar(ũ)
  calculate_residuals!(out, ũ, u₀, u₁, α, ρ, scalarnorm)
  out
end

"""
    calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)

Same as [`calculate_residuals`](@ref) but save result in `out`.
"""
@inline function calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)
  @. out = calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)
  out
end

@inline function calculate_residuals!(out::Array{<:Number}, E₁::Array{<:Number},
                                      E₂::Array{<:Number}, u₀::Array{<:Number},
                                      u₁::Array{<:Number}, α::Real, ρ::Real, δ::Number,
                                      scalarnorm)
  @tight_loop_macros @inbounds for i in eachindex(out)
      out[i] = calculate_residuals(E₁[i], E₂[i], u₀[i], u₁[i], α, ρ, δ, scalarnorm)
  end
  out
end

@inline function calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)
  @. calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)
end

@inline function calculate_residuals(E₁::Array{<:Number}, E₂::Array{<:Number},
                                     u₀::Array{<:Number}, u₁::Array{<:Number}, α::Real,
                                     ρ::Real, δ::Number, scalarnorm)
    out = similar(u₀)
    calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)
    out
end

macro cache(expr)
  name = expr.args[2].args[1].args[1]
  fields = expr.args[3].args[2:2:end]
  cache_vars = Expr[]
  rand_vars = Expr[]
  jac_vars = Pair{Symbol,Expr}[]
  ratenoise_vars = Expr[]
  for x in fields
    if x.args[2] == :uType || x.args[2] == :rateType ||
       x.args[2] == :kType || x.args[2] == :uNoUnitsType #|| x.args[2] == :possibleRateType
      push!(cache_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :JCType
      push!(cache_vars,:(c.$(x.args[1]).duals...))
    elseif x.args[2] == :GCType
      push!(cache_vars,:(c.$(x.args[1]).duals))
    elseif x.args[2] == :DiffCacheType
      push!(cache_vars,:(c.$(x.args[1]).du))
      push!(cache_vars,:(c.$(x.args[1]).dual_du))
    elseif x.args[2] == :JType || x.args[2] == :WType
      push!(jac_vars,x.args[1] => :(c.$(x.args[1])))
    elseif x.args[2] == :randType
      push!(rand_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :rateNoiseType || x.args[2] == :rateNoiseCollectionType
      # Should be a pair for handling non-diagonal
      push!(ratenoise_vars,:(c.$(x.args[1])))
    end
  end
  quote
    $expr
    $(esc(:full_cache))(c::$name) = tuple($(cache_vars...))
    $(esc(:jac_iter))($(esc(:c))::$name) = tuple($(jac_vars...))
    $(esc(:rand_cache))($(esc(:c))::$name) = tuple($(rand_vars...))
    $(esc(:ratenoise_cache))($(esc(:c))::$name) = tuple($(ratenoise_vars...))
  end
end

_reshape(v, siz) = reshape(v, siz)
_reshape(v::Number, siz) = v
_vec(v) = vec(v)
_vec(v::Number) = v

ext: jl


text:
@testset "DataStructures" begin
    @testset "BasicBinaryTree" begin
        @testset "BasicBinaryTree: BinaryTree" begin
            tree = BinaryTree{Int}(10, 0)
            insert!(tree, 1, 10)
            insert!(tree, 1, 20)
            insert!(tree, 2, 50)
            #=
            [1] 0
             |       \
            [2] 10  [3] 20
             |
            [4] 50
            =#
            @test tree.n == 4
            @test depth(tree, 4) == 2
            @test isleaf(tree, 2) == false
            @test ch(tree, 1, true) == 3
            @test tree.val[left(tree)] == 50
            @test height(tree) == 2
        end
        @testset "DisjointSet: DisjointSet" begin
            set = DisjointSet(10)
            set.par[4] = 1
            set.par[3] = 1
            set.par[2] = 3
            set.par[5] = 8
            @test find(set, 2) == 1
            merge!(set, 5, 4)
            @test find(set, 5) == find(set, 2)
        end

        @testset "LinkedList" begin
            list = LinkedList.create_list(3, 13)
            @test LinkedList.return_as_array(list) == [13, 13, 13]

            LinkedList.insert(list, 5, 2)
            LinkedList.push_back(list, 8)
            LinkedList.insert(list, 7, 1)
            LinkedList.insert(list, 7, 6)
            @test LinkedList.return_as_array(list) == [7, 13, 5, 13, 13, 7, 8]

            LinkedList.remove(list, 7)
            @test LinkedList.return_as_array(list) == [7, 13, 5, 13, 13, 7]

            LinkedList.remove_all(list, 13)
            @test LinkedList.return_as_array(list) == [7, 5, 7]

            @test LinkedList.indexOf(list, 5) == 2
            @test LinkedList.indexOf(list, 3) == -1

            LinkedList.remove_first(list, 7)
            @test LinkedList.return_as_array(list) == [5, 7]

            @test ismissing(LinkedList.get(list, 3)) == true

            LinkedList.clear(list)
            @test LinkedList.is_empty(list) == true
        end
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
export pcl_impose_Dirichlet_boundary_conditions, pcl_compute_fem_laplace_matrix1
@doc raw"""
    pcl_impose_Dirichlet_boundary_conditions(indices::Array{Int64, 2}, bdnode::Array{Int64, 1}, outdof::Int64)

Computes the Jacobian matrix for [`impose_Dirichlet_boundary_conditions`](@ref). Assume that `impose_Dirichlet_boundary_conditions`
transforms a sparse matrix `A` to `B`, and `v_A` and `v_B` are the nonzero entries, this function computes

$$J_{ij} = \frac{\partial (v_B)_j}{\partial (v_A)_i}$$

- `indices`: the **1-based** $n_A\times 2$ index array for `A`;
- `bdnode`: boundary DOF (the same as inputs to `impose_Dirichlet_boundary_conditions`)
- `outdof`: the number of nonzero entries in `B`, i.e., length of `v_B`
"""
function pcl_impose_Dirichlet_boundary_conditions(indices::Array{Int64, 2}, bdnode::Array{Int64, 1}, outdof::Int64)
    bdN = Int32(length(bdnode))
    sN = Int32(size(indices, 1))
    J = zeros(sN, outdof)
    @eval ccall((:pcl_ImposeDirichlet, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Clonglong}, Ptr{Clonglong}, Cint, Cint),
        $J, $indices, $bdnode, $bdN, $sN)
    return J
end

@doc raw"""
    pcl_compute_fem_laplace_matrix1(mmesh::Mesh)

Computes the Jacobian matrix for [`pcl_compute_fem_laplace_matrix1`](@ref). Assume we contruct a sparse matrix via
```julia
A = compute_fem_laplace_matrix1(κ, mmesh)
v = A.o.values
```
Then the function returns

$$J = \frac{\partial v_j}{\partial κ_i}$$
"""
function pcl_compute_fem_laplace_matrix1(mmesh::Mesh)
    J = zeros(get_ngauss(mmesh), mmesh.elem_ndof^2 * get_ngauss(mmesh))
    @eval ccall((:pcl_FemLaplaceScalar_Jacobian, $LIBMFEM), Cvoid, (Ptr{Cdouble},),$J)
    J
end


text:
###################################################
## geo/stats/stats.jl
## Data type to save geotrip stats
###################################################

"""
`GeoStats` : abstract type that stores statistics about an algorithm run
must implement attributes:
- `name` : some information about the method used to generate this
- `sdict` : dictionary mapping stat names to relevant values
must implement methods:
- `printStats` : print a summary of statistics on the algorithm run
"""

abstract type GeoStats end
function Base.show(io::IO, so::GeoStats)
	typeName = split(string(typeof(so)),".")[end]
	println(io, "$(typeName): $(so.name)")
end

"""
	`RealGeoStats`: instance of GeoStats, made for real trip data
"""
mutable struct RealGeoStats <: GeoStats
	name::AbstractString
	sdict::Dict{AbstractString, Union{Float64, Array{Float64,1}}}

	"array containing upper bounds for time breakdown of stats"
	timeBound::Array{Float64,1}
	function RealGeoStats(name::AbstractString, gt::GeoTimings, ds::DataSplit,
		timeBound::Vector{Float64} = [270., 450., 720., 900., 100_000.])
		obj = new()
		obj.name = name
		obj.timeBound = timeBound
		times = Float64[estimateTime(gt, ID) for ID in eachindex(gt.trips)]
		obj.sdict = Dict{AbstractString, Union{Float64, Array{Float64,1}}}(
			"testTripsLogError" => testTripsLogError(gt, ds, times),
			"trTripsLogError" => trTripsLogError(gt, ds, times),
			"testTripsLogBias" => testTripsLogBias(gt, ds, times),
			"trTripsLogBias" => trTripsLogBias(gt, ds, times),
			# "testTripsLogErrorbt" => 100 * testTripsLogErrorByTime(gt, ds, times, obj.timeBound),
			# "trTripsLogErrorbt" => 100 * trTripsLogErrorByTime(gt, ds, times, obj.timeBound),
			# "testTripsLogBiasbt" => 100 * testTripsLogBiasByTime(gt, ds, times, obj.timeBound),
			# "trTripsLogBiasbt" => 100 * trTripsLogBiasByTime(gt, ds, times, obj.timeBound)
		)
		return obj
	end
end

"""
	`printStats`: prints out all statistics identified by statName stored in given GeoStats
	Optional argument outputFileName if you want to write this to a file instead
"""
function printStats(so::RealGeoStats; outputFileName = "")
	if outputFileName == ""
		println(so)
		for statName in sort(collect(keys(so.sdict)))
			if occursin("bt", statName)
				println(statName, ":")
				for (i,time) in enumerate(so.timeBound)
					if occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
						@printf("\t\tTime < %.0f s:\t%.0fs\n", time, so.sdict[statName][i])
					else
						@printf("\t\tTime < %.0f s:\t%.5f\n", time, so.sdict[statName][i])
					end
				end
			elseif occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
				print(statName, ":\t")
				@printf("%.2f s\n", so.sdict[statName])
			else
				print(statName, ":\t")
				@printf("%.5f\n", so.sdict[statName])
			end
		end
	else
		f = open(outputFileName, "w")
		write(f, so)
		for statName in sort(collect(keys(so.sdict)))
			if occursin("bt", statName)
				for (i,time) in enumerate(so.timeBound)
					if occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
						@printf(f, "\t\tTime < %.0f s:\t%.0fs\n", time, so.sdict[statName][i])
					else
						@printf(f, "\t\tTime < %.0f s:\t%.5f\n", time, so.sdict[statName][i])
					end
				end
			elseif occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
				write(f, string(statName, ":\t"))
				@printf(f, "%.0fs\n", so.sdict[statName])
			else
				write(f, string(statName, ":\t"))
				@printf(f, "%.5f\n", so.sdict[statName])
			end
		end
		close(f)
	end
end

"""
	`printStats`: print relevant statistics, identified by statName, stored in given GeoStats object
"""
function printStats(stats::Vector{RealGeoStats}, statName::AbstractString)
	# check if stat is valid and if timebounds are same
	for so in stats
		if !(statName in collect(keys(so.sdict)))
			error("Statistic not found")
		end
		if occursin("bt", statName)
			if so.timeBound != stats[1].timeBound
				error("Time breakdown arrays don't match")
			end
		end
	end
	# set up top row
	println(statName)
	if occursin("bt", statName)
		print("\t\t")
		for time in stats[1].timeBound
			@printf("< %.0fs\t", time)
		end
		print("\n")
	end
	for so in stats
		if occursin("bt", statName)
			print(so.name, "\t")
			for (i,time) in enumerate(so.timeBound)
				if occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
					@printf("%.0fs\t", so.sdict[statName][i])
				else
					@printf("%.5f\t", so.sdict[statName][i])
				end
			end
			print("\n")
		elseif occursin("bias", lowercase(statName)) && !occursin("log", lowercase(statName))
			@printf("%s\t\t%.0fs\n", so.name, so.sdict[statName])
		else
			@printf("%s\t\t%.5f\n", so.name, so.sdict[statName])
		end
	end
end

"""
	`plotStats`: plot relevant statistic, identified by statName, stored in given GeoStats object
"""
function plotStats(stats::Vector{RealGeoStats}, statName::AbstractString; color::AbstractString="red")
	COLORS = ["red", "blue", "green", "orange", "black"]
	LABELS = ["<4min30s", "<7min30s", "<12min", "<15min", ">15 min"]
	# check if stat is valid and if timebounds are same
	for so in stats
		if !(statName in collect(keys(so.sdict)))
			error("Statistic not found")
		end
		if occursin("bt", statName)
			if so.timeBound != stats[1].timeBound
				error("Time breakdown arrays don't match")
			end
		end
	end
	x = collect(eachindex(stats))
	xlabels = [so.name for so in stats]
	if occursin("bt", statName)
		for (i, timeBound) in enumerate(stats[1].timeBound)
			stat = [stats[j].sdict[statName][i] for j=eachindex(stats)]
			plot(x, stat, "o", color = COLORS[i % 5 + 1], label = LABELS[i])
		end
		legend(loc="best")
	else
		plot(x, [so.sdict[statName] for so in stats], "o", color = color)
	end
	xticks(x, xlabels, rotation = "vertical")
	xlabel("Method used")
	ylabel(statName)
	title("$(statName)")
	margins(0.2)
end


text:
using StatsBase

function costMatrixjl(muhat, nuhat, p, d)
    mhsize = size(muhat)[1]
    nhsize = size(nuhat)[1]

    println("mhsize: ", mhsize)
    println("nhsize: ", nhsize)
    out = zeros(mhsize, nhsize)
    for i in 1:mhsize
        for j in 1:nhsize
            out[i,j] = (mean(abs.(muhat[:, i] - nuhat[:, j]).^q))^(p / q)
        end
    end
    return out
end

max_stars_repo_name: jacobusmmsmit/wasserstein-simstudy

max_issues_repo_name: jacobusmmsmit/wasserstein-simstudy

max_forks_repo_name: jacobusmmsmit/wasserstein-simstudy


text:
using MittagLeffler
using Test

@testset "mittleff" begin
    @test isapprox(mittleff(0.5, 0.5, 0.5),  1.5403698281390346)
    @test isapprox(mittleff(1.5, 0.5, 0.5), 1.1448466286155243)
    @test isapprox(mittleff(2.3, 0.7 + 2.0 * im), 1.201890136368392 + 0.7895394560075035 * im)
    @test isapprox(mittleff(2.3, 0.7 + 0.2 * im) , (1.268233154873853 + 0.07914994421659409im))
    @test mittleff(big".3",100.0) ==
        big"8.721285946907744692995882256235296113802695745418015206361825134909144332670706e+2015816"
    # allow alpha == beta == 1  --> exp(z)
    @test isapprox(mittleff(1, -2), 0.1353352832366127)
end

@testset "branches" begin
    # convert z to float or complex float to avoid integer to negative power.
    @test isapprox(mittleff(0.9, 0.5, 22 + 22im), -2.7808021618204008e13 - 2.8561425165239754e13im)
    # Test branch that calls `Pint`
    @test isapprox(mittleff(0.1, 1.05, 0.9 + 0.5im), 0.17617901349590603 + 2.063981943021305im)
    @test isapprox(mittleff(4.1,1), 1.0358176744122032)
end

@testset "derivative" begin
    myder(f,x,h) = (f(x+h/2)-f(x-h/2))/h
    @test isapprox(myder(z -> mittleff(.4,z),.4,1e-5), mittleffderiv(.4,.4))
    @test abs(myder(z -> mittleff(big".4",z),big".4",BigFloat(1//10^16)) - mittleffderiv(big".4",big".4")) < 1e-32
end


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

struct AppsPermissionsScopesApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_appsPermissionsScopesList(_api::AppsPermissionsScopesApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.scopes.list", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Returns list of scopes this app has on a team.
Return: Dict{String, Any}
"""
function appsPermissionsScopesList(_api::AppsPermissionsScopesApi; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsScopesList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsScopesList(_api::AppsPermissionsScopesApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsScopesList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export appsPermissionsScopesList


text:
# - - parameters - - - - - - - - - - -
#=
  Symanzik improved gauge actions debuged by
  https://link.springer.com/article/10.1140/epjc/s10052-017-5392-6
  Improved thermodynamics of SU(2) gauge theory Pietro Giudice & Stefano Piemonte
=#
# Physical setting
system["L"] = (8, 8, 8, 8)
system["β"] = 1.6
system["NC"] = 2
system["Nthermalization"] = 0
system["Nsteps"] = 201
system["initial"] = "cold"
system["initialtrj"] = 1
system["update_method"] = "Heatbath"
system["useOR"] = false
system["numOR"] = 1
system["Nwing"] = 1

# Physical setting(fermions)
system["quench"] = true
system["Dirac_operator"] = nothing
system["BoundaryCondition"] = [1, 1, 1, -1]

# System Control
system["log_dir"] = "./logs"
system["logfile"] = "Heatbath_L08080808_beta1.6_quenched.txt"
system["saveU_dir"] = ""
system["saveU_format"] = nothing
system["verboselevel"] = 2
system["randomseed"] = 111
measurement["measurement_basedir"] = "./measurements"
measurement["measurement_dir"] = "Heatbath_L08080808_beta1.6_quenched"

# HMC related

# Action parameter for SLMC
actions["use_autogeneratedstaples"] = false
actions["couplingcoeff"] = Any[]
actions["couplinglist"] = Any[]

# Measurement set
measurement["measurement_methods"] = Dict[
  Dict{Any,Any}("methodname" => "Polyakov_loop",
    "measure_every" => 1
  ),
  Dict{Any,Any}("methodname" => "Plaquette",
    "measure_every" => 1
  )
]

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

#-(\beta/2) Re Tr((5/3) plaquette - (1/12) r
beta = 1.6
actions["use_autogeneratedstaples"] = true
actions["couplingcoeff"] = [beta*5/3,-beta/12]
#actions["couplinglist"] = ["plaq","rect","polyx"]

#=
function makeplaq()
    loopset = []
    for μ=1:4
        for ν=μ:4
            if ν == μ
                continue
            end
            push!(loopset,[(μ,1),(ν,1),(μ,-1),(ν,-1)])
        end
    end
    return loopset
end

function makechair()
    loopset = []
    set1 = (1,2,3)
    set2 = (1,2,4)
    set3 = (2,3,4)
    set4 = (1,3,4)

    for set in (set1,set2,set3,set4)
        mu,nu,rho = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        mu,rho,nu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        nu,rho,mu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        nu,mu,rho = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        rho,mu,nu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        rho,nu,mu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

    end
    return loopset
end
=#

plaqloops = make_plaqloops()
rectloops = make_rectloops(2)

#=
chairloops = makechair()
L = system["L"]
poly1loops = make_polyakovloops(1,L[1])
poly2loops = make_polyakovloops(2,L[2])
poly3loops = make_polyakovloops(3,L[3])
poly4loops = make_polyakovloops(4,L[4])

=#

actions["coupling_loops"] = [plaqloops,rectloops]


text:
using KernelAbstractions
using Test

@testset "Unittests" begin
    include("test.jl")
end

@testset "Localmem" begin
    include("localmem.jl")
end

@testset "Private" begin
    include("private.jl")
end

@testset "Unroll" begin
    include("unroll.jl")
end

@testset "NDIteration" begin
    include("nditeration.jl")
end

@testset "async_copy!" begin
    include("async_copy.jl")
end

@testset "Events" begin
    include("events.jl")
end

@testset "Printing" begin
    include("print_test.jl")
end

@testset "Compiler" begin
    include("compiler.jl")
end

include("examples.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
	using Random, Distributions, Plots
	using HypothesisTests
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.05"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data1 = CSV.read(datadir("machine1.csv"), DataFrame; delim=",")[:, 1]
	data2 = CSV.read(datadir("machine2.csv"), DataFrame; delim=",")[:, 1]

	xBar1, xBar2 = mean(data1), mean(data2)
	s1, s2 = std(data1), std(data2)
	n1, n2 = length(data1), length(data2)
	alpha = 0.05

	v = (s1^2/n1 + s2^2/n2)^2 / ( (s1^2/n1)^2 / (n1-1) + (s2^2/n2)^2 / (n2-1) )

	t = quantile(TDist(v),1-alpha/2)

	Text("Calculating formula: $((xBar1 - xBar2 - t*sqrt(s1^2/n1 + s2^2/n2),
				  xBar1 - xBar2 + t*sqrt(s1^2/n1 + s2^2/n2)))")
end

# ╔═╡ 60cd0974-3b1d-11eb-0c62-0fe50e45fd7b
Text("Using confint(): $(confint(UnequalVarianceTTest(data1,data2),alpha))")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.05"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═60cd0974-3b1d-11eb-0c62-0fe50e45fd7b
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using Cointegration, Test, LinearAlgebra, StatsBase, Statistics, StableRNGs, CSV, DataFrames

@testset "Basic" begin
    k = 5
    n = 200

    rng = StableRNG(123)
    X = cumsum(randn(rng, n, k), dims = 1)

    @testset "I(1) model" begin
        f = civecmI1(X, lags=2)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), f) == """
Cointegration.CivecmI1

α:
5×5 Matrix{Float64}:
 -0.0947422   0.0761861   0.0257547  -0.0974986    0.0115806
  0.0632277   0.0336236   0.0637774   0.0303369    0.0670139
 -0.0428606  -0.0861965  -0.142879   -0.015907     0.032662
  0.165899    0.146969   -0.065865   -0.00883288  -0.0134935
  0.168351   -0.118823    0.0502608  -0.0596831    0.00882035

βᵀ:
5×5 Matrix{Float64}:
 0.150191    0.0384534   0.11206     0.0451019   -0.26955
 0.014997   -0.0768655   0.265054   -0.189269     0.268929
 0.0760236   0.145551    0.067286    0.00991004  -0.123557
 0.110725   -0.092395   -0.0527418   0.0568719    0.14021
 0.0712217  -0.0709114  -0.0881567  -0.0471977   -0.205761

Π:
5×5 Matrix{Float64}:
 -0.0210996   0.00243657   0.0154309   -0.024529     0.0267911
  0.0229809   0.00157467   0.0127809   -0.00431776  -0.0254161
 -0.0180272  -0.0166653   -0.0393039    0.010519    -0.0029248
  0.0201742  -0.0127313    0.0547689   -0.0208525    0.00448221
  0.0213436   0.0278115   -0.00687702   0.0267698   -0.093727"""

        normtest = normalitytest(f)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), normtest) == """
Cointegration.NormalityTest

Univarite tests:
Test values   df   p-values
       0.68    2       0.71
       4.39    2       0.11
       0.76    2       0.68
       3.00    2       0.22
       1.20    2       0.55

Multivariate test:
Test values   df   p-values
      10.03   10       0.44"""

        rt = ranktest(rng, f, 10000)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), rt) == """
Cointegration.TraceTest

 Rank    Value  p-value
    0   36.864    0.849
    1   21.054    0.865
    2   10.397    0.836
    3    4.067    0.703
    4    1.447    0.268"""

        f3 = setrank(f, 3)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), f3) == """
Cointegration.CivecmI1

α:
5×3 Matrix{Float64}:
 -0.0947422   0.0761861   0.0257547
  0.0632277   0.0336236   0.0637774
 -0.0428606  -0.0861965  -0.142879
  0.165899    0.146969   -0.065865
  0.168351   -0.118823    0.0502608

βᵀ:
3×5 Matrix{Float64}:
 0.150191    0.0384534  0.11206    0.0451019   -0.26955
 0.014997   -0.0768655  0.265054  -0.189269     0.268929
 0.0760236   0.145551   0.067286   0.00991004  -0.123557

Π:
5×5 Matrix{Float64}:
 -0.0111289  -0.0057506    0.0113096   -0.0184374    0.0428442
  0.0148491   0.00912971   0.0202887   -0.00288017  -0.0158808
 -0.0185921  -0.0158189   -0.0372635    0.0129653    0.00602609
  0.0221132  -0.0145042    0.0531135   -0.020987     0.00294423
  0.0273238   0.0229226   -0.00924725   0.0305804   -0.0835439"""
    end

    @testset "I(2) model" begin
        f = civecmI2(X, lags=2)
        rt = ranktest(f)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), rt) == """
5×6 Matrix{Float64}:
 768.781  591.41   430.957  280.984   149.057    36.8636
   0.0    455.037  297.547  148.5      24.8888   21.0544
   0.0      0.0    151.25    23.776    13.7991   10.397
   0.0      0.0      0.0     10.2462    4.07213   4.06734
   0.0      0.0      0.0      0.0       1.45113   1.44657"""

        rng = StableRNG(123)
        rtsim = ranktest(rng, f, 100)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), rtsim[2]) == """
5×6 Matrix{Float64}:
 0.0  0.0  0.0  0.0   0.0   0.88
 0.0  0.0  0.0  0.0   0.89  0.87
 0.0  0.0  0.0  0.79  0.89  0.81
 0.0  0.0  0.0  0.91  0.98  0.75
 0.0  0.0  0.0  0.0   0.56  0.29"""

        f11 = setrank(f, 1, 1)
        @test sprint((io, t) -> show(io, MIME"text/plain"(), t), f11) == """
β':
1×5 Matrix{Float64}:
 -0.0214383  -0.00947104  -0.0150298  -0.00571825  0.0427769

τ⊥δ:
5×1 Matrix{Float64}:
  0.0571834
 -0.0232661
 -0.611789
  0.761213
 -0.0896907

τ':
2×5 Matrix{Float64}:
  0.13321    0.0446713  0.103897  0.0415547  -0.282675
 -0.341134  -0.856249   0.283727  0.208753   -0.159012"""
    end
end

@testset "Illustrative examples from Johansen's book" begin
    danish_df  = CSV.read(joinpath(dirname(pathof(Cointegration)), "..", "data", "danish.csv"), DataFrame)
    danish_mat = [danish_df.log_real_money danish_df.log_real_income danish_df.bond_rate danish_df.deposit_rate]

    # Create constant term and seasonals.
    cnst = [ones(55) repeat(Matrix{Float64}(I, 4, 3), 14)[1:55, :]]
    # Johansen makes the seasonals orthogonal to the constant term
    proj(v, u) = (v'u)/(u'u)*u
    for k in 2:4
        cnst[:, k] .-= sum(proj(cnst[:, k], cnst[:, j]) for j in 1:k-1)
    end

    @testset "Chapter 2" begin
        f = civecmI1(danish_mat, unrestricted=cnst, lags=2)

        @testset "Table 2.1" begin
            @test f.Π[:,1:4] ≈ [-0.181  0.110 -1.042  0.638
                                 0.186 -0.309  0.658 -0.648
                                 0.014 -0.018  0.082 -0.167
                                -0.004  0.020  0.143 -0.314] atol=2e-3
        end

        @testset "Table 2.2" begin
            @test f.Γ[:,1:4] ≈ [0.195 -0.096 -0.138 -0.462
                                0.504 -0.045 -0.377  0.060
                                0.051  0.136  0.301  0.253
                                0.069 -0.022  0.227  0.265] atol=2e-3
        end

        @testset "Table 2.3" begin
            @testset "The constant" begin
                @test f.Γ[:, 5] ≈ [1.583, -0.390, -0.064, -0.071] atol=2e-3
            end

            @testset "The seasonal" begin
                @test_broken f.Γ[:, 6:8] ≈ [-0.023  0.016 -0.039
                                            -0.019 -0.007 -0.032
                                            -0.003 -0.007 -0.007
                                            -0.002  0.001 -0.003]
            end
        end

        @testset "Table 2.4" begin
            @test cor(residuals(f)) ≈ [ 1.00  0.53 -0.45 -0.31
                                        0.53  1.00 -0.08 -0.24
                                       -0.45 -0.08  1.00  0.25
                                       -0.31 -0.24  0.25  1.00] atol=1e-2
        end

        @testset "Table 2.5" begin
            # FIXME! ARCH and JB tests not implemented yet
            @test mapslices(skewness, residuals(f), dims=1) ≈ [0.552 0.524 -0.297 0.415] atol=1e-3
            @test mapslices(kurtosis, residuals(f), dims=1) ≈ [-0.075 -0.087 0.576 0.562] atol=1e-3
        end

        @testset "Table 2.6" begin
            @test sort(eigvals(f), by=real, rev=true) ≈ [ 0.9725,
                                                         0.7552 - 0.1571im,
                                                         0.7552 + 0.1571im,
                                                         0.6051,
                                                         0.5955 - 0.3143im,
                                                         0.5955 + 0.3143im,
                                                        -0.1425 - 0.2312im,
                                                        -0.1425 + 0.2312im] atol=1e-3
        end
    end

    @testset "Chapter 7" begin
        # Create centered seasonal dummies
        cseason = repeat(Matrix{Float64}(I, 4, 3), 14)[1:55,:] .- 1/4
        f = civecmI1(danish_mat, exogenous=ones(55, 1), unrestricted=cseason)

        @testset "Table 7.1" begin
            rt = ranktest(f)
            # Notice that there is a typo in the book where it's 8.89 instead of 8.69 (correct value in the 1990 paper)
            @test rt.values ≈ [49.14, 19.06, 8.69, 2.35] atol=1e-2
            @test all(rt.pvalues .> 0.05)
        end

        @testset "Table 7.2" begin
            # Signs flipped relative to book
            @test f.β ≈ [ 21.97  -14.66   7.95  -1.02
                         -22.70   20.05 -25.64   1.93
                         114.42   -3.56   4.28 -25.00
                         -92.64 -100.26 -44.88  14.65
                        -133.16   62.59  62.75   2.32] atol=2e-2
        end

        @testset "Single coitegration relation" begin
            f1 = setrank(f, 1)
            @test f1.β/f1.β[1,1] ≈ [1.00, -1.03, 5.21, -4.22, -6.06] atol=2e-2
            @test f1.α*f1.β[1,1] ≈ [-0.213, 0.115, 0.023, 0.029]     atol=2e-2

            @testset "restrictions on α and β" begin
                f1r1 = restrict(f1, Hβ=[1 0 0 0; -1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1])
                @test f1r1.β/f1r1.β[1,1]*-21.55      ≈ [-21.55, 21.55, -114.22, 92.45, 134.99] atol=2e-2
                # It appears that there are sign typos in the book and the magnitude is oof by 1000
                @test f1r1.α*f1r1.β[1,1]/-21.55*1000 ≈ -[-9.83, 4.99, 1.05, 1.38]              atol=2e-2
                @test lrtest(f1r1, f1, 1).value ≈ 0.043 atol=1e-2

                f1r2 = restrict(f1, Hβ=[1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1])
                @test lrtest(f1r2, f1r1, 1).value ≈ 0.89 atol=1e-2

                f1r3 = restrict(f1r2, Hα=[1 0; 0 1; 0 0; 0 0])
                @test lrtest(f1r3, f1r2, 2).value ≈ 5.81 atol=1e-2
            end
        end
    end
end


text:
using MCTS
using POMDPs
using POMDPModels
using Test
using NBInclude
using D3Trees
using Random
using POMDPPolicies

n_iter = 50
depth = 15
ec = 1.0

println("Testing vanilla MCTS solver.")

solver = MCTSSolver(n_iterations=n_iter, depth=depth, exploration_constant=ec, enable_tree_vis=true)
mdp = LegacyGridWorld()

struct A
    a::Vector{Int}
end

println("============== @requirements_info with only solver:")
@requirements_info solver
println("============== @requirements_info with solver and mdp:")
@requirements_info solver mdp
println("============== @requirements_info with solver, mdp, and state:")
@requirements_info solver mdp GridWorldState(1,1)
println("============== isequal and hash warnings:")
@requirements_info solver mdp A([1,2,3])

policy = solve(solver, mdp)

state = GridWorldState(1,1)

@testset "basic" begin
    a = @inferred action(policy, state)

    tree = policy.tree
    @test get_state_node(tree, state).id == 1
    @test get_state_node(tree, state, policy).id == 1

    clear_tree!(policy)
    @test policy.tree == nothing
end

@testset "value" begin
    include("value_test.jl")
end

@testset "options" begin
    include("options.jl")
end

@testset "dpw" begin
    include("dpw_test.jl")
end

@testset "visualization" begin
    include("visualization.jl")
end
@nbinclude("../notebooks/Test_Visualization.ipynb")

@testset "other" begin
    include("other.jl")
end

# # test the BeliefMCTSSolver docstring
# let
#     using ParticleFilters
#     using POMDPModels
#     using MCTS
#     using POMDPToolbox

#     pomdp = BabyPOMDP()
#     updater = SIRParticleFilter(pomdp, 1000)

#     solver = BeliefMCTSSolver(DPWSolver(), updater)
#     planner = solve(solver, pomdp)

#     @inferred action(planner, initialize_belief(updater, initialstate_distribution(pomdp)))

#     simulate(HistoryRecorder(max_steps=10), pomdp, planner, updater)
# end

@testset "timing" begin
    solver = DPWSolver(n_iterations=typemax(Int),
                       depth=depth,
                       max_time=1.0,
                       exploration_constant=ec)
    mdp = LegacyGridWorld()

    policy = solve(solver, mdp)
    state = GridWorldState(1,1)
    a = action(policy, state)
    t = @elapsed begin
        action(policy, state)
    end
    @test abs(t-1.0) < 0.5
end

@testset "timing" begin
    solver = DPWSolver(n_iterations=typemax(Int),
                       depth=depth,
                       max_time=1.0,
                       exploration_constant=ec)
    mdp = LegacyGridWorld()

    policy = solve(solver, mdp)
    state = GridWorldState(1,1,true)
    @test_throws ErrorException action(policy, state)
end

@testset "c=0" begin
    mdp = LegacyGridWorld()

    solver = DPWSolver(n_iterations=typemax(Int),
                       depth=depth,
                       max_time=1.0,
                       exploration_constant=0.0)

    policy = solve(solver, mdp)
    state = GridWorldState(1,1)
    action(policy, state)

    solver = MCTSSolver(exploration_constant=0.0)

    policy = solve(solver, mdp)
    state = GridWorldState(1,1)
    action(policy, state)
end

@nbinclude("../notebooks/Domain_Knowledge_Example.ipynb")


text:
module Sonums

# precompilation crashes, presumably due to large matrices being preallocated.
__precompile__(false)

export AbstractSonum, Sonum8, Sonum16, trainSonum8, trainSonum16,
    fillSonum8Tables, fillSonum16Tables, fillTable,
    returnSonum8List, returnSonum16List,
    setSonum8, setSonum16, setOverflow8, setOverflow16, setUnderflow8, setUnderflow16,
    maxentropy

import Base: Float64, Float32, Float16, Int32, Int64, Int,
    UInt8, UInt16,
    (+), (-), (*), (/), sqrt, bitstring, round,
    one, zero, promote_rule, floatmin, floatmax,
    signbit, sign, isfinite, nextfloat, prevfloat, abs, inv

using Printf, Random, LinearAlgebra

include("typedef.jl")
include("helpers.jl")
include("feedback.jl")
include("conversionsInt.jl")
include("sign.jl")
include("print.jl")
include("constants.jl")
include("preallocate.jl")

# set-up zeroed number lists and bounds as constants
const sonum8 = createSonumList(Float64,8)
const sonum16 = createSonumList(Float64,16)
const bounds8 = createSonumList(Float64,8)
const bounds16 = createSonumList(Float64,16)

# set-up empty (undef) look up tables and their symmetric views
const TableAdd8S, TableAdd8 = createSonumTable(Sonum8,8)
const TableSub8S, TableSub8 = createSonumTable(Sonum8,8)
const TableMul8S, TableMul8 = createSonumTable(Sonum8,8)
const TableDiv8S, TableDiv8 = createSonumTable(Sonum8,8)

const TableAdd16S, TableAdd16 = createSonumTable(Sonum16,16)
const TableSub16S, TableSub16 = createSonumTable(Sonum16,16)
const TableMul16S, TableMul16 = createSonumTable(Sonum16,16)
const TableDiv16S, TableDiv16 = createSonumTable(Sonum16,16)

# set-up zeroed look up lists
const ListSqrt8 = createSonumList(Sonum8,8)
const ListInv8 = createSonumList(Sonum8,8)

const ListSqrt16 = createSonumList(Sonum16,16)
const ListInv16 = createSonumList(Sonum16,16)

# include functions that rely on sonum8, sonum16, bounds8, bounds16
include("returnSonumList.jl")
include("conversionsFloat.jl")

# define arithmetics based on lookup tables and lists
include("arithmetics8.jl")
include("arithmetics16.jl")

# functions for training and filling the lookup tables/lists
include("clustering.jl")
include("training.jl")
include("bounds.jl")
include("lookup_tables.jl")

end


text:
@testset "ruleset_loading.jl" begin
    @testset "on_new_rule" begin
        frule_history = []
        rrule_history = []
        on_new_rule(frule) do sig
            op = sig.parameters[1]
            push!(frule_history, op)
        end
        on_new_rule(rrule) do sig
            op = sig.parameters[1]
            push!(rrule_history, op)
        end

        @testset "new rules hit the hooks" begin
            # Now define some rules
            @scalar_rule x + y (1, 1)
            @scalar_rule x - y (1, -1)
            refresh_rules()

            @test Set(frule_history[end-1:end]) == Set((typeof(+), typeof(-)))
            @test Set(rrule_history[end-1:end]) == Set((typeof(+), typeof(-)))
        end

        @testset "# Make sure nothing happens anymore once we clear the hooks" begin
            ChainRulesOverloadGeneration.clear_new_rule_hooks!(frule)
            ChainRulesOverloadGeneration.clear_new_rule_hooks!(rrule)

            old_frule_history = copy(frule_history)
            old_rrule_history = copy(rrule_history)

            @scalar_rule sin(x) cos(x)
            refresh_rules()

            @test old_rrule_history == rrule_history
            @test old_frule_history == frule_history
        end
    end

    @testset "_primal_sig" begin
        _primal_sig = ChainRulesOverloadGeneration._primal_sig
        @testset "frule" begin
            @test isequal(  # DataType without shared type but with constraint
                _primal_sig(frule, Tuple{typeof(frule), Any, typeof(*), Int, Vector{Int}}),
                Tuple{typeof(*), Int, Vector{Int}}
            )
            @test isequal(  # UnionAall without shared type but with constraint
                _primal_sig(frule, Tuple{typeof(frule), Any, typeof(*),  T, Int} where T<:Real),
                Tuple{typeof(*), T, Int} where T<:Real
            )
            @test isequal(  # UnionAall with share type
                _primal_sig(frule, Tuple{typeof(frule), Any, typeof(*), T, Vector{T}} where T),
                Tuple{typeof(*), T, Vector{T}} where T
            )
        end

        @testset "rrule" begin
            @test isequal(  # DataType without shared type but with constraint
                _primal_sig(rrule, Tuple{typeof(rrule), typeof(*), Int, Vector{Int}}),
                Tuple{typeof(*), Int, Vector{Int}}
            )
            @test isequal(  # UnionAall without shared type but with constraint
                _primal_sig(rrule, Tuple{typeof(rrule), typeof(*),  T, Int} where T<:Real),
                Tuple{typeof(*), T, Int} where T<:Real
            )
            @test isequal(  # UnionAall with share type
                _primal_sig(rrule, Tuple{typeof(rrule), typeof(*), T, Vector{T}} where T),
                Tuple{typeof(*), T, Vector{T}} where T
            )
        end
    end

    @testset "_rule_list" begin
        _rule_list = ChainRulesOverloadGeneration._rule_list
        @testset "should not have frules that need RuleConfig" begin
            old_frule_list = collect(_rule_list(frule))
            function ChainRulesCore.frule(
                ::RuleConfig{>:Union{HasForwardsMode,HasReverseMode}}, dargs, sum, f, xs
            )
                return 1.0, 1.0 # this will not be call so return doesn't matter
            end
            # New rule should not have appeared
            @test collect(_rule_list(frule)) == old_frule_list
        end

        @testset "should not have rrules that need RuleConfig" begin
            @testset "normal type sigs" begin
                old_rrule_list = collect(_rule_list(rrule))
                function ChainRulesCore.rrule(
                    ::RuleConfig{>:Union{HasForwardsMode,HasReverseMode}}, sum, f, xs
                )
                    return 1.0, x->(x,x,x) # this will not be call so return doesn't matter
                end
                # New rule should not have appeared
                @test collect(_rule_list(rrule)) == old_rrule_list
            end
            @testset "UnionAll type sigs" begin
                old_rrule_list = collect(_rule_list(rrule))
                function ChainRulesCore.rrule(
                    ::RuleConfig{>:Union{HasForwardsMode,HasReverseMode}}, sum, f::F, xs
                ) where F <: Function
                    return 1.0, x->(x,x,x) # this will not be call so return doesn't matter
                end
                # New rule should not have appeared
                @test collect(_rule_list(rrule)) == old_rrule_list
                # Above would error if we were not handling UnionAll's right
            end
        end

        @testset "opting out" begin
            oa_id(x, y) = x
            @scalar_rule(oa_id(x::Number), 1)
            @opt_out ChainRulesCore.rrule(::typeof(oa_id), x::Float32)
            @opt_out ChainRulesCore.frule(::Any, ::typeof(oa_id), x::Float32)

            # In theses tests we `@assert` the behavour that `methods` has
            # and then `@test` that `_rule_list` differs from that, in the way we want

            @test !isempty([m for m in _rule_list(rrule) if m.sig <: Tuple{Any,typeof(oa_id),Number}])
            # Opted out
            @assert !isempty([m for m in methods(rrule) if m.sig <: Tuple{Any,typeof(oa_id),Float32}])
            @test isempty([m for m in _rule_list(rrule) if m.sig <: Tuple{Any,typeof(oa_id),Float32}])
            # fallback
            @test !isempty([m for m in methods(rrule) if m.sig == Tuple{typeof(rrule),Any,Vararg{Any}}])
            @test isempty([m for m in _rule_list(rrule) if m.sig == Tuple{typeof(rrule),Any,Vararg{Any}}])

            @test !isempty([m for m in _rule_list(frule) if m.sig <: Tuple{Any,Any,typeof(oa_id),Number}])
            # Opted out
            @assert !isempty([m for m in methods(frule) if m.sig <: Tuple{Any,Any,typeof(oa_id),Float32}])
            @test isempty([m for m in _rule_list(frule) if m.sig <: Tuple{Any,Any,typeof(oa_id),Float32}])
            # fallback
            @assert !isempty([m for m in methods(frule) if m.sig == Tuple{typeof(frule),Any,Any,Vararg{Any}}])
            @test isempty([m for m in _rule_list(frule) if m.sig == Tuple{typeof(frule),Any,Any,Vararg{Any}}])
        end
    end
end


text:
using ModelingToolkit, Unitful
using Test

t = Variable{u"s"}(:t)()
x = Variable{u"kg"}(:x)(t)
y = Variable{u"kg"}(:y)(t)
D = Differential(t)

eq1 = x ~ y*t
eq2 = x*10u"s" ~ y*t

@test ModelingToolkit.instantiate(t) == 1u"s"
@test ModelingToolkit.instantiate(x) == 1u"kg"
@test ModelingToolkit.instantiate(y) == 1u"kg"

@test !ModelingToolkit.validate(eq1)
@test ModelingToolkit.validate(eq2)

eqs = [
        D(x) ~ y/t
        D(y) ~ (x*y)/(t*10u"kg")
]

sys = ODESystem(eqs,t,[x,y],[])
@test ModelingToolkit.validate(sys)


text:
const opt_C = "**C** : **color** : **cmap** : -- Str --		Flags = [cpt |master[+izinc] |color1,color2[,*color3*,…]]

    Name of the CPT (for grd_z only). Alternatively, supply the name of a GMT color
    master dynamic CPT.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#c)"

const opt_J = "**J** : **proj** : **projection** : -- Str --

    Select map projection. Defaults to 12x8 cm with linear (non-projected) maps.
    [`-J`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#j-full)"

const opt_Jz = "**Jz** : **z_axis** : -- Str --

    Set z-axis scaling.
    [`-Jz`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#jz-full)"

const opt_R = "**R** : **region** : **limits** : -- Str or list or GMTgrid|image --		Flags = xmin/xmax/ymin/ymax[+r][+uunit]

    Specify the region of interest. Set to data minimum BoundinBox if not provided.
    [`-R`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#r-full)"

const opt_B = "**B** : **frame** : **axes** : -- Str --

    Set map boundary frame and axes attributes.
    [`-B`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#b-full)"

const opt_P = "**P** : **portrait** : --- Bool or [] --

    Tell GMT to **NOT** draw in portriat mode (that is, make a Landscape plot)"

const opt_U = "**U** : **stamp** : -- Str or Bool or [] --	Flags = [[just]/dx/dy/][c|label]

    Draw GMT time stamp logo on plot.
    [`-U`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#u-full)"

const opt_V = "**V** : **verbose** : -- Bool or Str --		Flags = [level]

    Select verbosity level, which will send progress reports to stderr.
    [`-V`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#v-full)"

const opt_X = "**X** : **x_offset** : -- Str --Flags = [a|c|f|r][x-shift[u]]"

const opt_Y = "**Y** : **y_offset** : -- Str --Flags = [a|c|f|r][y-shift[u]]

	Shift plot origin relative to the current origin by (x-shift,y-shift) and optionally
	append the length unit (c, i, or p).
    [`-Y`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#xy-full)"

const opt_a = "**a** : **aspatial** : -- Str --			Flags = [col=]name[…]

      Control how aspatial data are handled in GMT during input and output.
      [`-a`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#aspatial-full)"

const opt_b = "**b** : **binary** : -- Str --

      [`-b`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#b-full)"

const opt_bi = "**bi** : **binary_in** : -- Str --			Flags = [ncols][type][w][+L|+B]

      Select native binary format for primary input (secondary inputs are always ASCII).
      [`-bi`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#bi-full)"

const opt_bo = "**bo** : **binary_out** : -- Str --			Flags = [ncols][type][w][+L|+B]

      Select native binary output.
      [`-bo`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#bo-full)"

const opt_c = "**c** : **nodata** : -- Str --				Flags = row,col

      Used to advance to the selected subplot panel. Only allowed when in subplot mode.
      [`-c`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#c-full)"

const opt_d = "**d** : **nodata** : -- Str or Number --		Flags = [i|o]nodata

      Control how user-coded missing data values are translated to official NaN values in GMT.
      [`-d`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#d-full)"

const opt_di = "**di** : **nodata_in** : -- Str or Number --      Flags = nodata

      Examine all input columns and if any item equals nodata we interpret this value as a
      missing data item and substitute the value NaN.
      [`-di`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#di-full)"

const opt_do = "**do** : **nodata_out** : -- Str or Number --     Flags = nodata

      Examine all output columns and if any item equals NAN substitute it with
      the chosen missing data value.
      [`-do`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#do-full)"

const opt_e = "**e** : **pattern** : -- Str --        Flags = [~]”pattern” | -e[~]/regexp/[i]

      Only accept ASCII data records that contains the specified pattern.
      [`-e`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#e-full)"

const opt_f = "**f** : **colinfo** : -- Str --        Flags = [i|o]colinfo

      Specify the data types of input and/or output columns (time or geographical data).
      [`-f`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#f-full)"

const opt_g = "**g** : **gaps** : -- Str --           Flags = [a]x|y|d|X|Y|D|[col]z[+|-]gap[u]

      Examine the spacing between consecutive data points in order to impose breaks in the line.
      [`-g`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#g-full)"

const opt_h = "**h** : **headers** : -- Str --        Flags = [i|o][n][+c][+d][+rremark][+ttitle]

      Primary input file(s) has header record(s).
      [`-h`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#h-full)"

const opt_i = "**i** : **input_col** : -- Str --      Flags = cols[+l][+sscale][+ooffset][,…]

      Select specific data columns for primary input, in arbitrary order.
      [`-i`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#icols-full)"

const opt_n = "**n** : **interp** : -- Str --         Flags = [b|c|l|n][+a][+bBC][+c][+tthreshold]

      Select grid interpolation mode by adding b for B-spline smoothing, c for bicubic interpolation,
      l for bilinear interpolation, or n for nearest-neighbor value.
      [`-n`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#n-full)"

const opt_o = "**o** : **output_col** : -- Str --     Flags = cols[,…]

      Select specific data columns for primary output, in arbitrary order.
      [`-o`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#ocols-full)"

const opt_p = "**p** : **view** : **perspective** : -- Str or List --   Flags = [x|y|z]azim[/elev[/zlevel]][+wlon0/lat0[/z0]][+vx0/y0]

      Selects perspective view and sets the azimuth and elevation of the viewpoint [180/90].
      [`-p`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#perspective-full)"

const opt_r = "**r** : **reg** : **registration** : -- Bool or [] --

      Force pixel node registration [Default is gridline registration].
      [`-r`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#r-full)"

const opt_s = "**s** : **skip_col** : -- Str --       Flags = [cols][a|r]

      Suppress output for records whose z-value equals NaN.
      [`-s`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#s-full)"

const opt_t = "**t** : **alpha** : **transparency** : -- Str --   Flags = [transp]

      Set PDF transparency level for an overlay, in (0-100] percent range. [Default is 0, i.e., opaque].
      [`-t`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#t-full)"

const opt_x = "**x** : **n_threads** : -- Str or Number --  Flags = [[-]n]

      Limit the number of cores to be used in any OpenMP-enabled multi-threaded algorithms.
      [`-x`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#x-full)"

const opt_swap_xy = "**swap_xy** : Str or Bool or [] --     Flags = [i|o]

      Swap 1st and 2nd column on input and/or output.
      [`-:`](http://gmt.soest.hawaii.edu/doc/latest/gmt.html#colon-full)"

const opt_write = "**write** : **|>** : Str --     Flags = ``fname``

      Save result to ASCII file instead of returning to a Julia variable. Give file name as argument.
      Use the bo option to save as a binary file."

const opt_append = "**append** : Str --     Flags = ``fname``

      Append result to an existing file named ``fname`` instead of returning to a Julia variable.
      Use the bo option to save as a binary file."


text:
abstract type FlatProj <: FieldMetadata end

struct ProjLambert{T, V<:AbstractVector{T}, M<:AbstractMatrix{T}} <: FlatProj
    # these must be the same to broadcast together
    Ny        :: Int
    Nx        :: Int
    θpix      :: Float64
    center    :: Tuple{Float64,Float64}
    # these can be different and still broadcast (including different types)
    storage
    Δx        :: T
    Ωpix      :: T
    nyquist   :: T
    Δℓx       :: T
    Δℓy       :: T
    ℓy        :: V
    ℓx        :: V
    ℓmag      :: M
    sin2ϕ     :: M
    cos2ϕ     :: M
end

real_type(T) = promote_type(real(T), Float32)
@init @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" real_type(::Type{<:Unitful.Quantity{T}}) where {T} = real_type(T)

ProjLambert(;Ny, Nx, θpix=1, center=(0,0), T=Float32, storage=Array) =
    ProjLambert(Ny, Nx, Float64(θpix), Float64.(center), real_type(T), storage)

@memoize function ProjLambert(Ny, Nx, θpix, center, ::Type{T}, storage) where {T}

    Δx           = T(deg2rad(θpix/60))
    Δℓx          = T(2π/(Nx*Δx))
    Δℓy          = T(2π/(Ny*Δx))
    nyquist      = T(2π/(2Δx))
    Ωpix         = T(Δx^2)
    ℓy           = adapt(storage, (ifftshift(-Ny÷2:(Ny-1)÷2) .* Δℓy)[1:Ny÷2+1])
    ℓx           = adapt(storage, (ifftshift(-Nx÷2:(Nx-1)÷2) .* Δℓx))
    ℓmag         = @. sqrt(ℓx'^2 + ℓy^2)
    ϕ            = @. angle(ℓx' + im*ℓy)
    sin2ϕ, cos2ϕ = @. sin(2ϕ), cos(2ϕ)
    if iseven(Ny)
        sin2ϕ[end, end:-1:(Nx÷2+2)] .= sin2ϕ[end, 2:Nx÷2]
    end

    ProjLambert(Ny,Nx,θpix,center,storage,Δx,Ωpix,nyquist,Δℓx,Δℓy,ℓy,ℓx,ℓmag,sin2ϕ,cos2ϕ)

end

typealias_def(::Type{<:ProjLambert{T}}) where {T} = "ProjLambert{$T}"

### promotion

# used in broadcasting to decide the resulting metadata when
# broadcasting over two fields
function promote_metadata_strict(metadata₁::ProjLambert{T₁}, metadata₂::ProjLambert{T₂} ) where {T₁,T₂}

    if (
        metadata₁.θpix === metadata₂.θpix &&
        metadata₁.Ny   === metadata₂.Ny   &&
        metadata₁.Nx   === metadata₂.Nx
    )

        # always returning the "wider" metadata even if T₁==T₂ helps
        # inference and is optimized away anyway
        promote_type(T₁,T₂) == T₁ ? metadata₁ : metadata₂

    else
        error("""Can't broadcast two fields with the following differing metadata:
        1: $(select(fields(metadata₁),(:θpix,:Ny,:Nx)))
        2: $(select(fields(metadata₂),(:θpix,:Ny,:Nx)))
        """)
    end

end

# used in non-broadcasted algebra to decide the resulting metadata
# when performing some operation across two fields. this is free to do
# more generic promotion than promote_metadata_strict (although this
# is currently not used, but in the future could include promoting
# resolution, etc...). the result should be a common metadata which we
# can convert both fields to then do a succesful broadcast
promote_metadata_generic(metadata₁::ProjLambert, metadata₂::ProjLambert) =
    promote_metadata_strict(metadata₁, metadata₂)

### preprocessing
# defines how ImplicitFields and BatchedReals behave when broadcasted
# with ProjLambert fields. these can return arrays, but can also
# return `Broadcasted` objects which are spliced into the final
# broadcast, thus avoiding allocating any temporary arrays.

function preprocess((_,proj)::Tuple{<:Any,<:ProjLambert{T,V}}, r::Real) where {T,V}
    r isa BatchedReal ? adapt(V, reshape(r.vals, 1, 1, 1, :)) : r
end
# need custom adjoint here bc Δ can come back batched from the
# backward pass even though r was not batched on the forward pass
@adjoint function preprocess(m::Tuple{<:Any,<:ProjLambert{T,V}}, r::Real) where {T,V}
    preprocess(m, r), Δ -> (nothing, Δ isa AbstractArray ? batch(real.(Δ[:])) : Δ)
end

function preprocess((_,proj)::Tuple{BaseFieldStyle{S,B},<:ProjLambert}, ∇d::∇diag) where {S,B}

    (B <: Union{Fourier,QUFourier,IQUFourier}) ||
        error("Can't broadcast ∇[$(∇d.coord)] as a $(typealias(B)), its not diagonal in this basis.")

    # turn both into 2D matrices so this function is type-stable
    # (reshape doesnt actually make a copy here, so this doesn't
    # impact performance)
    if ∇d.coord == 1
        broadcasted(*, ∇d.prefactor * im, reshape(proj.ℓx, 1, :))
    else
        broadcasted(*, ∇d.prefactor * im, reshape(proj.ℓy, :, 1))
    end
end

function preprocess((_,proj)::Tuple{BaseFieldStyle{S,B},<:ProjLambert}, ::∇²diag) where {S,B}

    (B <: Union{Fourier,<:Basis2Prod{<:Any,Fourier},<:Basis3Prod{<:Any,<:Any,Fourier}}) ||
        error("Can't broadcast a BandPass as a $(typealias(B)), its not diagonal in this basis.")

    broadcasted(+, broadcasted(^, proj.ℓx', 2), broadcasted(^, proj.ℓy, 2))
end

function preprocess((_,proj)::Tuple{<:Any,<:ProjLambert}, bp::BandPass)
    Cℓ_to_2D(bp.Wℓ, proj)
end

function Cℓ_to_2D(Cℓ, proj::ProjLambert{T}) where {T}
    Complex{T}.(nan2zero.(Cℓ.(proj.ℓmag)))
end

### adapting

# dont adapt the fields in proj, instead re-call into the memoized
# ProjLambert so we always get back the singleton ProjLambert object
# for the given set of parameters (helps reduce memory usage and
# speed-up subsequent broadcasts which would otherwise not hit the
# "===" branch of the "promote_*" methods)
function adapt_structure(storage, proj::ProjLambert{T}) where {T}
    @unpack Ny, Nx, θpix = proj
    T′ = eltype(storage)
    ProjLambert(;Ny, Nx, θpix, T=(T′==Any ? T : real(T′)), storage)
end
adapt_structure(::Nothing, proj::ProjLambert{T}) where {T} = proj

@doc doc"""
    pixwin(θpix, ℓ)

Returns the pixel window function for square flat-sky pixels of width `θpix` (in
arcmin) evaluated at some `ℓ`s. This is the scaling of k-modes, the scaling of
the power spectrum will be pixwin^2.
"""
pixwin(θpix, ℓ) = @. sinc(ℓ*deg2rad(θpix/60)/2π)

# ### serialization
# makes it so the arrays in ProjLambert objects aren't actually
# serialized, instead just (Ny, Nx, θpix, center, T) are stored, and
# deserializing just recreates the ProjLambert object, possibly from
# the memoized cache if it already exists.

function _serialization_key(proj::ProjLambert{T}) where {T}
    @unpack Ny, Nx, θpix, center, storage = proj
    (;Ny, Nx, θpix, center, T, storage)
end

# Julia serialization
function Serialization.serialize(s::AbstractSerializer, proj::ProjLambert)
    @unpack Ny, Nx, θpix, center = proj
    Serialization.writetag(s.io, Serialization.OBJECT_TAG)
    Serialization.serialize(s, ProjLambert)
    Serialization.serialize(s, _serialization_key(proj))
end
function Serialization.deserialize(s::AbstractSerializer, ::Type{ProjLambert})
    ProjLambert(; Serialization.deserialize(s)...)
end

# JLD2 serialization
# (always deserialize as Array)
function JLD2.writeas(::Type{<:ProjLambert})
    Tuple{Val{ProjLambert},NamedTuple{(:Ny,:Nx,:θpix,:center,:T),Tuple{Int,Int,Float64,Tuple{Float64,Float64},DataType}}}
end
function JLD2.wconvert(::Type{<:Tuple{Val{ProjLambert},NamedTuple}}, proj::ProjLambert)
    (Val(ProjLambert), delete(_serialization_key(proj), :storage))
end
function JLD2.rconvert(::Type{<:ProjLambert}, (_,s)::Tuple{Val{ProjLambert},NamedTuple})
    ProjLambert(; storage=Array, s...)
end


text:
"""
    haltonsequence(b, n)

Generate the first `n` numbers in Halton's sequence with base `b`.
"""
haltonsequence(b, n) = haltonnumber.(b, 1:n)

"""
    haltonnumber(base, index)

Generate the `n`-th Halton number in the sequence with base `b`.

# Note
    Implementation is based on the psudo code in:
        http://en.wikipedia.org/wiki/Halton_sequence
"""
function haltonnumber(base::Integer, index::Integer)
    res = 0
    f = 1 / base
    i = index

    while (i > 0)
        res += f * (i % base)
        i = floor(Integer, i / base)
        f = f / base
    end

    return res
end


text:
module RulesAdditionMarginalsTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_marginalrules

@testset "marginalrules:typeof(+)" begin

    @testset ":in1_in2 (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(3.0, 4.0), m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(5/4, 3/4), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(3, 4), m_in1 = NormalMeanPrecision(1, 4), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(12.0, 8.0), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(3.0, 4.0), m_in1 = NormalWeightedMeanPrecision(1.0, 4.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(0.0, 8.0), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(2.0, 4.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(9/2, 9/2), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanVariance(2.0, 4.0), m_in1 = NormalMeanPrecision(1.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(9/4, 9/4), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(3.0, 4.0), m_in1 = NormalWeightedMeanPrecision(1.0, 2.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(5.0, 6.0), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 2.0), m_in1 = NormalMeanPrecision(1.0, 2.0), m_in2 = PointMass(-1.0)),
                output = (in1 = NormalWeightedMeanPrecision(6.0, 4.0), in2 = PointMass(-1.0))
            ),
            (
                input = (m_out = NormalMeanVariance(3.0, 3.0), m_in1 = NormalWeightedMeanPrecision(2.0, 1.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(7/3, 4/3), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 4.0), m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(-1.0, 9/2), in2 = PointMass(1.0))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(3.0, 4.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(2.0, 2.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(5/4, 3/4))
            ),
            (
                input = (m_out = NormalMeanPrecision(3, 4), m_in1 = PointMass(1.0), m_in2 = NormalMeanPrecision(1, 4)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(12.0, 8.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(3.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalWeightedMeanPrecision(1.0, 4.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(0.0, 8.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(1.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(9/2, 9/2))
            ),
            (
                input = (m_out = NormalMeanVariance(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(9/4, 9/4))
            ),
            (
                input = (m_out = NormalMeanPrecision(3.0, 4.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(5.0, 6.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 2.0), m_in1 = PointMass(-1.0), m_in2 = NormalMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(-1.0), in2 = NormalWeightedMeanPrecision(6.0, 4.0))
            ),
            (
                input = (m_out = NormalMeanVariance(3.0, 3.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(2.0, 1.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(7/3, 4/3))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(-1.0, 9/2))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanCovariance, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanCovariance([2.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([1/4, 11/8], [1.0 -1/4; -1/4 7/8]), in2 = PointMass([1.0, 1.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([4.0, 7/2], [5/2 1.0; 1.0 5/2]), in2 = PointMass([1.0, 1.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in2 = PointMass([2.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([7/6, 13/6], [7/3 -1/6; -1/6 7/3]), in2 = PointMass([2.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, -1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]), in2 = PointMass([1.0, -1.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 3.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([2.0, 2.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([-1.0, -1.0], [5.0 2.0; 2.0 5.0]), in2 = PointMass([2.0, 2.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in2 = PointMass([-2.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([23.0, 16.0], [7.0  3.0; 3.0 5.0]), in2 = PointMass([-2.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalWeightedMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalWeightedMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, -1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]), in2 = PointMass([1.0, -1.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([2.0, 2.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([-4.0, -2.0], [5.0 2.0; 2.0 4.0]), in2 = PointMass([2.0, 2.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [4.0 1.0; 1.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 0.0; 0.0 2.0]), m_in2 = PointMass([-1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([5.0, 0.0], [7.0  1.0; 1.0 6.0]), in2 = PointMass([-1.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanCovariance, m_in1::PointMass, m_in2::PointMassMultivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanCovariance([2.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = PointMass([1.0, 1.0]), m_in2 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([1/4, 11/8], [1.0 -1/4; -1/4 7/8]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = PointMass([1.0, 1.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([4.0, 7/2], [5/2 1.0; 1.0 5/2]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = PointMass([2.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [2.0 0.0; 0.0 2.0])),
                output = (in1 = PointMass([2.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([7/6, 13/6], [7/3 -1/6; -1/6 7/3]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanPrecision, m_in1::, m_in2::MultivariateNormalDistributionsFamilyPointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = PointMass([1.0, -1.0]), m_in2 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, -1.0]), in2 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 3.0]), m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([2.0, 2.0]), in2 = MvNormalWeightedMeanPrecision([-1.0, -1.0], [5.0 2.0; 2.0 5.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = PointMass([-2.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0])),
                output = (in1 = PointMass([-2.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([23.0, 16.0], [7.0  3.0; 3.0 5.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalWeightedMeanPrecision, m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = MvNormalWeightedMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = PointMass([1.0, -1.0]), m_in2 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, -1.0]), in2 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([2.0, 2.0]), in2 = MvNormalWeightedMeanPrecision([-4.0, -2.0], [5.0 2.0; 2.0 4.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [4.0 1.0; 1.0 4.0]), m_in1 = PointMass([-1.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 0.0; 0.0 2.0])),
                output = (in1 = PointMass([-1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([5.0, 0.0], [7.0  1.0; 1.0 6.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::NormalDistributionsFamily, m_in1::NormalDistributionsFamily, m_in2::NormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(+)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = NormalMeanVariance(3.0, 4.0), m_in2 = NormalMeanVariance(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([5/4, 4/3], [3/4 1/2; 1/2 2/3]))
            ),
            (
                input = (m_out = NormalMeanPrecision(1.0, 2.0), m_in1 = NormalMeanPrecision(3.0, 4.0), m_in2 = NormalMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([14.0, 32.0], [6.0 2.0; 2.0 8.0]))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(1.0, 2.0), m_in1 = NormalWeightedMeanPrecision(3.0, 4.0), m_in2 = NormalWeightedMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([4.0, 6.0], [6.0 2.0; 2.0 8.0]))
            ),
            (
                input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = NormalMeanPrecision(3.0, 4.0), m_in2 = NormalWeightedMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([25/2, 11/2], [9/2 1/2; 1/2 13/2]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalMeanCovariance([1.0, 2.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([1/5, 12/5, 0.0, 2.0], [[0.8 -0.4 0.4 -0.2; -0.4 1.2 -0.2 0.6]; [0.4 -0.2 0.8 -0.4; -0.2 0.6 -0.4 1.2]]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([10.0, 10.0, 10.0, 10.0], [6.0 2.0 3.0 1.0; 2.0 4.0 1.0 2.0; 3.0 1.0 6.0 2.0; 1.0 2.0 2.0 4.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0])),
                output = (MvNormalWeightedMeanPrecision([2.0, 4.0, 2.0, 4.0], [6.0 2.0 3.0 1.0; 2.0 2.0 1.0 1.0; 3.0 1.0 6.0 2.0; 1.0 1.0 2.0 2.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([21/5, 17/5, 6/5, 7/5], [17/5 4/5 2/5 -1/5; 4/5 13/5 -1/5 3/5; 2/5 -1/5 17/5 4/5; -1/5 3/5 4/5 13/5]))
            ),

        ]

    end

end
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--dim"
            help = "number of dimensions of the problem"
            default = 100
            arg_type = Int

        "--phi1"
            help = "effect of the local best for the velocity update"
            default = 2.0
            arg_type = Float64

        "--phi2"
            help = "effect of the global best for the velocity update"
            default = 2.0
            arg_type = Float64

        "--num"
            help = "number of the particles"
            default = 1024
            arg_type = Int

        "--epoch"
            help = "number of epochs"
            default = 1000
            arg_type = Int

        "--l"
            help = "lower bound of x"
            default = -100.
            arg_type = Float64

        "--h"
            help = "upper bound of x"
            default = 100.
            arg_type = Float64

        "--w"
            help = "inertia value"
            default = 0.8
            arg_type = Float64

        "--f"
            help = "function that will be optimized"
            default = "qing"

        "--vmin"
            help = "lower bound for the velocity"
            default = -1.0
            arg_type = Float64

        "--vmax"
            help = "upper bound for the velocity"
            default = 1.0
            arg_type = Float64

        "--vnorm"
            help = "norm constraint for the velocity"
            default = 1.0
            arg_type = Float64

        "--batchsize"
            help = "batch size"
            arg_type=Int
            default=100
        "--gpu"
            action = :store_true
    end
    return parse_args(s)
end


text:
println("===== Testing Yeppp!  =====")

## length of test array
N = 1000
Random.seed!(7)

## YepCore Signed Integer Tests
for T in (Int8, Int16, Int32, Int64)
    @testset "Yeppp: Integer Operations::$T" begin
        X = round.(T, randn(N))
        Y = round.(T, randn(N))
        @testset "Testing $fy::$T" for (f, fy) in [(+, :add), (-, :sub), (*, :mul)]
            @eval fb = $f
            @eval fyp = Vectorize.Yeppp.$fy
            @test fb.(X, Y) ≈ fyp(X, Y)
        end
    end
end

## YepCore Floating Point Tests
for T in (Float32, Float64)
    @testset "Yeppp: Floating Point Operations::$T" begin
        X = randn(T, N)
        Y = randn(T, N)
        @testset "Testing $fy::$T" for (f, fy) in [(+, :add), (-, :sub), (*, :mul)]
            @eval fb = $f
            @eval fyp = Vectorize.Yeppp.$fy
            @test fb.(X, Y) ≈ fyp(X, Y)
        end
    end
end

# YeppCore Max/Min Integer
for T in (Int8, Int16, Int32)
    @testset "Yeppp: Integer Max/Min::$T" begin
        X = round.(T, randn(N))
        Y = round.(T, randn(N))
        @testset "Testing $fy::$T" for (f, fy) in [(max, :max), (min, :min)]
            @eval fb = $f
            @eval fyp = Vectorize.Yeppp.$fy
            @test fb.(X, Y) ≈ fyp(X, Y)
        end
    end
end

# YeppCore Max/Min Floating Point
for T in (Float32, Float64)
    @testset "Yeppp: Floating Point Max/Min::$T" begin
        X = randn(T, N)
        Y = randn(T, N)
        @testset "Testing $fy::$T" for (f, fy) in [(max, :max), (min, :min)]
            @eval fb = $f
            @eval fyp = Vectorize.Yeppp.$fy
            @test fb.(X, Y) ≈ fyp(X, Y)
        end
    end
end

## YepCore Unsigned Integer Tests
# for (Tin, Tout) in ((UInt8, UInt16), (UInt16, UInt32), (UInt32, UInt64))
#     @testset "Yeppp: YepppCore::$Tin" begin
#         X = trunc(Tin, 10*abs(randn(N)) .+ 1)
#         Y = trunc(Tin, 10*abs(randn(N)) .+ 1)
#         @testset "Testing $fy::$Tin" for (f, fy) in [(.+, :add), (.-, :sub), (./, :div), (.*, :mul)]
#             @eval fb = $f
#             @eval fyp = Vectorize.Yeppp.$fy
#             println(trunc(Tout, abs(fb(X, Y))))
# #            println(trunc(T, fyp(X, Y)))
#             println("fyp: ",  trunc(Tout, fyp(X, Y)))
#             @test trunc(Tout, abs(fb(X, Y))) ≈ trunc(Tout, fyp(X, Y))
#         end
#     end
# end

# YepppMath Floating Point
for T in [Float64]
    @testset "Yeppp: Floating Point Math::$T" begin
        X = convert(Vector{T}, randn(N))
        @testset "Testing $f::$T" for f in [:sin, :cos, :tan]
            @eval fb = $f
            @eval fy = Vectorize.Yeppp.$f
            @test fb.(X) ≈ fy(X)
        end

        @testset "Testing $f::$T" for f in [:log, :exp]
            @eval fb = $f
            @eval fy = Vectorize.Yeppp.$f
            @test fb.(abs.(X)) ≈ fy(abs.(X))

        end
    end
end

# YepppMath Floating Point - return scalar
for T in (Float32, Float64)
    @testset "Yeppp: Floating Point Reductions::$T" begin
        X = convert(Vector{T}, randn(N))
        @testset "Testing $f::$T" for f in [:sum]
            @eval fb = $f
            @eval fy = Vectorize.Yeppp.$f
            @test fb(X) ≈ fy(X)
        end

        @testset "Testing $f::$T" for f in [:sumsqr]
            @eval fy = Vectorize.Yeppp.$f
            @test sum(X .* X) ≈ fy(X)
        end
    end
end

println("===== Yeppp! Tests Successful =====\n\n")


text:
using Plots
using Printf
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity

grid = RectilinearGrid(size=32, z=(-64, 0), topology=(Flat, Flat, Bounded))

closure = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1,
                                                  background_κz = 1e-5,
                                                  convective_νz = 1e-3,
                                                  background_νz = 1e-4)

Qᵇ = +1e-8
Qᵘ = +1e-5

b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = :b,
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = FPlane(f=1e-4),
                                    boundary_conditions = (; b=b_bcs, u=u_bcs),
                                    closure = closure)

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

z = znodes(model.tracers.b)

b = view(interior(model.tracers.b), 1, 1, :)
u = view(interior(model.velocities.u), 1, 1, :)
v = view(interior(model.velocities.v), 1, 1, :)

b_plot = plot(b, z, linewidth = 2, label = "t = 0", xlabel = "Buoyancy", ylabel = "z", legend=:bottomright)

u_plot = plot(u, z, linewidth = 2, label = "u at t = 0",
              linestyle = :solid, color = :black, xlabel = "velocity", ylabel = "z", legend = :bottomright)

plot!(u_plot, v, z, linewidth = 2, label = "v at t = 0",
      linestyle = :dash, color = :black)

simulation = Simulation(model, Δt = 20.0, stop_time = 48hours)

run!(simulation)

plot!(b_plot, b, z, linewidth = 2, label = @sprintf("t = %s", prettytime(model.clock.time)))

plot!(u_plot, u, z, label = @sprintf("u at t = %s", prettytime(model.clock.time)),
      linewidth = 2, linestyle = :dash, color = :blue)

plot!(u_plot, v, z, label = @sprintf("v at t = %s", prettytime(model.clock.time)),
      linewidth = 2, linestyle = :solid, color = :blue)

ub_plot = plot(b_plot, u_plot, layout=(1, 2))

display(ub_plot)


text:
using ITensors
using ITensorGLMakie
using Test

starts_and_ends_with(file, st, en) = startswith(file, st) && endswith(file, en)
starts_and_ends_with(st, en) = file -> starts_and_ends_with(file, st, en)

test_path = joinpath(@__DIR__)
test_files = filter(starts_and_ends_with("test_", ".jl"), readdir(test_path))
@testset "ITensorGLMakie.jl" for file in test_files
  file_path = joinpath(test_path, file)
  println("Running test $(file_path)")
  include(file_path)
end


text:
using ParametricMachinesDemos
using Documenter

makedocs(;
    modules=[ParametricMachinesDemos],
    authors="Pietro Vertechi <pietro.vertechi@protonmail.com> & Mattia G. Bergomi <mattiagbergomi@gmail.com>",
    repo="https://github.com/BeaverResearch/ParametricMachinesDemos.jl/blob/{commit}{path}#{line}",
    sitename="Parametric Machines",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
    ),
    pages=Any[
        "Home" => "index.md",
        "Machines" => "machines.md",
        "Internals" => "internals.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/BeaverResearch/ParametricMachinesDemos.jl",
    push_preview=true,
)


text:
using Basics
using Documenter

makedocs(;
    modules=[Basics],
    authors="Alexandre Gravier",
    repo="https://github.com/agravier/Basics.jl/blob/{commit}{path}#L{line}",
    sitename="Basics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://agravier.github.io/Basics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/agravier/Basics.jl",
)


text:
"""
    StateSpaceProblem{isinplace, dType, uType,tType,noiseType,oType, lType} <: AbstractStateSpaceProblem{isinplace}

A state space problem contains all the information necessary to simulate a state space
process.
"""
struct StateSpaceProblem{isinplace, dType, uType,tType,noiseType,oType} <: AbstractStateSpaceProblem{isinplace}
    f::Function # Evolution
    g::Function # Noise function
    h::Function # Observation function
    latent_noise::dType # Noise distribution
    u0::uType # Initial condition
    tspan::tType # Timespan to use
    noise::noiseType # The observation noise to use. No noise → perfect obs
    observables::oType # Data to use, if any
end

noisetype(prob::StateSpaceProblem{isinplace, dType}) where {isinplace, dType} = dType
hasobservables(prob::StateSpaceProblem{a,b,c,d,e,f}) where {a,b,c,d,e,f} = ~(f <: Nothing)

"""
    StateSpaceProblem(f, g, h, D, u0, tspan, noise, observables, likelihood=((x...) -> 0.0))
"""
function StateSpaceProblem(f, g, h, latent_noise, u0, tspan, noise=nothing, observables=nothing)
    return StateSpaceProblem{
        Val(false), # Default to out-of-place
        typeof(latent_noise),
        typeof(u0),
        typeof(tspan),
        typeof(noise),
        typeof(observables)
    }(f, g, h, latent_noise, u0, SciMLBase.promote_tspan(tspan), noise, observables)
end

"""
    index_or_nothing(x::AbstractArray, i::Int)
    index_or_nothing(x::Nothing, i::Int)

Extracts index `i` from `x`, if `x` is a subtype of `AbstractArray`. Otherwise,
this function returns `nothing`.
"""
index_or_nothing(x::AbstractArray, i::Int) = x[i]
index_or_nothing(x::Nothing, i::Int) = nothing

latent_noise(prob::StateSpaceProblem{inplace, <:AbstractArray}, t) where {inplace} = prob.latent_noise[t]
latent_noise(prob::StateSpaceProblem{inplace, <:Distribution}, t) where {inplace} = rand(prob.latent_noise)
latent_noise(prob::StateSpaceProblem{inplace, <:AbstractMatrix}, t) where {inplace} = prob.latent_noise * randn(size(prob.u0))
latent_noise(prob::StateSpaceProblem, t) = randn(size(prob.u0))

# No observations, return zero likelihood
function likelihood(
    prob::StateSpaceProblem,
    alg::ConditionalGaussian,
    z_t,
    obs_t::Nothing
)
    return 0
end

# Have observables, dispatch on inner likelihood
function likelihood(
    prob::StateSpaceProblem{isinplace, dType, uType,tType,noiseType,oType},
    alg::ConditionalGaussian,
    z_t,
    obs_t
) where {isinplace, dType, uType,tType,noiseType,oType}
    return _likelihood(prob, alg, z_t, obs_t, noiseType)
end

function _likelihood(
    prob::StateSpaceProblem,
    alg::ConditionalGaussian,
    z_t,
    obs_t,
    noiseType::Type{Nothing}
)
    # Unclear whether this should be default behavior
    return 0
end

function _likelihood(
    prob::StateSpaceProblem,
    alg::ConditionalGaussian,
    z_t,
    obs_t::Number,
    noiseType::Type{<:UnivariateDistribution}
)
    # Unclear whether this should be default behavior
    return logpdf(prob.noise, obs_t - only(z_t))
end

function _likelihood(
    prob::StateSpaceProblem,
    alg::ConditionalGaussian,
    z_t,
    obs_t::Number,
    noiseType::Type{<:MultivariateDistribution}
)
    # Unclear whether this should be default behavior
    return logpdf(prob.noise, [obs_t] .- z_t)
end

function _likelihood(
    prob::StateSpaceProblem,
    alg::ConditionalGaussian,
    z_t,
    obs_t,
    noiseType
)
    # If given a noise distribution, provide the log pdf w.r.t. the
    # distribution.
    throw(ArgumentError("No implementation for prob = StateSpaceProblem, alg = ConditionalGaussian, and noiseType = $noiseType"))
end

"""
    _solve(prob::StateSpaceProblem, params)

Returns a `StateSpaceSolution` using the provided parameters `params`.
Calculates the likelihood if the likelihood function is available.
"""
function _solve(
    prob::StateSpaceProblem{isinplace, dType, uType,tType,noiseType,oType},
    alg::ConditionalGaussian,
    params
) where {isinplace, dType, uType,tType,noiseType,oType <: Nothing} # No observables
    t0, T = prob.tspan

    # Initial noise
    z0 = prob.h(prob.u0, params, 1)

    # Preallocate arrays
    u = Vector{typeof(prob.u0)}(undef, T+1)
    z = Vector{typeof(z0)}(undef, T+1)
    n = noisetype(prob) <: Vector ?
        Vector{eltype(prob.noise)}(undef, T+1) :
        Vector{typeof(prob.u0)}(undef, T+1) # Store sampled transition noise

    u[1] = prob.u0
    z[1] = z0

    for t in 2:T+1
        n[t-1] = latent_noise(prob, t-1)
        u[t] = prob.f(u[t - 1], params, t-1) + prob.g(u[t - 1], params, t-1) .* n[t-1] # Should this be n[t] or n[t-1]?
        z[t] = prob.h(u[t], params, t)
    end

    return StateSpaceSolution(copy(z),copy(u),prob.noise,nothing,nothing,n)
end