text:
module StaticNumbers

using Requires
import Base.Broadcast: broadcasted, DefaultArrayStyle

export Static, static,
       StaticBool, StaticInteger, StaticReal, StaticNumber, StaticOrInt, StaticOrBool,
       @staticnumbers, @generate_static_methods, unstatic

function __init__()
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" include("StaticArrays_glue.jl")
    @require SIMD="fdea26ae-647d-5447-a871-4b548cad5224" include("SIMD_glue.jl")
end

const StaticError = ErrorException("Illegal type parameter for Static.")

"""
A `StaticInteger` is an `Integer` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticInteger{X} <: Integer
    function StaticInteger{X}() where {X}
        X isa Integer && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
@inline StaticInteger{X}(x::Number) where {X} = x==X ? StaticInteger{X}() : throw(InexactError(:StaticInteger, StaticInteger{X}, x))
@inline StaticInteger(x::Number) = StaticInteger{Integer(x)}()

"""
A `StaticReal` is a `Real` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticReal{X} <: Real
    function StaticReal{X}() where {X}
        X isa Real && !(X isa Integer) && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
@inline StaticReal{X}(x::Number) where {X} = x==X ? StaticReal{X}() : throw(InexactError(:StaticReal, StaticReal{X}, x))
@inline StaticReal(x::Number) = StaticReal{Real(x)}()

"""
A `StaticNumber` is a `Number` whose value is stored in the type, and which
contains no runtime data.
"""
struct StaticNumber{X} <: Number
    function StaticNumber{X}() where {X}
        X isa Number && !(X isa Real) && !(X isa Static) && isimmutable(X) || throw(StaticError)
        new{X}()
    end
end
@inline StaticNumber{X}(x::Number) where {X} = x==X ? StaticNumber{X}() : throw(InexactError(:StaticNumber, StaticReal{X}, x))
@inline StaticNumber(x::Number) = StaticNumber{x}()

"""
`Static{X}` is short-hand for the `Union` of `StaticInteger{X}`, `StaticReal{X}`
and `StaticNumber{X}`.
"""
const Static{X} = Union{StaticInteger{X}, StaticReal{X}, StaticNumber{X}}

# We'll define this constructor, but not recommend it.
Static{X}() where X = static(X)

# This is the recommended constructor.
"""
`static(X)` is shorthand for `StaticInteger{X}()`, `StaticReal{X}()` or `StaticNumber{X}()`,
depending on the type of `X`.
"""
@inline static(x::Integer) = StaticInteger{x}()
@inline static(x::Real) = StaticReal{x}()
@inline static(x::Number) = StaticNumber{x}()
@inline static(x::Bool) = x ? StaticInteger{true}() : StaticInteger{false}() # help inference
@inline static(x::StaticInteger) = x
@inline static(x::StaticReal) = x
@inline static(x::StaticNumber) = x

Base.Irrational{X}(::StaticReal{Y}) where {X, Y} = Irrational{X}(Y)

# There's no point crating a Val{Static{X}} since any function that would accept
# it should treat it as equivalent to Val{X}.
# But maybe it is a bad idea to add a method to a @pure function?
# @inline Base.Val(::Static{X}) where X = Val(X)

"""
`StaticNumbers.map(f, t::Tuple)` works the same as `Base.map(f, t)`, but does not stop
in-lining after a certain tuple length. Long tuples will cause code blowup.
"""
# Code copied from julia/base/tuples.jl and slightly modified
@inline map(f::F, t::Tuple{}) where {F} = ()
@inline map(f::F, t::Tuple) where {F} = (f(t[1]), map(f, Base.tail(t))...)
@inline map(f::F, t::Tuple{}, s::Tuple{}) where {F} = ()
@inline map(f::F, t::Tuple, s::Tuple) where {F} = (f(t[1],s[1]), map(f, Base.tail(t), Base.tail(s))...)
@inline heads(ts::Tuple...) = map(first, ts)
@inline tails(ts::Tuple...) = map(Base.tail, ts)
@inline map(f, ::Tuple{}...) = ()
@inline map(f::F, t1::Tuple{Any,Vararg{Any,N}}, t2::Tuple{Any,Vararg{Any,N}}, ts::Tuple{Any,Vararg{Any,N}}...) where {F,N} = (f(heads(t1, t2, ts...)...), map(f, tails(t1, t2, ts...)...)...)
@inline map(args...) = Base.map(args...) # fallback for non-tuples

"""
`ntuple(f, static(n))` yields the tuple `(f(static(1)), f(static(2)), ..., f(static(n)))`.
This will ususally yield the same result as `ntuple(f, n)` but may make the compiler try
even harder to constant-propagate the integer input into `f`.
"""
# Recursive version only worked for tuple lengths up to 32, so resorting to @generated.
@generated function Base.ntuple(f::F, ::StaticInteger{N}) where {F,N}
    N::Int
    (N >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", N)))
    v = Expr(:tuple, [ :( f(StaticInteger{$i}()) ) for i in 1:N ]... )
    :( Base.@_inline_meta; $v )
end

"""
`StaticOrInt` is the type union of `Int` and `StaticInteger`

(Functions that take only `Int` may be too restrictive.)
"""
const StaticOrInt = Union{StaticInteger, Int}

# Promotion
# We need to override promote and promote_typeof because they don't even call
# promote_rule for all-same types.
Base.promote(::ST, ys::ST...) where {ST <: Static{X}} where {X} = ntuple(i->X, static(1+length(ys)))
Base.promote_type(::Type{ST}, ::Type{ST})  where {ST <: Static{X}} where {X} = typeof(X)
Base.promote_typeof(::ST, ::ST...) where {ST <: Static{X}} where {X} = typeof(X)

# To avoid infinite recursion, we need this:
Base.promote_type(::Type{<:Static{X}}, T::Type...) where {X} = promote_type(typeof(X), promote_type(T...))

# Loop over all three types specifically, instead of dispatching on the Union.
for ST in (StaticInteger, StaticReal, StaticNumber)
    Base.promote_rule(::Type{ST{X}}, ::Type{T}) where {X,T<:Number} = promote_type(typeof(X), T)

    # Constructors
    (::Type{Complex{T}})(::ST{X}) where {T<:Real, X} = Complex{T}(X)
    (::Type{Rational{T}})(::ST{X}) where {T<:Integer, X} = Rational{T}(X)
end

# TODO: Constructors to avoid Static{Static}

# Some of the more common constructors that do not default to `convert`
# Note:  We cannot have a (::Type{T})(x::Static) where {T<:Number} constructor
# instead of all of these, because of ambiguities with user-defined types.
for T in (:Bool, :Integer, :AbstractFloat, :Unsigned, :Signed,
          :BigInt, :Int128, :Int16, :Int32, :Int64, :Int8,
          :UInt128, :UInt16, :UInt32, :UInt64, :UInt8,
          :BigFloat, :Float16, :Float32, :Float64)
    @eval Base.$T(::Static{X}) where {X} = $T(X)
end

# big(x) still defaults to convert.

# Single-argument functions that do not already work.
# Note: We're not attempting to support every function in Base.
# TODO: Should have a macro for easily extending support.
for fun in (:-, :zero, :one, :oneunit, :trailing_zeros, :decompose)
    @eval Base.$fun(::Static{X}) where X = Base.$fun(X)
end
for fun in (:trunc, :floor, :ceil, :round, :isnan)
    @eval Base.$fun(::Union{StaticReal{X}, StaticNumber{X}}) where {X} = Base.$fun(X)
end
for fun in (:zero, :one, :oneunit)
    @eval Base.$fun(::Type{<:Static{X}}) where {X} = Base.$fun(typeof(X))
end
Base.widen(::Static{X}) where {X} = X isa Bool ? X : widen(X)

# It's a pity there's no AbstractBool supertype.
"StaticBool is a shorthand for Union{StaticInteger{false}, StaticInteger{true}}"
const StaticBool = Union{StaticInteger{false}, StaticInteger{true}}
"StaticOrBool can be either a `Bool` or a `StaticBool`"
const StaticOrBool = Union{StaticBool, Bool}

Base.:!(x::StaticBool) = !Bool(x)

# false is a strong zero
for T in (Integer, Rational, Real, Number, Complex{<:Real}, Complex{Bool}, StaticInteger, StaticReal, StaticNumber)
    Base.:*(::StaticInteger{false}, y::T) = zero(y)
    Base.:*(x::T, ::StaticInteger{false}) = zero(x)
end
Base.:*(::StaticInteger{false}, ::StaticInteger{false}) = false # disambig

# Handle static(Inf)*false
Base.:*(x::Bool, ::StaticReal{Y}) where Y = x*Y
Base.:*(::StaticReal{X}, y::Bool) where X = X*y

# Until https://github.com/JuliaLang/julia/pull/32117 is merged
Base.:*(::StaticInteger{false}, ::AbstractIrrational) = 0.0
Base.:*(::AbstractIrrational, ::StaticInteger{false}) = 0.0

# For complex-valued inputs, there's no auto-convert to floating-point.
# We only support a limited subset of functions, which the user can extend
# as needed.
# TODO: Should have a macro for making functions accept Static input.
for fun in (:abs, :abs2, :cos, :sin, :exp, :log, :isinf, :isfinite, :isnan)
    @eval Base.$fun(::StaticNumber{X}) where {X} = Base.$fun(X)
end
Base.sign(::StaticInteger{X}) where {X} = Base.sign(X) # work around problem with Bool

# Other functions that do not already work
Base.:(<<)(::StaticInteger{X}, y::UInt64) where {X} = X << y
Base.:(>>)(::StaticInteger{X}, y::UInt64) where {X} = X >> y

# Two-argument functions that have methods in promotion.jl that give no_op_err:
for f in (:+, :-, :*, :/, :^)
    @eval Base.$f(::ST, ::ST) where {ST <: Static{X}} where {X} = $f(X,X)
end
# ...where simplifications are possible:
# Note: We allow creation of specific static numbers, like 1 and 0 (as an exception)
# since this cannot lead to the set of static numbers growing uncontrollably.
Base.:&(x::ST, ::ST) where {ST<:StaticInteger} = x
Base.:|(x::ST, ::ST) where {ST<:StaticInteger} = x
Base.xor(::ST, ::ST) where {ST<:StaticInteger} = static(zero(ST))
Base.rem(::ST, ::ST) where {ST<:Static{X}} where {X} = (X==0 || isinf(X)) ? X isa AbstractFloat ? static(oftype(X, NaN)) : throw(DivideError()) : static(zero(X))
Base.mod(::ST, ::ST) where {ST<:Static{X}} where {X} = (X==0 || isinf(X)) ? X isa AbstractFloat ? static(oftype(X, NaN)) : throw(DivideError()) : static(zero(X))
Base.div(::ST, ::ST) where {ST<:Static{X}} where {X} = (X == 0 || isinf(X)) ? ( X isa Integer ? throw(DivideError()) : oftype(X, NaN) ) : static(one(X)) # Needed for Julia > 1.3
Base.:<(::ST, ::ST) where {ST<:StaticReal{X}} where {X} = false
Base.:<=(::ST, ::ST) where {ST<:StaticReal{X}} where {X} = true
# Bypass promotion in comparisons involving static unsigned integers
for fun in (:(<), :(<=))
    @eval Base.$fun(::StaticInteger{X}, y::Integer) where {X} = $fun(X, y)
    @eval Base.$fun(x::Integer, ::StaticInteger{Y}) where {Y} = $fun(x, Y)
    @eval Base.$fun(::StaticInteger{X}, y::BigInt) where {X} = $fun(X, y)
    @eval Base.$fun(x::BigInt, ::StaticInteger{Y}) where {Y} = $fun(x, Y)
    @eval Base.$fun(::StaticInteger{X}, ::StaticInteger{Y}) where {X,Y} = $fun(X, Y)
end

Base.:(==)(x::AbstractIrrational, ::StaticReal{Y}) where {Y} = x == Y
Base.:(==)(::StaticReal{X}, y::AbstractIrrational) where {X} = X == y

# Three-argument function that gives no_op_err
Base.fma(::ST, ::ST, ::ST) where {ST<:Static{X}} where {X} = fma(X,X,X)

# Static powers using Base.literal_pow.
# This avoids DomainError in some cases?
for T in (Bool, Int32, Int64, Float32, Float64, ComplexF32, ComplexF64, Irrational)
    Base.:^(x::T, ::StaticInteger{p}) where {p} = Base.literal_pow(^, x, Val(p))
end
Base.:^(x::Static{X}, ::StaticInteger{p}) where {X,p} = Base.literal_pow(^, X, Val(p))
Base.:^(x::ST, ::ST) where {ST<:StaticInteger{p}} where {p} = Base.literal_pow(^, x, Val(p)) #disambig
Base.:^(::Irrational{:ℯ}, p::StaticInteger) = exp(p) #disambig

# For brevity, all `Static` numbers are displayed as `static(X)`, rather than, for
# example, `StaticInteger{X}()`. It is possible to discern between the different
# types of `Static` by looking at `X`.
# To get the default behaviour back, run:
#   methods(Base.show, (IO, Static{X} where X)) |> first |> Base.delete_method
function Base.show(io::IO, x::Static{X}) where X
    print(io, "static(")
    show(io, X)
    print(io, ")")
end

for f in (:isfinite, :isnan, :isinf)
    @eval Base.$f(::StaticReal{X}) where {X} = $f(X)
end

# Dont promote when it's better to treat real and imaginary parts separately
Base.:/(::StaticNumber{X}, y::Real) where {X} = X/y
Base.:*(::StaticNumber{X}, y::Real) where {X} = X*y
Base.:*(x::Real, ::StaticNumber{Y}) where {Y} = x*Y

@inline function Base.setindex(x::Tuple, v, i::StaticInteger)
    @boundscheck 1 <= i <= length(x) || throw(BoundsError(x, i))
    return ntuple(n -> n == i ? v : x[n], static(length(x)))
end

Base.oftype(::Static{X}, y) where {X} = oftype(X, y)

"""
`unstatic(x)` returns a non-static version of `x`.
This function is rarely needed, as most operations on a static number (e.g. `x+0`)
will yield a non-static result.
"""
unstatic(x) = x
unstatic(::Static{X}) where {X} = X

include("generate_static_methods.jl")

include("stat_macro.jl")

include("LengthRanges.jl")

include("trystatic.jl")

include("deprecated.jl")

end # module


text:
using Random

module Why
QUERY_COUNT = 0

function why()
    return make_ans()
end

const begins = ("because ", "")
const fins = ("wanted it", "decided so", "made it happen", "caused it", "could not prevent it", "permited it")
const mess_persistent = ("Why do you keep asking these questions?", "Please stop.")
const single = ("I do not know.", "42", "You must find that answer yourself.", "It is axiomatic.")
function make_ans()
    if (global QUERY_COUNT+=1) > 5 + rand(1:15)
        global QUERY_COUNT -= rand(1:10)
        return message = mess_persistent[rand(1:length(mess_persistent))]
    elseif rand() < 0.15
        return message = single[rand(1:length(single))]
    end
    return uppercasefirst(rand_element(begins) * noun(1.0) * rand_element(fins) * ".")
end

#Nouns
const preps = ("the", "a", "some")
const nouns = ("compiler", "person", "research software engineer", "type pirate", "cockatoo", "manager", "programmer", "software", "penguin", "config file", "divine being", "struct")
function noun(p)
    (rand()>p) && return ""
    return rand_element(preps)*" "*adj(0.80)*rand_element(nouns)*" "*and(noun(p/2))
end

#Ajectives
const adjs = ("fast", "type-safe", "harmless", "deceitful", "cool", "self aware", "scary", "insecure", "paranoid", "cunning", "irrelevant", "undocumented", "rudimentary", "primitive")
function adj(p)
    (rand()>p) && return ""
    return adv(0.5)*rand_element(adjs)*" "*and(adj(p/2.5))
end

#Adverbs
const advs = ("vaguely", "extraordinarily", "properly", "very", "not at all", "probably")
function adv(p)
    (rand()>p) && return ""
    return rand_element(advs)*" "
end

#Extra functions
function and(s::String)
    (s=="") && return s
    return "and " * s
end

rand_element(tup::Tuple) = tup[rand(1:length(tup))]

export why

#Preparing for expanding.

#Nouns
#const proper_nouns = ("Bob", "Alice", "Kim")
#const concrete_nouns = ("compiler", "person", "research software engineer", "type pirate", "cockatoo", "manager", "programmer", "software", "penguin", "config file", "divine being")
#const abstract_nouns = ("truth", "joy")
#const collective_nouns = ("online community")

#Verbs

#Adjective
#const descriptive_adjectives = ("scary", "cool", "backwards")
#const quantitative_adjectives = ("one", "two", "seventeen", "many", "few", "several")
#const possessive_adjectives = ("their")
#const distributive_adjectives = ("every", "any")
#const articles = ("a", "the", "some") # "An"?

#Adverbs
#const time_adverbs = ("never", "soon", "recentyl")
#const space_adverbs = ("there", "somewhere", "everywhere")
#const manner_adverbs = ("polietly", "neatly", "promptly")
#const degree_adverbs = ("quite", "too", "so")
#const ferquency_adverbs = ("seldom", "usually", "again")

#Pronouns

#Prepositions

#Conjuctions

#Determiners

#Exclamations

end # module

ext: jl


text:
import TexasHoldem
const TH = TexasHoldem
using TexasHoldem
using BenchmarkTools
using Logging

struct BotCheckCall <: AbstractAI end

TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CheckRaiseFold) = check!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallRaiseFold) = call!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallAllInFold) = call!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallFold) = call!(game, player)

