text:
include("pi.jl")
include("md.jl")

Operator = Union{Matrix, Diagonal, Symmetric, UniformScaling}

abstract type IsomorphicClassicalSystem <: ClassicalSystem end

kinetic(system::IsomorphicClassicalSystem, p::AbstractArray) = .5p' * inv(system.M) * p
potentialInternal(system::IsomorphicClassicalSystem, q::AbstractArray) = .5q' * system.A * q
potentialExternal(system::IsomorphicClassicalSystem, q::AbstractArray) = sum(system.quantum.potential, system.S * q) / system.n
potential(system::IsomorphicClassicalSystem, q::AbstractArray) = potentialInternal(system, q) + potentialExternal(system, q)

"""
    Primitive(quantum, n)

Creates an primitive isomorphic PIMD system from a quantum system, with beads number `n`.
"""
struct Primitive <: IsomorphicClassicalSystem
    n::Integer
    β::Real
    M::Operator
    A::Operator
    S::Operator
    quantum::OneDimensionalQuantumSystem
    Primitive(quantum::OneDimensionalQuantumSystem, n::Int) = begin
        M = Diagonal(quantum.m / n * I, n)
        A = stiffnessMatrix(quantum, n)
        S = I
        new(n, quantum.β, M, A, S, quantum)
    end
end

"""
    Staging(quantum, n)

Creates an staging isomorphic PIMD system from a quantum system, with beads number `n`.
"""
struct Staging <: IsomorphicClassicalSystem
    n::Integer
    β::Real
    M::Operator
    A::Operator
    S::Operator
    quantum::OneDimensionalQuantumSystem
    Staging(quantum::OneDimensionalQuantumSystem, n::Int) = begin
        M = Diagonal([quantum.m * i / (i - 1) for i in 1:n])
        M[1,1] = quantum.m
        K = stiffnessMatrix(quantum, n)
        S = zeros(n,n)
        S[:,1] .= 1.
        for i in 2:n
            S[i,i:n] = (i-1) ./ ((i-1):(n-1))
        end
        A = Diagonal(S' * K * S)
        new(n, quantum.β, M, A, S, quantum)
    end
end

"""
    Normal(quantum, n)

Creates an normal isomorphic PIMD system from a quantum system, with beads number `n`.
"""
struct Normal <: IsomorphicClassicalSystem
    n::Integer
    β::Real
    M::Diagonal
    A::Operator
    S::Operator
    quantum::OneDimensionalQuantumSystem
    Normal(quantum::OneDimensionalQuantumSystem, n::Int) = begin
        a = 0:n-1
        M = Diagonal([quantum.m * 4 * sin(i * π / n)^2 for i in a])
        M[1,1] = quantum.m
        K = stiffnessMatrix(quantum, n)
        S = zeros(n,n)
        S[:,1] .= sqrt(1 / n)
        for i in 1:(n-1)÷2
            c1, c2 = i+1, n+1-i
            @. S[:,c1] = sqrt(2 / n) * cos(2π * i * a / n)
            @. S[:,c2] = sqrt(2 / n) * sin(2π * i * a / n)
        end
        if n % 2 == 0
            @. S[:,n÷2+1] = sqrt(1 / n) * (-1) ^ a
        end
        A = Diagonal(S' * K * S)
        new(n, quantum.β, M, A, S, quantum)
    end
end

@doc raw"""
    estimatorPotential(phase, system)

Estimates the potential energy for an isomorphic classical system under a phase.

``\hat V(\boldsymbol x)=\varphi(\boldsymbol x)``
"""
const estimatorPotential = (phase::HamiltonianPhase, system::IsomorphicClassicalSystem) -> sum(system.quantum.potential, system.S * phase.q) / system.n
@doc raw"""
    estimatorKinetic(phase, system)

Estimates the kinetic energy (using the primitive estimator) for an isomorphic classical system under a phase.

``\hat K_{\rm p}(\boldsymbol x)=\frac n{2\beta}-\frac12\boldsymbol x^T\mathsf K\boldsymbol x``
"""
const estimatorKinetic = (phase::HamiltonianPhase, system::IsomorphicClassicalSystem) -> system.n / 2system.β - phase.q' * system.A * phase.q / 2

@doc raw"""
    estimatorKinetic(phase, system)

Estimates the kinetic energy (using the virial estimator) for an isomorphic classical system under a phase.

``\hat K_{\rm v}(\boldsymbol x)=\frac 1{2\beta}+\frac1{2n}(\bar x-\boldsymbol x)\partial{\varphi(x)}{\boldsymbol x}``
"""
const estimatorKineticVirial = (phase::HamiltonianPhase, system::IsomorphicClassicalSystem) -> begin
    x = system.S * phase.q
    centroid = sum(x) / length(x)
    1 / 2system.β + (centroid .- x)' * system.quantum.force.(x) / system.n / 2.
end

@doc raw"""
    estimatorHessian(phase, system)

Estimates the hessian for an isomorphic classical system under a phase.

``\hat H(\boldsymbol x)=\frac1n\sum_iV''(x_i)``
"""
const estimatorHessian = (phase::HamiltonianPhase, system::IsomorphicClassicalSystem) -> sum(system.quantum.hessian, system.S * phase.q) / system.n


text:
using DataFrames

using Dates

using CSV: File

function tabix_regions_from_file(
        bed_file_path::String,
        vcf_file_path::String,
        output_dir::String
    )

    start_time = now()

    println("Checking variants in this file: $bed_file_path\n")

    bed_io = open(bed_file_path)

    bed = DataFrame(File(bed_io; delim='\t', header=false))

    println("bed file: \n")

    println(bed)

    println("\n")

    close(bed_io)

    bed_split_file_path = Array(split(bed_file_path, "/"))

    bed_name = split(bed_split_file_path[size(bed_split_file_path)[1]], ".")[1]

    output_file_path = joinpath(output_dir, bed_name)

    print_and_run_cmd(pipeline(`tabix -R $bed_file_path $vcf_file_path`, stdout="$output_file_path.tsv"))

    println("\nSaved output at: $output_file_path.tsv")

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export tabix_regions_from_file

ext: jl


text:
module ObservableInterfaces

using QueryOperators, DataFrames

export Foo, TextFile, MapObservable, ConsoleSink, Subject

abstract type Observable end

function onnext end
function oncompleted end

function subscribe end

mutable struct Subject{T} <: Observable
    observers::Vector{Any}

    function Subject{T}()
        return new{T}(Vector{Any}())
    end
end

function subscribe(s::Subject, observer)
    push!(s.observers, observer)
end

function onnext(s::Subject, el)
    for o in s.observers
        onnext(o, el)
    end
end

function oncompleted(s::Subject)
    for o in s.observers
        oncompleted(o)
    end
end

mutable struct Foo <: Observable
    range::Int
    observers::Vector{Any}

    function Foo(range)
        return new(range, Vector{Any}(0))
    end
end

function subscribe(f::Foo, observer)
    push!(f.observers, observer)
end

function process(f::Foo)
    for i in 1:f.range
        for s in f.observers
            onnext(s,i)
        end
    end

    for s in f.observers
        oncompleted(s)
    end
end

mutable struct MapObservable <: Observable
    f
    observers::Vector{Any}
    function MapObservable(f, observable)
        new_instance = new(f, Vector{Any}(0))
        subscribe(observable, new_instance)
        return new_instance
    end
end

function subscribe(f::MapObservable, observer)
    push!(f.observers, observer)
end

function onnext(mo::MapObservable, el)
    for o in mo.observers
        onnext(o, mo.f(el))
    end
end

function oncompleted(mo::MapObservable)
    for o in mo.observers
        oncompleted(o)
    end
end

mutable struct TextFile
    handle

    function TextFile(filename::AbstractString, observable)
        h = open(filename, "w")
        new_instance = new(h)
        subscribe(observable, new_instance)
        return new_instance
    end
end

function onnext(tx::TextFile, el)
    println(tx.handle, el)
end

function oncompleted(tx::TextFile)
    close(tx.handle)
end

mutable struct ConsoleSink
    function ConsoleSink(observable)
        new_instance = new()
        subscribe(observable, new_instance)
        return new_instance
    end
end

function onnext(tx::ConsoleSink, el)
    println(el)
end

function oncompleted(tx::ConsoleSink)
end

include("query.jl")

end # module


text:
##
# Represent a binary hierarchical Domain, Space, and Fun
##

for (HDSF,DSF) in ((:HierarchicalDomain,:(ApproxFun.UnivariateDomain)),(:HierarchicalSpace,:(ApproxFun.UnivariateSpace)))
    @eval begin
        type $HDSF{S,T,HS} <: $DSF{T}
            data::HS
            $HDSF(data::HS) = new(data)
        end
    end
end

type HierarchicalFun{S,T,HS}
    data::HS
    HierarchicalFun(data::HS) = new(data)
end

for HDSF in (:HierarchicalDomain,:HierarchicalSpace,:HierarchicalFun)
    @eval begin
        export $HDSF

        $HDSF{S1,S2}(data::Tuple{S1,S2}) = $HDSF{promote_type(S1,S2),promote_type(eltype(S1),eltype(S2)),Tuple{S1,S2}}(data)
        $HDSF{S1,S2,T1,T2,HS1,HS2}(data::Tuple{$HDSF{S1,T1,HS1},$HDSF{S2,T2,HS2}}) = $HDSF{promote_type(S1,S2),promote_type(eltype(S1),eltype(S2),T1,T2),Tuple{$HDSF{S1,T1,HS1},$HDSF{S2,T2,HS2}}}(data)
        $HDSF{S,S1,T,HS}(data::Tuple{S1,$HDSF{S,T,HS}}) = $HDSF{promote_type(S,S1),promote_type(eltype(S),eltype(S1),T),Tuple{S1,$HDSF{S,T,HS}}}(data)
        $HDSF{S,S1,T,HS}(data::Tuple{$HDSF{S,T,HS},S1}) = $HDSF{promote_type(S,S1),promote_type(eltype(S),eltype(S1),T),Tuple{$HDSF{S,T,HS},S1}}(data)

        $HDSF(data::Vector) = $HDSF(data,round(Int,log2(length(data))))

        function $HDSF(data::Vector,n::Int)
            @assert length(data) == 2^n
            if n == 1
                return $HDSF(tuple(data...))
            elseif n ≥ 2
                return $HDSF(($HDSF(data[1:2^(n-1)],n-1),$HDSF(data[1+2^(n-1):end],n-1)))
            end
        end

        Base.similar(H::$HDSF) = $HDSF(map(similar,data(H)))
        Base.similar{SS,V,T}(H::$HDSF{SS,V,T}, S) = $HDSF(map(A->similar(A,S),data(H)))

        data(H::$HDSF) = H.data

        nlevels(H::$HDSF) = 1+mapreduce(nlevels,max,data(H))

        partition(H::$HDSF) = data(H)

        function partition{H<:$HDSF}(VH::Vector{H})
            n = length(VH)
            H11,H12 = partition(VH[1])
            H1,H2 = fill(H11,n),fill(H12,n)
            for i=1:n
                H1[i],H2[i] = partition(VH[i])
            end
            H1,H2
        end

        collectdata{S,T}(H::$HDSF{S,T,NTuple{2,S}}) = collect(data(H))
        collectdata{S,T,HS}(H::$HDSF{S,T,Tuple{S,$HDSF{S,T,HS}}}) = vcat(H.data[1],collectdata(H.data[2]))
        collectdata{S,T,HS}(H::$HDSF{S,T,Tuple{$HDSF{S,T,HS},S}}) = vcat(collectdata(H.data[1]),H.data[2])
        function collectdata{S}(H::$HDSF{S})
            ret = S[]
            append!(ret,mapreduce(collectdata,vcat,data(H)))
            ret
        end

        Base.convert{S,T,HS}(::Type{$HDSF{S,T,HS}},M::$HDSF) = $HDSF(convert(Vector{S},collectdata(M)))
        Base.promote_rule{S,T,HS,SS,TT,HSS}(::Type{$HDSF{S,T,HS}},::Type{$HDSF{SS,TT,HSS}})=$HDSF{promote_type(S,SS),promote_type(T,TT),promote_type(HS,HSS)}
        Base.eltype{S,T}(::$HDSF{S,T})=T
        Base.eltype{S,T,HS}(::Type{$HDSF{S,T,HS}})=T
    end
end

domain(H::HierarchicalSpace) = HierarchicalDomain(map(domain,data(H)))
space(H::HierarchicalFun) = HierarchicalSpace(map(space,data(H)))

ApproxFun.basistype(H::HierarchicalSpace) = mapreduce(typeof,promote_type,data(H))
ApproxFun.domaindimension(H::HierarchicalSpace) = 1
Space(H::HierarchicalDomain) = HierarchicalSpace(map(Space,data(H)))

PiecewiseSpace(H::HierarchicalSpace) = PiecewiseSpace(map(PiecewiseSpace,data(H))...)

Fun(f::Function,H::HierarchicalDomain) = HierarchicalFun((Fun(f,H.data[1]),Fun(f,H.data[2])))
Fun(f::Function,H::HierarchicalSpace) = HierarchicalFun((Fun(f,H.data[1]),Fun(f,H.data[2])))

ApproxFun.depiece(f::Fun) = f
ApproxFun.depiece(H::HierarchicalFun) = depiece(map(depiece,data(H)))

# algebra

for op in (:+,:-,:.+,:.-,:.*)
    @eval begin
        $op(H::HierarchicalFun) = HierarchicalFun(map($op,data(H)))
        $op(H::HierarchicalFun,a::Number) = HierarchicalFun(($op(H.data[1],a),$op(H.data[2],a)))
        $op(a::Number,H::HierarchicalFun) = HierarchicalFun(($op(a,H.data[1]),$op(a,H.data[2])))

        $op(H::HierarchicalFun,J::HierarchicalFun) = HierarchicalFun(map($op,data(H),data(J)))
        $op(H::HierarchicalFun,J::Fun) = $op(J,full(H))
    end
end

*(H::HierarchicalFun,a::Number) = HierarchicalFun((H.data[1]*a,H.data[2]*a))
*(a::Number,H::HierarchicalFun) = H*a

for op in (:linebilinearform,:bilinearform)
    @eval begin
        $op(H::HierarchicalFun,J::HierarchicalFun) = $op(H.data[1],J.data[1])+$op(H.data[2],J.data[2])
        $op{S<:PiecewiseSpace}(H::HierarchicalFun,J::Fun{S}) = sum(map($op,collectdata(H),pieces(J)))
        $op{S<:PiecewiseSpace}(J::Fun{S},H::HierarchicalFun) = sum(map($op,pieces(J),collectdata(H)))
    end
end

Base.cumsum(H::HierarchicalFun) = HierarchicalFun((cumsum(H.data[1]),sum(H.data[1])+cumsum(H.data[2])))
Base.conj!(H::HierarchicalFun) = (map(conj!,data(H));H)
Base.copy!(H::HierarchicalFun,J::HierarchicalFun) = (map(copy!,data(H),data(J));H)

for op in (:(Base.zero),:(Base.ones),:(Base.abs),:(Base.abs2),:(Base.conj),:(Base.copy),:.^)
    @eval begin
        $op(H::HierarchicalFun) = HierarchicalFun(map($op,data(H)))
    end
end

# Formatted operations

for (op,opformatted) in ((:+,:⊕),(:-,:⊖))
    @eval begin
        function $opformatted{S,T,HS1<:HierarchicalFun,HS2<:HierarchicalFun,P<:PiecewiseSpace}(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P})
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($opformatted(H1,depiece(pieces(J)[1:n1])),$opformatted(H2,depiece(pieces(J)[1+n1:n1+n2]))))
        end
        function $opformatted{S,T,HS1<:Fun,HS2<:HierarchicalFun,P<:PiecewiseSpace}(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P})
            H1,H2 = partition(H)
            n1,n2 = 1,length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($op(H1,depiece(pieces(J)[1:n1])),$opformatted(H2,depiece(pieces(J)[1+n1:n1+n2]))))
        end
        function $opformatted{S,T,HS1<:HierarchicalFun,HS2<:Fun,P<:PiecewiseSpace}(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P})
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),1
            return HierarchicalFun(($opformatted(H1,depiece(pieces(J)[1:n1])),$op(H2,depiece(pieces(J)[1+n1:n1+n2]))))
        end
        function $opformatted{S,T,HS1<:Fun,HS2<:Fun,P<:PiecewiseSpace}(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P})
            H1,H2 = partition(H)
            n1,n2 = 1,1
            return HierarchicalFun(($op(H1,depiece(pieces(J)[1:n1])),$op(H2,depiece(pieces(J)[1+n1:n1+n2]))))
        end

        function $opformatted{S,T,HS1<:HierarchicalFun,HS2<:HierarchicalFun,P<:PiecewiseSpace}(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}})
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($opformatted(depiece(pieces(J)[1:n1]),H1),$opformatted(depiece(pieces(J)[1+n1:n1+n2]),H2)))
        end
        function $opformatted{S,T,HS1<:Fun,HS2<:HierarchicalFun,P<:PiecewiseSpace}(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}})
            H1,H2 = partition(H)
            n1,n2 = 1,length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($op(depiece(pieces(J)[1:n1]),H1),$opformatted(depiece(pieces(J)[1+n1:n1+n2]),H2)))
        end
        function $opformatted{S,T,HS1<:HierarchicalFun,HS2<:Fun,P<:PiecewiseSpace}(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}})
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),1
            return HierarchicalFun(($opformatted(depiece(pieces(J)[1:n1]),H1),$op(depiece(pieces(J)[1+n1:n1+n2]),H2)))
        end
        function $opformatted{S,T,HS1<:Fun,HS2<:Fun,P<:PiecewiseSpace}(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}})
            H1,H2 = partition(H)
            n1,n2 = 1,1
            return HierarchicalFun(($op(depiece(pieces(J)[1:n1]),H1),$op(depiece(pieces(J)[1+n1:n1+n2]),H2)))
        end
    end