players() = ntuple(i->(Player(BotCheckCall(), i)), 4)

with_logger(NullLogger()) do
    play!(Game(players()))
end

using Profile
@profile begin
    for k in 1:100
        with_logger(NullLogger()) do
            play!(Game(players()))
        end
    end
end

Profile.print(format=:flat, sortedby=:count)


text:
# This file is part of Kpax3. License is MIT.

function test_plotk()
  (k, pk) = Kpax3.readposteriork("../build/mcmc_6")
  p = Kpax3.kpax3plotk(k, pk, xticks=[1; 2; 3; 4; 5; 6])

  Plots.png(p, "../build/MCMC_posterior_k.png")

  nothing
end

test_plotk()

function test_plotc()
  (site, aa, freq, C) = Kpax3.readposteriorC("../build/mcmc_6")

  p = Kpax3.kpax3plotc(site, freq, C)

  Plots.png(p, "../build/MCMC_posterior_column_classifier.png")

  nothing
end

test_plotc()

function test_plotd()
  settings = Kpax3.KSettings("data/mcmc_6.fasta", "../build/tmp")
  x = Kpax3.AminoAcidData(settings)
  state = Kpax3.optimumstate(x, "data/mcmc_6.csv", settings)

  p = Kpax3.kpax3plotd(x, state, clusterorder=[4; 2; 1; 3],
                       clusterlabel=["d"; "b"; "a"; "c"])

  Plots.png(p, "../build/MCMC_posterior_dataset.png")

  nothing
end

test_plotd()

function test_plotp()
  (id, P) = Kpax3.readposteriorP("../build/mcmc_6")
  R = Kpax3.normalizepartition("data/mcmc_6.csv", id)

  p = Kpax3.kpax3plotp(R, P, clusterorder=[4; 2; 1; 3],
                       clusterlabel=["d"; "b"; "a"; "c"])

  Plots.png(p, "../build/MCMC_posterior_probability_matrix.png")

  nothing
end

test_plotp()

function test_plottrace()
  (entropy_R, avgd_R) = Kpax3.traceR("../build/mcmc_6", maxlag=50)
  (entropy_C, avgd_C) = Kpax3.traceC("../build/mcmc_6", maxlag=50)

  p = Kpax3.kpax3plottrace(entropy_R, maxlag=50, main="Trace")
  q = Kpax3.kpax3plottrace(entropy_C, maxlag=50, main="Trace")

  Plots.png(p, "../build/MCMC_posterior_trace_R.png")
  Plots.png(q, "../build/MCMC_posterior_trace_C.png")

  nothing
end

test_plottrace()

function test_plotdensity()
  (entropy_R, avgd_R) = Kpax3.traceR("../build/mcmc_6", maxlag=50)
  (entropy_C, avgd_C) = Kpax3.traceC("../build/mcmc_6", maxlag=50)

  p = Kpax3.kpax3plotdensity(entropy_R, maxlag=50, main="Density")
  q = Kpax3.kpax3plotdensity(entropy_C, maxlag=50, main="Density")

  Plots.png(p, "../build/MCMC_posterior_density_R.png")
  Plots.png(q, "../build/MCMC_posterior_density_C.png")

  nothing
end

test_plotdensity()

function test_plotjump()
  (entropy_R, avgd_R) = Kpax3.traceR("../build/mcmc_6", maxlag=50)
  (entropy_C, avgd_C) = Kpax3.traceC("../build/mcmc_6", maxlag=50)

  p = Kpax3.kpax3plotjump(avgd_R, main="Jump distance")
  q = Kpax3.kpax3plotjump(avgd_C, main="Jump distance")

  Plots.png(p, "../build/MCMC_posterior_jump_R.png")
  Plots.png(q, "../build/MCMC_posterior_jump_C.png")

  nothing
end

test_plotjump()


text:
@doc raw"""
    PSFModels.Moffat([T=Float64]; fwhm, x=0, y=0, amp=1, alpha=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.Moffat([T=Float64]; fwhm, pos, amp=1, alpha=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.Moffat([T=Float64]; fwhm, r, theta, amp=1, alpha=1, maxsize=3, extent=maxsize .* fwhm)

Two dimensional Moffat model. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. By default the model is placed at the origin. The position can also be given as a polar coordinate using `r`/`ρ` and `theta`/`θ`, optionally centered around `origin`.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). For efficient calculations, we recommend using [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl). Here, `maxsize` is a multiple of the fwhm, and can be given as a scalar or as a tuple for each axis. The `extent` defines the bounding box for the model and is used for the default rendering size.

# Functional form
```
f(x | x̂, FWHM, α) = A / (1 + ||x - x̂|| / (FWHM / 2)^2)^α
```
where `x̂` and `x` are position vectors (indices) `||⋅||` represents the square-distance, and `FWHM` is the full width at half-maximum. If `FWHM` is a vector or tuple, the weighting is applied along each axis.
"""
struct Moffat{T,FT,VT<:AbstractVector,IT<:Tuple} <: PSFModel{T}
    fwhm::FT
    pos::VT
    amp::T
    alpha::Int
    indices::IT

    function Moffat(::Type{T}, fwhm::FT, pos::VT, amp::T, alpha, indices::IT) where {T,FT,VT<:AbstractVector,IT<:Tuple}
        new{T,FT,VT,IT}(fwhm, pos, amp, alpha, indices)
    end
end

## constructors
# default type is Float64
Moffat(; kwargs...) = Moffat(Float64; kwargs...)
function Moffat(T; fwhm, amp=one(T), alpha=1, maxsize=3, extent = maxsize .* fwhm, position...)
    # get the position from keyword distpatch
    pos = _position(; position...)
    return Moffat(T, fwhm, pos, convert(T, amp), alpha, indices_from_extent(pos, extent))
end

Base.size(m::Moffat) = map(length, m.indices)
Base.axes(m::Moffat) = m.indices

# scalar case
function (m::Moffat{T})(point::AbstractVector) where T
    hwhm = m.fwhm / 2
    Δ = sqeuclidean(point, m.pos)
    val = m.amp / (1 + Δ / hwhm^2)^m.alpha
    return convert(T, val)
end

# vector case
function (m::Moffat{T,<:Union{AbstractVector,Tuple}})(point::AbstractVector) where T
    weights = SA[(2 / first(m.fwhm))^2, (2 / last(m.fwhm))^2]
    Δ = wsqeuclidean(point, m.pos, weights)
    val = m.amp  / (1 + Δ)^m.alpha
    return convert(T, val)
end


text:
using Documenter
using SpinMonteCarlo

makedocs( format=Documenter.HTML(),
          sitename="SpinMonteCarlo.jl",
          modules=[SpinMonteCarlo],
          pages = [
                   "Home" => "index.md",
                   "Manual" => Any[
                                   "Run Monte Carlo" => "runmc.md",
                                   "Generate lattice" => "lattice.md",
                                   "Develop Monte Carlo" => "develop.md",
                                  ],
                   "Library" => Any[
                                    "Public" => "lib/public.md",
                                    "Internals" => "lib/internals.md",
                                   ]
                  ]
        )

deploydocs(
           repo = "github.com/yomichi/SpinMonteCarlo.jl.git",
           target = "build"
          )


text:
import ACE: OneParticleBasis, AbstractState, Scal1pBasis, VectorPool
import ACE.OrthPolys: TransformedPolys, transformed_jacobi
import ACE.Transforms: IdTransform
import NamedTupleTools
using NamedTupleTools: namedtuple

export Pop1PBasis

@doc raw"""
`struct Pop1PBasis <: OneParticleBasis`

One-particle basis similar to Scal1pBasis, but it has always length 1 such that P(x) = x
"""
mutable struct Pop1PBasis <: OneParticleBasis{Float64}
  # P::TransformedPolys
end

symbols(basis::Pop1PBasis) = [:P]

indexrange(basis::Pop1PBasis) = Dict( :P => 1:1 ) # Dict( :P => 1:length(basis) )

isadmissible(b, basis::Pop1PBasis) = (b.P == 1) # (1 <= b.P <= length(basis))

get_index(basis::Pop1PBasis, b) = 1

#function Pop1PBasis() #(num_basis::Int)
  #P = transformed_jacobi(num_basis, IdTransform(), 120.0, 0.0)
#  return Pop1PBasis(P)
#end

Base.length(basis::Pop1PBasis) = 1  # length(basis.P)

_val(X::AbstractState, basis::Pop1PBasis) = X.pop

# Dict(:C => (4.0, 7.0), :H => (0.0, 2.0))

# (x-(a+b)/2)/(b-a)

ACE.valtype(basis::Pop1PBasis, args...) = Float64

ACE.gradtype(basis::Pop1PBasis, X::AbstractState) = ACE.dstate_type(valtype(basis, X), X)

function ACE.evaluate!(B, basis::Pop1PBasis, x::Number)
  B[1] = x
  return B
  #return evaluate!(B, basis.P, x)
end

function ACE.evaluate!(B, basis::Pop1PBasis, X::AbstractState)
  return evaluate!(B, basis, _val(X, basis))
end

function ACE.evaluate_ed!(B, dB, basis::Pop1PBasis, X::AbstractState)
   evaluate!(B, basis, X)
   dB[1] = zero(eltype(dB))
   return B, dB
end

degree(b, basis::Pop1PBasis, args...) = 0 #degree(basis.P)

write_dict(V::Pop1PBasis) =
      Dict( "__id__" => "ACE_Pop1PBasis")

function read_dict(::Val{:ACE_Pop1PBasis}, D::Dict)
   return Pop1PBasis()
end


text:
const CIENV = get(ENV, "TRAVIS", "") == "true"  || get(ENV, "CI", "") == "true"

using DocStringExtensions, LinearAlgebra, LogDensityProblems, OffsetArrays, Parameters,
    Random, Test, TransformVariables, StaticArrays
import Tracker, ForwardDiff
using LogDensityProblems: logdensity, logdensity_and_gradient
using TransformVariables:
    AbstractTransform, ScalarTransform, VectorTransform, ArrayTransform,
    unit_triangular_dimension, logistic, logistic_logjac, logit, inverse_and_logjac

include("utilities.jl")

Random.seed!(1)

####
#### utilities
####

@testset "misc utilities" begin
    @test unit_triangular_dimension(1) == 0
    @test unit_triangular_dimension(2) == 1
    @test unit_triangular_dimension(5) == 10
end

@testset "logistic and logit" begin
    for _ in 1:1000
        x = randn(Float64) * 50
        bx = BigFloat(x)
        lbx = 1/(1+exp(-bx))
        @test logistic(x) ≈ lbx
        ljx = logistic_logjac(x)
        ljbx = -(log(1+exp(-bx))+log(1+exp(bx)))
        @test ljx ≈ ljbx rtol = eps(Float64)
    end
    for _ in 1:1000
        y = rand(Float64)
        @test logistic(logit(y)) ≈ y
    end
end

####
#### scalar transformations correctness checks
####

@testset "scalar transformations consistency" begin
    for _ in 1:100
        a = randn() * 100
        test_transformation(as(Real, -∞, a), y -> y < a)
        test_transformation(as(Real, a, ∞), y -> y > a)
        b = a + 0.5 + rand(Float64) + exp(randn() * 10)
        test_transformation(as(Real, a, b), y -> a < y < b)
    end
    test_transformation(as(Real, -∞, ∞), _ -> true)
end

@testset "scalar transformation corner cases" begin
    @test_throws ArgumentError as(Real, "a fish", 9)
    @test as(Real, 1, 4.0) == as(Real, 1.0, 4.0)
    @test_throws ArgumentError as(Real, 3.0, -4.0)

    t = as(Real, 1.0, ∞)
    @test_throws DomainError inverse(t, 0.5)

    t = as(Real, -∞, 10.0)
    @test_throws DomainError inverse(t, 11.0)

    t = as(Real, 1.0, 10.0)
    @test_throws DomainError inverse(t, 0.5)
    @test_throws DomainError inverse(t, 11.0)
    @test_throws DomainError inverse_and_logjac(t, 0.5)
    @test_throws DomainError inverse_and_logjac(t, 11.0)
end

@testset "scalar alternatives" begin
    @test as_real ≡ asℝ
    @test as_positive_real ≡ asℝ₊
    @test as_negative_real ≡ asℝ₋
    @test as_unit_interval ≡ as𝕀
end

####
#### special array transformation correctness checks
####

@testset "to unit vector" begin
    @testset "dimension checks" begin
        U = UnitVector(3)
        x = zeros(3)               # incorrect
        @test_throws ArgumentError U(x)
        @test_throws ArgumentError transform(U, x)
        @test_throws ArgumentError transform_and_logjac(U, x)
    end

    @testset "consistency checks" begin
        for K in 1:10
            t = UnitVector(K)
            @test dimension(t) == K - 1
            if K > 1
                test_transformation(t, y -> sum(abs2, y) ≈ 1,
                                    vec_y = y -> y[1:(end-1)])
            end
        end
    end
end

@testset "to unit simplex" begin
    @testset "dimension checks" begin
        S = UnitSimplex(3)
        x = zeros(3)               # incorrect
        @test_throws ArgumentError S(x)
        @test_throws ArgumentError transform(S, x)
        @test_throws ArgumentError transform_and_logjac(S, x)
    end

    @testset "consistency checks" begin
        for K in 1:10
            t = UnitSimplex(K)
            @test dimension(t) == K - 1
            if K > 1
                test_transformation(t, y -> (sum(y) ≈ 1) & (all(y.>=0)),
                                    vec_y = y -> y[1:(end-1)])
            end
            x = zeros(dimension(t))
            @test transform(t, x) ≈ 1 ./ fill(K, K)
        end
    end
end

@testset "to correlation cholesky factor" begin
    @testset "dimension checks" begin
        C = CorrCholeskyFactor(3)
        wrong_x = zeros(dimension(C) + 1)

        @test_throws ArgumentError C(wrong_x)
        @test_throws ArgumentError transform(C, wrong_x)
        @test_throws ArgumentError transform_and_logjac(C, wrong_x)
    end

    @testset "consistency checks" begin
        for K in 1:8
            t = CorrCholeskyFactor(K)
            @test dimension(t) == (K - 1)*K/2
            CIENV && @info "testing correlation cholesky K = $(K)"
            if K > 1
                test_transformation(t, is_valid_corr_cholesky;
                                    vec_y = vec_above_diagonal, N = 100)
            end
        end
    end
end

####
#### aggregation
####

###
### array correctness checks
###

@testset "to array scalar" begin
    dims = (3, 4, 5)
    t = as𝕀
    ta = as(Array, t, dims...)
    @test dimension(ta) == prod(dims)
    x = random_arg(ta)
    y = @inferred transform(ta, x)
    @test typeof(y) == Array{Float64, length(dims)}
    @test size(y) == dims
    @test inverse(ta, y) ≈ x
    ℓacc = 0.0
    for i in 1:length(x)
        yi, ℓi = @inferred transform_and_logjac(t, x[i])
        @test yi == y[i]
        ℓacc += ℓi
    end
    y2, ℓ2 = transform_and_logjac(ta, x)
    @test y == y2
    @test ℓ2 ≈ ℓacc
end

@testset "as array fallback" begin
    is_expected(t, dims) = t isa ArrayTransform && t.transformation == asℝ && t.dims == dims
    @test is_expected(as(Array, 2, 3), (2, 3))
    @test is_expected(as(Array, (2, 3)), (2, 3))
    @test is_expected(as(Matrix, 2, 3), (2, 3))
    @test is_expected(as(Matrix, (2, 3)), (2, 3))
    @test_throws ArgumentError as(Vector, 2, 3)
    @test_throws ArgumentError as(Vector, (2, 3))
    @test is_expected(as(Vector, 2), (2, ))
end

@testset "as array w/ Identity" begin
    d = (2, 3, 4)
    t = as(Array, d)
    v = random_arg(t)
    y1 = transform(t, v)
    y2, lj = transform_and_logjac(t, v)
    @test y1 == y2 == reshape(v, d)
    @test lj == 0
end

###
### tuple correctness checks
###

@testset "to tuple" begin
    t1 = asℝ
    t2 = as𝕀
    t3 = CorrCholeskyFactor(7)
    tt = as((t1, t2, t3))
    @test dimension(tt) == dimension(t1) + dimension(t2) + dimension(t3)
    x = random_arg(tt)
    y = @inferred transform(tt, x)
    @test inverse(tt, y) ≈ x
    TransformVariables.inverse_eltype(tt, y)
    index = 0
    ljacc = 0.0
    for (i, t) in enumerate((t1, t2, t3))
        d = dimension(t)
        xpart = t isa ScalarTransform ? x[index + 1] : x[index .+ (1:d)]
        @test y[i] == transform(t, xpart)
        ypart, ljpart = transform_and_logjac(t, xpart)
        @test ypart == y[i]
        ljacc += ljpart
        index += d
    end
    y2, lj2 = transform_and_logjac(tt, x)
    @test y == y2
    @test lj2 ≈ ljacc
end

###
### named tuple correctness checks
###

@testset "to named tuple" begin
    t1 = asℝ
    t2 = CorrCholeskyFactor(7)
    t3 = UnitVector(3)
    tn = as((a = t1, b = t2, c = t3))
    @test dimension(tn) == dimension(t1) + dimension(t2) + dimension(t3)
    x = randn(dimension(tn))
    y = @inferred transform(tn, x)
    @test y isa NamedTuple{(:a,:b,:c)}
    @test inverse(tn, y) ≈ x
    index = 0
    ljacc = 0.0
    for (i, t) in enumerate((t1, t2, t3))
        d = dimension(t)
        xpart = t isa ScalarTransform ? x[index + 1] : x[index .+ (1:d)]
        @test y[i] == transform(t, xpart)
        ypart, ljpart = transform_and_logjac(t, xpart)
        @test ypart == y[i]
        ljacc += ljpart
        index += d
    end
    y2, lj2 = transform_and_logjac(tn, x)
    @test y == y2
    @test lj2 ≈ ljacc
end

@testset "empty Tuple, NamedTuple aggregators" begin
    zt = as(())
    znt = as(NamedTuple())
    za = as(Array, asℝ₊, 0)
    @test dimension(zt) == dimension(znt) == 0
    @test @inferred(transform(zt, Float64[])) == ()
    @test_skip inverse(zt, ()) == []
    @test @inferred(transform_and_logjac(zt, Float64[])) == ((), 0.0)
    @test @inferred(transform(znt, Float64[])) == NamedTuple()
    @test @inferred(transform_and_logjac(znt, Float64[])) == (NamedTuple(), 0.0)
    @test_skip inverse(znt, ()) == []
    @test @inferred(transform(za, Float64[])) == Float64[]
    @test @inferred(transform_and_logjac(za, Float64[])) == (Float64[], 0.0)
    @test_skip inverse(za, []) == []
end

@testset "nested combinations" begin
    # for https://github.com/tpapp/TransformVariables.jl/issues/57
    for _ in 1:10
        N = rand(3:7)
        tt = as((a = as(Tuple(as(Vector, asℝ₊, 2) for _ in 1:N)),
                 b = as(Tuple(UnitVector(n) for n in 1:N))))
        x = randn(dimension(tt))
        y = tt(x)
        x′ = inverse(tt, y)
        @test x ≈ x′
    end
end

####
#### log density correctness checks
####

@testset "transform logdensity: correctness" begin
    # the density is p(σ) = σ⁻³
    # let z = log(σ), so σ = exp(z)
    # the transformed density is q(z) = -3z + z = -2z
    f(σ) = -3*log(σ)
    q(z) = -2*z
    for _ in 1:1000
        z = randn()
        qz = @inferred transform_logdensity(asℝ₊, f, z)
        @test q(z) ≈ qz
    end
end

@testset "transform logdensity: type inference" begin
    t = as((a = asℝ₋, b = as𝕀, c = as((d = UnitVector(7), e = CorrCholeskyFactor(3))),
            f = as(Array, 9)))
    z = zeros(dimension(t))
    f(θ) = θ.a + θ.b + sum(abs2, θ.c.d) + sum(abs2, θ.c.e)
    @test (@inferred f(t(z))) isa Float64
    @test (@inferred transform_logdensity(t, f, z)) isa Float64
end

####
#### custom transformations
####

@testset "custom transformation: triangle below diagonal in [0,1]²" begin
    tfun(y) = y[1], y[1]*y[2]   # triangle below diagonal in unit square
    t = CustomTransform(as(Array, as𝕀, 2), tfun, collect;)
    test_transformation(t, ((y1, y2),) -> 0 ≤ y2 ≤ y1 ≤ 1;
                        vec_y = collect, test_inverse = false)

    # test inference w/ manually specified chunk
    function f(x)
        CustomTransform(as(Array, as𝕀, 2), tfun, collect;
                            chunk = ForwardDiff.Chunk{2}())
        transform_and_logjac(t, x)
    end
    y, lj = @inferred f(zeros(2))
    @test y == (0.5, 0.25)
    @test lj ≈ -3.465735902799726
end