end


text:
@testset "1598.crawler-log-folder.jl" begin
    @test min_operations_1198(["d1/","d2/","../","d21/","./"]) == 2
    @test min_operations_1198(["d1/","d2/","./","d3/","../","d31/"]) == 3
    @test min_operations_1198(["d1/","../","../","../"]) == 0
end


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables

function PS4()

url="https://raw.githubusercontent.com/OU-PhD-EConometrics/fall-2021/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df=CSV.read(HTTP.get(url).body, DataFrame)
X= [df.age df.white df.collgrad]
Z=hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4, df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y=df.occ_code

function mlogit(Beta, X, Z,y)

        K = size(X,2)+1
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigBeta = [reshape(Beta,K,J-1) zeros(K)]

        bigZ = zeros(N,J)
        for j=1:J
            bigZ[:,j] = Z[:,j]-Z[:,J]
        end

        num = zeros(N,J)
        dem = zeros(N)
        for j=1:J
            XZ=cat(X,bigZ[:,j],dims=2)
            num[:,j] = exp.(XZ*bigBeta[:,j])
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    mlogit_hat_optim = optimize(b-> mlogit(b,X,Z, y), rand(7*(size(X,2)+1)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    println(mlogit_hat_optim.minimizer)

2. #The coefficient is positive which is closer to our expecation.

3.
using Distributions
include ("lgwt.jl)

#a

d=Normal(0,1)

nodes,weights = lgwt(7,-4,4)

sum(weights.*pdf.(d,nodes))

sum(weights.*nodes.*pdf.(d,nodes))

#b

d=Normal(0,2)
nstd = 5*sqrt(2)
nodes, weights = lgwt (7, -nstd, nstd)

sum(weights.*nodes.*nodes,*pdf.(d,nodes))

nodes, weights = lgwt (10, -nstd, 5*nstd)

sum(weights.*nodes.*nodes,*pdf.(d,nodes))

#the 10 grid setup provides a closer approximation

#c

NumObs = 10^6

x= rand(Normal(0,2),NumObs)
x_sq= x.*X

function f(y)
1/(2*sqrt(2*pi)) * exp(-.5*(y/2)^2)
end

w=(nstd)/NumObs
func1=map(y->f(y),x) .*x_sq
MC1= w * sum(func1)
println

#close to 4

func2= map(y->f(y),x) .* X
MC1= w * sum(func2)
println

#close to 0

func3= map(y->f(y),x)
MC3= w * sum(func3)
println

# Is not close to one

#d

w .* (x_sq .*pdf.(Normal(0,2),x))

w .* (x.* pdf.(Normal(0,2),x))

w .* pdf(Normal(0,2),x)

10^6 is the better approximation

#d Not sure there was really a question there

#4 ??????????????????????????????

#5 ??????????????????????????????

#6 PS(4)


text:
# Simple Send and Accept with Arguments

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

@process lock_process() begin
    while true
        @accept caller my_lock(x::String) begin
            println(x)
        end
        @accept caller my_unlock()
    end
end

the_lock = Nothing

@process p1(i::Integer) begin
    println("$(current_time()): Process p1($i) started")
    @send the_lock my_lock("message from $i")
    println("$(current_time()): Process p1($i) acquired lock")
    work(rand(Uniform(1.0, 10.0)))
    println("$(current_time()): Process p1($i) released lock")
    @send the_lock my_unlock()
    println("$(current_time()): Process p1($i) ended")
end

function run_simulation(n::Integer)
    @simulation begin
        # current_trace!(true)
        global the_lock = @schedule now lock_process()
        for i = 1:n
            @schedule at rand(Uniform(0.0, 10.0)) p1(i)
        end
        start_simulation()
    end
end

run_simulation(10)


text:
"""
    Optimizer

An Optim.jl wrapper to make the choice of the optimizer straightforward in StateSpaceModels.jl
Users can choose among all suitable Optimizers in Optim.jl using very similar syntax.

# Example
```@jldoctest
julia> using Optim

# use a semicolon to avoid displaying the big log
julia> opt = Optimizer(Optim.LBFGS(), Optim.Options(show_trace = true));
```
"""
struct Optimizer
    method::Optim.AbstractOptimizer
    options::Optim.Options
end

function Optimizer(
    method::Optim.AbstractOptimizer;
    options=Optim.Options(; f_tol=1e-6, g_tol=1e-6, iterations=10^5, show_trace=false),
)
    return Optimizer(method, options)
end

# General to every StateSpaceModel, some of them haave trouble to converge
# or have numerical errors with LBFGS
function default_optimizer(::StateSpaceModel)
    return Optimizer(Optim.LBFGS())
end


text:
# OptAnalyzer
# ===========

"""
Every [entry point of optimization analysis](@ref optanalysis-entry) can accept
any of [general JET configurations](@ref JET-configurations) as well as
the following additional configurations that are specific to the optimization analysis.

---
- `skip_noncompileable_calls::Bool = true`:\\
  Julia's runtime dispatch is "powerful" because it can always compile code with concrete
  runtime arguments so that [a "kernel" function](https://docs.julialang.org/en/v1/manual/performance-tips/#kernel-functions)
  runs very effectively even if it's called from a type-instable call site.
  This means, we (really) often accept that some parts of our code are not inferred statically,
  and rather we want to just rely on information that is only available at runtime.
  To model this programming style, the optimization analyzer by default does NOT report any
  optimization failures or runtime dispatches detected within non-concrete calls
  (more correctly, "non-compileable" calls are ignored: see also the note below).
  We can turn off this `skip_noncompileable_calls` configuration to get type-instabilities
  within those calls.
  ```julia-repl
  # the following examples are adapted from https://docs.julialang.org/en/v1/manual/performance-tips/#kernel-functions
  julia> function fill_twos!(a)
             for i = eachindex(a)
                 a[i] = 2
             end
         end;

  julia> function strange_twos(n)
             a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
             fill_twos!(a)
             return a
         end;

  # by default, only type-instabilities within concrete call (i.e. `strange_twos(3)`) are reported
  # and those within non-concrete calls (`fill_twos!(a)`) are not reported
  julia> @report_opt strange_twos(3)
  ═════ 2 possible errors found ═════
  ┌ @ REPL[2]:2 %45(Main.undef, n)
  │ runtime dispatch detected: %45::Type{Vector{_A}} where _A(Main.undef, n::Int64)
  └─────────────
  ┌ @ REPL[2]:3 Main.fill_twos!(%46)
  │ runtime dispatch detected: Main.fill_twos!(%46::Vector)
  └─────────────
  Vector (alias for Array{_A, 1} where _A)

  # we can get reports from non-concrete calls with `skip_noncompileable_calls=false`
  julia> @report_opt skip_noncompileable_calls=false strange_twos(3)
  ═════ 4 possible errors found ═════
  ┌ @ REPL[2]:3 Main.fill_twos!(a)
  │┌ @ REPL[1]:3 Base.setindex!(a, 2, %14)
  ││ runtime dispatch detected: Base.setindex!(a::Vector, 2, %14::Int64)
  │└─────────────
  │┌ @ REPL[1]:3 Base.setindex!(a, 2, i)
  ││┌ @ array.jl:877 Base.convert(_, x)
  │││ runtime dispatch detected: Base.convert(_::Any, x::Int64)
  ││└────────────────
  ┌ @ REPL[2]:2 %45(Main.undef, n)
  │ runtime dispatch detected: %45::Type{Vector{_A}} where _A(Main.undef, n::Int64)
  └─────────────
  ┌ @ REPL[2]:3 Main.fill_twos!(%46)
  │ runtime dispatch detected: Main.fill_twos!(%46::Vector)
  └─────────────
  Vector (alias for Array{_A, 1} where _A)
  ```

  !!! note "Non-compileable calls"
      Julia runtime system sometimes generate and execute native code of an abstract call.
      More technically, when some of call arguments are annotated as `@nospecialize`,
      Julia compiles the call even if those `@nospecialize`d arguments aren't fully concrete.
      `skip_noncompileable_calls = true` also respects this behavior, i.e. doesn't skip
      compileable abstract calls:

          julia> function maybesin(@nospecialize x)
               if isa(x, Number)
                   return sin(x) # this call is dynamically dispatched
               else
                   return 0
               end
           end
           maybesin (generic function with 1 method)

           julia> report_opt((Vector{Any},)) do xs
               for x in xs
                   s = maybesin(x) # this call is resolved statically and compiled
                   s !== 0 && return s
               end
           end
           ═════ 1 possible error found ═════
           ┌ @ none:3 s = Main.maybesin(x)
           │┌ @ none:3 Main.sin(%3)
           ││ runtime dispatch detected: Main.sin(%3::Number)
           │└──────────

           julia> function maybesin(x)  # now `maybesin` is always called with concrete `x`
                      if isa(x, Number)
                          return sin(x) # this call is dynamically dispatched
                      else
                          return 0
                      end
                  end
                  maybesin (generic function with 1 method)

           julia> report_opt((Vector{Any},)) do xs
                      for x in xs
                          s = maybesin(x) # this call is dynamically dispatched
                          s !== 0 && return s
                      end
                  end
           ═════ 1 possible error found ═════
           ┌ @ none:3 Main.maybesin(%21)
           │ runtime dispatch detected: Main.maybesin(%21::Any)
           └──────────

---
- `function_filter = @nospecialize(ft)->true`:\\
  A predicate which takes a function type and returns `false` to skip runtime dispatch
  analysis on the function call. This configuration is particularly useful when your program
  uses a function that is intentionally written to use runtime dispatch.

  ```julia-repl
  # ignores `Core.Compiler.widenconst` calls (since it's designed to be runtime-dispatched):
  julia> function_filter(@nospecialize(ft)) = ft !== typeof(Core.Compiler.widenconst)

  julia> @test_opt function_filter=function_filter f(args...)
  ...
  ```

---
- `skip_unoptimized_throw_blocks::Bool = true`:\\
  By default, Julia's native compilation pipeline intentionally disables inference (and so
  succeeding optimizations too) on "throw blocks", which are code blocks that will eventually
  lead to `throw` calls, in order to ease [the compilation latency problem, a.k.a. "first-time-to-plot"](https://julialang.org/blog/2020/08/invalidations/).
  Accordingly, the optimization analyzer also ignores any performance pitfalls detected
  within those blocks since we _usually_ don't mind if code involved with error handling
  isn't optimized.
  If `skip_unoptimized_throw_blocks` is set to `false`, it doesn't ignore them and will
  report type instabilities detected within "throw blocks".

  See also <https://github.com/JuliaLang/julia/pull/35982>.

  ```julia-repl
  # by default, unoptimized "throw blocks" are not analyzed
  julia> @test_opt sin(10)
  Test Passed
    Expression: #= none:1 =# JET.@test_opt sin(10)

  # we can turn on the analysis on unoptimized "throw blocks" with `skip_unoptimized_throw_blocks=false`
  julia> @test_opt skip_unoptimized_throw_blocks=false sin(10)
  JET-test failed at none:1
    Expression: #= REPL[6]:1 =# JET.@test_call analyzer = JET.OptAnalyzer skip_unoptimized_throw_blocks = false sin(10)
    ═════ 1 possible error found ═════
    ┌ @ math.jl:1221 Base.Math.sin(xf)
    │┌ @ special/trig.jl:39 Base.Math.sin_domain_error(x)
    ││┌ @ special/trig.jl:28 Base.Math.DomainError(x, "sin(x) is only defined for finite x.")
    │││ runtime dispatch detected: Base.Math.DomainError(x::Float64, "sin(x) is only defined for finite x.")
    ││└──────────────────────

  ERROR: There was an error during testing

  # we can also turns off the heuristic itself
  julia> @test_opt unoptimize_throw_blocks=false skip_unoptimized_throw_blocks=false sin(10)
  Test Passed
    Expression: #= REPL[7]:1 =# JET.@test_call analyzer = JET.OptAnalyzer unoptimize_throw_blocks = false skip_unoptimized_throw_blocks = false sin(10)
  ```

---
"""
struct OptAnalyzer{RP,FF} <: AbstractAnalyzer
    state::AnalyzerState
    report_pass::RP
    skip_noncompileable_calls::Bool
    function_filter::FF
    skip_unoptimized_throw_blocks::Bool
    __frame_checks::BitVector # temporary stash to keep per-frame analysis-skip configuration
    __cache_key::UInt
end
function OptAnalyzer(;
    report_pass = OptAnalysisPass(),
    skip_noncompileable_calls::Bool = true,
    function_filter = @nospecialize(ft)->true,
    skip_unoptimized_throw_blocks::Bool = true,
    jetconfigs...)
    state = AnalyzerState(; jetconfigs...)
    # we want to run different analysis with a different filter, so include its hash into the cache key
    cache_key = state.param_key
    cache_key = hash(skip_noncompileable_calls, cache_key)
    cache_key = @invoke hash(function_filter::Any, cache_key::UInt) # HACK avoid dynamic dispatch
    cache_key = hash(skip_unoptimized_throw_blocks, cache_key)
    return OptAnalyzer(
        state,
        report_pass,
        skip_noncompileable_calls,
        function_filter,
        skip_unoptimized_throw_blocks,
        #=__frame_checks=# BitVector(),
        cache_key,
        )
end

# AbstractAnalyzer API requirements
JETInterface.AnalyzerState(analyzer::OptAnalyzer) = analyzer.state
function JETInterface.AbstractAnalyzer(analyzer::OptAnalyzer, state::AnalyzerState)
    return OptAnalyzer(
        state,
        analyzer.report_pass,
        analyzer.skip_noncompileable_calls,
        analyzer.function_filter,
        analyzer.skip_unoptimized_throw_blocks,
        analyzer.__frame_checks,
        analyzer.__cache_key,
        )
end
JETInterface.ReportPass(analyzer::OptAnalyzer) = analyzer.report_pass
JETInterface.get_cache_key(analyzer::OptAnalyzer) = analyzer.__cache_key

JETInterface.vscode_diagnostics_order(analyzer::OptAnalyzer) = false

struct OptAnalysisPass <: ReportPass end

function CC.finish(frame::InferenceState, analyzer::OptAnalyzer)
    ret = @invoke CC.finish(frame::InferenceState, analyzer::AbstractAnalyzer)

    skip = false
    if !skip && analyzer.skip_noncompileable_calls
        if !(is_compileable_frame(frame) || get_entry(analyzer) === get_linfo(frame))
            skip |= true
        end
    end

    push!(analyzer.__frame_checks, skip)
    if !skip
        # report pass for captured variables
        ReportPass(analyzer)(CapturedVariableReport, analyzer, frame)

        # report pass for optimization failure due to recursive calls, etc.
        ReportPass(analyzer)(OptimizationFailureReport, analyzer, frame.result)
    end

    return ret
end

@reportdef struct CapturedVariableReport <: InferenceErrorReport
    name::Union{Nothing,Symbol}
end
get_msg(::Type{CapturedVariableReport}, _, name::Union{Nothing,Symbol}) =
    isnothing(name) ? "captured variable detected" : "captured variable `$name` detected"
print_error_report(io, report::CapturedVariableReport) = printlnstyled(io, "│ ", report.msg; color = ERROR_COLOR)
function (::OptAnalysisPass)(::Type{CapturedVariableReport}, analyzer::OptAnalyzer, frame::InferenceState)
    local reported = false
    for (pc, typ) in enumerate(frame.src.ssavaluetypes::Vector{Any})
        if typ === Core.Box
            stmt = frame.src.code[pc]
            if isexpr(stmt, :(=))
                lhs = first(stmt.args)
                if isa(lhs, SlotNumber)
                    name = frame.src.slotnames[slot_id(lhs)]
                else
                    name = nothing
                end
                add_new_report!(analyzer, frame.result, CapturedVariableReport((frame, pc), name))
                reported |= true
            end
        end
    end
    return reported
end

@reportdef struct OptimizationFailureReport <: InferenceErrorReport end
get_msg(::Type{OptimizationFailureReport}, args...) = "failed to optimize"
function (::OptAnalysisPass)(::Type{OptimizationFailureReport}, analyzer::OptAnalyzer, caller::InferenceResult)
    if !isa(caller.src, OptimizationState)
        add_new_report!(analyzer, caller, OptimizationFailureReport(caller.linfo))
        return true
    end
    return false
end

function CC.finish!(analyzer::OptAnalyzer, frame::InferenceState)
    caller = frame.result

    # get the source before running `finish!` to keep the reference to `OptimizationState`
    src = caller.src

    ret = @invoke CC.finish!(analyzer::AbstractAnalyzer, frame::InferenceState)

    popfirst!(analyzer.__frame_checks) || ReportPass(analyzer)(RuntimeDispatchReport, analyzer, caller, src)

    return ret
end

@reportdef struct RuntimeDispatchReport <: InferenceErrorReport end
get_msg(::Type{RuntimeDispatchReport}, _) = "runtime dispatch detected"
function (::OptAnalysisPass)(::Type{RuntimeDispatchReport}, analyzer::OptAnalyzer, caller::InferenceResult, @nospecialize(src))
    if isa(src, Const) # the optimization was very successful, nothing to report
        return false
    elseif isnothing(src) # the optimization didn't happen, `OptAnalyzer` should have reported `OptimizationFailureReport` for this frame
        return false
    elseif isa(src, OptimizationState) # compiler cached the optimized IR, just analyze it
        opt = src
    else
        # we should already report `OptimizationFailureReport` for this case,
        # and thus this pass should never happen
        Core.eval(@__MODULE__, :(src = $src))
        throw("unexpected state happened, inspect `$(@__MODULE__).src`") # this pass should never happen
    end

    (; src, sptypes, slottypes) = opt

    # branch on https://github.com/JuliaLang/julia/pull/42149
    @static if !isdefined(CC, :mark_throw_blocks!)
        throw_blocks =
            analyzer.skip_unoptimized_throw_blocks && opt.inlining.params.unoptimize_throw_blocks ?
            CC.find_throw_blocks(src.code) : nothing
    end

    local reported = false
    for (pc, x) in enumerate(src.code)
        # branch on https://github.com/JuliaLang/julia/pull/42149
        @static if isdefined(CC, :mark_throw_blocks!)
            if analyzer.skip_unoptimized_throw_blocks
                CC.is_stmt_throw_block(src.ssaflags[pc]) && continue
            end
        else
            if !isnothing(throw_blocks)
                # optimization is intentionally turned off for this block, let's ignore anything here
                CC.in(pc, throw_blocks) && continue
            end
        end
        if isexpr(x, :call)
            ft = widenconst(argextype(first(x.args), src, sptypes, slottypes))
            ft <: Builtin && continue # ignore `:call`s of language intrinsics
            if analyzer.function_filter(ft)
                add_new_report!(analyzer, caller, RuntimeDispatchReport((opt, pc)))
                reported |= true
            end
        end
    end
    return reported
end

# entries
# -------

"""
    report_opt(f, types = Tuple{}; jetconfigs...) -> JETCallResult
    report_opt(tt::Type{<:Tuple}; jetconfigs...) -> JETCallResult

Analyzes the generic function call with the given type signature with [the optimization analyzer](@ref optanalysis),
which collects optimization failures and runtime dispatches involved within the call stack.
"""
function report_opt(@nospecialize(args...);
                    analyzer = OptAnalyzer,
                    jetconfigs...)
    if !(analyzer === OptAnalyzer)
        throw(ArgumentError("`analyzer` is fixed to $OptAnalyzer"))
    end
    return report_call(args...; analyzer, jetconfigs...)
end

"""
    @report_opt [jetconfigs...] f(args...)

Evaluates the arguments to the function call, determines its types, and then calls
[`report_opt`](@ref) on the resulting expression.
As with `@code_typed` and its family, any of [JET configurations](@ref JET-configurations)
or [optimization analysis specific configurations](@ref optanalysis-config) can be given
as the optional arguments like this:
```julia-repl
# reports `rand(::Type{Bool})` with `unoptimize_throw_blocks` configuration turned on
julia> @report_opt unoptimize_throw_blocks=true rand(Bool)
```
"""
macro report_opt(ex0...)
    return var"@report_call"(__source__, __module__, :(analyzer=$OptAnalyzer), ex0...)
end

"""
    @test_opt [jetconfigs...] [broken=false] [skip=false] f(args...)

Tests the generic function call `f(args...)` is free from runtime dispatch.
Returns a `Pass` result if it is, a `Fail` result if if contains any location where runtime
dispatch or optimization failure happens, or an `Error` result if this macro encounters an
unexpected error. When the test `Fail`s, abstract call stack to each problem location will
also be printed to `stdout`.

```julia-repl
julia> @test_opt sincos(10)
Test Passed
  Expression: #= none:1 =# JET.@test_opt sincos(10)
```

As with [`@report_opt`](@ref), any of [JET configurations](@ref JET-configurations) or
[optimization analysis specific configurations](@ref optanalysis-config) can be given
as the optional arguments like this:
```julia-repl
julia> function f(n)
            r = sincos(n)
            # `println` is full of runtime dispatches,
            # but we can ignore the corresponding reports from `Base`
            # with the `target_modules` configuration
            println(r)
            return r
       end;

julia> @test_opt target_modules=(@__MODULE__,) f(10)
Test Passed
  Expression: #= REPL[3]:1 =# JET.@test_call analyzer = JET.OptAnalyzer target_modules = (#= REPL[3]:1 =# @__MODULE__(),) f(10)
```

Like [`@test_call`](@ref), `@test_opt` is fully integrated with [`Test` standard library](https://docs.julialang.org/en/v1/stdlib/Test/).
See [`@test_call`](@ref) for the details.
"""
macro test_opt(ex0...)
    return var"@test_call"(__source__, __module__, :(analyzer=$OptAnalyzer), ex0...)
end

"""
    test_opt(f, types = Tuple{}; broken::Bool = false, skip::Bool = false, jetconfigs...)
    test_opt(tt::Type{<:Tuple}; broken::Bool = false, skip::Bool = false, jetconfigs...)

Tests the generic function call with the given type signature is free from runtime dispatch.
Except that it takes a type signature rather than a call expression, this function works
in the same way as [`@test_opt`](@ref).
"""
function test_opt(@nospecialize(args...);
                  analyzer = OptAnalyzer,
                  kwargs...)
    if !(analyzer === OptAnalyzer)
        throw(ArgumentError("`analyzer` is fixed to $OptAnalyzer"))
    end
    return test_call(args...; analyzer, kwargs...)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_datepickerrange

"""
    ''_datepickerrange(;kwargs...)

A DatepickerRange component.

Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `endDate` (String; optional): A endDate that will be printed when this component is rendered.
- `startDate` (String; optional): A startDate that will be printed when this component is rendered.
"""
function ''_datepickerrange(; kwargs...)
        available_props = Symbol[:id, :endDate, :startDate]
        wild_props = Symbol[]
        return Component("''_datepickerrange", "DatepickerRange", "datepicker_range_component", available_props, wild_props; kwargs...)
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

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

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test

source_terms = source_terms_convergence_test

# BCs must be passed as Dict
boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict(
  :all => boundary_condition
)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

###############################################################################
# Get the uncurved mesh from a file (downloads the file if not available locally)

# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                              mesh_file)

mesh = P4estMesh(mesh_file, initial_refinement_level=0)

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

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

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

summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

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

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback)
###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module StochasticDelayDiffEq
#########################################################
#########################################################
#########################################################
#########################################################
# using DifferentialEquations
using Reexport
@reexport using StochasticDiffEq
# import StochasticDiffEq: calc_J, calc_J!, calc_tderivative!, calc_tderivative
using LinearAlgebra, StaticArrays
using UnPack, DataStructures
using Logging
using RecursiveArrayTools
using DiffEqBase: AbstractSDDEProblem, AbstractSDDEAlgorithm, AbstractRODESolution, AbstractRODEFunction, AbstractSDEIntegrator, AbstractSDDEIntegrator, DEIntegrator, DEAlgorithm, AbstractRODEAlgorithm, AbstractSDEAlgorithm

import DelayDiffEq: constant_extrapolant!, constant_extrapolant, AbstractMethodOfStepsAlgorithm, iscomposite, MethodOfSteps
using DiffEqNoiseProcess

using DelayDiffEq: Discontinuity, HistoryFunction

import RandomNumbers: Xorshifts
using Random
import Base: convert

# #########################################################
# #########################################################
# #########################################################
# #########################################################

include("integrators/type.jl")
include("integrators/interface.jl")
include("integrators/utils.jl")
include("functionwrapper.jl")
include("utils.jl")
include("solve.jl")

end # module


text:
using CreditMetrics
using Statistics

# Consider three counterparties
# Please be sure to understand the single counterparty example first.

# Probabilities of ending up in each rating class after a period
pvs = ProbabilityVector{Float64}[]
push!(pvs, ProbabilityVector([0.01, 0.04, 0.10, 0.20, 0.40, 1 - (0.01 + 0.04 + 0.10 + 0.20 + 0.40)]))
push!(pvs, ProbabilityVector([0.03, 0.07, 0.20, 0.40, 0.20, 1 - (0.03 + 0.07 + 0.20 + 0.40 + 0.20)]))
push!(pvs, ProbabilityVector([0.10, 0.20, 0.35, 0.20, 0.10, 1 - (0.10 + 0.20 + 0.35 + 0.20 + 0.10)]))

# Value of the exposure to each counterparty in each rating class after a period
vals = [10 * [-1.0, -0.2, -0.1, -0.05, 0.0, 0.03], 5 * [-1.0, -0.2, -0.05, 0.0, 0.03, 0.04], 8 * [-1.0, -0.15, 0.0, 0.05, 0.08, 0.1]]

# Consider these in a -2 sigma macroeconomic scenario

# Correlations to a single macro factor
rho = [0.3, 0.2, 0.4]

# Convert to z (Gauss) based thresholds
link = NormalThresholdLink{Float64}()
thrs = [RealLineThresholds(UnitIntervalThresholds(pv), link) for pv in pvs]

# Deform thresholds according to macroeconomic environment
dthrs = [deform(thrs[i], 2.0 * rho[i], sqrt(1 - rho[i]^2)) for i in 1:3]

# Now simulate idiosyncratic effects in this environment
my_stressed_outcomes = [outcome_value(dthrs[i], vals[i], randn(1_000_000)) for i in 1:3]

# Mean outcome is depressed
[mean(my_stressed_outcomes[i]) for i in 1:3]

# The three counterparties are independent given the macro environment
[cov(my_stressed_outcomes[i], my_stressed_outcomes[j]) for i in 1:3, j in 1:3]

# Oppositely, do a lot of macroeconomic scenarios, but mean over idiosyncratic effect.
zmacro = randn(1_000_000)
my_macro_outcomes = [deformation_value(thrs[i], vals[i], link, -rho[i] * zmacro, sqrt(1 - rho[i]^2)) for i in 1:3]

my_portfolio_macro_outcomes = [sum([my_macro_outcomes[i][j] for i in 1:3]) for j in 1:1_000_000]

# Expected loss
expeced_loss_macro = -mean(my_portfolio_macro_outcomes)

# Unexpected loss at 0.5%
unexpected_loss_macro = -quantile(my_portfolio_macro_outcomes, 0.005) - expeced_loss_macro

# Combine both, to do both macro and idiosyncratic
my_full_outcomes = [outcome_value(thrs[i], vals[i], rho[i] * zmacro + sqrt(1 - rho[i]^2) * randn(1_000_000)) for i in 1:3]
my_portfolio_full_outcomes = [sum([my_full_outcomes[i][j] for i in 1:3]) for j in 1:1_000_000]

# Expected loss
expeced_loss_full = -mean(my_portfolio_full_outcomes)

# Unexpected loss at 0.5%
unexpected_loss_full = -quantile(my_portfolio_full_outcomes, 0.005) - expeced_loss_full

# No difference in expected loss
expeced_loss_full - expeced_loss_macro

# Difference in unexpected loss due to idiosyncratic effects
unexpected_loss_full - unexpected_loss_macro


text:
function format_codeblock(meta)
    netString =  """
    register(DataDep(
        \"$(meta.shortname)\",
        \"\"\""""
    netString *= format_meta(meta.fullname, "Dataset")
    netString *= format_meta(meta.website, "Website")
    ## Start of the message
    netString *= format_meta(format_authors(meta.author), "Author")
    netString *= format_meta(meta.maintainer, "Maintainer")
    netString *= format_meta(format_dates(meta.published_date), "Date of Publication")
    netString *= format_meta(format_dates(meta.create_date), "Date of Creation")
    netString *= format_meta(format_dates(meta.modified_date), "Date of Last Modification")
    netString *= format_meta(meta.license, "License")
    netString *= "\n"
    netString *= format_meta(meta.description)
    netString *= format_meta(meta.paper_cite, "\nPlease cite this paper")
    netString *= format_meta(meta.dataset_cite, "\nPlease cite this dataset")
    netString *= "\n\t\"\"\","
    netString = netString |> strip
    ## End of the message
    netString *= format_meta(ismissing(meta.dataurls) ? "missing" : meta.dataurls)
    netString *= ","
    netString *= format_meta(format_checksums(meta.datachecksums))
    netString *= "\n))"

    netString
end

format_meta(::Missing, args...; kwargs...) = ""

function format_meta(data::Any, label="")
    if label != ""
        label*=": "
    end
    return "\n" * indent(label * string(data))
end
format_dates(rawdate) = rawdate #Catches the missing and the string
format_dates(rawdate::Dates.Date) = Dates.format(rawdate, "U d, yyyy")

function format_checksums(csums::Vector)
    csumvec = join(format_checksums.(csums), ", ")
    "[$csumvec]"
end

function format_checksums(csum::Tuple{T,<:AbstractString}) where T<:Symbol
    func = string(csum[1])
    hashstring = format_checksums(csum[2])
    "($func, $hashstring)"
end

format_checksums(::Missing) = missing

function format_checksums(csum::AbstractString)
    if length(csum)>0 "\"$csum\"" else "" end
end

format_authors(authors::Missing) = "Authors not specified"

format_authors(authors::AbstractString) = authors

function format_authors(authors::Vector)
    if length(authors) == 1
        authors[1]
    elseif length(authors) == 2
        authors[1] * ", " * authors[2]
    elseif length(authors) >2
        authors[1] * " et al."
    else
        @warn("Not able to retrieve any authors")
        format_authors(missing)
    end
end

function format_papers(authors::Vector, year::String, name::String, link::String)
    #APA format. Other formats can be included here later.
    join(authors, ", ") * " ($year). " * name * " " * link
end


text:
using Revise
includet(joinpath("..", "init.jl"))

# Load Spine SOC evolution
# NOTE: To do this 1) open the Spine DB out databases 2) copy paste the values
# into the excel file "data/Spine_battery_storage_soc.csv"
df = CSV.read(joinpath(dataPath, "Spine_output_battery_storage_soc.csv"))
eSpineRepDays = df.Spine_rep_days
eSpineRolling = df.Spine_rolling_horizon

# Run GEPPR model as a check:
configFile = joinpath.(dataPath, "GEPPR_input_data", [
        "booleans.yaml",
        "parameters.yaml",
        "storage.yaml",
        "technologies_and_fuels_no_storage.yaml"
    ]
)
param = Dict(
    "modelType" => "Bram storage model",
)
gep = GenerationExpansionPlanningModel(configFile, param)
make_JuMP_model!(gep)
optimize_GEP_model!(gep)
# NOTE: "Spine storage model" should also work above, but I (Seb) had some
# troubles with it (coding bugs on my end, not intrinsic problems
# with the formulation I should say)

if GEPPR.get_model_type(gep) == "Bram storage model"
    eGEPPR = permutedims(
        dropdims(
            GEPPR.calculate_energy_content_evolution(gep).data,
            dims=(2,3)
        )[1,:,:],
        [2,1]
    )[:]
elseif GEPPR.get_model_type(gep) == "Spine storage model"
    eGEPPR = permutedims(
        dropdims(
           gep[:ltseap].data,
           dims=(2,3)
        )[1,:,:],
        [2,1]
    )[:]
end

# Plot simulations side by side
x = 1:7*2*24
plot(
    x ./ 24,
    hcat(eSpineRepDays, eGEPPR, eSpineRolling)[x,:],
    label = ["Spine Rep Days" "Validation model" "Spine rolling horizon"],
    ylabel = "Battery State of charge [GWh]",
    xlabel = "Days [-]"
)

# Save image
savefig(joinpath(dataPath, "battery_soc.png"))


text:
## PARAMETER RANGES

abstract type ParamRange{T} end

Base.isempty(::ParamRange) = false

abstract type Boundedness end

abstract type   Bounded <: Boundedness end
abstract type Unbounded <: Boundedness end

abstract type   LeftUnbounded <: Unbounded end
abstract type  RightUnbounded <: Unbounded end
abstract type DoublyUnbounded <: Unbounded end

struct NumericRange{T,B<:Boundedness,D} <: ParamRange{T}
    field::Union{Symbol,Expr}
    lower::Union{T,Float64}     # Float64 to allow for -Inf
    upper::Union{T,Float64}     # Float64 to allow for Inf
    origin::Float64
    unit::Float64
    scale::D
end

struct NominalRange{T,N} <: ParamRange{T}
    field::Union{Symbol,Expr}
    values::NTuple{N,T}
end

function Base.show(stream::IO,
#                   ::MIME"text/plain",
                   r::NumericRange{T}) where T
    fstr = string(r.field)
    repr = "NumericRange($(r.lower) ≤ $fstr ≤ $(r.upper); "*
        "origin=$(r.origin), unit=$(r.unit))"
    if r.scale isa Symbol
        r.scale !== :linear && (repr *= " on $(r.scale) scale")
    else
        repr = "transformed "*repr
    end
    print(stream, repr)
    return nothing
end

function Base.show(stream::IO,
#                   ::MIME"text/plain",
                   r::NominalRange{T}) where T
    fstr = string(r.field)
    seqstr = sequence_string(collect(r.values))
    repr = "NominalRange($fstr = $seqstr)"
    print(stream, repr)
    return nothing
end

"""
    r = range(model, :hyper; values=nothing)

Define a one-dimensional `NominalRange` object for a field `hyper` of
`model`. Note that `r` is not directly iterable but `iterator(r)` is.

A nested hyperparameter is specified using dot notation. For example,
`:(atom.max_depth)` specifies the `max_depth` hyperparameter of
the submodel `model.atom`.

    r = range(model, :hyper; upper=nothing, lower=nothing,
              scale=nothing, values=nothing)

Assuming `values` is not specified, define a one-dimensional
`NumericRange` object for a `Real` field `hyper` of `model`.  Note
that `r` is not directly iteratable but `iterator(r, n)`is an iterator
of length `n`. To generate random elements from `r`, instead apply
`rand` methods to `sampler(r)`. The supported scales are `:linear`,`
:log`, `:logminus`, `:log10`, `:log2`, or a callable object.

Note that `r` is not directly iterable, but `iterator(r, n)` is, for
given resolution (length) `n`.

By default, the behaviour of the constructed object depends on the
type of the value of the hyperparameter `:hyper` at `model` *at the
time of construction.* To override this behaviour (for instance if
`model` is not available) specify a type in place of `model` so the
behaviour is determined by the value of the specified type.

A nested hyperparameter is specified using dot notation (see above).

If `scale` is unspecified, it is set to `:linear`, `:log`,
`:logminus`, or `:linear`, according to whether the interval `(lower,
upper)` is bounded, right-unbounded, left-unbounded, or doubly
unbounded, respectively.  Note `upper=Inf` and `lower=-Inf` are
allowed.

If `values` is specified, the other keyword arguments are ignored and
a `NominalRange` object is returned (see above).

See also: [`iterator`](@ref), [`sampler`](@ref)

"""
function Base.range(model::Union{Model, Type}, field::Union{Symbol,Expr};
                    values=nothing, lower=nothing, upper=nothing,
                    origin=nothing, unit=nothing, scale::D=nothing) where D
    all(==(nothing), [values, lower, upper, origin, unit]) &&
        throw(ArgumentError("You must specify at least one of these: "*
                            "values=..., lower=..., upper=..., origin=..., "*
                            "unit=..."))

    if model isa Model
        value = recursive_getproperty(model, field)
        T = typeof(value)
    else
        T = model
    end
    if T <: Real && values === nothing
        return numeric_range(T, D, field, lower, upper, origin, unit, scale)
    else
        return nominal_range(T, field, values)
    end
end

function numeric_range(T, D, field, lower, upper, origin, unit, scale)
    lower === Inf &&
        throw(ArgumentError("`lower` must be finite or `-Inf`."))
    upper === -Inf &&
        throw(ArgumentError("`upper` must be finite or `Inf`."))

    lower === nothing && (lower = -Inf)
    upper === nothing && (upper = Inf)

    lower < upper ||
        throw(ArgumentError("`lower` must be strictly less than `upper`."))

    is_unbounded = (lower === -Inf || upper === Inf)

    if origin === nothing
        is_unbounded &&
            throw(DomainError("For an unbounded range you must specify " *
                              "`origin=...` to define a centre.\nTo make " *
                              "the range bounded, specify finite " *
                              "`upper=...` and `lower=...`."))
        origin = (upper + lower)/2
    end
    if unit === nothing
        is_unbounded &&
            throw(DomainError("For an unbounded range you must specify " *
                              "`unit=...` to define a unit of scale.\nTo " *
                              "make the range bounded, specify finite " *
                              "`upper=...` and `lower=...`."))
        unit = (upper - lower)/2
    end
    unit > 0 || throw(DomainError("`unit` must be positive."))
    origin < upper && origin > lower ||
        throw(DomainError("`origin` must lie strictly between `lower` and " *
                          "`upper`."))
    if lower === -Inf
        if upper === Inf
            B = DoublyUnbounded
            scale === nothing && (scale = :linear)
        else
            B = LeftUnbounded
            scale === nothing && (scale = :logminus)
        end
    else
        if upper === Inf
            B = RightUnbounded
            scale === nothing && (scale = :log)
        else
            B = Bounded
            scale === nothing && (scale = :linear)
        end
    end
    lower isa Union{T, Float64} || (lower = convert(T, lower) )
    upper isa Union{T, Float64} || (upper = convert(T, upper) )
    scale isa Symbol && (D = Symbol)
    return NumericRange{T,B,D}(field, lower, upper, origin, unit, scale)
end

nominal_range(T, field, values) = throw(ArgumentError(
   "`$values` must be an instance of type `AbstractVector{<:$T}`."
    * (T <: Model ? "\n Perharps you forgot to instantiate model"
     * "as `$(T)()`" : "") ))

nominal_range(T, field, ::Nothing) = throw(ArgumentError(
"The inferred hyper-parameter type is $T, which is nominal. "*
"If this is true, you must specify values=... "*
"If this is false, specify the correct type as "*
"first argument of `range`, as in "*
"the example, "*
"`range(Int, :dummy, lower=1, upper=10)`. "  ))

function nominal_range(::Type{T}, field, values::AbstractVector{<:T}) where T
    return NominalRange{T,length(values)}(field, Tuple(values))
end

#specific def for T<:AbstractFloat(Allows conversion btw AbstractFloats and Signed types)
function nominal_range(::Type{T}, field,
        values::AbstractVector{<:Union{AbstractFloat,Signed}}) where T<: AbstractFloat
    return NominalRange{T,length(values)}(field, Tuple(values))
end

#specific def for T<:Signed (Allows conversion btw Signed types)
function nominal_range(::Type{T}, field,
               values::AbstractVector{<:Signed}) where T<: Signed
    return NominalRange{T,length(values)}(field, Tuple(values))
end


text:
# pRF fitting in julia - using LsqFit.jl
#
# a worked example of doing the pRF fitting in julia.
# also trying to establish how fast it is + how easy to extend.
#
# ds 2020-10

# dependencies
using NIfTI, MAT
using Plots
using Colors
using Images # mosaicview
using LinearAlgebra
using ImageFiltering # convolution (more choice in algos)
using StatsBase # mean, etc
using Distributions # Normal, for quick bootstrapping
using LsqFit # curve_fit lmfit

# helper functions
include("pRF_nls_helpers.jl")

# set gr backend for plotting
gr()

# for now: data in parent folder
tseries_file = "../tSeries_10_31_18.mat"
stim_image_file = "../pRFStimImage_example.mat"

# unpack the dict that is being read in:
ts = matread(tseries_file)["tSeries"]
stimImage = matread(stim_image_file)["pRFStimImage"]

# display an image from stim_image, 3d video is easy
# just one line of code:
img = permutedims(stimImage["im"], (2, 1, 3));
heatmap(img[:,:,39],c=:hot)

## make a montage of stim image
stimim_ = heatmap(mosaicview(img[:,:,1:5:end]; ncol=6, npad=3, fillvalue=1, rowmajor=true), c=:grays, yflip=true, aspect_ratio=1)
xaxis!(showaxis=:hide, title="Stimulus image moasic", xlabel="frames _across_ then down")
savefig(stimim_, "fig-stim-image.png")
# vectors / extent of x, y
x = stimImage["x"][:,1,1];
y = stimImage["y"][1,:,1];

# mesh versions
mx = stimImage["x"];
my = stimImage["y"];

# Parameters that control the pRF
# (the "region of visual space that drives the response")
x0 = -2.0;
y0 = +2.0;
sigma_xy = 3.0;
amp = 100.0; # an additional scale factor.

# calculate
p = [x0; y0; sigma_xy; amp];
pRFhat = mgauss(p, mx, my);

# heatmap the pRF
# display - a transpose is needed on the Z image / pRFhat aspect_ratio
# as per plotting conventions for heatmap.
# the following should be equivalent to what happens in mrTools.
hm_ = heatmap(x,y,transpose(pRFhat),
        c=:hot,
        yflip=true,
        aspect_ratio=1,
        xaxis="y (deg)", yaxis="x (deg)") # use jet1 for matlab look!
hline!([0.0], c=:white, label="")
vline!([0.0], c=:white, label="")
savefig(hm_, "fig-pRFimage.png") # really best to use SVG/PDF

# get a haemodynamic response function
TR = 1.5; # acq TR
t_hrf = 0:TR:30;
hrf = fmribHRF(t_hrf);
# msw - markerstrokewidth, mc - markercolor, ms = markersize
hrf_ = plot(t_hrf, hrf,
            lw=3, msw=3, m=:circle, ms=8, mc=:white, c=:black,
            xlabel="Time (s)", ylabel="response (a.u.)",
            label="",
            title="haemodynamic response")
savefig(hrf_, "fig-hrf-plot.png") # really best to use SVG/PDF

# details of the stim image sequence
t = transpose(stimImage["t"])
szStimIm = size(stimImage["im"])
nPixels = prod(szStimIm[1:2])
nT = szStimIm[3]

# unwrap the stim image into a nx * ny by t matrix.
# only needed once, so saves time by being available in this form.
@time stimImageUW = reshape(stimImage["im"], nPixels, nT)

# pRF response without HRF effect
@time r = calculate_prf_response(stimImageUW, pRFhat);
# try BenchmarkTools.jl and @btime macro (instead of @time)
# to measure performance numbers a bit more accurately (including warmup)
# on my (slightly tired laptop... ~1.25ms w/o )

# ... and with HRF - doesn't take longer (@btime says :)
@time r_hrf = calculate_prf_response(stimImageUW, pRFhat, hrf);

# make a plot of stimulated pRF response
prf1_ = plot(t, r, c=:black, lw=2, label="pRF", title="simple, simulated pRF response")
plot!(t, r_hrf, c=:blue, lw=2, label="pRF w/ HRF", xlabel="Time (s)", ylabel="response")
savefig(prf1_, "fig-prf-example-plot.png") # really best to use SVG/PDF

## now set up the optimisation / fitting
#
# 1. COST FUNCTION
#
# define the cost function for nonlinear least squares
# nb. data, mx and my and hrf are "baked" into this closure
function cost_func(x)
    measured = [ts...;]
    model = x[4] .* calculate_prf_response(stimImageUW, mgauss(x[1:3], mx, my), hrf)
    # residual = abs.(measured .- model)
    residual = (measured .- model)
    return residual
end

# 2. OPTIMIZE / run NLS
#
# starting parameters and bounds
p0 = [-3.0, -3.0, 1.0, 10000.0]; # starting vals
lb = [-30.0, -30.0, 0.5, 1.0];
ub = [+30.0, +30.0, 8.0, +Inf]

# Levenberg-Marquardt as implemented in LsqFit.
@time res = LsqFit.lmfit(cost_func, p0, Float64[];
                         autodiff=:forward,
                         lower=lb,
                         upper=ub)
res.converged ? println("✓ lmfit converged...") :  println("X lmfit didn't converge...")

# 2b. ELABORATED VERSION / using gradient w/ autodiff
#
# there may be a bit more accuracy and speed if we also try to use
# gradient (forwardiff / numerical, but could think about explicitly defining
# jacobian?)
@time resad = LsqFit.lmfit(cost_func, p0, Float64[],
            autodiff=:forward)
resad.param

# 3. CHECK and PLOT results
#
# the solution parameters are here
res.param
# LsqFit.jl includes fumctions for CI and SE.
ci = confidence_interval(res, 0.01)
se = standard_errors(res)

# this function takes a param vector and turns it into a timeseries prediction
# cf the cost_func... this makes plotting and bootstrapping responses easier.
function calcFit(x)
    x[4] .* calculate_prf_response(stimImageUW, mgauss(x[1:3], mx, my), hrf)
end

# this is taking the SE of each parameter and making the
# extreme prediction. NB! ignores any interactions / covariance in the
# parameters in the model... but easy to visualise
seMinus = calcFit(res.param .- se)
sePlus = calcFit(res.param .+ se)

perror_ = plot(ts, c=:black, lw=2, label="", title="Naive approach to error bars, pHat ± 1 SE ")
plot!(calcFit(res.param), c=:pink, ribbon=(seMinus, sePlus), label="+/- 1 SE")
plot!(calcFit(res.param), c=:red, lw=2, label="NLS fit")
savefig(perror_, "fig-prf-tc+error-plot.png") # really best to use SVG/PDF

# 4. ELABORATED version of error bars / bootstrapping
#
# make a multivariate distribution of params
d = Distributions.MultivariateNormal(res.param, se)

println("get a 1000 samples of params")
@time d_boostrapped = rand(d, 1000)

println("get a 1000 boostrap replicates of timeseries")
@time boots = mapslices(calcFit, d_boostrapped; dims=1)

# and plot those boostrap replicates
pboot_ = plot(t, boots, la=0.1, c=:gray, label="")
plot!(t, ts, c=:black, lw=2, label="", title="Boostrapped error bars, pHat, 1000 boostraps ")
plot!(t, calcFit(res.param), c=:red, lw=2, label="NLS fit")
savefig(pboot_, "fig-prf-tc+boot-plot.png") # really best to use SVG/PDF

## a way to plot pRF shapes (parametric plots)
# r1 = res.param .+ se
# r2 = res.param .- se
# xₜ(t) = r1[3] .* cos(t) .+ r1[1] # , r2[3] .* cos(t) .+ r2[1])
# yₜ(t) = r1[3] .* sin(t) .+ r1[2] # , r2[3] .* sin(t) .+ r2[2])

# plot(xₜ, yₜ, 0, 2π, leg=false, fill=(0, :orange),        xlim=(-15, +15), ylim=(-15, +15), aspect_ratio=1)

# https://stackoverflow.com/a/56174228/4961292

## display

max_stars_repo_name: schluppeck/prf-explainer

max_issues_repo_name: schluppeck/prf-explainer

max_forks_repo_name: schluppeck/prf-explainer


text:
function linear_scheme(offset_min::Int64, offset_max::Int64)
    float_coeff = float_flux_coefficients(
                      [i for i in offset_min - 1: offset_max])
    return Linear(offset_min,
                  offset_max,
                  float_coeff)
end

function flux_positive_linear(scheme::Linear,
                              vector_in::Array{<:AbstractFloat, 1})
    return dot(scheme.flux_coefficients, vector_in)
end # flux_positive_linear


text:
## Linear Algebra
Base.pointer(x::ComponentArray{T,N,A,Axes}) where {T,N,A<:DenseArray,Axes} = pointer(getdata(x))

Base.unsafe_convert(::Type{Ptr{T}}, x::ComponentArray{T,N,A,Axes}) where {T,N,A,Axes} = Base.unsafe_convert(Ptr{T}, getdata(x))

Base.strides(x::ComponentArray) = strides(getdata(x))

Base.stride(x::ComponentArray, k) = stride(getdata(x), k)
Base.stride(x::ComponentArray, k::Int64) = stride(getdata(x), k)

ArrayInterface.lu_instance(jac_prototype::ComponentArray) = ArrayInterface.lu_instance(getdata(jac_prototype))

## Vector to matrix concatenation
Base.hcat(x::CV...) where {CV<:ComponentVector} = ComponentArray(reduce(hcat, getdata.(x)), getaxes(x[1])[1], FlatAxis())

Base.vcat(x::ComponentVector, y::AbstractVector) = vcat(getdata(x), y)
Base.vcat(x::AbstractVector, y::ComponentVector) = vcat(x, getdata(y))
function Base.vcat(x::ComponentVector, y::ComponentVector)
    if reduce((accum, key) -> accum || (key in keys(x)), keys(y); init=false)
        return vcat(getdata(x), getdata(y))
    else
        return ComponentArray(x; NamedTuple(y)...)
    end
end
Base.vcat(x::CV...) where {CV<:AdjOrTransComponentArray} = ComponentArray(reduce(vcat, map(y->getdata(y.parent)', x)), getaxes(x[1]))
Base.vcat(x::ComponentVector...) = reduce((x1, x2) -> vcat(x1, x2), x)
Base.vcat(x::ComponentVector, args...) = vcat(getdata(x), getdata.(args)...)

# While there are some cases where these were faster, it is going to be almost impossible to
# to keep up with method ambiguity errors due to other array types overloading *, /, and \.
# Leaving these here and commented out for now, but will delete them later.

# # Avoid slower fallback
# for f in [:(*), :(/), :(\)]
#     @eval begin
#         # The normal stuff
#         Base.$f(x::ComponentArray, y::AbstractArray) = $f(getdata(x), y)
#         Base.$f(x::AbstractArray, y::ComponentArray) = $f(x, getdata(y))
#         Base.$f(x::ComponentArray, y::ComponentArray) = $f(getdata(x), getdata(y))

#         # A bunch of special cases to avoid ambiguous method errors
#         Base.$f(x::ComponentArray, y::AbstractMatrix) = $f(getdata(x), y)
#         Base.$f(x::AbstractMatrix, y::ComponentArray) = $f(x, getdata(y))

#         Base.$f(x::ComponentArray, y::AbstractVector) = $f(getdata(x), y)
#         Base.$f(x::AbstractVector, y::ComponentArray) = $f(x, getdata(y))
#     end
# end

# # Adjoint/transpose special cases
# for f in [:(*), :(/)]
#     @eval begin
#         Base.$f(x::Adjoint, y::ComponentArray) = $f(getdata(x), getdata(y))
#         Base.$f(x::Transpose, y::ComponentArray) = $f(getdata(x), getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector) where T = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector) where T = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentMatrix{T,A,Axes}) where {T,A,Axes} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentMatrix{T,A,Axes}) where {T,A,Axes} = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractMatrix{T}}, y::ComponentVector) where {T} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractMatrix{T}}, y::ComponentVector) where {T} = $f(x, getdata(y))

#         Base.$f(x::ComponentArray, y::Adjoint{T,<:AbstractVector{T}}) where T = $f(getdata(x), y)
#         Base.$f(x::ComponentArray, y::Transpose{T,<:AbstractVector{T}}) where T = $f(getdata(x), y)

#         Base.$f(x::ComponentArray, y::Adjoint{T,<:ComponentVector}) where T = $f(getdata(x), getdata(y))
#         Base.$f(x::ComponentArray, y::Transpose{T,<:ComponentVector}) where T = $f(getdata(x), getdata(y))

#         # There seems to be a new method in Julia > v.1.4 that specializes on this
#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Number,A,Axes} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Number,A,Axes} = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Real,A,Axes} = $f(getdata(x), getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Real,A,Axes} = $f(getdata(x), getdata(y))
#     end
# end


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: csv2h5.jl
#
#  Decription:
#       Transform CSV files into HDF5 stores
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-24   wm              Initial version
#
################################################################################

using HDF5

const THIS_DIR = dirname(Base.source_path())

function main()
    const DATA_DIR = "$(THIS_DIR)/../../data"
    const CSV_DIR = DATA_DIR
    const TRAIN_CSV = "$(CSV_DIR)/trainingData.csv"
    const TEST_CSV = "$(CSV_DIR)/testingData.csv"

    const TRAIN = readdlm(TRAIN_CSV, ',', ASCIIString, header=true)[1]
    const TEST = readdlm(TEST_CSV, ASCIIString)

    h5write("trainXmp3.h5", "lid/train/X/mp3", TRAIN[:, 1])
    h5write("trainYlang.h5", "lid/train/y/lang", TRAIN[:, 2])
    h5write("testXmp3.h5", "lid/test/X/mp3", TEST[:, 1])

    const LANGUAGES = Set(TRAIN[:, 2])
    ids = [0:length(LANGUAGES) - 1]
    lang2id = [x::ASCIIString => y::Int for (x, y) in zip(LANGUAGES, ids)]

    labels = [lang2id[lang] for lang in TRAIN[:, 2]]
    h5write("trainYlabels.h5", "lid/train/y/labels", labels)

    h5write("ydict.h5", "lid/data/y/lang", vcat(keys(lang2id)...))
    fid = h5open("ydict.h5", "r+")
    lang_dset = d_create(fid, "lid/data/y/labels", datatype(Int), dataspace(size(ids)))
    lang_dset[:] = vcat(values(lang2id)...)
    close(fid)
end

if ~isinteractive()
    main()
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -6.148474362391325 <= q <= 10.677081718106185)

                     add_NL_constraint(m, :(log(1 + exp(-0.2518902526786948 + 0.9847866884384731*log(1 + exp(0.2752793536861313 + -0.1568397657479923*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1062303426643516*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + 0.5642457669318222*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + -0.7649756572019379*log(1 + exp(0.9918790053549298 + 0.85981905523253*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + 0.9639515148457134*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.08816238542180388*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + 0.10087195663512549*log(1 + exp(-0.8202416682813327 + -0.17766510212211006*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1793020696087071*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.41892665263312834*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))))) + log(1 + exp(0.638435553115452 + 0.09453389081965424*log(1 + exp(0.2752793536861313 + -0.1568397657479923*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1062303426643516*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + 0.5642457669318222*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + -0.07927014075141203*log(1 + exp(0.9918790053549298 + 0.85981905523253*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + 0.9639515148457134*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.08816238542180388*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + 0.2624391124261729*log(1 + exp(-0.8202416682813327 + -0.17766510212211006*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1793020696087071*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.41892665263312834*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))))) + log(1 + exp(0.5273095272255199 + 0.4025366346817978*log(1 + exp(0.2752793536861313 + -0.1568397657479923*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1062303426643516*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + 0.5642457669318222*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + -0.585290488444234*log(1 + exp(0.9918790053549298 + 0.85981905523253*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + 0.9639515148457134*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.08816238542180388*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))) + 0.1647823489958915*log(1 + exp(-0.8202416682813327 + -0.17766510212211006*log(1 + exp(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5]))) + -0.1793020696087071*log(1 + exp(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5]))) + -0.41892665263312834*log(1 + exp(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@inline function _abs(x::T) where T
    CUDA.abs(x)
end

@inline function _asin(x::T) where T
    CUDA.asin(x)
end

@inline function _cathethus(hyp::T, cos_theta::T) where T
    CUDA.sqrt(hyp * hyp * (1 - cos_theta * cos_theta))
end

@inline function _clamp(x::T, lo::T, hi::T) where T
    CUDA.max(lo, CUDA.min(hi, x))
end

@inline function _cos(u::CuPosition{T}, unorm::T, v::CuPosition{T}, vnorm::T) where T
    _dot(u, v) / (unorm * vnorm)
end

@inline function _cross(u::CuPosition{T}, v::CuPosition{T}) where T
    (
        x = u.y * v.z - u.z * v.y,
        y = u.z * v.x - u.x * v.z,
        z = u.x * v.y - u.y * v.x
    )
end

@inline function _dot(u::CuPosition{T}, v::CuPosition{T}) where T
    CUDA.fma(u.x, v.x, CUDA.fma(u.y, v.y, u.z * v.z))
end

@inline function _log(x::T) where T
    CUDA.log(x)
end

@inline function _norm(v::CuPosition{T}) where T
    CUDA.sqrt(_dot(v, v))
end

@inline function _sign(x::T) where T
    x == 0 ? zero(x) : x / _abs(x)
end

@inline function _sqrt(x::T) where T
    CUDA.sqrt(x)
end

@inline function _smul(v::CuPosition{T}, s::T) where T
    (x = s * v.x, y = s * v.y, z = s * v.z)
end

@inline function _sub(u::CuPosition{T}, v::CuPosition{T}) where T
    (x = u.x - v.x, y = u.y - v.y, z = u.z - v.z)
end


text:
#First Level Abstraction

const AbstractACPowerModel = Union{PM.AbstractACPForm, PM.AbstractACRForm, PM.AbstractWRForm}

const AbstractDCPowerModel = Union{PM.AbstractDCPForm, PM.AbstractDCPLLForm} #Adopt from PowerModels, LL -> Line Losses

abstract type CopperPlatePowerModel <: PM.AbstractPowerFormulation end

#Second Level Abstraction AC

#adopted from PowerModels

const StandardAC = PM.StandardACPForm

#Second Level Abstraction DC

##This line is from PowerModels, needs to be removed later

abstract type DCAngleForm <: PM.DCPlosslessForm end

abstract type DCAngleLLForm <: PM.AbstractDCPLLForm end

#abstract type StandardDCPLLForm<: PM.AbstractDCPLLForm end

abstract type AbstractFlowForm <: PM.NFAForm end

abstract type AbstractFlowLLForm <: PM.AbstractDCPLLForm end

#Third Level of Abstraction.

abstract type StandardPTDF <: AbstractFlowForm end

abstract type StandardPTDFLL <: AbstractFlowLLForm end

abstract type StandardNetFlow <: AbstractFlowForm end

abstract type StandardNetFlowLL <: AbstractFlowLLForm end


text:
Try.Ok() = Try.Ok(nothing)
Try.Ok(::Type{T}) where {T} = Try.Ok{Type{T}}(T)

Try.Err(value) = Try.Err(value, maybe_backtrace())
Try.Err{E}(value) where {E} = Try.Err{E}(value, maybe_backtrace())

Try.unwrap(ok::Ok) = ok.value
Try.unwrap(err::Err) = _throw(err)

Try.unwrap_err(ok::Ok) = throw(Try.IsOkError(ok))
Try.unwrap_err(err::Err) = err.value

function _throw(err::Err)
    if err.backtrace === nothing
        throw(err.value)
    else
        throw(ErrorTrace(err.value, err.backtrace))
    end
end

function Try.IsOkError(ok)
    if Try.iserr(ok)
        error("unexpected error value: ", ok)
    end
    return _IsOkError(ok)
end

Base.convert(::Type{Ok{T}}, ok::Ok) where {T} = Ok{T}(ok.value)
Base.convert(::Type{Err{E}}, err::Err) where {E} = Err{E}(err.value)
# An interesting approach may be to simply throw the `err.value` if it is not a
# subtype of `E`.  It makes the error value propagation pretty close to the
# chain-of-custody Julep.  Maybe this should be done only when the destination
# type is `AbstractResult{<:Any,E′}` s.t. `!(err.value isa E′)`.

Try.oktype(::Type{R}) where {T,R<:AbstractResult{T}} = T
Try.oktype(result::AbstractResult) = Try.oktype(typeof(result))

Try.errtype(::Type{R}) where {E,R<:AbstractResult{<:Any,E}} = E
Try.errtype(result::AbstractResult) = Try.errtype(typeof(result))

Try.isok(::Ok) = true
Try.isok(::Err) = false

Try.iserr(::Ok) = false
Try.iserr(::Err) = true


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mc01pd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 10
    REZ = Array{Float64,1}(undef, KMAX)
    IMZ = Array{Float64,1}(undef, KMAX)
    P = Array{Float64,1}(undef, KMAX+1)
    DWORK = Array{Float64,1}(undef, KMAX+1)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    K = parse(BlasInt, vs[1])
    if ( K<0 || K>KMAX )
        @error "Illegal K=$K"
    end

    vs = String[]
    _isz,_jsz = (K,2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
        REZ[i],IMZ[i] = parsex.(Float64, vs[2i-1:2i])
    end
    close(f)
# interp call 1

    INFO = SLICOT.mc01pd!(K, REZ, IMZ, P)
    @test INFO == 0
    INFO == 0 || return

    println(io, "P:")
    show(io, "text/plain", P[1:K+1])
    println(io)

end # run_mc01pd()


text:
using Flux, Functors

"""
radial function
"""
struct Radfunc
    p::AbstractVector
    f
    rmin
    rmax
end
@functor Radfunc
Flux.trainable(m::Radfunc) = [m.f]

function Radfunc(; rmin = 0, rmax = 1e16, n = 32)
    p = ones(4)
    f = Chain(Dense(1, n, leakyrelu), Dense(n, 1))
    Radfunc(p, f, rmin, rmax)
end

function (m::Radfunc)(r)
    @unpack f, p, rmin, rmax = m
    c, k = p[1:2]
    # rmin < r < rmax ? c * exp(-abs(k) * r) : 0.0
     f([r])[1]
end


text:
using GridapDistributedPETScWrappers
using GridapDistributedPETScWrappers.C
using MPI
using LinearAlgebra
using SparseArrays
using Test

# determine scalar type of current run
global ST = Float64  # scalar type

function RC(x::Number)
# used to test real, complex
  if ST == Float64
    return Float64(real(x))
  elseif ST == Float32
    return Float32(real(x))
  else  # scalar_type == 3
    return complex(x)
  end
end

function RC(x::AbstractArray)
# used to test real, complex
  if ST == Float64
    tmp = similar(x, ST)
    for i=1:length(x)
      tmp[i] = Float64(real(x[i]))
    end
    return tmp
  elseif ST == Float32
    tmp = similar(x, ST)
    for i=1:length(x)
      tmp[i] = ST(real(x[i]))
    end
    return tmp
  else  # scalar_type == 3
    return x
  end
end

# convert to PetscReal
function RT(x::Number)
  if ST == Float64 || ST == Complex128
    return Float64(x)
  else
    return Float32(x)
  end

end

function mymult(A::GridapDistributedPETScWrappers.C.Mat{T}, x::GridapDistributedPETScWrappers.C.Vec, b::GridapDistributedPETScWrappers.C.Vec) where {T}
# matrix multiplication function for the shell matrix A
# A performs the action of A = diagm(1:sys_size)

  bigx = Vec{T}(x, first_instance=false)
  bigb = Vec{T}(b, first_instance=false)
  localx = LocalVector_readonly(bigx)
  localb = LocalVector(bigb)
  for i=1:length(localx)
    localb[i] = i*localx[i]
  end

  restore(localx)
  restore(localb)
  return GridapDistributedPETScWrappers.C.PetscErrorCode(0)
end


text:
# Autogenerated wrapper script for rtmidi_jll for aarch64-apple-darwin
export rtmidi

JLLWrappers.@generate_wrapper_header("rtmidi")
JLLWrappers.@declare_library_product(rtmidi, "@rpath/librtmidi.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        rtmidi,
        "lib/librtmidi.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset verbose = true "Test Init Opts" begin
    using JavaCall

    @testset "init defaultopts" begin
        opts = defaultopts()
        @test_false opts.fromcurrentvm
        @test_nothing opts.javahome
        @test opts.jvmclasspath == []
        @test opts.jvmopts == []
    end

    @testset "init fromcurrentvm" begin
        opts = fromcurrentvm()
        @test opts.fromcurrentvm
        @test_nothing opts.javahome
        @test opts.jvmclasspath == []
        @test opts.jvmopts == []
    end

    @testset "init forjavahome" begin
        javahome = "JavaHome"
        opts = forjavahome(javahome)
        @test_false opts.fromcurrentvm
        @test opts.javahome == javahome
        @test opts.jvmclasspath == []
        @test opts.jvmopts == []
    end

    @testset "set/unset fromcurrentvm" begin
        opts = defaultopts()
        @test_false opts.fromcurrentvm
        setfromcurrentvm!(opts)
        @test opts.fromcurrentvm
        unsetfromcurrentvm!(opts)
        @test_false opts.fromcurrentvm
    end

    @testset "set/unset javahome" begin
        opts = defaultopts()
        @test_false opts.fromcurrentvm
        setfromcurrentvm!(opts)
        @test opts.fromcurrentvm
        unsetfromcurrentvm!(opts)
        @test_false opts.fromcurrentvm
    end

    @testset "push classpath" begin
        opts = defaultopts()
        firstentry = "First Entry"
        secondentry = "Second Entry"
        thirdentry = "Third Entry"
        fourthentry = "Fourth Entry"

        @test opts.jvmclasspath == []
        pushclasspath!(opts, firstentry)
        @test opts.jvmclasspath == [firstentry]
        pushclasspath!(opts, secondentry)
        @test opts.jvmclasspath == [firstentry, secondentry]
        pushclasspath!(opts, thirdentry, fourthentry)
        @test opts.jvmclasspath == [firstentry, secondentry, thirdentry, fourthentry]
    end

    @testset "push options" begin
        opts = defaultopts()
        firstopt = "First Option"
        secondopt = "Second Option"
        thirdopt = "Third Option"
        fourthopt = "Fourth Option"

        @test opts.jvmopts == []
        pushoptions!(opts, firstopt)
        @test opts.jvmopts == [firstopt]
        pushoptions!(opts, secondopt)
        @test opts.jvmopts == [firstopt, secondopt]
        pushoptions!(opts, thirdopt, fourthopt)
        @test opts.jvmopts == [firstopt, secondopt, thirdopt, fourthopt]
    end
end


text:
using FactCheck: facts, context, @fact, not
using LatBo.geometry: is_in_pipe, is_in_half_space
using LatBo.Playground: SOLID, FLUID, INLET, OUTLET, NOTHING
import LatBo.Playground.initialize

facts("Constructing playground with a single pipe") do

    context("2d") do
        # Create playground a single pipe
        playpen = initialize((40, 40)) do i, j
            coords = Float64[i, j]
            if !is_in_pipe(coords, Float64[0, 1], Float64[20.5, 20.5], 18.)
                return SOLID
            elseif is_in_half_space(coords, Float64[0, -1], Float64[0, 1])
                return INLET
            elseif is_in_half_space(coords, Float64[0, 1], Float64[0, 40])
                return OUTLET
            else
                return FLUID
            end
        end

        # Check we got what we asked for
        @fact playpen .!= NOTHING --> all
        for indices in ([3, 20], [20, 5], [38, 35])
            @fact playpen[indices...] --> FLUID
        end
        for indices in ([2, 20], [39, 35], [2, 1], [2, 40])
            @fact playpen[indices...] --> SOLID
        end
        @fact playpen[3:end-3, 1] .== INLET --> all
        @fact playpen[3:end-3, 40] .== OUTLET --> all
    end

    context("3d") do
        # Create playground a single pipe
        playpen = initialize((10, 10, 20)) do i, j, k
            coords = Float64[i, j, k]
            if !is_in_pipe(coords, [0., 0, 1], [5.5, 5.5, 10.5], 4.)
                return SOLID
            elseif is_in_half_space(coords, [0., 0, -1], [0., 0, 1])
                return INLET
            elseif is_in_half_space(coords, [0., 0, 1], [0., 0, 20])
                return OUTLET
            else
                return FLUID
            end
        end

        # Check we got what we asked for
        @fact playpen .!= NOTHING --> all
        in_sphere(i, j) =  (i - 5.5)^2 + (j - 5.5)^2 <= 16
        for i = 1:size(playpen, 1), j = 1:size(playpen, 2)
            @fact playpen[i, j, 10] --> (in_sphere(i, j) ? FLUID: SOLID)
        end
        for i = 1:size(playpen, 1), j = 1:size(playpen, 2)
            @fact playpen[i, j, 1] --> (in_sphere(i, j) ? INLET: SOLID)
        end
        for i = 1:size(playpen, 1), j = 1:size(playpen, 2)
            @fact playpen[i, j, 20] --> (in_sphere(i, j) ? OUTLET: SOLID)
        end
    end
end


text:
module CompoundDistributions

using StatsFuns.RFunctions: betarand, gammarand
using LogProbs

# distribution types
export BetaBern, DirMul, DirCat, UniCat, CatDist, ChineseRest, DirCatVar, DirCatVarArith

# conditional distribution type
export SimpleCond

# functions
export support, sample, logscore, add_obs!, rm_obs!

# helper function
export categorical_sample

abstract type Distribution{T} end

function categorical_sample(tokens, weights)
    x = rand() * sum(weights)
    cum_weights = zero(eltype(weights))
    for (t, w) in zip(tokens, weights)
        cum_weights += w
        if cum_weights > x
            return t
        end
    end
end

categorical_sample(d::Dict) = categorical_sample(keys(d), values(d))
categorical_sample(v::Vector) = categorical_sample(1:length(v), v)

############################
### Beta Bernoulli Class ###
############################

mutable struct BetaBern{T, C} <: Distribution{T}
    heads :: T
    tails :: T
    alpha :: C
    beta  :: C
end

BetaBern(a, b) = BetaBern(true, false, a, b)
BetaBern(support) = BetaBern(support..., 1, 1)

support(bb::BetaBern) = [bb.heads, bb.tails]

function sample(bb::BetaBern)
    p = betarand(bb.alpha, bb.beta)
    categorical_sample((bb.heads, bb.tails), (p, 1-p))
end

function logscore(bb::BetaBern, obs)
    k = obs == bb.heads
    LogProb(lbeta(k+bb.alpha, 1-k+bb.beta) - lbeta(bb.alpha, bb.beta))
end

add_obs!(bb::BetaBern, obs) = obs == bb.heads ? bb.alpha += 1 : bb.beta += 1
rm_obs!(bb::BetaBern, obs)  = obs == bb.heads ? bb.alpha -= 1 : bb.beta -= 1

###################################
### Dirichlet Multinomial Class ###
###################################

mutable struct DirMul{T, C} <: Distribution{T}
    counts :: Dict{T, C}
end

DirMul(support) = DirMul(Dict(x => 1.0 for x in support))
support(dm::DirMul) = keys(dm.counts)

function sample(dm::DirMul, n)
    weights = [gammarand(c, 1) for c in values(dm.counts)]
    d = Dict(k=>0 for k in keys(dm.counts))
    for i in 1:n
        d[categorical_sample(keys(dm.counts), weights)] += 1
    end
    d
end

function logscore(dm::DirMul, obs::AbstractDict)
    n = sum(values(obs))
    LogProb(
        log(n) + lbeta(sum(values(dm.counts)), n) -
        sum(
            log(obs[x]) + lbeta(dm.counts[x], obs[x])
            for x in keys(dm.counts) if obs[x] > 0
        )
    )
end

function add_obs!(dm::DirMul, obs::AbstractDict)
    for x in keys(obs)
        dm.counts[x] += obs[x]
    end
end

function rm_obs!(dm::DirMul, obs::AbstractDict)
    for x in keys(obs)
        dm.counts[x] -= obs[x]
    end
end

###################################
### Dirichlet Categorical Class ###
###################################

mutable struct DirCat{T, C} <: Distribution{T}
    counts :: Dict{T, C}
end

DirCat(support) = DirCat(Dict(x => 1.0 for x in support))
support(dc::DirCat) = keys(dc.counts)

function sample(dc::DirCat)
    weights = [gammarand(c, 1) for c in values(dc.counts)]
    categorical_sample(keys(dc.counts), weights)
end

function logscore(dc::DirCat, obs)
    LogProb(lbeta(sum(values(dc.counts)), 1) - lbeta(dc.counts[obs], 1))
end

function add_obs!(dc::DirCat, obs)
    dc.counts[obs] += 1
end

function rm_obs!(dc::DirCat, obs)
    dc.counts[obs] -= 1
end

#########################################################
### Dirichlet Categorical for Variational Bayes Class ###
#########################################################

mutable struct DirCatVar{T, C} <: Distribution{T}
    counts :: Dict{T, C}
end

DirCatVar(support) = DirCatVar(Dict(x => 1.0 for x in support))
support(dc::DirCatVar) = keys(dc.counts)

function sample(dc::DirCatVar)
    weights = [gammarand(c, 1) for c in values(dc.counts)]
    categorical_sample(keys(dc.counts), weights)
end

function logscore(dc::DirCatVar, obs)
    LogProb(digamma(dc.counts[obs]) - digamma(sum(values(dc.counts))))
end

add_obs!(dc::DirCatVar, obs) = dc.counts[obs] += 1
rm_obs!(dc::DirCatVar, obs) = dc.counts[obs] -= 1

###########################################################################
### Dirichlet Categorical for Variational Bayes Class - Arithmetic mean ###
###########################################################################

mutable struct DirCatVarArith{T, C} <: Distribution{T}
    counts :: Dict{T, C}
end

DirCatVarArith(support) = DirCatVarArith(Dict(x => 1.0 for x in support))
support(dc::DirCatVarArith) = keys(dc.counts)

function sample(dc::DirCatVarArith)
    weights = [gammarand(c, 1) for c in values(dc.counts)]
    categorical_sample(keys(dc.counts), weights)
end

function logscore(dc::DirCatVarArith, obs)
    LogProb((dc.counts[obs])/(sum(values(dc.counts))))
end

add_obs!(dc::DirCatVarArith, obs) = dc.counts[obs] += 1
rm_obs!(dc::DirCatVarArith, obs) = dc.counts[obs] -= 1

######################################
### Categorical Distribution Class ###
######################################

mutable struct CatDist{T} <: Distribution{T}
    probs :: Dict{T, LogProb}
end

support(cd::CatDist) = keys(cd.probs)

logscore(cd::CatDist, x) = cd.probs[x]
sample(cd::CatDist) = categorical_sample(cd.probs)
add_obs!(cd::CatDist, obs, context) = nothing
remove_obs!(cd::CatDist, obs, context) = nothing

#################################
### Uniform Categorical Class ###
#################################

struct UniCat{T} <: Distribution{T}
    support :: Vector{T}
end

UniCat(support) = UniCat(collect(support))
support(uc::UniCat) = uc.support

sample(uc::UniCat) = uc.support[rand(1:length(uc.support))]
logscore(uc::UniCat, obs) = obs in uc.support ? LogProb(1 / length(uc.support)) : zero(LogProb)
add_obs!(uc::UniCat, obs, context) = nothing
rm_obs!(uc::UniCat, obs, context) = nothing

################################
### Chinese Restaurant Class ###
################################

mutable struct ChineseRest{Dish, Dist <: Distribution{Dish}} <: Distribution{Dish}
    a :: Float64 # discount parameter
    b :: Float64 # crp parameter
    basedist :: Dist
    tables :: Dict{Dish, Vector{Int}}
    num_tables :: Int
    num_customers :: Int
end

ChineseRest(a, b, basedist::Distribution{T}) where T =
    ChineseRest(a, b, basedist, Dict{T, Vector{Int}}(), 0, 0)
ChineseRest(support) = ChineseRest(0.0, 0.1, UniCat(support))
support(r::ChineseRest) = support(r.basedist)

function sample(r::ChineseRest)
    if flip((r.num_tables * r.a + r.b) / (r.num_customers + r.b))
        sample(r.basedist)
    else
        categorical_sample(
            Dict(dish => sum(r.tables[dish]) for dish in keys(r.tables))
        )
    end
end

function new_table_logscore(r::ChineseRest, dish)
    new_table_prob = LogProb(
        log(r.num_tables * r.a + r.b) - log(r.num_customers + r.b)
    )
    base_dist_prob = logscore(r.basedist, dish)
    new_table_prob * base_dist_prob
end

function logscore(r::ChineseRest, dish)
    if haskey(r.tables, dish)
        numerator = sum(n - r.a for n in r.tables[dish])
        LogProb(log(numerator) - log(r.num_customers + r.b)) + new_table_logscore(r, dish)
    else
        new_table_logscore(r, dish)
    end
end

function add_obs!(r::ChineseRest, dish)
    r.num_customers += 1
    if !haskey(r.tables, dish)
        r.num_tables += 1
        r.tables[dish] = [1]
    else
        r.tables[dish][1] += 1
    end
    nothing
end

function rm_obs!(r::ChineseRest, dish)
    r.num_customers -= 1
    r.tables[dish][1] -= 1
    if r.tables[dish][1] == 0
        r.num_tables -= 1
        delete!(r.tables, dish)
    end
    nothing
end

################################
### Simple Conditional Class ###
################################

mutable struct SimpleCond{C, D, S} # context, distribution, support
    dists   :: Dict{C, D}
    support :: S
end

SimpleCond(dists::AbstractDict) = SimpleCond(
    dists,
    vcat([collect(support(dist)) for dist in values(dists)]...)
)

SimpleCond(g::Base.Generator) = SimpleCond(Dict(g))

sample(sc::SimpleCond, context, args...) = sample(sc.dists[context], args...)
logscore(sc::SimpleCond, obs, context) = logscore(sc.dists[context], obs)
rm_obs!(sc::SimpleCond, obs, context) = rm_obs!(sc.dists[context], obs)

function add_obs!(cond::SimpleCond{C,D,S}, obs, context) where {C,D,S}
    if !haskey(cond.dists, context)
        cond.dists[context] = D(cond.support)
    end
    add_obs!(cond.dists[context], obs)
end

end # module