@testset "custom transformation: covariance matrix" begin
    "Transform to a `n×n` covariance matrix."
    to_covariance(n) = CustomTransform(
        # pre-transform to standard deviations and correlation Cholesky factor
        as((as(Array, asℝ₊, n), CorrCholeskyFactor(n))),
        # use these to construct a covariance matrix
        (((σ, Ω),) -> (Ω*Diagonal(σ) |> x -> Symmetric(x'*x))),
        # flatten to elements above the diagonal
        A -> A[axes(A, 1) .≤ axes(A, 2)'])
    C5 = to_covariance(5)
    test_transformation(C5, A -> all(diag(cholesky(A).U) .> 0);
                        vec_y = C5.flatten, test_inverse = false)
end

@testset "offset arrays" begin
    t = as((λ = asℝ₊, a = CorrCholeskyFactor(4),
            θ = as((as(Array, as𝕀, 2, 3), UnitVector(4)))))
    x = random_arg(t)
    xo = OffsetVector(x, axes(x, 1) .- 7)
    @test transform(t, x) == transform(t, xo)
    @test transform_and_logjac(t, x) == transform_and_logjac(t, xo)
end

@testset "random value" begin
    t1 = asℝ
    t2 = CorrCholeskyFactor(7)
    t3 = UnitVector(3)
    tn = as((a = t1, b = t2, c = t3))
    y = random_value(tn)
    @test y isa NamedTuple{(:a, :b, :c), <: Tuple{Float64, AbstractMatrix, Vector}}
    @test size(y.b) == (7, 7)
    @test size(y.c) == (3, )
end

@testset "random arg" begin
    t = as(Array, 5)
    for _ in 1:1000
        @test sum(abs2, random_arg(t; cauchy = false, scale = 1.0)) ≤ 100
    end
end

####
#### AD compatibility tests
####

@testset "AD tests" begin
    t = as((μ = asℝ, σ = asℝ₊, β = asℝ₋, α = as(Real, 0.0, 1.0),
            u = UnitVector(3), L = CorrCholeskyFactor(4),
            δ = as((asℝ₋, as𝕀))))
    function f(θ)
        @unpack μ, σ, β, α, δ = θ
        -(abs2(μ) + abs2(σ) + abs2(β) + α + δ[1] + δ[2])
    end
    P = TransformedLogDensity(t, f)
    x = zeros(dimension(t))
    v = logdensity(P, x)
    g = ForwardDiff.gradient(x -> logdensity(P, x), x)

    # ForwardDiff
    P1 = ADgradient(:ForwardDiff, P)
    @test v == logdensity(P1, x)
    v1, g1 = @inferred logdensity_and_gradient(P1, x)
    @test v1 == v
    @test g1 ≈ g

    # Tracker
    P2 = ADgradient(:Tracker, P)
    v2, g2 = @inferred logdensity_and_gradient(P2, x)
    @test v2 == v
    @test g2 ≈ g

end

# if VERSION ≥ v"1.1"
#     if CIENV
#         @info "installing Zygote"
#         import Pkg
#         Pkg.API.add(Pkg.PackageSpec(; name = "Zygote"))
#     end

#     import Zygote

#     @testset "Zygote AD" begin
#         # Zygote
#         # NOTE @inferred removed as it currently fails
#         # NOTE tests simplified disabled as they currently fail
#         t = as((μ = asℝ, ))
#         function f(θ)
#             @unpack μ = θ
#             -(abs2(μ))
#         end
#         P = TransformedLogDensity(t, f)
#         x = zeros(dimension(t))
#         PF = ADgradient(:ForwardDiff, P)
#         PZ = ADgradient(:Zygote, P)
#         @test @inferred(logdensity(PZ, x)) == logdensity(P, x)
#         vZ, gZ = logdensity_and_gradient(PZ, x)
#         @test vZ == logdensity(P, x)
#         @test gZ ≈ last(logdensity_and_gradient(PF, x))
#     end
# end

@testset "inverse_and_logjac" begin
    # WIP, test separately until integrated
    for _ in 1:100
        x = randn()
        a = randn()
        t = as(Real, a, a + abs(randn()) + 0.1)
        y, lj = transform_and_logjac(t, x)
        x2, lj2 = TransformVariables.inverse_and_logjac(t, y)
        @test x2 ≈ x
        @test lj2 ≈ -lj
    end
end

@testset "inference of nested tuples" begin
    # An MWE adapted from a real-life problem
    ABOVE1 = as(Real, 1, ∞)   # transformation for μ ≥ 1

    trans_β̃s = as((asℝ, asℝ))     # a tuple of 2 elements, otherwise identity

    PARAMS_TRANSFORMATIONS =
        (EE = as((β̃s = trans_β̃s, μs = as((as𝕀, as𝕀)))),
         EN = as((w̃₂ = asℝ, β̃s = trans_β̃s, μs = as((as𝕀, ABOVE1)))),
         NE = as((w̃₁ = asℝ, β̃s = trans_β̃s, μs = as((ABOVE1, as𝕀)))),
         NN = as((w̃s = as((asℝ, asℝ)), β̃s = trans_β̃s, μs = as((ABOVE1, ABOVE1)))))

    function make_transformation(ls)
        as((hyper_parameters = as((μ = as(Array, 6),
                                   σ = as(Array, asℝ₊, 6),
                                   LΩ = CorrCholeskyFactor(6))),
            couple_parameters = as(map((t, l) -> as(Array, t, l),
                                       PARAMS_TRANSFORMATIONS, ls))))
    end
    t = make_transformation((EE = 1, EN = 2 , NE = 3, NN = 4,))
    x = zeros(dimension(t))
    @test_nowarn @inferred transform(t, x)
    @test_nowarn @inferred transform_and_logjac(t, x)
end

@testset "support abstract array inverses in ArrayTransform" begin
    t = as(Array, 2, 3)
    @test inverse(t, ones(SMatrix{2,3})) == ones(6)
end

####
#### broadcasting
####

@testset "broadcasting" begin
    @test as𝕀.([0, 0]) == [0.5, 0.5]

    t = UnitVector(3)
    d = dimension(t)
    x = [zeros(d), zeros(d)]
    @test t.(x) == map(t, x)
end

####
#### show
####

@testset "scalar show" begin
    @test string(asℝ) == "asℝ"
    @test string(asℝ₊) == "asℝ₊"
    @test string(asℝ₋) == "asℝ₋"
    @test string(as𝕀) == "as𝕀"
    @test string(as(Real, 0.0, 2.0)) == "as(Real, 0.0, 2.0)"
    @test string(as(Real, 1.0, ∞)) == "as(Real, 1.0, ∞)"
    @test string(as(Real, -∞, 1.0)) == "as(Real, -∞, 1.0)"
end


text:
using Documenter
using SeaPearl

makedocs(
    sitename = "SeaPearl",
    doctest = VERSION >= v"1.4",
    format = Documenter.HTML(),
    modules = [SeaPearl],
    pages = ["Home" => "index.md",
    "Constraint Programming solver" => [
        "Variables" => "CP/int_variable.md",
        "Trailer" => "CP/trailer.md",
    ],
    "Building Models" =>
      ["Basics" => "models/basics.md"],
    "Community" => "community.md"],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
##########################
# Geomagnetic Index Data #
##########################

let
    @test length(KpIndices(53941)) == 8
    @test isapprox(KpDailyIndex(53941), 12.3333, atol=1e-4)
    @test length(ApIndices(53941)) == 8
    @test isapprox(ApDailyIndex(53941), 6.0, atol=1e-4)
end

###################
# Solar Flux Data #
###################

let
    @test f107Observed(53941) == 75.9
    @test f107Adjusted(53941) == 78.3
    @test isapprox(f107ObservedAvg(53941), 72.2951, atol=1e-4)
    @test isapprox(f107AdjustedAvg(53941), 74.5593, atol=1e-4)

    data = f107Data(53941)
    @test data[1] == 75.9
    @test data[2] == 78.3
    @test isapprox(data[3], 72.2951, atol=1e-4)
    @test isapprox(data[4], 74.5593, atol=1e-4)
end


text:
"""
    Dependent Bounds

Specify lower and upper bounds that are not independent.
###
It can take as input a function, `f(pomdp, belief)`, that returns both lower and upper bounds.
It can take as input any object for which a `bounds` function is implemented

"""

init_bounds(bds, pomdp, sol, rng) = bds

# Used to initialize lower and upper bounds with a single function
bounds(f::Function, pomdp::P, s::S, max_depth::Int) where {S,P<:POMDP{S}} = f(pomdp, s)

# Used to initialize both the lower and upper bound with an object for which a `bounds` function is implemented
function bounds!(L::Vector{Float64}, U::Vector{Float64}, bd::B, pomdp::P, particles::Vector{S}, max_depth::Int) where {S,P<:POMDP{S},B}
    @inbounds for (i, s) in enumerate(particles)
        if isterminal(pomdp, s)
            L[i], U[i] = 0.0, 0.0
        else
            L[i], U[i] = bounds(bd, pomdp, s, max_depth)
        end
    end
    return L, U
end

"""
    IndependentBounds(lower, upper)

Specify lower and upper bounds that are independent of each other (the most common case).
A lower or upper bound can be a Number, a Function, `f(pomdp, belief)`, that returns a bound, an object for which a `bound` function
is implemented. Specifically, for FOValue, POValue, FORollout, SemiPORollout, and PORollout, a `bound` function is already implemented.
You can also implement a `bound!` function to initialize sibling beliefs simultaneously, if it could provide a further performace gain.
"""

mutable struct IndependentBounds{L, U}
    lower::L
    upper::U
end

function IndependentBounds(l, u)
    return IndependentBounds(l, u)
end

function init_bounds(bds::IndependentBounds, pomdp::POMDP, sol::DESPOTAlphaSolver, rng::AbstractRNG)
    return IndependentBounds(convert_estimator(bds.lower, sol, pomdp),
                             convert_estimator(bds.upper, sol, pomdp)
                             )
end

function bounds!(L::Vector{Float64}, U::Vector{Float64}, bds::IndependentBounds{LB,UB}, pomdp::P, particles::Vector{S}, max_depth::Int) where {LB,UB,S,P<:POMDP{S}}
    bound!(L, bds.lower, pomdp, particles, max_depth)
    bound!(U, bds.upper, pomdp, particles, max_depth)
    for (i, s) in enumerate(particles)
        if isterminal(pomdp, s)
            L[i], U[i] = 0.0, 0.0
        end
    end
    return L, U
end

# Used when the lower or upper bound is a fixed number
bound!(V::Vector{Float64}, n::Float64, pomdp::P, particles::Vector{S}, max_depth::Int) where {S,P<:POMDP{S}} = fill!(V, n)
bound!(V::Vector{Float64}, n::Int, pomdp::P, particles::Vector{S}, max_depth::Int) where {S,P<:POMDP{S}} = fill!(V, convert(Float64, n))

# Used when the lower or upper bound is an object for which a `bound` function is implemented
function bound!(V::Vector{Float64}, bd::B, pomdp::P, particles::Vector{S}, max_depth::Int) where {S,P<:POMDP{S},B}
    @inbounds for s in particles
        V[i] = bound(bd, pomdp, s, max_depth)
    end
    return V
end

function bound!(V::Vector{Float64}, bd::SolvedFORollout, pomdp::M, particles::Vector{S}, max_depth::Int) where {S,M<:POMDP{S}}
    sim = RolloutSimulator(bd.rng, max_depth)
    broadcast!((s)->simulate(sim, pomdp, bd.policy, s), V, particles)
end

function bound!(V::Vector{Float64}, bd::SolvedFOValue{P}, pomdp::M, particles::Vector{S}, max_depth::Int) where {P,S,M<:POMDP{S}}
    broadcast!((s)->value(bd.policy, s), V, particles)
end


text:
using Flux
using ConditionalDists
using GenerativeModels
import GenerativeModels: GAN
using ValueHistories
using MLDataPattern: RandomBatches
using Distributions
using DistributionsAD
using StatsBase
using Random

"""
	gan_constructor(;idim::Int=1, zdim::Int=1, activation="relu", hdim=128, nlayers::Int=3,
		init_seed=nothing, kwargs...)

Constructs a generative adversarial net.

# Arguments
	- `idim::Int`: input dimension.
	- `zdim::Int`: latent space dimension.
	- `activation::String="relu"`: activation function.
	- `hdim::Int=128`: size of hidden dimension.
	- `nlayers::Int=3`: number of generator/discriminator layers, must be >= 2.
	- `init_seed=nothing`: seed to initialize weights.
	- `last_linear=false`: is output of discriminator linear?
"""
function gan_constructor(;idim::Int=1, zdim::Int=1, activation = "relu", hdim=128, nlayers::Int=2,
	init_seed=nothing, last_linear=false, kwargs...)
	(nlayers < 2) ? error("Less than 3 layers are not supported") : nothing

	# if seed is given, set it
	(init_seed != nothing) ? Random.seed!(init_seed) : nothing

	# construct the model
	# generator
	generator_map = build_mlp(zdim, hdim, idim, nlayers, activation=activation)
	generator = ConditionalMvNormal(generator_map)

	# discriminator
	lastlayer = last_linear ? "linear" : "σ"
	discriminator_map = build_mlp(idim, hdim, 1, nlayers, activation=activation, lastlayer=lastlayer)
	discriminator = ConditionalMvNormal(discriminator_map)

	# reset seed
	(init_seed != nothing) ? Random.seed!() : nothing

	# constructor form GenerativeModels.jl
	model = GAN(zdim, generator, discriminator)
end

"""
	conv_gan_constructor(idim=(2,2,1), zdim::Int=1, activation="relu", hdim=1024, kernelsizes=(1,1),
		channels=(1,1), scalings=(1,1), init_seed=nothing, batchnorm=false, kwargs...)

Constructs a convolutional GAN.

# Arguments
	- `idim::Int`: input dimension.
	- `zdim::Int`: latent space dimension.
	- `activation::String="relu"`: activation function.
	- `hdim::Int=1024`: size of hidden dimension.
	- `kernelsizes=(1,1)`: kernelsizes in consequent layers.
	- `channels=(1,1)`: number of channels.
	- `scalings=(1,1)`: scalings in subsequent layers.
	- `init_seed=nothing`: seed to initialize weights.
	- `batchnorm=false`: use batchnorm (discouraged).
"""
function conv_gan_constructor(;idim=(2,2,1), zdim::Int=1, activation="relu", hdim=1024,
	kernelsizes=(1,1), channels=(1,1), scalings=(1,1),
	init_seed=nothing,	batchnorm=false, kwargs...)
	# if seed is given, set it
	(init_seed != nothing) ? Random.seed!(init_seed) : nothing

	# construct the model
	# generator
	generator_map = conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels),
		reverse(scalings); activation=activation, batchnorm=batchnorm)
	generator = ConditionalMvNormal(generator_map)

	# decoder - we will optimize only a shared scalar variance for all dimensions
	# also, the decoder output will be vectorized so the usual logpdfs vcan be used
	vecdim = reduce(*,idim[1:3]) # size of vectorized data
	discriminator_map = Chain(conv_encoder(idim, 1, kernelsizes, channels, scalings,
			activation=activation, batchnorm=batchnorm)..., x->σ.(x))
	discriminator = ConditionalMvNormal(discriminator_map)

	# reset seed
	(init_seed != nothing) ? Random.seed!() : nothing

	# get the vanilla VAE
	model = GAN(zdim, generator, discriminator)
end

"""
	StatsBase.fit!(model::GenerativeModels.GAN, data::Tuple, gloss:Function, dloss::Function;
		max_iter=10000, max_train_time=82800, lr=0.001, batchsize=64, patience=30, check_interval::Int=10,
		weight_clip=nothing, discriminator_advantage::Int=1, stop_threshold=0.01, usegpu=false,
		kwargs...)
"""
function StatsBase.fit!(model::GenerativeModels.GAN, data::Tuple, gloss::Function, dloss::Function;
	max_iters=10000, max_train_time=82800, lr=0.001, batchsize=64, patience=30, check_interval::Int=10,
	weight_clip=nothing, discriminator_advantage::Int=1, stop_threshold=0.01, usegpu=false,
	kwargs...)
	history = MVHistory()
	dopt = RMSProp(lr)
	gopt = RMSProp(lr)

	tr_model = deepcopy(model)
	dps = Flux.params(tr_model.discriminator)
	gps = Flux.params(tr_model.generator)
	_patience = patience

	tr_x = data[1][1]
	if ndims(tr_x) == 2
		val_x = data[2][1][:,data[2][2] .== 0]
	elseif ndims(tr_x) == 4
		val_x = data[2][1][:,:,:,data[2][2] .== 0]
	else
		error("not implemented for other than 2D and 4D data")
	end
	val_N = size(val_x,2)

	# on large datasets, batching loss is faster
	best_val_dloss = Inf
	best_val_gloss = Inf
	i = 1
	start_time = time() # end the training loop after 23hrs
	for xbatch in RandomBatches(tr_x, batchsize)
		xbatch = usegpu ? gpu(Array(xbatch)) : xbatch
		# disc loss
		batch_dloss = 0f0
		# some recommend training the discriminator mupltiple times
		for n in 1:discriminator_advantage
			gs = gradient(() -> begin
				batch_dloss = dloss(tr_model,xbatch)
			end, dps)
		 	Flux.update!(dopt, dps, gs)

		 	# clip weights of discriminator - needed for wgan
		 	(weight_clip != nothing) ? clip_weights!(dps, weight_clip) : nothing
		end

	 	# gen loss
		batch_gloss = 0f0
		gs = gradient(() -> begin
			batch_gloss = gloss(tr_model,xbatch)
		end, gps)
		Flux.update!(gopt, gps, gs)

		push!(history, :training_dloss, i, batch_dloss)
		push!(history, :training_gloss, i, batch_gloss)

		if mod(i, check_interval) == 0

			# validation/early stopping
			val_dloss = (val_N > 5000) ? dloss(tr_model, val_x, 256) : dloss(tr_model, val_x)
			val_gloss = (val_N > 5000) ? gloss(tr_model, val_x, 256) : gloss(tr_model, val_x)

			@info "$i - dloss: $(batch_dloss) | $(val_dloss) (validation)"
			@info "$i - gloss: $(batch_gloss) | $(val_gloss) (validation)"

			if isnan(val_gloss) || isnan(batch_dloss) || isnan(val_dloss) || isnan(batch_gloss)
				error("Encountered invalid values in loss function.")
			end

			push!(history, :validation_dloss, i, val_dloss)
			push!(history, :validation_gloss, i, val_gloss)

			if val_dloss > stop_threshold
				best_val_dloss = val_dloss
				_patience = patience

				# this should save the model at least once
				# when the validation loss is decreasing
				model = deepcopy(tr_model)
			else # else stop if the model has not improved for `patience` iterations
				_patience -= 1
				if _patience == 0
					@info "Stopped training after $(i) iterations."
					break
				end
			end
		end
		if (time() - start_time > max_train_time) | (i > max_iters) # stop early if time is running out
			model = deepcopy(tr_model)
			@info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours."
			break
		end
		i += 1
	end
	# again, this is not optimal, the model should be passed by reference and only the reference should be edited
	(history=history, iterations=i, model=model, npars=sum(map(p->length(p), Flux.params(model))))
end

"""
	dloss(model::GenerativeModels.GAN,x[,batchsize])

Classical discriminator loss of the GAN model.
"""
dloss(m::GenerativeModels.GAN,x) =
	dloss(m.discriminator.mapping,m.generator.mapping,x,rand(m.prior,size(x,ndims(x))))
dloss(m::GenerativeModels.GAN,x,batchsize::Int) =
	mean(map(y->dloss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))

"""
	gloss(model::GenerativeModels.GAN,x[,batchsize])

Classical generator loss of the GAN model.
"""
gloss(m::GenerativeModels.GAN,x) =
	gloss(m.discriminator.mapping,m.generator.mapping,rand(m.prior,size(x,ndims(x))))
gloss(m::GenerativeModels.GAN,x,batchsize::Int) =
	mean(map(y->gloss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))

"""
	fmloss(model,x[,batchsize])

Feature-matching loss matches the output of the penultimate layer of the discriminator on
real and fake data.
"""
function fmloss(m::GenerativeModels.GAN,x)
	z = rand(m.prior, size(x,ndims(x)))
	h = m.discriminator.mapping[1:end-1]
	hx = h(x)
	hz = h(m.generator.mapping(z))
	Flux.mse(hx, hz)
end
fmloss(m::GenerativeModels.GAN,x,batchsize::Int) =
	mean(map(y->fmloss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))

"""
	generate(model::GenerativeModels.GAN, N::Int)

Generate novel samples.
"""
generate(m::GenerativeModels.GAN, N::Int) = m.generator.mapping(rand(m.prior, N))

"""
	discriminate(model::GenerativeModels.GAN, x)

Discriminate the input - lower score belongs to samples not coming from training distribution.
"""
discriminate(m::GenerativeModels.GAN, x) = m.discriminator.mapping(x)


text:
using Test
using LinearAlgebra
using DynamicPolynomials
using SetProg, SetProg.Sets

@testset "Sets" begin
    @testset "ConvexPolySet" begin
        @polyvar x y
        P = SetProg.SumOfSquares.SymMatrix(Float64[1, 2, 3], 2)
        Q = SetProg.SumOfSquares.SymMatrix(BigInt[2, 3, 4], 2)
        basis = SetProg.Sets.MonoBasis(monovec([x, y]))
        q = SetProg.Sets.ConvexPolySet(2, SetProg.SumOfSquares.GramMatrix(P, basis), Q)
        @test q isa SetProg.Sets.ConvexPolySet{BigFloat}
    end
    @testset "zero_eliminate" begin
        @polyvar x y z
        p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i + j),
                                           monovec([x, y, z]))
        set = Sets.ConvexPolySet(2, p, nothing)
        el = Sets.zero_eliminate(set, 1:2)
        @test el.p.Q == 6ones(1, 1)
        el = Sets.zero_eliminate(set, 3:3)
        @test el.p.Q == [2 3; 3 4]
        el = Sets.zero_eliminate(set, 2:2)
        @test el.p.Q == [2 4; 4 6]

        @testset "Householder" begin
            p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i + j),
                                            monovec([x, y]))
            set = SetProg.perspective_dual_polyset(2, p, SetProg.InteriorPoint(zeros(2)), z, [x, y])
            @test set.set.p == 2x^2 + 6x*y + 4y^2 - z^2
            @test set.set.h == zeros(2)
            @test set.set.x == [x, y]
            @test Sets.gauge1(set.set.set) == 2x^2 + 6x*y + 4y^2
            set2 = Sets.project(set, [2])
            @test set2.set.p == 4y^2 - z^2
            @test set2.set.h == zeros(1)
            @test set2.set.x == [y]
        end
    end
end


text:
using AbstractGPs
using MLJAbstractGPsGlue
using MLJBase
using MLJ
using ParameterHandling
using Test

@testset "MLJAbstractGPsGlue.jl" begin

    X, y = @load_boston

    model = MLJAbstractGP()
    regressor = machine(model, X, y)

    # Just ensure that this runs as intended.
    evaluate!(regressor, resampling=CV(), measure=logpdf_loss, verbosity=0)

    # Train a single model and make some predictions.
    train_rows, test_rows = partition(eachindex(y), 0.7, shuffle=true)
    fit!(regressor; rows=train_rows, verbosity=0)

    @test ==(
        mean(predict(regressor; rows=test_rows)),
        predict_mean(regressor; rows=test_rows),
    )
end


text:
"""
    StopWatch(start, interval, callback)

Initialize a stopwatch.

# Arguments
- `start::Float64`: initial time (in seconds)
- `interval::Float64` : interval to click (in seconds)
- `callback` : callback function after each click (interval seconds)
"""
mutable struct StopWatch
    start::Float64
    interval::Float64
    f::Function
    StopWatch(_interval, callback) = new(time(), _interval, callback)
end

"""
    check(stopwatch, parameter...)

Check stopwatch. If it clicks, call the callback function with the unpacked parameter
"""
function check(watch::StopWatch, parameter...)
    now = time()
    if now - watch.start > watch.interval
        watch.f(parameter...)
        watch.start = now
    end
end


text:
"""
Bayesian 1dim Poisson Hidden Markov Model
"""
module PoissonHMM
using StatsFuns.logsumexp, SpecialFunctions.digamma
using Distributions

export Gam, GHMM, Poi, HMM
export sample_HMM, sample_data, winner_takes_all
export learn_VI

####################
## Types
struct Gam
    # Parameters of Gamma distribution
    # 1dim
    a::Float64
    b::Float64
end

struct BHMM
    # Parameters of Bayesian Bernoulli Mixture Model
    K::Int
    alpha_phi::Vector{Float64}
    alpha_A::Matrix{Float64}
    cmp::Vector{Gam}
end

struct Poi
    # Parameters of Poisson Distribution
    # 1 dim
    lambda::Float64
end

struct HMM
    # Parameters of Bernoulli Mixture Model
    K::Int
    phi::Vector{Float64}
    A::Matrix{Float64}
    cmp::Vector{Poi}
end

####################
## Common functions
"""
Sample an HMM from prior
"""
function sample_HMM(bhmm::BHMM)
    cmp = Vector{Poi}()
    for c in bhmm.cmp
        lambda = rand(Gamma(c.a, 1.0/c.b))
        push!(cmp, Poi(lambda))
    end
    phi = rand(Dirichlet(bhmm.alpha_phi))
    A = zeros(size(bhmm.alpha_A))
    for k in 1 : bhmm.K
        A[:,k] = rand(Dirichlet(bhmm.alpha_A[:,k]))
    end
    return HMM(bhmm.K, phi, A, cmp)
end

"""
Sample data from a specific Poisson HMM
"""
function sample_data(hmm::HMM, N::Int)
    X = zeros(N)
    Z = zeros(hmm.K, N)

    # sample (n=1)
    Z[:,1] = categorical_sample(hmm.phi)
    k = indmax(Z[:, 1])
    X[1] = rand(Poisson(hmm.cmp[k].lambda))

    # sample (n>1)
    for n in 2 : N
        Z[:,n] = categorical_sample(hmm.A[:,k])
        k = indmax(Z[:, n])
        X[n] = rand(Poisson(hmm.cmp[k].lambda))
    end
    return X, Z
end

categorical_sample(p::Vector{Float64}) = categorical_sample(p, 1)[:,1]
function categorical_sample(p::Vector{Float64}, N::Int)
    K = length(p)
    S = zeros(K, N)
    S_tmp = rand(Categorical(p), N)
    for k in 1 : K
        S[k,find(S_tmp.==k)] = 1
    end
    return S
end

function init_Z(X::Vector{Float64}, bhmm::BHMM)
    N = size(X, 1)
    K = bhmm.K

    Z = rand(Dirichlet(ones(K)/K), N)
    ZZ = [zeros(K,K) for _ in 1 : N - 1]
    for n in 1 : N - 1
        ZZ[n] = Z[:,n+1] * Z[:,n]'
    end

    return Z, ZZ
end

"""
Not implemented yet.
"""
function calc_ELBO(X::Matrix{Float64}, pri::BHMM, pos::BHMM)
end

function add_stats(bhmm::BHMM, X::Vector{Float64},
                   Z::Matrix{Float64}, ZZ::Vector{Matrix{Float64}})
    K = bhmm.K
    sum_Z = sum(Z, 2)
    alpha_phi = [bhmm.alpha_phi[k] + Z[k,1] for k in 1 : K]
    alpha_A = bhmm.alpha_A + sum(ZZ)
    cmp = Vector{Gam}()

    ZX = Z*X # (KxN) x (Nx1) = Kx1
    for k in 1 : K
        a = bhmm.cmp[k].a + ZX[k]
        b = bhmm.cmp[k].b + sum_Z[k]
        push!(cmp, Gam(a, b))
    end
    return BHMM(K, alpha_phi, alpha_A, cmp)
end

remove_stats(bhmm::BHMM, X::Vector{Float64}, Z::Matrix{Float64}) = add_stats(bhmm, X, -Z)

####################
## used for Variational Inference
function update_Z(bhmm::BHMM, X::Vector{Float64}, Z::Matrix{Float64})
    N = size(X, 1)
    K = bhmm.K
    ln_expt_Z = zeros(K, N)

    ln_lkh = zeros(K, N)
    for k in 1 : K
        ln_lambda = digamma.(bhmm.cmp[k].a) - log.(bhmm.cmp[k].b)
        lambda = bhmm.cmp[k].a / bhmm.cmp[k].b
        for n in 1 : N
            ln_lkh[k,n] = X[n]'*(ln_lambda) - lambda
        end
    end

    expt_ln_A = zeros(size(bhmm.alpha_A))
    for k in 1 : K
        expt_ln_A[:,k] = digamma.(bhmm.alpha_A[:,k]) - digamma.(sum(bhmm.alpha_A[:,k]))
    end

    # copy
    ln_expt_Z = log.(Z)

    # n = 1
    ln_expt_Z[:,1] = (digamma.(bhmm.alpha_phi) - digamma.(sum(bhmm.alpha_phi))
                      + expt_ln_A' * exp.(ln_expt_Z[:,2])
                      + ln_lkh[:,1]
                      )
    ln_expt_Z[:,1] = ln_expt_Z[:,1] - logsumexp(ln_expt_Z[:,1])

    # 2 <= n <= N - 1
    for n in 2 : N - 1
        ln_expt_Z[:,n] =( expt_ln_A * exp.(ln_expt_Z[:,n-1])
                          + expt_ln_A' * exp.(ln_expt_Z[:,n+1])
                          + ln_lkh[:,n]
                          )
        ln_expt_Z[:,n] = ln_expt_Z[:,n] - logsumexp(ln_expt_Z[:,n])
    end

    # n = N
    ln_expt_Z[:,N] =( expt_ln_A * exp.(ln_expt_Z[:,N-1])
                      + ln_lkh[:,N]
                      )
    ln_expt_Z[:,N] = ln_expt_Z[:,N] - logsumexp(ln_expt_Z[:,N])

    # calc output
    Z_ret = exp.(ln_expt_Z)
    ZZ_ret = [zeros(K,K) for _ in 1 : N - 1]
    for n in 1 : N - 1
        ZZ_ret[n] = Z_ret[:,n+1] * Z_ret[:,n]'
    end
    return Z_ret, ZZ_ret
end

"""
Pick single states having a max probability.
"""
function winner_takes_all(Z::Matrix{Float64})
    Z_ret = zeros(size(Z))
    for n in 1 : size(Z_ret, 2)
        idx = indmax(Z[:,n])
        Z_ret[idx,n] = 1
    end
    return Z_ret
end

function logmatprod(ln_A::Array{Float64}, ln_B::Array{Float64})
    I = size(ln_A, 1)
    J = size(ln_B, 2)
    ln_C = zeros(I, J)
    for i in 1 : I
        for j in 1 : J
            ln_C[i, j] = logsumexp(ln_A[i, :] + ln_B[:, j])
        end
    end
    return ln_C
end

function update_Z_fb(bhmm::BHMM, X::Vector{Float64})
    K = bhmm.K
    N = length(X)

    # calc likelihood
    ln_lik = zeros(K, N)
    for k in 1 : K
        ln_lambda = digamma.(bhmm.cmp[k].a) - log.(bhmm.cmp[k].b)
        lambda = bhmm.cmp[k].a / bhmm.cmp[k].b
        for n in 1 : N
            ln_lik[k,n] =X[n]'*(ln_lambda) - lambda
        end
    end
    expt_ln_phi = digamma.(bhmm.alpha_phi) - digamma.(sum(bhmm.alpha_phi))
    expt_ln_A = zeros(K,K)
    for k in 1 : K
        expt_ln_A[:,k] = digamma.(bhmm.alpha_A[:,k]) - digamma.(sum(bhmm.alpha_A[:,k]))
    end

    Z, ZZ = fb_alg(ln_lik, expt_ln_phi, expt_ln_A)

    # different notation
    ZZ_ret = [ZZ[:,:,n] for n in 1:size(ZZ, 3)]

    return Z, ZZ_ret
end

function fb_alg(ln_lik::Matrix{Float64}, ln_phi::Vector{Float64}, ln_A::Matrix{Float64})
    K, T = size(ln_lik)
    ln_Z = zeros(K, T)
    ln_ZZ = zeros(K, K, T)
    ln_alpha = zeros(K, T)
    ln_beta = zeros(K, T)
    ln_st = zeros(T)

    for t in 1 : T
        if t == 1
            ln_alpha[:, 1] = ln_phi + ln_lik[:, 1]
        else
            ln_alpha[:, t] = logmatprod(ln_A, ln_alpha[:, t-1]) + ln_lik[:, t]
        end
        ln_st[t] = logsumexp(ln_alpha[:, t])
        ln_alpha[:,t] = ln_alpha[:,t] - ln_st[t]
    end

    for t in T-1 : -1 : 1
        ln_beta[:, t] = logmatprod(ln_A', ln_beta[:, t+1] + ln_lik[:,t+1])
        ln_beta[:, t] = ln_beta[:, t] - ln_st[t+1]
    end

    ln_Z = ln_alpha + ln_beta
    for t in 1 : T
        if t < T
            ln_ZZ[:,:,t] = (repmat(ln_alpha[:, t]', K, 1) + ln_A
            + repmat(ln_lik[:, t+1] + ln_beta[:,t+1], 1, K))
            ln_ZZ[:,:,t] = ln_ZZ[:,:,t] - ln_st[t+1]
        end
    end
    return exp.(ln_Z), exp.(ln_ZZ)
end

"""
Compute approximate posterior distributions via variational inference.
"""
function learn_VI(X::Vector{Float64}, prior_bhmm::BHMM, max_iter::Int)
    # initialisation
    expt_Z, expt_ZZ = init_Z(X, prior_bhmm)
    bhmm = add_stats(prior_bhmm, X, expt_Z, expt_ZZ)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # E-step
        #expt_Z, expt_ZZ = update_Z(bhmm, X, expt_Z)
        expt_Z, expt_ZZ = update_Z_fb(bhmm, X)
        # M-step
        bhmm = add_stats(prior_bhmm, X, expt_Z, expt_ZZ)
    end

    return expt_Z, bhmm
end
end


text:
const SEED    = 5357
const MODELS  = ["kendall", "mmek", "autoreg", "dimer-decay", "yeast"]

const ParamSet = Tuple{Float64,Int,Int,Int,Int}

parameters = Dict{String,ParamSet}()

##### Model 1 #####

# simulation parameters
t_final  = 4.0
n_saves  = 1000
n_trial  = 10

# benchmark parameters
n_sample = 1000
t_limit  = 30 * 60

parameters[MODELS[1]] = (t_final, n_saves, n_trial, n_sample, t_limit)

##### Model 2 #####

# simulation parameters
t_final  = 50.0
n_saves  = 1000
n_trial  = 10

# benchmark parameters
n_sample = 1000
t_limit  = 30 * 60

parameters[MODELS[2]] = (t_final, n_saves, n_trial, n_sample, t_limit)

##### Model 3 #####

# simulation parameters
t_final  = 500.0
n_saves  = 1000
n_trial  = 10

# benchmark parameters
n_sample = 1000
t_limit  = 30 * 60

parameters[MODELS[3]] = (t_final, n_saves, n_trial, n_sample, t_limit)

##### Model 4 #####

# simulation parameters
t_final  = 100.0
n_saves  = 1000
n_trial  = 10

# benchmark parameters
n_sample = 1000
t_limit  = 30 * 60

parameters[MODELS[4]] = (t_final, n_saves, n_trial, n_sample, t_limit)

##### Model 5 #####

# simulation parameters
t_final  = 100.0
n_saves  = 1000
n_trial  = 10

# benchmark parameters
n_sample = 1000
t_limit  = 30 * 60

parameters[MODELS[5]] = (t_final, n_saves, n_trial, n_sample, t_limit)


text:
"""
TestControlLogicCore.jl defines the TestControlLogicCore type and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

using Boid: Boid
using Boid: AbstractControlLogicCore, InputChannel

# --- Type definitions

mutable struct TestControlLogicCore <: AbstractControlLogicCore
    #=
      Control state
    =#
    count::Int

    # Default constructor
    function TestControlLogicCore()
        count = 0
        new(count)
    end
end

@enum TestControlSignal begin
    INCREMENT
end

TEST_CONTROL_SIGNALS = instances(TestControlSignal)

# --- Method definitions

Boid.wait_for_input_ready(logic_core::TestControlLogicCore,
                          input_channels::Vector{InputChannel}) = nothing

function Boid.decode_control_signal(::Type{TestControlLogicCore},
                                    bytes::Vector{UInt8})

    signal = reinterpret(TestControlSignal, bytes)[1]

    if !(signal in TEST_CONTROL_SIGNALS)
        throw(ArgumentError("Unknown signal $(Int(signal))"))
    end

    return signal
end

function Boid.process_control_signal!(logic_core::TestControlLogicCore, signal)
    if signal == INCREMENT
        logic_core.count += 1
    else
        throw(ArgumentError("Unknown signal $(Int(signal))"))
    end

    return Boid.SUCCESS
end


text:
# Benchmark on map

using NumericFuns
using NumericExtensions

macro bench_map1(name, rp, f1, f2)
    quote
        println("Benchmark on $($name) ...")

        ($f1)(a)
        ($f2)(a)

        t1 = @elapsed for i in 1 : ($rp)
            ($f1)(a)
        end

        t2 = @elapsed for i in 1 : ($rp)
            ($f2)(a)
        end

        add_row!(perftable, $name, [t1, t2, t1/t2])
    end
end

macro bench_map2(name, rp, f1, f2)
    quote
        println("Benchmark on $($name) ...")

        ($f1)(a, b)
        ($f2)(a, b)

        t1 = @elapsed for i in 1 : ($rp)
            ($f1)(a, b)
        end

        t2 = @elapsed for i in 1 : ($rp)
            ($f2)(a, b)
        end

        add_row!(perftable, $name, [t1, t2, t1/t2])
    end
end

macro bench_map3(name, rp, f1, f2)
    quote
        println("Benchmark on $($name) ...")

        ($f1)($a, $b, $c)
        ($f2)($a, $b, $c)

        t1 = @elapsed for i in 1 : ($rp)
            ($f1)($a, $b, $c)
        end

        t2 = @elapsed for i in 1 : ($rp)
            ($f2)($a, $b, $c)
        end

        add_row!(perftable, $name, [t1, t2, t1/t2])
    end
end

# data

a = rand(1000, 1000)
b = rand(1000, 1000) .+ 0.5
c = rand(1000, 1000)

const perftable = BenchmarkTable("Comparison of element-wise map", ["Julia-expr", "Functor-map", "gain"])

add_map(a::Array, b::Array) = map(Add(), a, b)
@bench_map2("add", 10, .+, add_map)

mul_map(a::Array, b::Array) = map(Multiply(), a, b)
@bench_map2("multiply", 10, .*, mul_map)

my_inv(a::Array) = 1.0 ./ a
inv_map(a::Array) = map(Divide(), 1.0, a)
@bench_map1("inv", 10, my_inv, inv_map)

pow_map(a::Array, b::Array) = map(Pow(), a, b)
@bench_map2("pow", 10, .^, pow_map)

max_map(a::Array, b::Array) = map(MaxFun(), a, b)
@bench_map2("max", 10, max, max_map)

min_map(a::Array, b::Array) = map(MinFun(), a, b)
@bench_map2("min", 10, min, min_map)

abs_map(a::Array) = map(AbsFun(), a)
@bench_map1("abs", 10, abs, abs_map)

abs2_map(a::Array) = map(Abs2Fun(), a)
@bench_map1("abs2", 10, abs2, abs2_map)

sqrt_map(a::Array) = map(SqrtFun(), a)
@bench_map1("sqrt", 10, sqrt, sqrt_map)

exp_map(a::Array) = map(ExpFun(), a)
@bench_map1("exp", 10, exp, exp_map)

log_map(a::Array) = map(LogFun(), a)
@bench_map1("log", 10, log, log_map)

ju_absdiff(a::Array, b::Array) = abs(a .- b)
@bench_map2("absdiff", 10, ju_absdiff, absdiff)

ju_sqrdiff(a::Array, b::Array) = abs2(a .- b)
@bench_map2("sqrdiff", 10, ju_sqrdiff, sqrdiff)

ju_fma(a::Array, b::Array, c::Array) = a .+ b .* c
@bench_map3("fma", 10, ju_fma, fma)

println()
println(perftable)


text:
####
#### Coverage summary, printed as "(percentage) covered".
####
#### Useful for CI environments that just want a summary (eg a Gitlab setup).
####

using Coverage
cd(joinpath(@__DIR__, "..", "..")) do
    covered_lines, total_lines = get_summary(process_folder())
    percentage = covered_lines / total_lines * 100
    println("($(percentage)%) covered")
end


text:
# this file implements a shared buffer structure.

"""
Mutable buffer, sharing content with token and ['SubString']@ref instances.

content is guaranteed to be immutable when returned in a SubString or Token reference.
implementation is currently NOT threadsafe.
"""
mutable struct SharedIO <: IO
    buffer :: String # PRIVATE!! memory with token/substring text data.
    first :: UInt32 # read position offset of first not consumed byte
    free :: UInt32 # write position offset: offset of first unused byte
    shared :: UInt32 # last index in buffer shared with other objects (0 is valid)
    flags :: UInt32 # bit flags which define behavior
end

"""
Mutable string implementation.

"""
struct MutableString <: AbstractString
    io:: SharedIO
end

# verhaltensvarianten: über spezielle functions, automatische Auswahl
# über smart variante die flags auswertet ==> overhead durch if...
# dann kann anwendung entscheiden ob sie explizit steuert oder nicht
#
#oder doch typvarianten?

# tokenvector: load store in sharedIO. wann shared hochsetzen?
# wenn vector getrennt, dann share beim speichern. bedeutet shared==free

# Base._string_n
# allocates and returns a String witn n code units

#= alloc und Schreiben in String

# allocate Vector{UInt8}s for IOBuffer storage that can efficiently become Strings
StringVector(n::Integer) = unsafe_wrap(Vector{UInt8}, _string_n(n))
# das hat dann aber

unsafe_convert(::Type{Ptr{UInt8}}, s::String) = convert(Ptr{UInt8}, pointer_from_objref(s)+sizeof(Int))

aus strings.jl 253ff
function getindex(s::String, r::UnitRange{Int})
    isempty(r) && return ""
    i, j = first(r), last(r)
    @boundscheck begin
        checkbounds(s, r)
        @inbounds isvalid(s, i) || string_index_err(s, i)
        @inbounds isvalid(s, j) || string_index_err(s, j)
    end
    j = nextind(s, j) - 1
    n = j - i + 1
    ss = _string_n(n) # allocate garbage collected string
    p = pointer(ss) # ist wohl nötig, aber kann ja variable sein!
    for k = 1:n
        unsafe_store!(p, codeunit(s, i + k - 1), k) # schreibt in ss[k]
        # unsafe_store!(p, someUInt8, k) #
    end
    return ss # und fertich!
end

# out == IObuffer transferiert String ptr
String(take!(out))

so liest julia selbst aus ...

@inline function codeunit(s::String, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(pointer(s, i))
end

doc dazu
GC.@preserve x1 x2 ... xn expr
Temporarily protect the given objects from being garbage collected,
even if they would otherwise be unreferenced.

heißt für mich: da wir string im SharedIO gespeichert haben,
kann nix passieren. -> in Forum fragen!!

=#


text:
# # Making an Image of a Black Hole

# `Comrade` has been designed to work with the EHT and ngEHT.
# In this tutorial we will show how to reproduce some of the results
# from [EHTC VI 2019](https://iopscience.iop.org/article/10.3847/2041-8213/ab1141).

# In EHTC VI, they considered fitting simple geometric models to the data
# to estimate the image size, shape, brightness profile etc of the black hole.
# In this page we will construct a similar model and fit it to the data in under
# 50 lines of code (sans comments). To start we load some packages we will need

using Comrade
using Plots
using Distributions
using Pathfinder
using AdvancedHMC

# Pathfinder and AdvancedHMC will be required for sampling the posterior.
# Distributions is a general purpose probability distribution package in Julia.
# The next step is to load the data. For this we will use the publically
# available M 87 data which can be downloaded
# from [cyverse](https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/EHTC_FirstM87Results_Apr2019).
# For an introduction to data loading see [Loading Data into Comrade](@ref).

load_ehtim()
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "../assets/SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true)
# grab data products we want to fit: log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs; count="min")
dcphase = extract_cphase(obs, count="min")

# For this demo only consider closure products since these are invariant to station specific
# gain systematics. Given these data products we can then form our visibility likelihood:

lklhd = RadioLikelihood(dlcamp, dcphase)

# The `lklhd` constructs a measure using `MeasureTheory.jl`. To evaluate the likelihood
# we need to pass it a type that implements the `Comrade` model interface which is described
# in [Model Interface](@ref).

# To finish the construction of our posterior we need to specify an image model and a prior.
# For the image model we will be using a modified `MRing`, which is a
# infinitely thin delta ring with an azimuthal structure given by a Fourier expansion.
# To give the ring some width we will convolve the ring with a gaussian, and add an
# additional gaussian to the image to model any non-ring flux. The resulting model function
# is given by

function model(θ)
  (;radius, width, α, β, f, σG, τG, ξG, xG, yG) = θ
  ring = f*smoothed(stretched(MRing((α,), (β,)), radius, radius), width)
  g = (1-f)*shifted(rotated(stretched(Gaussian(), σG, σG*(1+τG)), ξG), xG, yG)
  return ring + g
end

# We now need to specify the priors for our model. The easiest way to do this is to
# specify a NamedTuple of distributions.

prior = (
          radius = Uniform(μas2rad(10.0), μas2rad(30.0)),
          width = Uniform(μas2rad(1.0), μas2rad(10.0)),
          α = Uniform(-0.5, 0.5),
          β = Uniform(-0.5, 0.5),
          f = Uniform(0.0, 1.0),
          σG = Uniform(μas2rad(1.0), μas2rad(40.0)),
          τG = Uniform(0.0, 0.75),
          ξG = Uniform(0.0, 1π),
          xG = Uniform(-μas2rad(80.0), μas2rad(80.0)),
          yG = Uniform(-μas2rad(80.0), μas2rad(80.0))
        )

# To form the posterior we now call

post = Posterior(lklhd, prior, model)

# This constructs a posterior density that can be evaluated by calling `logdensity`.
# For example

logdensity(post, (radius = μas2rad(20.0),
                  width = μas2rad(10.0),
                  α = 0.3,
                  β = 0.3,
                  f = 0.6,
                  σG = μas2rad(20.0),
                  τG = 0.1,
                  ξG = 0.5,
                  xG = 0.0,
                  yG = 0.0))

# We can now try to sample from our posterior `post` so that we can make probabilstic
# inferences about our data.

# Our strategy here will be to use Hamiltonian Monte Carlo. However, to lower burn-in time
# we will first use the excellent Pathfinder algorithm[^Zhang2021], which uses variational
# inference to find a starting point that is roughly drawn from the posterior. The benefit
# this is that we will greatly minimize burn-in time, and we get a variational approximation
# for free.

q, ϕ, _ = multipathfinder(post, 100)

# The q's are a uniformly weighted mixture of multivariate normal distributions, while
# the `ϕ`'s are our approximate posterior draws. We can then start an HMC run by calling

ndim = dimension(post)
chain, stats = sample(post, HMC(metric=DiagEuclideanMetric(ndim)), 2000; nadapts=1000, init_params=ϕ[1])

# That's it! We now have a converged chain. To finish it up we can then plot some diagnostics.

# First to plot the image we call

plot(model(chain[end]))

# That looks similar to the EHTC VI, and it took us no time at all!. To see how well the
# model is fitting the data we can plot the model and data products

plot(model(chain[end]), dlcamp)

# or even better the residuals

residual(model(chain[end]), dlcamp)

# [^Zhang2021]: Lu Zhang, Bob Carpenter, Andrew Gelman, Aki Vehtari (2021).
# Pathfinder: Parallel quasi-Newton variational inference.
# arXiv: [2108.03782](https://arxiv.org/abs/2108.03782) [stat.ML].
# [Code](https://github.com/LuZhangstat/Pathfinder))


text:
using Mft-3-1-2
using Test

@testset "Mft-3-1-2.jl" begin
    # Write your tests here.
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler

equations = CompressibleEulerEquations2D(1.4)

@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D)
  # set the freestream flow parameters
  rho_freestream = 1.0
  u_freestream = 0.3
  p_freestream = inv(equations.gamma)

  theta = pi / 90.0 # analogous with a two degree angle of attack
  si, co = sincos(theta)
  v1 = u_freestream * co
  v2 = u_freestream * si

  prim = SVector(rho_freestream, v1, v2, p_freestream)
  return prim2cons(prim, equations)
end

initial_condition = uniform_flow_state

boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
boundary_conditions = Dict( :Body    => boundary_condition_uniform_flow,
                            :Button1 => boundary_condition_slip_wall,
                            :Button2 => boundary_condition_slip_wall,
                            :Eye1    => boundary_condition_slip_wall,
                            :Eye2    => boundary_condition_slip_wall,
                            :Smile   => boundary_condition_slip_wall,
                            :Bowtie  => boundary_condition_slip_wall )

volume_flux = flux_ranocha
solver = DGSEM(polydeg=5, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# Get the unstructured quad mesh from a file (downloads the file if not available locally)
default_mesh_file = joinpath(@__DIR__, "abaqus_gingerbread_man.inp")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/0e9e990a04b5105d1d2e3096a6e41272/raw/0d924b1d7e7d3cc1070a6cc22fe1d501687aa6dd/abaqus_gingerbread_man.inp",
                                      default_mesh_file)
mesh_file = default_mesh_file

mesh = P4estMesh{2}(mesh_file)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 3.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorLöhner(semi, variable=density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=0,
                                      med_level=1, med_threshold=0.05,
                                      max_level=3, max_threshold=0.1)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback)

###############################################################################
# run the simulation
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7,
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using MPI

function MPISendrecvFloat(sendbuf::Float64, dest::Integer, comm::MPI.Comm)::Float64
    sbuf = [sendbuf]
    rbuf = [sendbuf]
    MPI.Sendrecv!(sbuf, dest, 0, rbuf, dest, 0, comm)
    return rbuf[1]
end

function MPISendBool(sendbuf::Bool, dest::Integer, comm::MPI.Comm)
    sbuf = [sendbuf]
    MPI.Send(sbuf, dest, 0, comm)
    return nothing
end

function MPIRecvBool(dest::Integer, comm::MPI.Comm)::Bool
    rbuf = [ false ]
    MPI.Recv!(rbuf, dest, 0, comm)
    return rbuf[1]
end

function MPIBcastBool(sendbuf::Bool, root::Integer, comm::MPI.Comm)::Bool
    sbuf = [ sendbuf ]
    MPI.Bcast!(sbuf, root, comm)
    return sbuf[1]
end


text:
push!(LOAD_PATH,"../src/")

# include package
info("Include all...")
try
  include("../src/JuliaOpenGL.jl")
catch e # do not exit this run!
  warn(e)
end
info("Include done.")

info("Create Docs...")
using Documenter, App

makedocs(
  build     = joinpath(@__DIR__, "../build/docs"),
  modules   = [App],
  clean     = true,
  doctest   = true, # :fix
  #linkcheck = true,
  strict    = false,
  checkdocs = :none,
  format    = :html, #:latex
  sitename  = "JuliaOpenGL",
  authors   = "Gilga",
  #analytics = "UA-89508993-1",
  #html_prettyurls = true,
  #html_canonical = "https://gilga.github.io/JuliaOpenGL/",
  pages = Any[ # Compat: `Any` for 0.4 compat
      "Home" => "index.md",
      "Manual" => Any[
          "manual/start.md",
          "manual/algorithm.md",
          "manual/optimization.md",
          "manual/references.md",
      ],
      "Source Files" => Any[
          "files/JuliaOpenGL.md",
          "files/build.md",
          "files/camera.md",
          "files/chunk.md",
          "files/compileAndLink.md",
          "files/cubeData.md",
          "files/frustum.md",
          "files/lib_math.md",
          "files/lib_opengl.md",
          "files/lib_time.md",
          "files/lib_window.md",
          "files/libs.md",
          "files/matrix.md",
          "files/shader.md",
          "files/test.md",
          "files/texture.md",
          "files/vector.md",
      ],
  ],
)

deploydocs(
  deps   = Deps.pip("mkdocs", "python-markdown-math"), #, "curl"
  repo = "https://github.com/Gilga/JuliaOpenGL.git",
  branch = "gh-pages",
  julia  = "0.6.2",
)

info("Docs done.")

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#=

  From https://reference.wolfram.com/language/ref/BinomialDistribution.html
  """
  A basketball player has a free-throw percentage of 0.75.

  Find the probability that the player hits 2 out of 3 free throws in a game.
  [free_throwIs2]
    Answer: 0.421875

  Find the probability that the player hits the last 2 of 5 free throws.
  [q1]
     Answer: 0.00878906

  Find the expected number of hits in a game with n free throws.
  [q2]
     Answer: 0.75 n
  """

  This model don't find the answer on q1 since the probability is too small...

  q2:  0.7504

  Distributions of variable free_throwIs2 (num:0)
  0.00000 =>   57552  (0.575520)
  1.00000 =>   42448  (0.424480)
  Distributions of variable q1 (num:0)
  0.00000 =>  100000  (1.000000)
  Distributions of variable q2 (num:0)
  0.80000 =>   27864  (0.278640)
  0.70000 =>   25049  (0.250490)
  0.90000 =>   18829  (0.188290)
  0.60000 =>   14636  (0.146360)
  0.50000 =>    5906  (0.059060)
  1.00000 =>    5803  (0.058030)
  0.40000 =>    1606  (0.016060)
  0.30000 =>     269  (0.002690)
  0.20000 =>      37  (0.000370)
  0.10000 =>       1  (0.000010)

  Cf ~/webppl/binomial_basketball.wppl

=#

using Turing
include("jl_utils.jl")

@model function binomial_basketball()

    p = 0.75
    n = 2

    # Find the probability that the player hits 2 out of 3 free throws in a game.
    free_throw ~ Binomial(3,p)

    free_throwIs2 ~ Dirac(free_throw == 2)

    # Find the probability that the player hits the last 2 of 5 free throws.
    free_throw2a ~ Binomial(3,p) # This should be 0
    free_throw2b ~ Binomial(2,p) # This should be 2

    # Find the expected number of hits in a game with n free throws.
    free_throw3 ~ Binomial(10,p)

    q1 ~ Dirac(free_throw2a == 0 & free_throw2b == 2)
    q2 ~ Dirac(free_throw3 / 10.0)

end

model = binomial_basketball()

num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 40_000)
# chs = sample(model, IS(), 50_000)
chs = sample(model, SMC(), 100_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chains)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

show_var_dist_pct(chs,:free_throwIs2)
show_var_dist_pct(chs,:q1)
show_var_dist_pct(chs,:q2)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Timescale for diffusion across one cell
min_Δxyz(grid) = min(grid.Δx, grid.Δy, grid.Δz)
min_Δxy(grid) = min(grid.Δx, grid.Δy)

min_Δx(grid) = grid.Δx
min_Δy(grid) = grid.Δy
min_Δz(grid) = grid.Δz

cell_diffusion_timescale(model) = cell_diffusion_timescale(model.closure, model.diffusivities, model.grid)

function cell_diffusion_timescale(closure::IsotropicDiffusivity{N, <:NamedTuple{()}},
                                  diffusivities, grid) where N
    return min_Δxyz(grid)^2 / closure.ν
end

function cell_diffusion_timescale(closure::IsotropicDiffusivity, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_κ = maximum(closure.κ)
    return min(Δ^2 / closure.ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::AnisotropicDiffusivity{NX, NY, NZ, <:NamedTuple{()}, <:NamedTuple{()}, <:NamedTuple{()}},
                                  diffusivities, grid) where {NX, NY, NZ}
    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)

    return min(Δx^2 / closure.νx,
               Δy^2 / closure.νy,
               Δz^2 / closure.νz)
end

function cell_diffusion_timescale(closure::AnisotropicDiffusivity, diffusivities, grid)

    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)

    max_κx = maximum(closure.κx)
    max_κy = maximum(closure.κy)
    max_κz = maximum(closure.κz)

    return min(Δx^2 / closure.νx,
               Δy^2 / closure.νy,
               Δz^2 / closure.νz,
               Δx^2 / max_κx,
               Δy^2 / max_κy,
               Δz^2 / max_κz)
end

function cell_diffusion_timescale(closure::AnisotropicBiharmonicDiffusivity{V, <:NamedTuple{()}, <:NamedTuple{()}},
                                  diffusivities, grid) where V
    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)

    return min(Δx^4 / closure.νx,
               Δy^4 / closure.νy,
               Δz^4 / closure.νz)
end

function cell_diffusion_timescale(closure::AnisotropicBiharmonicDiffusivity, diffusivities, grid)
    Δx = min_Δx(grid)
    Δy = min_Δy(grid)
    Δz = min_Δz(grid)

    max_κx = maximum(closure.κx)
    max_κy = maximum(closure.κy)
    max_κz = maximum(closure.κz)

    return min(Δx^4 / closure.νx,
               Δy^4 / closure.νy,
               Δz^4 / closure.νz,
               Δx^4 / max_κx,
               Δy^4 / max_κy,
               Δz^4 / max_κz)
end

function cell_diffusion_timescale(closure::SmagorinskyLilly{FT, P, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, P}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::BlasiusSmagorinsky{ML, FT, P, <:NamedTuple{()}},
                                  diffusivities, grid) where {ML, FT, P}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::AbstractSmagorinsky, diffusivities, grid)
    Δ = min_Δxyz(grid)
    min_Pr = minimum(closure.Pr)
    max_κ = maximum(closure.κ)
    max_νκ = maximum(diffusivities.νₑ.data.parent) * max(1, 1/min_Pr)
    return min(Δ^2 / max_νκ, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::RozemaAnisotropicMinimumDissipation{FT, <:NamedTuple{()}},
                                  diffusivities, grid) where FT
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::VerstappenAnisotropicMinimumDissipation{FT, PK, PN, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, PK, PN}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::AbstractAnisotropicMinimumDissipation, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    max_κ = max(Tuple(maximum(κₑ.data.parent) for κₑ in diffusivities.κₑ)...)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::AbstractLeith, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

cell_diffusion_timescale(closure::Tuple, diffusivities, grid) =
    min(Tuple(cell_diffusion_timescale(c, diffusivities, grid) for c in closure)...)


text:
#TYPES
"""
    AbstractModel

abstract type for models
"""
abstract type AbstractModel end

"""
    AbstractDependenceStructure

Types inheriting from abstract type `AbstractDependenceStructure`
"""
abstract type AbstractDependenceStructure end

"""
    FullIndependence <: AbstractDependenceStructure

Type acting as a Flag for Full independent Zig-Zag sampler
"""
struct FullIndependence <: AbstractDependenceStructure end

"""

Type acting as a Flag for partial independent Zig-Zag sampler
"""
struct PartialIndependence <: AbstractDependenceStructure end

"""
    SamplingScheme

Abstact type for Sampling scheme
"""
abstract type SamplingScheme end

"""
    SubSampling <: SamplingScheme

If you cannot sample from the inhomogeneous Poisson rate, sample it by
subsampling from Poisson rate with higher intensity (bound) and accept reject the event.
"""
struct SubSampling <: SamplingScheme end

"""
    Regular <: SamplingScheme

For linear sdes, where the actual imhogeneous Poisson rate can be sampled directly.
"""
struct Regular <: SamplingScheme end

"""
    System

CHANGE IN TUPLE WITH FREE NUMBER OF PARAMETERS
contains all the information needed for the ZigZag sampler
    ξ::Vector{Float64} := vector for the position of the coefficients
    θ::Vector{Float64} := vector containing the velocities (need to be changes in float64)
    ϕ::Vector{Fs} := Faber Schauder functions information (see Fs)
    τ::Vector{Float64} := vector containing waiting time (could be sorted)
    L::Int64 := Number of Levels
    T::Float64 := time length of the bridge
    b1::Vector{Float64} := free vector needed for linear growth sde
    b2::Vector{Float64} := free vector needed for linear growth sde
    tt::Vector{Float64} := free vector needed for linear growth sde
"""
struct System
    ξ::Vector{Float64} #
    θ::Vector{Float64}
    ϕ::Vector{Fs}
    τ::Vector{Float64}
    L::Int64
    T::Float64
    function System(L::Int64, T::Float64, ξ = fill(0.0, 2<<L - 1), θ = fill(1.0, 2<<L - 1), τ = fill(0.0, 2<<L - 1))
        new(ξ, θ, generate(L, T), τ, L, T)
    end
end

"""
    Skeleton
container for output zigzag sampler
"""
struct Skeleton
    ξ::Vector{Float64}
    t::Float64
end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


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

g = load(joinpath(testdir,"testdata","graph-50-500.jgz"), "graph-50-500")

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

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

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

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

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

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


text:
module AWSInfinidash
using AWS: AWS, @service

export I8h, i8h

# low level SDK
const i8h = try
    AWS.AWSServices.infinidash
catch err
    # defer any errors til the user time
    (args...; kwargs...) -> throw(err)
end

# High level SDK
function __init__()
    # TODO This should just be `const I8h = @service Infinidash` but that doesn't support
    # checking if AWS.jl has support for it. So we unwreap the macro by hand.
    aws_package_loc = Base.find_package("AWS")
    if aws_package_loc !== nothing
        service_path = joinpath(dirname(aws_package_loc), "services", "infinidash.jl")
        if isfile(service_path)
            include(service_path)
            return nothing  # No
        end
    end

    # Would have returned before now if it loaded successfully
    @warn "SDK for Infinidash not found. Is your version of AWS.jl up to date? Functionality may be limitted."
    @eval I8h = Module(:Infinidash)
end

# Very High level API
# TODO implement a more idiomatic julia interface over the SDK methods

end  # module


text:
# example of passing a pointer to a julia method to Cxx

using Cxx
__current_compiler__ = Cxx.new_clang_instance()

cxx"""
void evaluate(void (*f)() ){
    f();
}
"""

function evaluate(f)
    # void pointer on f
    f_void_ptr = cfunction(f, Void, ())

    # cast the void pointer into a function pointer
    icxx"""
    typedef void (*fptr)();
    fptr f_fun_ptr = reinterpret_cast<fptr>(reinterpret_cast<long>($f_void_ptr));
    """

    # finally we can do the desired call
    @cxx evaluate(f_fun_ptr);

end

f() = println("hi")
evaluate(f)

g() = println("ho")
evaluate(g)


text:
@with_kw struct fitDTOF

    #required inputs
    t::Union{AbstractRange{Float64}, AbstractVector{Float64}}
    DTOF::Array{Float64,1}
    IRF::Array{Float64,1}

    nmed::Float64 = 1.5
	ndet::Float64 = 1.51
	lambda::Float64 = 750.0

    #semi-infinite and slab
    ρ::Union{Float64, Missing} = missing

    #slab
    s::Union{Float64, Missing} = missing

    #parallelepiped
    rd::Union{Array{Float64,1}, Missing} = missing
    rs::Union{Array{Float64,1}, Missing} = missing
    L::Union{Array{Float64,1}, Missing} = missing
end

@with_kw struct fitg2

    #required inputs
    τ::Union{AbstractRange{Float64}, AbstractVector{Float64}}
    g2::Array{Float64,1}

    #semi-infinite
    ρ::Union{Float64, Missing}
end

@with_kw struct DTOF_fitparams
    model::Function = TPSF_DA_semiinf_refl          # model to be fitted against
	initparams::Array{Float64,1} = [0.2, 10.0]      # initial starting point of LsqFit
	lb::Array{Float64,1} = [0.001, 1.0]             # lowerbound of fit for [mua, musp]
	ub::Array{Float64,1} = [1.0, 20.0]              # upperbound of fit for [mua, musp]
	risefactor::Float64  = 0.9			            # percent of peak to fit on rising edge
	tailfactor::Float64	 = 0.1			            # percent of peak to fit on falling tail
end

@with_kw struct g2_fitparams
    model::Function = g2_DA_SI_R                    # model to be fitted against
	initparams::Array{Float64,1} = [2.e-8, 0.5]     # initial starting point of LsqFit
	lb::Array{Float64,1} = [1.e-12, 0.0]            # lowerbound of fit for [BFI, β]
	ub::Array{Float64,1} = [5.e-6, 1.0]             # upperbound of fit for [BFI, β]
end

struct fitresult
    #Fit inputs
    xdata::Array{Float64,1}     #independent variable
    ydata::Array{Float64,1}     #dependent variable
    yerrors::Array{Float64,1}   #uncertainties in dependent variable
	model::Function             #model to be fitted against

    #Fit results
    residuals::Array{Float64,1}  	 #weighted residuals
    param::Array{Float64,1}    		#best fit parameters for model
	marginerror::Array{Float64,1}        #product of standard error and the critical value of each parameter at a certain significance level (default is 5%) from t-distribution.
    perrors::Array{Float64,1}            #uncertainties in best fit parameters
    dof::Int64                  		#number of degrees of freedom in fit
    chisq::Float64             		 #chi^2 value of fit
    probability::Float64      		  #probability of chi^2 value
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
module VulkanCore

include("LibVulkan.jl")
using .LibVulkan

const api = LibVulkan
const vk = api
export vk

end # module