text:
using CUDA

function D!(DA::AbstractArray{T,3}, A::AbstractArray{T,2},
            dx::Real, dz::Real) where T
    @assert size(DA,1) == size(A,1)
    @assert size(DA,2) == size(A,2)
    @assert size(DA,3) == 2

    @views DA[1:end-1,:,1] .= (A[2:end,:] .- A[1:end-1,:])./dx
    @views DA[:,1:end-1,2] .= (A[:,2:end] .- A[:,1:end-1])./dz
    # Periodic boundary conditions
    @views DA[end,:,1] .= (A[1,:] .- A[end,:])./dx
    @views DA[:,end,2] .= (A[:,1] .- A[:,end])./dz

    DA
end

function D!(DA::AbstractArray{T}, A::AbstractArray{T}) where T
    D!(DA, A, 1, 1)
end

function Dᵀ!(A::AbstractArray{T,2}, DA::AbstractArray{T,3},
             dx::Real, dz::Real) where T
    @assert size(DA,1) == size(A,1)
    @assert size(DA,2) == size(A,2)
    @assert size(DA,3) == 2

    @views A[2:end,:] .= (DA[1:end-1,:,1] .- DA[2:end,:,1])./dx
    @views A[1,:] .= (DA[end,:,1] .- DA[1,:,1])./dx
    @views A[:,2:end] .+= (DA[:,1:end-1,2] .- DA[:,2:end,2])./dz
    @views A[:,1] .+= (DA[:,end,2] .- DA[:,1,2])./dx
    A
end

function Dᵀ!(A::AbstractArray{T}, DA::AbstractArray{T}) where T
    Dᵀ!(A, DA, 1, 1)
end

function proj_G1!(G::AbstractArray{T,3}) where T
    map!(x -> x / max(1, abs(x)), G, G)
    G
end

function proj_G2!(G::AbstractArray{T,3}) where T
    n, m, d = size(G)
    for j in 1:m
        for i in 1:n
            @views nl2 = sqrt.(sum(abs.(G[i,j,:]).^2))
            for k in 1:d
                G[i,j,k] /= max(1, nl2)
            end
        end
    end
    G
end

function proj_G2_kernel(G::CuDeviceArray{T,3}) where T
    n, m, d = size(G)
    ix = (blockIdx().x - 1)*blockDim().x + threadIdx().x
    jx = (blockIdx().y - 1)*blockDim().y + threadIdx().y

    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    @inbounds for j in jx:stride_y:m
        for i in ix:stride_x:n
            nl2 = real(T)(0)
            for k in 1:d
                nl2 += abs(G[i,j,k])^2
            end
            nl2 = sqrt(nl2)
            if nl2 > 1
                for k in 1:d
                    G[i,j,k] /= nl2
                end
            end
        end
    end
    nothing
end

function proj_G2!(G::CuArray{T,3}) where T
    n, m, d = size(G)
    @cuda threads=(16,16) blocks=(Int(ceil(n/16)),Int(ceil(m/16))) (
        proj_G2_kernel(G))
    G
end

struct TV_FGP{T}
    Z
    Xt
    Dt
    Gt
    Gt_tmp
    qt :: Ref{T}
    τ :: T
    γ :: T
    projG! :: Function
    prox! :: Function

    function TV_FGP(Z::AbstractArray{T}, τ::Real,
                    q0::Real=1, projG!::Function=proj_G2!;
                    prox!::Function=identity) where T
        RT = real(T)
        @assert τ > 0
        γ = 1/(12τ)
        Xt = similar(Z)
        Dt = similar(Z, (size(Z)..., ndims(Z)))
        Gt = similar(Dt)
        Gt_tmp = similar(Dt)
        Xt .= Z
        prox!(Xt)
        D!(Dt, Xt)
        Dt .*= γ
        projG!(Dt)
        new{RT}(copy(Z), Xt, Dt, Gt, Gt_tmp, Ref(RT(q0)), RT(τ), RT(γ), projG!, prox!)
    end
end

function TV1_FGP(Z::AbstractArray{T}, τ::Real, q0::Real=1;
                 prox!::Function=identity) where T
    TV_FGP(Z, τ, q0, proj_G1!; prox! = prox!)
end

struct TV_GP{T}
    Z
    Xt
    Gt
    Gt_tmp
    τ :: T
    γ :: T
    projG! :: Function
    prox! :: Function

    function TV_GP(Z::AbstractArray{T}, τ::Real,
                   projG!::Function=proj_G2!;
                   prox!::Function=identity) where T
        RT = real(T)
        @assert τ > 0
        γ = 1/(12τ)
        Xt = similar(Z)
        Gt = similar(Z, (size(Z)..., ndims(Z)))
        Gt_tmp = similar(Gt)
        Xt .= Z
        prox!(Xt)
        D!(Gt_tmp, Xt)
        Gt_tmp .*= γ
        projG!(Gt_tmp)
        new{RT}(copy(Z), Xt, Gt, Gt_tmp, RT(τ), RT(γ), projG!, prox!)
    end
end

function TV1_GP(Z::AbstractArray{T}, τ::Real;
                 prox!::Function=identity) where T
    TV_GP(Z, τ, proj_G1!; prox! = prox!)
end

function TV_step(tv::TV_FGP)
    Dᵀ!(tv.Xt, tv.Dt)
    tv.Xt .= tv.Z .- tv.τ.*tv.Xt
    tv.prox!(tv.Xt)
    D!(tv.Gt, tv.Xt)
    tv.Gt .= tv.Dt .+ tv.γ.*tv.Gt
    tv.projG!(tv.Gt)

    Dᵀ!(tv.Xt, tv.Gt)
    tv.Xt .= tv.Z .- tv.τ.*tv.Xt
    tv.prox!(tv.Xt)

    qt_tmp = tv.qt[]
    tv.qt[] = (1 + sqrt(1+4*tv.qt[]^2))/2
    β = (qt_tmp - 1)/tv.qt[]
    tv.Dt .= tv.Gt .+ β.*(tv.Gt .- tv.Gt_tmp)
    tv.Gt_tmp .= tv.Gt
    tv.Xt
end

function TV_step(tv::TV_GP)
    Dᵀ!(tv.Xt, tv.Gt_tmp)
    tv.Xt .= tv.Z .- tv.τ.*tv.Xt
    tv.prox!(tv.Xt)
    D!(tv.Gt, tv.Xt)
    tv.Gt .= tv.Gt_tmp .+ tv.γ.*tv.Gt
    tv.projG!(tv.Gt)

    Dᵀ!(tv.Xt, tv.Gt)
    tv.Xt .= tv.Z .- tv.τ.*tv.Xt
    tv.prox!(tv.Xt)

    tv.Gt_tmp .= tv.Gt
    tv.Xt
end

function initialize!(tv::TV_FGP, Z::AbstractArray{T}) where T
    tv.qt[] = 1
    tv.Z .= Z
    tv.Xt .= Z
    tv.prox!(tv.Xt)
    D!(tv.Dt, tv.Xt)
    tv.Dt .*= tv.γ
    tv.projG!(tv.Dt)
end

function initialize!(tv::TV_GP, Z::AbstractArray{T}) where T
    tv.Z .= Z
    tv.Xt .= Z
    tv.prox!(tv.Xt)
    D!(tv.Gt_tmp, tv.Xt)
    tv.Gt_tmp .*= tv.γ
    tv.projG!(tv.Gt_tmp)
end

function iter_TV_prox!(tv::Union{TV_GP,TV_FGP},
                       X::AbstractArray{T}, n_iter::Int) where T
    initialize!(tv, X)
    for i in 1:n_iter
        TV_step(tv)
    end
    X .= tv.Xt
end

function get_FTV_prox!(Z::AbstractArray{T}, τ::Real;
                       n_iter::Int=50,
                       prox!::Function=identity) where T
    tv = TV_FGP(Z, τ, prox! = prox!)
    X -> iter_TV_prox!(tv, X, n_iter)
end

function get_FTV1_prox!(Z::AbstractArray{T}, τ::Real;
                       n_iter::Int=50,
                       prox!::Function=identity) where T
    tv = TV1_FGP(Z, τ, prox! = prox!)
    X -> iter_TV_prox!(tv, X, n_iter)
end

function get_TV_prox!(Z::AbstractArray{T}, τ::Real;
                      n_iter::Int=50,
                      prox!::Function=identity) where T
    tv = TV_GP(Z, τ, prox! = prox!)
    X -> iter_TV_prox!(tv, X, n_iter)
end

function get_TV1_prox!(Z::AbstractArray{T}, τ::Real;
                       n_iter::Int=50,
                       prox!::Function=identity) where T
    tv = TV1_GP(Z, τ, prox! = prox!)
    X -> iter_TV_prox!(tv, X, n_iter)
end


text:
push!(LOAD_PATH, pwd() )
using Documenter, BumpFuncs

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = String[],
)

makedocs(
    modules = [BumpFuncs],
    sitename = "BumpFuncs.jl",
    authors = "Joris Pinkse",
    format = format,
    pages = [
        "Home" => "index.md"
    ]
)

deploydocs(
    repo="github.com/NittanyLion/BumpFuncs.jl.git"
#    versions = ["stable" => "v^", "v#.#", devurl => devurl]
)

ext: jl


text:
function egarchfit(data::Array{Float64, 1})
    ## fit egarch function to historic return data
    ##
    ## log(s^2) = kappa + gamma log(s^2) + xi*Z + alpha(|Z| -
    ## sqrt(2/pi))

    nObs = length(data)
    mu = mean(data)
    centeredData = data .- mu

    ## subsample to estimate initial sigma
    nInit = 50
    if nInit > nObs
        nInit = nObs
    end

    ## preallocate values
    volaT = zeros(nObs, 1)
    logVarT = zeros(nObs, 1)

    function egarchNllh(params::Array{Float64, 1}, grad::Vector)
        kappa, gamma, xi, alpha = params;

        initSigma = std(centeredData[1:nInit])
        volaT[1] = initSigma

        ## retrieve volas
        for ii=2:nObs
            lastZ = centeredData[ii-1]/volaT[ii-1]

            logVarT[ii] = kappa + gamma*logVarT[ii-1] +
            xi*lastZ + alpha*(abs(lastZ) - sqrt(2/pi))

            volaT[ii] = sqrt(exp(logVarT[ii]))
        end

        zVals = centeredData ./ volaT
        nllh = nObs * log(sqrt(2*pi)) +
        sum(log(volaT) + 0.5*(centeredData./volaT).^2)

        return nllh
    end

    opt = Opt(:LN_SBPLX, 4)
    min_objective!(opt, egarchNllh)
    (minf, minx, ret) = optimize(opt, [std(centeredData), 0.8, 0.1, 0.1])

end

function getVola(data::Array{Float64, 1}, params::Array{Float64, 1})

    kappa, gamma, xi, alpha = params;
    nObs = length(data)
    mu = mean(data)
    centeredData = data .- mu

    ## subsample to estimate initial sigma
    nInit = 50
    if nInit > nObs
        nInit = nObs
    end

    ## preallocate values
    volaT = zeros(nObs, 1)
    logVarT = zeros(nObs, 1)

    initSigma = std(centeredData[1:nInit])
    volaT[1] = initSigma

    ## retrieve volas
    for ii=2:nObs
        lastZ = centeredData[ii-1]/volaT[ii-1]

        logVarT[ii] = kappa + gamma*logVarT[ii-1] +
        xi*lastZ + alpha*(abs(lastZ) - sqrt(2/pi))

        volaT[ii] = sqrt(exp(logVarT[ii]))
    end

    return volaT
end


text:
# METIS4 dependency for HSL.jl

libmetis4 = library_dependency("libmetis")  # METIS 5 not ok!

metis4_tap = "brewsci/num/brewsci-metis@4"
@static if Sys.isapple()
  provides(Homebrew.HB, metis4_tap, libmetis4, os=:Darwin)
end

provides(Sources,
         URI("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/OLD/metis-4.0.3.tar.gz"),
         libmetis4,
         SHA="5efa35de80703c1b2c4d0de080fafbcf4e0d363a21149a1ad2f96e0144841a55",
         unpacked_dir="metis-4.0.3", os=:Unix)

metis_depsdir = BinDeps.depsdir(libmetis4)
metis_srcdir = joinpath(metis_depsdir, "src", "metis-4.0.3")
metis_prefix = joinpath(metis_depsdir, "usr")

# The METIS Makefile is quite rudimentary and will only build a static library.
provides(SimpleBuild,
         (@build_steps begin
            GetSources(libmetis4)
            (@build_steps begin
               ChangeDirectory(metis_srcdir)
               CreateDirectory(metis_prefix)
               CreateDirectory(joinpath(metis_prefix, "lib"))
               CreateDirectory(joinpath(metis_prefix, "include"))
               (@build_steps begin
                  `make COPTIONS=-fPIC`
                  `cc -fPIC -shared -Wl,$all_load $metis_srcdir/libmetis.a -Wl,$noall_load -o $metis_prefix/lib/libmetis.$so`
                end)
             end)
          end), [libmetis4], os=:Unix)

# try to guess what BinDeps installed
@static if Sys.isunix()
  # Set metis_libpath for all Unix platforms.
  # TODO: update this when there are more providers.
  metis_libdir = joinpath(metis_prefix, "lib")
  @static if Sys.isapple()
    # It's different on OSX because we used Homebrew.
    metis_libdir = joinpath(Homebrew.prefix(metis4_tap), "lib")
  end
  # This is run before metis_libdir is created, so checking if isdir fails
  metis_libpath = "-L$metis_libdir -lmetis"
end


text:
using Test, MPI

using ClimateMachine
using ClimateMachine.MPIStateArrays

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD

@testset "MPIStateArray broadcasting" begin
    let
        localsize = (4, 6, 8)
        A = rand(Float32, localsize)
        B = rand(Float32, localsize)

        QA = MPIStateArray{Float32}(mpicomm, ArrayType, localsize...)
        QB = similar(QA)

        QA .= A
        QB .= B

        @test Array(QA) == A
        @test Array(QB) == B

        QC = QA .+ QB
        @test typeof(QC) == typeof(QA)
        C = Array(QC)
        @test C == A .+ B

        QC = QA .+ sqrt.(QB)
        C = Array(QC)
        @test C ≈ A .+ sqrt.(B)

        QC = QA .+ sqrt.(QB) .* exp.(QA .- QB .^ 2)
        C = Array(QC)
        @test C ≈ A .+ sqrt.(B) .* exp.(A .- B .^ 2)

        # writing to an existing array instead of creating a new one
        fill!(QC, 0)
        QC .= QA .+ sqrt.(QB) .* exp.(QA .- QB .^ 2)
        C = Array(QC)
        @test C ≈ A .+ sqrt.(B) .* exp.(A .- B .^ 2)
    end

    let
        numelems = 12
        realelems = 1:7
        ghostelems = 8:12

        QA = MPIStateArray{Int}(
            mpicomm,
            ArrayType,
            1,
            1,
            numelems,
            realelems = realelems,
            ghostelems = ghostelems,
        )
        QB = similar(QA)

        fill!(QA, 1)
        fill!(QB, 3)

        QB .= QA .+ QB

        @test all(Array(QB)[realelems] .== 4)
        @test all(Array(QB)[ghostelems] .== 3)
    end
end


text:
abstract type AocDay end
abstract type InputDay end


text:
using Navigation
using Printf
using Roots
# const GridColor = RGB(107.0/255.0, 191.0/255.0, 91.0/255.0)
# const CoastColor = RGB(117.0/255.0, 207.0/255.0, 224.0/255.0)
# const SeaColor = RGB(205.0/255.0, 235.0/255.0,240.0/255.0)
# const LandColor = RGB(1.0, 1.0, 1.0)
# using Printf

struct MapSettings
    map_filename::String
    lake_filename::String
    projection::Union{Equirectangular, Mercator, Lambert, Gnomonic}
    lonticks::StepRangeLen
    latticks::StepRangeLen
    gridcolor::RGB{Float64}
    coastcolor::RGB{Float64}
    seacolor::RGB{Float64}
    landcolor::RGB{Float64}
    gridlinewidth::Float64
    gridalpha::Float64
    function MapSettings(;map_filename::String="data\\GSHHS_c_L1.shp",
        lake_filename::String="data\\GSHHS_c_L2.shp",
        projection::Union{Equirectangular, Mercator, Lambert, Gnomonic}=
        Equirectangular(), lonticks::StepRangeLen=-180.0:10.0:180.0,
        latticks::StepRangeLen=-85.0:2.5:85.0,
        gridcolor::RGB{Float64} = RGB(107.0/255.0, 191.0/255.0, 91.0/255.0),
        coastcolor::RGB{Float64} = RGB(117.0/255.0, 207.0/255.0, 224.0/255.0),
        seacolor::RGB{Float64} = RGB(205.0/255.0, 235.0/255.0,240.0/255.0),
        landcolor::RGB{Float64} = RGB(1.0, 1.0, 1.0),
        gridlinewidth::Float64 = 0.01, gridalpha::Float64 = 0.33)
        new(map_filename, lake_filename, projection, lonticks, latticks, gridcolor,
        coastcolor, seacolor, landcolor, gridlinewidth, gridalpha)
    end
end

# function add_map!(shp::String, projparam=Equirectangular(); kwargs...)
#     add_polygon!(shp, projparam; kwargs...)
# end

function add_map!(shapefile::String, projparam=Equirectangular(); kwargs...)
    shapes = get_shapefile(shapefile)
    add_polygon!(shapes, projparam; kwargs...)
end

add_map!(mapsettings::MapSettings) = add_map!(mapsettings.map_filename,
mapsettings.projection; linecolor=mapsettings.coastcolor,
fillcolor=mapsettings.landcolor, background_color = mapsettings.seacolor,
fillalpha=1.0, fillrange=0)

add_lakes!(mapsettings::MapSettings) = add_map!(
mapsettings.lake_filename, mapsettings.projection;
linecolor=mapsettings.coastcolor, fillcolor=mapsettings.seacolor,
background_color = mapsettings.seacolor, fillalpha=1.0, fillrange=0)

function add_grid!(mapsettings::MapSettings)
    plot!(mapsettings.lonticks, mapsettings.latticks, mapsettings.projection;
    linewidth=mapsettings.gridlinewidth, linecolor=mapsettings.gridcolor,
    linealpha = mapsettings.gridalpha)
end

function boundaries_map(mapsettings::MapSettings;
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))

    if mapsettings.projection isa Gnomonic
        left = projection(lonlims[1], (latlims[1] + latlims[2])/2.0, mapsettings.projection)[1]
        right = projection(lonlims[2], (latlims[1] + latlims[2])/2.0, mapsettings.projection)[1]
        top = projection((lonlims[1] + lonlims[2])/2.0, latlims[2], mapsettings.projection)[2]
        bottom = projection((lonlims[1] + lonlims[2])/2.0, latlims[1], mapsettings.projection)[2]
    else
        left = projection(lonlims[1], latlims[1], mapsettings.projection)[1]
        right = projection(lonlims[2], latlims[1], mapsettings.projection)[1]
        top = projection(lonlims[1], latlims[2], mapsettings.projection)[2]
        bottom = projection(lonlims[1], latlims[1], mapsettings.projection)[2]
    end
    left, right, bottom, top
end

function lims_map!(mapsettings::MapSettings;
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))
    left, right, bottom, top = boundaries_map(mapsettings, lonlims=lonlims,
    latlims=latlims)
    plot!(xlims=(left, right), ylims=(bottom, top))
end

function axis_map!(mapsettings::MapSettings;
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))
    left, right, bottom, top = boundaries_map(mapsettings, lonlims=lonlims,
    latlims=latlims)
    x = [left, right, right, left, left]
    y = [bottom, bottom, top, top, bottom]
    plot!(x,y, linewidth=1.0, linecolor=:black)
end

function ticks_map!(mapsettings::MapSettings;
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))
    if mapsettings.projection isa Gnomonic
        return ticks_map_gnomonic!(mapsettings, lonlims, latlims)
    else
        return ticks_map_straight!(mapsettings, lonlims, latlims)
    end
end

function ticks_map_straight!(mapsettings::MapSettings,
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))

    x = Vector{Float64}()
    lontext = Vector{String}()
    for lontick in mapsettings.lonticks
        if lonlims[1] ≤ lontick ≤ lonlims[2]
            loc = projection(lontick, 0.0, mapsettings.projection)[1]
            append!(x, loc)
            lontex = Printf.@sprintf("%4.0f°", lontick)
            lontext = vcat(lontext, lontex)
        end
    end

    y = Vector{Float64}()
    lattext = Vector{String}()
    for lattick in mapsettings.latticks
        if latlims[1] ≤ lattick ≤ latlims[2]
            loc = projection(0.0, lattick, mapsettings.projection)[2]
            append!(y, loc)
            lattex = Printf.@sprintf("%5.1f°", lattick)
            lattext = vcat(lattext, lattex)
        end
    end

    plot!(showaxis=true, xticks = (x, lontext),
    yticks = (y, lattext), ticks = true)
end

function ticks_map_gnomonic!(mapsettings::MapSettings,
    lonlims::Tuple{Float64, Float64} = (-180.0, 180.0),
    latlims::Tuple{Float64, Float64} = (-85.0, 85.0))

    left, right, bottom, top = boundaries_map(mapsettings, lonlims=lonlims,
    latlims=latlims)

    #xtick Gnomonic
    λ₃, ϕ₃ = projection_inv(left, bottom, mapsettings.projection)
    λ₄, ϕ₄ = projection_inv(right, bottom, mapsettings.projection)
    λ₅, ϕ₅ = projection_inv(left, top, mapsettings.projection)
    pos₃ = Navigation.Point_deg(ϕ₃, λ₃)
    pos₄ = Navigation.Point_deg(ϕ₄, λ₄)
    pos₅ = Navigation.Point_deg(ϕ₅, λ₅)
    bearing₃₄ = Navigation.bearing(pos₃, pos₄)

    x = Vector{Float64}()
    t = Vector{String}()
    for lontick in mapsettings.lonticks
        pos₁ = Navigation.Point_deg(latlims[2], lontick)
        bearing₁₂ = 180.0
        pos_tick = Navigation.intersection_point(pos₁, pos₃, bearing₁₂, bearing₃₄)
        if !isinf(pos_tick.ϕ)
            loc = projection(lontick, pos_tick.ϕ, mapsettings.projection)[1]
            txt = Printf.@sprintf("%4.0f°", lontick)
            append!(x, loc)
            t = vcat(t, txt)
        end
    end
    plot!(showaxis=true, xticks = (x,t), ticks = true)

    #ytick Gnomonic
    y = Vector{Float64}()
    t = Vector{String}()

    for lattick in mapsettings.latticks
        if ϕ₃ ≤ lattick ≤ ϕ₅
            f(yloc) = projection_inv(left, yloc, mapsettings.projection)[2] - lattick
            loc = find_zero(f, (bottom, top))
            txt = Printf.@sprintf("%5.1f°", lattick)
            append!(y, loc)
            t = vcat(t, txt)
        end

    end
    plot!(yticks = (y,t))
end


text:
platform_not_supported() = error("platform not yet supported. File Pull Request or issue at Toxcore.jl")

zip_eding = @windows? ".zip" : @unix? ".tar.xz" : platform_not_supported()

const base 				= "https://jenkins.libtoxcore.so/"
const lastbuild 		= "lastSuccessfulBuild/artifact/"
const name_unix 		= "libtoxcore$zip_eding"
base_win(wordsize) 		= "view/Libs/job/toxcore_win$(wordsize)_dll/"
name_win(wordsize) 		= "tox_win$(wordsize)_dll$zip_eding"
base_linux(wordsize) 	= "job/libtoxcore-linux-$(wordsize == 32 ? "i686" : "x86_64")/"
const base_macos 		= "job/libtoxcore-osx/"

@windows? begin
	const liburl = base*base_win(WORD_SIZE)*lastbuild*name_win(WORD_SIZE)
end : @linux? begin
	const liburl = base*base_linux(WORD_SIZE)*lastbuild*name_unix
end : @macos? begin
	const liburl = base*base_macos*lastbuild*name_unix
end : platform_not_supported()

println(base*base_linux(WORD_SIZE)*lastbuild*name_unix)

path = download(liburl)

!(isdir("lib") || isfile("lib")) && mkdir("lib")
!(isdir("builds") || isfile("builds")) && mkdir("builds")

lib_archive = Pkg.dir("Toxcore", "deps", "lib", "lib$zip_eding")
isfile(lib_archive) && rm(lib_archive)
mv(path, lib_archive)

const unzip = @windows? `7z x -obuilds -y $lib_archive` : @unix? `tar xf $lib_archive -C builds` : platform_not_supported()

run(unzip)

!(isdir("bin") || isfile("bin")) && mkdir("bin")
const file = @windows? "libtox" : @linux? "libtoxcore" : @macos? "NO_IDEA" : error("platform not support")
const ending = @windows? ".dll" : @linux? ".a" : @macos? ".dylib" : error("platform not support")
const folder = @windows? "bin" : @linux? "lib" : error("platform not support")

inzip_location = Pkg.dir("Toxcore", "deps", "builds", folder, "$file$ending")
# platfrom specific extraction
mv(inzip_location, joinpath("bin", "libtoxcore$ending"))


text:
println("Hello Word!")

function words()
    println("Nine")
end


text:
export error_if_empty

import Base: show

"""
    error_if_empty(err)

Creates a `error_if_empty` operator, which emits a given error if the source Observable completes
without emitting any next value, otherwise mirrors the source Observable.

```jldoctest
using Rocket

source = completed(Int) |> error_if_empty("Empty")

subscribe!(source, logger())
;

# output
[LogActor] Error: Empty
```
See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`logger`](@ref), [`map`](@ref)
"""
error_if_empty(err) = ErrorIfEmptyOperator(err)

struct ErrorIfEmptyOperator <: InferableOperator
    err
end

operator_right(operator::ErrorIfEmptyOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::ErrorIfEmptyOperator, source) where L
    return proxy(L, source, ErrorIfEmptyProxy(operator.err))
end

struct ErrorIfEmptyProxy <: ActorProxy
    err
end

actor_proxy!(::Type{L}, proxy::ErrorIfEmptyProxy, actor::A) where { L, A } = ErrorIfEmptyActor{L, A}(actor, false, proxy.err)

mutable struct ErrorIfEmptyActor{L, A} <: Actor{L}
    actor      :: A
    is_emitted :: Bool
    err
end

function on_next!(actor::ErrorIfEmptyActor{L}, data::L) where L
    actor.is_emitted = true
    next!(actor.actor, data)
end

function on_error!(actor::ErrorIfEmptyActor, err)
    actor.is_emitted = true
    error!(actor.actor, err)
end

function on_complete!(actor::ErrorIfEmptyActor)
    if !actor.is_emitted
        error!(actor.actor, actor.err)
    else
        complete!(actor.actor)
    end
end

Base.show(io::IO, ::ErrorIfEmptyOperator)         = print(io, "ErrorIfEmptyOperator()")
Base.show(io::IO, ::ErrorIfEmptyProxy)            = print(io, "ErrorIfEmptyProxy()")
Base.show(io::IO, ::ErrorIfEmptyActor{L}) where L = print(io, "ErrorIfEmptyActor($L)")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

Makie.plottype(::CartesianGrid) = Viz{<:Tuple{CartesianGrid}}

function Makie.plot!(plot::Viz{<:Tuple{CartesianGrid}})
  # retrieve grid object
  grid = plot[:object][]

  color        = plot[:color][]
  alpha        = plot[:alpha][]
  colorscheme  = plot[:colorscheme][]
  facetcolor   = plot[:facetcolor][]
  showfacets   = plot[:showfacets][]

  # process color spec into colorant
  colorant = process(color, colorscheme, alpha)

  # relevant settings
  nd = embeddim(grid)
  or = coordinates(minimum(grid))
  sp = spacing(grid)
  sz = size(grid)

  if color isa AbstractVector
    # create a full heatmap or volume
    xyz = cartesiancenters(or, sp, sz, nd)
    if nd == 2
      xs, ys = xyz
      xs′ = xs .- sp[1] / 2
      ys′ = ys .- sp[2] / 2
      C   = reshape(colorant, sz)
      Makie.heatmap!(plot, xs′, ys′, C)
    elseif nd == 3
      xs, ys, zs = xyz
      coords = [(x,y,z) for x in xs for y in ys for z in zs]
      Makie.meshscatter!(plot, coords,
        marker = Makie.Rect3(-sp, sp),
        markersize = 1,
        color = colorant,
      )
    end
  else
    # create the smallest mesh of simplices
    mesh = cartesianmesh(or, sp, sz, nd)
    viz!(plot, mesh,
      color = color,
      alpha = alpha,
      showfacets = false
    )

    if showfacets
      # create a minimum number of segments
      xyz  = cartesiansegments(or, sp, sz, nd)
      Makie.lines!(plot, xyz...,
        color = facetcolor
      )
    end
  end
end

# helper function to create the smallest mesh
# of simplices covering the Cartesian grid
function cartesianmesh(or, sp, sz, nd)
  if nd == 2
    A = Point2(or) + Vec2(0, 0)
    B = Point2(or) + Vec2(sp[1]*sz[1], 0)
    C = Point2(or) + Vec2(sp[1]*sz[1], sp[2]*sz[2])
    D = Point2(or) + Vec2(0, sp[2]*sz[2])
    points = [A, B, C, D]
    connec = connect.([(1,2,3),(3,4,1)])
    SimpleMesh(points, connec)
  elseif nd == 3
    A = Point3(or) + Vec3(0, 0, 0)
    B = Point3(or) + Vec3(sp[1]*sz[1], 0, 0)
    C = Point3(or) + Vec3(sp[1]*sz[1], sp[2]*sz[2], 0)
    D = Point3(or) + Vec3(0, sp[2]*sz[2], 0)
    E = Point3(or) + Vec3(0, 0, sp[3]*sz[3])
    F = Point3(or) + Vec3(sp[1]*sz[1], 0, sp[3]*sz[3])
    G = Point3(or) + Vec3(sp[1]*sz[1], sp[2]*sz[2], sp[3]*sz[3])
    H = Point3(or) + Vec3(0, sp[2]*sz[2], sp[3]*sz[3])
    points = [A, B, C, D, E, F, G, H]
    connec = connect.([(4,3,2,1),(5,6,7,8),(6,5,1,2),
                       (3,7,6,2),(4,8,7,3),(1,5,8,4)])
    SimpleMesh(points, connec)
  else
    throw(ErrorException("not implemented"))
  end
end

# helper function to create a minimum number
# of line segments within Cartesian grid
function cartesiansegments(or, sp, sz, nd)
  if nd == 2
    xs = range(or[1], step=sp[1], length=sz[1]+1)
    ys = range(or[2], step=sp[2], length=sz[2]+1)
    coords = []
    for x in xs
      push!(coords, (x, first(ys)))
      push!(coords, (x, last(ys)))
      push!(coords, (NaN, NaN))
    end
    for y in ys
      push!(coords, (first(xs), y))
      push!(coords, (last(xs), y))
      push!(coords, (NaN, NaN))
    end
    x = getindex.(coords, 1)
    y = getindex.(coords, 2)
    x, y
  elseif nd == 3
    xs = range(or[1], step=sp[1], length=sz[1]+1)
    ys = range(or[2], step=sp[2], length=sz[2]+1)
    zs = range(or[3], step=sp[3], length=sz[3]+1)
    coords = []
    for y in ys, z in zs
      push!(coords, (first(xs), y, z))
      push!(coords, (last(xs), y, z))
      push!(coords, (NaN, NaN, NaN))
    end
    for x in xs, z in zs
      push!(coords, (x, first(ys), z))
      push!(coords, (x, last(ys), z))
      push!(coords, (NaN, NaN, NaN))
    end
    for x in xs, y in ys
      push!(coords, (x, y, first(zs)))
      push!(coords, (x, y, last(zs)))
      push!(coords, (NaN, NaN, NaN))
    end
    x = getindex.(coords, 1)
    y = getindex.(coords, 2)
    z = getindex.(coords, 3)
    x, y, z
  else
    throw(ErrorException("not implemented"))
  end
end

# helper function to create the center of
# the elements of the Cartesian grid
function cartesiancenters(or, sp, sz, nd)
  if nd == 2
    xs = range(or[1]+sp[1], step=sp[1], length=sz[1])
    ys = range(or[2]+sp[2], step=sp[2], length=sz[2])
    xs, ys
  elseif nd == 3
    xs = range(or[1]+sp[1], step=sp[1], length=sz[1])
    ys = range(or[2]+sp[2], step=sp[2], length=sz[2])
    zs = range(or[3]+sp[3], step=sp[3], length=sz[3])
    xs, ys, zs
  else
    throw(ErrorException("not implemented"))
  end
end


text:
fp = dirname(@__FILE__)
@testset "Impulse responses of a VAR using a DSGE as a prior (wrapper function)" begin
    jlddata = load(joinpath(fp, "../../reference/test_dsgevar_lambda_irfs.jld2"))
    m = Model1002("ss10", custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    m <= Setting(:impulse_response_horizons, 10)
    dsgevar = DSGEVAR(m, collect(keys(m.exogenous_shocks)), "ss11")
    DSGE.update!(dsgevar, λ = 1.)
    DSGE.update!(dsgevar, jlddata["modal_param"])

    Random.seed!(1793)
    out = impulse_responses(dsgevar, jlddata["modal_param"],
                            jlddata["data"],
                            :mode, :cholesky; parallel = false,
                            create_meansbands = false, flip_shocks = false,
                            n_obs_shock = 1)
    out_lr = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                               :mode, :choleskyLR; parallel = false,
                               create_meansbands = false, flip_shocks = false,
                               n_obs_shock = 1)
    out_maxbc = impulse_responses(dsgevar,  jlddata["modal_param"], jlddata["data"],
                                  :mode, :maxBC; parallel = false,
                                  create_meansbands = false, flip_shocks = false,
                                  n_obs_shock = 1)

    Random.seed!(1793)
    _ = impulse_responses(dsgevar, jlddata["modal_param"],
                          jlddata["data"],
                          :mode, :cholesky; parallel = false,
                          create_meansbands = false, flip_shocks = false,
                          n_obs_shock = 1)

    out_lr2 = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                :mode, :cholesky_long_run; parallel = false,
                                create_meansbands = false, flip_shocks = false,
                                n_obs_shock = 1)
    out_maxbc2 = impulse_responses(dsgevar,  jlddata["modal_param"], jlddata["data"],
                                   :mode, :maximum_business_cycle_variance; parallel = false,
                                   create_meansbands = false, flip_shocks = false,
                                   n_obs_shock = 1)

    Random.seed!(1793)
    out_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                 :mode, :cholesky; parallel = false,
                                 create_meansbands = false, flip_shocks = true,
                                 n_obs_shock = 1)
    out_lr_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                    :mode, :choleskyLR; parallel = false,
                                    create_meansbands = false, flip_shocks = true,
                                    n_obs_shock = 1)
    out_maxbc_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                       :mode, :maxBC; parallel = false,
                                       create_meansbands = false, flip_shocks = true,
                                       n_obs_shock = 1)

    # Not testing but checking no errors are caused when creating a MeansBands
    mb = impulse_responses(dsgevar, reshape(jlddata["modal_param"], 1,
                                            length(jlddata["modal_param"])),
                            jlddata["data"],
                            :mode, :cholesky; parallel = false,
                            create_meansbands = true, test_meansbands = true,
                            flip_shocks = false, n_obs_shock = 1)
    mb = impulse_responses(dsgevar, reshape(jlddata["modal_param"], 1,
                                            length(jlddata["modal_param"])),
                            jlddata["data"],
                            :mode, :choleskyLR; parallel = false,
                            create_meansbands = true, test_meansbands = true,
                            flip_shocks = false, n_obs_shock = 1)
    mb = impulse_responses(dsgevar, reshape(jlddata["modal_param"], 1,
                                            length(jlddata["modal_param"])),
                            jlddata["data"],
                            :mode, :maxBC; parallel = false,
                            create_meansbands = true, test_meansbands = true,
                            flip_shocks = false, n_obs_shock = 1)

    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] dropdims(out, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] dropdims(out_lr, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] dropdims(out_maxbc, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] dropdims(out_lr2, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] dropdims(out_maxbc2, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] -dropdims(out_flip, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] -dropdims(out_lr_flip, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] -dropdims(out_maxbc_flip, dims = 3)
end

@testset "Impulse responses of a VAR using parallel (1 worker) and using a DSGE as a prior (wrapper function)" begin
    jlddata = load(joinpath(fp, "../../reference/test_dsgevar_lambda_irfs.jld2"))
    m = Model1002("ss10", custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    m <= Setting(:impulse_response_horizons, 10)
    dsgevar = DSGEVAR(m, collect(keys(m.exogenous_shocks)), "ss11")
    DSGE.update!(dsgevar, λ = 1.)
    DSGE.update!(dsgevar, jlddata["modal_param"])

    Random.seed!(1793)
    out = impulse_responses(dsgevar, jlddata["modal_param"],
                            jlddata["data"],
                            :mode, :cholesky; parallel = true,
                            create_meansbands = false, flip_shocks = false,
                            n_obs_shock = 1)
    out_lr = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                               :mode, :choleskyLR; parallel = true,
                               create_meansbands = false, flip_shocks = false,
                               n_obs_shock = 1)
    out_maxbc = impulse_responses(dsgevar,  jlddata["modal_param"], jlddata["data"],
                                  :mode, :maxBC; parallel = true,
                                  create_meansbands = false, flip_shocks = false,
                                  n_obs_shock = 1)

    Random.seed!(1793)
    out_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                 :mode, :cholesky; parallel = true,
                                 create_meansbands = false, flip_shocks = true,
                                 n_obs_shock = 1)
    out_lr_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                    :mode, :choleskyLR; parallel = true,
                                    create_meansbands = false, flip_shocks = true,
                                    n_obs_shock = 1)
    out_maxbc_flip = impulse_responses(dsgevar, jlddata["modal_param"], jlddata["data"],
                                       :mode, :maxBC; parallel = true,
                                       create_meansbands = false, flip_shocks = true,
                                       n_obs_shock = 1)

    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] dropdims(out, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] dropdims(out_lr, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] dropdims(out_maxbc, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_cholesky_irf"] -dropdims(out_flip, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_choleskyLR_irf"] -dropdims(out_lr_flip, dims = 3)
    @test @test_matrix_approx_eq jlddata["exp_modal_maxBC_irf"] -dropdims(out_maxbc_flip, dims = 3)
end

@testset "Impulse responses of a VAR using a DSGE as a prior (wrapper function)" begin
    jlddata = load(joinpath(fp, "../../reference/test_dsgevar_lambda_irfs.jld2"))
    m = Model1002("ss10", custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true)))
    m <= Setting(:impulse_response_horizons, 10)
    dsgevar = DSGEVAR(m, collect(keys(m.exogenous_shocks)), "ss11")
    DSGE.update!(dsgevar, λ = 1.)
    DSGE.update!(dsgevar, jlddata["modal_param"])
    Random.seed!(1793)
    params = convert(Matrix{Float64},
                     hcat(jlddata["modal_param"],
                          jlddata["modal_param"] + vcat(randn(3) ./ 50, zeros(92)))')
    Random.seed!(1793)
    out = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation, normalize_rotation = false)
    Random.seed!(1793)
    out_draw = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation, draw_shocks = true, normalize_rotation = false)
    Random.seed!(1793)
    out_dev = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation,
                                deviations = true, normalize_rotation = false)
    Random.seed!(1793)
    out_dev_draw = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation, draw_shocks = true,
                                     deviations = true, normalize_rotation = false)

    Random.seed!(1793)
    out_parallel = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation, parallel = true, normalize_rotation = false)
    Random.seed!(1793)
    out_draw_parallel = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation,
                                          draw_shocks = true, parallel = true, normalize_rotation = false)
    Random.seed!(1793)
    out_dev_parallel = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation,
                                         deviations = true, parallel = true, normalize_rotation = false)
    Random.seed!(1793)
    out_dev_draw_parallel = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation, draw_shocks = true,
                                              deviations = true, parallel = true, normalize_rotation = false)

    # Not testing but just checking no error when creating MeansBands
    mb = impulse_responses(dsgevar, params, jlddata["data"], :full, :rotation,
                           create_meansbands = true, test_meansbands = true, normalize_rotation = false)

    @test @test_matrix_approx_eq out out_parallel
    @test @test_matrix_approx_eq out_dev out_dev_parallel
    @test @test_matrix_approx_eq out[:, :, :, 1] jlddata["rotation_irf_by_shock"]
    @test @test_matrix_approx_eq out_dev[:, :, :, 1] jlddata["deviations_rotation_irf_by_shock"]
    @test @test_matrix_approx_eq out_draw out_draw_parallel
    @test @test_matrix_approx_eq out_dev_draw out_dev_draw_parallel
    @test @test_matrix_approx_eq out_draw[:, :, 1] jlddata["rotation_irf_draw_shock"]
    @test @test_matrix_approx_eq out_dev_draw[:, :, 1] jlddata["deviations_rotation_irf_draw_shock"]
    @test size(out) == (4, 10, 24, 2) # 4 observables, horizon is 10, 24 shocks, 2 parameter draws
    @test size(out_draw) == (4, 10, 2) # 4 observables, horizon is 10, 24 shocks, 2 parameter draws
    @test size(out_dev) == (4, 10, 24, 2) # 4 observables, horizon is 10, 24 shocks, 2 parameter draws
    @test size(out_dev_draw) == (4, 10, 2) # 4 observables, horizon is 10, 24 shocks, 2 parameter draws
end


text:
using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic

using Oceananigans.Architectures: AbstractArchitecture, CPU
using Oceananigans.Advection: CenteredSecondOrder
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions,
                                       VVelocityBoundaryConditions,
                                       TracerBoundaryConditions

using Oceananigans.Fields: Field, tracernames, TracerFields, XFaceField, YFaceField, CenterField
using Oceananigans.Forcings: model_forcing
using Oceananigans.Grids: with_halo, topology, inflate_halo_size, halo_size
using Oceananigans.TimeSteppers: Clock, TimeStepper
using Oceananigans.TurbulenceClosures: with_tracers, DiffusivityFields
using Oceananigans.Utils: tupleit

function ShallowWaterTendencyFields(arch, grid, tracer_names)

    uh = XFaceField(arch, grid, UVelocityBoundaryConditions(grid))
    vh = YFaceField(arch, grid, VVelocityBoundaryConditions(grid))
    h  = CenterField(arch,  grid, TracerBoundaryConditions(grid))
    tracers = TracerFields(tracer_names, arch, grid)

    return merge((uh=uh, vh=vh, h=h), tracers)
end

function ShallowWaterSolutionFields(arch, grid, bcs)

    uh_bcs = :uh ∈ keys(bcs) ? bcs.uh : UVelocityBoundaryConditions(grid)
    vh_bcs = :vh ∈ keys(bcs) ? bcs.vh : VVelocityBoundaryConditions(grid)
    h_bcs  = :h  ∈ keys(bcs) ? bcs.h  : TracerBoundaryConditions(grid)

    uh = XFaceField(arch, grid, uh_bcs)
    vh = YFaceField(arch, grid, vh_bcs)
    h = CenterField(arch, grid, h_bcs)

    return (uh=uh, vh=vh, h=h)
end

struct ShallowWaterModel{G, A<:AbstractArchitecture, T, V, R, F, E, B, Q, C, K, TS} <: AbstractModel{TS}

                          grid :: G         # Grid of physical points on which `Model` is solved
                  architecture :: A         # Computer `Architecture` on which `Model` is run
                         clock :: Clock{T}  # Tracks iteration number and simulation time of `Model`
    gravitational_acceleration :: T         # Gravitational acceleration, full, or reduced
                     advection :: V         # Advection scheme for velocities _and_ tracers
                      coriolis :: R         # Set of parameters for the background rotation rate of `Model`
                       forcing :: F         # Container for forcing functions defined by the user
                       closure :: E         # Diffusive 'turbulence closure' for all model fields
                    bathymetry :: B         # Bathymetry/Topography for the model
                      solution :: Q         # Container for transports `uh`, `vh`, and height `h`
                       tracers :: C         # Container for tracer fields
                 diffusivities :: K         # Container for turbulent diffusivities
                   timestepper :: TS        # Object containing timestepper fields and parameters

end

function ShallowWaterModel(;
                           grid,
                           gravitational_acceleration,
  architecture::AbstractArchitecture = CPU(),
                               clock = Clock{eltype(grid)}(0, 0, 1),
                           advection = UpwindBiasedFifthOrder(),
                            coriolis = nothing,
                 forcing::NamedTuple = NamedTuple(),
                             closure = nothing,
                          bathymetry = nothing,
                            solution = nothing,
                             tracers = (),
                       diffusivities = nothing,
     boundary_conditions::NamedTuple = NamedTuple(),
                 timestepper::Symbol = :RungeKutta3)

    tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example)

    Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, grid.Hz, topology(grid), advection, closure)
    grid = with_halo((Hx, Hy, Hz), grid)

    model_field_names = (:uh, :vh, :h, tracers...)
    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, model_field_names)

    solution = ShallowWaterSolutionFields(architecture, grid, boundary_conditions)
    tracers  = TracerFields(tracers, architecture, grid, boundary_conditions)
    diffusivities = DiffusivityFields(diffusivities, architecture, grid,
                                      tracernames(tracers), boundary_conditions, closure)

    # Instantiate timestepper if not already instantiated
    timestepper = TimeStepper(timestepper, architecture, grid, tracernames(tracers);
                              Gⁿ = ShallowWaterTendencyFields(architecture, grid, tracernames(tracers)),
                              G⁻ = ShallowWaterTendencyFields(architecture, grid, tracernames(tracers)))

    # Regularize forcing and closure for model tracer and velocity fields.
    model_fields = merge(solution, tracers)
    forcing = model_forcing(model_fields; forcing...)
    closure = with_tracers(tracernames(tracers), closure)

    return ShallowWaterModel(grid,
                             architecture,
                             clock,
                             eltype(grid)(gravitational_acceleration),
                             advection,
                             coriolis,
                             forcing,
                             closure,
                             bathymetry,
                             solution,
                             tracers,
                             diffusivities,
                             timestepper)
end


text:
using ITensors, Test

@testset "Basic DMRG" begin
  @testset "Spin-one Heisenberg" begin
    N = 100
    sites = spinOneSites(N)

    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    H = toMPO(ampo,sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10,10)
    cutoff!(sweeps,1E-11)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy,psi = dmrg(H,psi,sweeps;quiet=true)
    @test energy < -120.0
    # test with SVD too!
    psi = randomMPS(sites)
    energy,psi = dmrg(H,psi,sweeps;quiet=true,which_factorization="svd")
    @test energy < -120.0
  end

  @testset "Transverse field Ising" begin
    N = 32
    sites = spinHalfSites(N)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"Sz",j,"Sz",j+1)
      add!(ampo,-0.5,"Sx",j)
    end
    H = toMPO(ampo,sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    energy,psi = dmrg(H,psi0,sweeps,quiet=true)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 0.25 - 0.25/sin(π/(2*(2*N+1)))
    @test abs(energy - energy_exact) < 0.05
  end

  @testset "DMRGObserver" begin
    N = 10
    sites = spinHalfSites(N)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"Sz",j,"Sz",j+1)
      add!(ampo,-0.2,"Sx",j)
    end
    H = toMPO(ampo,sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10)
    cutoff!(sweeps,1E-12)

    observer = DMRGObserver(["Sz","Sx"], sites)

    E,psi = dmrg(H,psi0,sweeps,observer=observer,quiet=true)
    @test length(measurements(observer)["Sz"])==3
    @test length(measurements(observer)["Sx"])==3
    @test all(length.(measurements(observer)["Sz"]) .== N)
    @test all(length.(measurements(observer)["Sx"]) .== N)
    @test length(energies(observer))==3
    @test energies(observer)[end]==E

    orthogonalize!(psi,1)
    m = scalar(dag(psi[1])*noprime(op(sites, "Sz", 1)*psi[1]))
    @test measurements(observer)["Sz"][end][1] ≈ m
  end
end


text:
const ReverseOrdering = Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}

mutable struct CoveringFitter
    covering::Covering
    queue::PriorityQueue{Int, Float64, ReverseOrdering}
    approximated_global_error::Float64
    iteration::Int

    desired_maximal_error::Float64
    vertices_accuracy::Float64
    maximal_iterations::Int

    F::AmoebaFiber2D
    startvalue_generator::TorusStartValueGenerator{2}
    membership_test_options::MembershipTestOptions
end

function CoveringFitter(C::Covering, F::AmoebaFiber2D;
    desired_maximal_error=5e-2,
    vertices_accuracy = desired_maximal_error * 1e-3,
    maximal_iterations = 2000,
    options=MembershipTestOptions())

    covering = deepcopy(C)
    startvalue_generator = TorusStartValueGenerator{2}()
    initial_refinement!(covering, F, startvalue_generator, vertices_accuracy, options)

    queue = PriorityQueue{Int, Float64}(Base.Reverse)
    approximated_global_error = 0.0
    for (id, poly) in polygons(covering)
        loc_err = approximated_local_error(covering, poly, vertices_accuracy)
        approximated_global_error += loc_err
        enqueue!(queue, id => approximated_local_error(covering, poly, vertices_accuracy))
    end

    CoveringFitter(covering, queue, approximated_global_error, 0,
        desired_maximal_error, vertices_accuracy, maximal_iterations,
        F, startvalue_generator, options)
end

function update_error!(M::CoveringFitter, id::Int)
    if id != 0
        old_err = M.queue[id]
        new_err = approximated_local_error(M.covering, id, M.vertices_accuracy)
        M.queue[id] = new_err
        M.approximated_global_error -= (old_err - new_err)
    end
    nothing
end

function add2queue!(M::CoveringFitter, id::Int)
    if id != 0
        err = approximated_local_error(M.covering, id, M.vertices_accuracy)
        M.queue[id] = err
        M.approximated_global_error += err
    end
    nothing
end

"""
    fit!(covering_fitter)

Improve the covering until the `approximated_global_error` is less
than the `desired_maximal_error` or `maximal_iterations` were executed.
"""
function fit!(M::CoveringFitter; report_progress=false)
    k = 0
    while !done(M, k)
        k += 1
        if report_progress && mod(k, 5) == 0
            println("Iteration: $k - Error: $(M.approximated_global_error)")
        end
        next_step!(M)
    end
    k
end

function dequeue_pairs!(M::CoveringFitter, N::Int)
    ids = Int[]
    for k=1:N
        id, err = dequeue_pair!(M.queue)
        M.approximated_global_error -= err
        push!(ids, id)
    end
    ids
end

@inline function next_step!(M::CoveringFitter)
    M.iteration += 1
    id, err = dequeue_pair!(M.queue)
    M.approximated_global_error -= err
    pillar = central_pillar(M, id)
    poly = polygon(M.covering, id)

    id_left, id_right = split!(M.covering, id, pillar)
    update_error!(M, poly.left_id)
    update_error!(M, poly.right_id)
    add2queue!(M, id_left)
    add2queue!(M, id_right)
end

function Base.iterate(M::CoveringFitter, state=nothing)
    if state === nothing
        next_step!(M)
        M, state+1
    elseif done(M, state)
        return nothing
    else
        next_step!(M)
        M, state + 1
    end
end
done(M::CoveringFitter, k) = k > M.maximal_iterations || M.approximated_global_error < M.desired_maximal_error

Base.IteratorSize(::Type{<:CoveringFitter}) = Base.SizeUnknown()
Base.eltype(::Type{<:CoveringFitter}) = M

function initial_refinement!(covering::Covering, F::AbstractAmoebaFiber, startvalue_generator, ε, o::MembershipTestOptions)
    scalepillars!(covering) do origin, dir
        boundary_approximation(F, startvalue_generator, origin, dir, ε, o)
    end
    covering
end

"""
    boundary_approximation(F::AmoebaFiber, startvalue_generator, origin, direction, ε, membershiptest_options)

Starting from `origin` track the amoeba along `direction` with an initial step length.

Returns `t` s.t. `origin` + `t` * `direction` is *not* in the amoeba and the distance
A point is considered in the amoeba if we find a root (withi tolerance `membership_tol`) within `ntries` tries
and maximal `iterations` iterations. If we cannot find a root the step length is halfed.
"""
function boundary_approximation(F, startvalue_generator, origin::SVector{2}, direction::SVector{2}, ε::Float64, o::MembershipTestOptions, s0=0.0, Δs = 0.5)
    s = s0
    tol = ε / norm(direction)
    if tol ≥ 1.0
        return 1.0
    end
    while s ≤ 1.0
        w = origin + (s + Δs) * direction
        if membershiptest(F, w, startvalue_generator, o).successfull
            s += Δs
        else
            Δs *= 0.5
        end
        if Δs < tol
            return s + 2Δs
        end
    end
    return 1.0
end

function central_pillar(M::CoveringFitter, polygon_id::Int)
    poly = polygon(M.covering, polygon_id)
    orig, dir = average_pillar(M.covering, poly)
    pillar_line = (orig, orig + dir)
    l1, l2 = approximation_lines(M.covering, poly, M.vertices_accuracy)

    p1 = lineintersection(pillar_line, l1)
    p2 = lineintersection(pillar_line, l2)

    t1 = (p1 - orig)[1] / dir[1]
    t2 = (p2 - orig)[1] / dir[1]

    t1 = t1 > 1 ? 0.0 : t1
    t2 = t2 > 1 ? 0.0 : t2
    s0 = max(0.0, t1, t2)
    Δs = 0.5 * (1.0 - s0)

    t = boundary_approximation(M.F, M.startvalue_generator, orig, dir, M.vertices_accuracy, M.membership_test_options, s0, Δs)

    orig, t * dir
end

function average_pillar(covering::Covering, poly::CoveringPolygon)
    p1 = pillar(covering, poly.p1_id)
    p2 = pillar(covering, poly.p2_id)
    0.5 .* (p1 .+ p2)
end


text:
# please instantiate ./test/Manifest.toml for this

using Combinatorics, JLD2, DataStructures
import Statistics

# ===============================
# Roll Fuzz Fixtures
# ===============================

# format: [(initial_values, roll_values, expected_medians)]
# window_size == length(initial_values)
# length(roll_values) == length(expected_medians)
# check each expected_median after pushing the corresponding push_value

roll_fixtures = [([1.], [2., 0.], [2., 0.])]

function generate_roll_fixture(initial_values, N)
    roll_values = rand(N)
    cb = CircularBuffer{Float64}(length(initial_values))
    for i in eachindex(initial_values)
        push!(cb, initial_values[i])
    end
    expected_medians = Array{Float64}(undef, N)
    for i in 1:N
        push!(cb, roll_values[i])
        expected_medians[i] = Statistics.median(cb)
    end
    return (initial_values, roll_values, expected_medians)
end

function append_roll_fixtures(window_size, push_N)
    initial_values = rand(window_size)
    push!(roll_fixtures, generate_roll_fixture(initial_values, push_N))
end

append_roll_fixtures(2, 20)
append_roll_fixtures(3, 20)
append_roll_fixtures(9, 5000)
append_roll_fixtures(10, 5000)
append_roll_fixtures(1_000, 10_000)

roll_filename = dirname(@__FILE__) * "/roll.jld2"
@save roll_filename roll_fixtures

max_stars_repo_name: Firionus/FastRunningMedian

max_issues_repo_name: Firionus/FastRunningMedian


text:
module OpenQuantumTools

using DocStringExtensions
using Reexport
using RecipesBase

import SciMLBase:
    ODEProblem,
    EnsembleProblem,
    ODEFunction,
    DiscreteCallback,
    ContinuousCallback,
    u_modified!,
    full_cache,
    solve,
    ODESolution,
    DiffEqArrayOperator,
    CallbackSet,
    terminate!
import DiffEqCallbacks:
    FunctionCallingCallback, PresetTimeCallback, IterativeCallback

@reexport using LinearAlgebra
@reexport using SparseArrays
@reexport using OpenQuantumBase

include("math_util.jl")
include("QControl/callback_lib.jl")

include("QSolver/util.jl")
include("QSolver/closed_system_solvers.jl")
include("QSolver/lindblad_solver.jl")
include("QSolver/stochastic_schrodinger_solver.jl")
include("QSolver/ame_solver.jl")
include("QSolver/redfield_solver.jl")
include("QSolver/ensemble_builder.jl")

include("plot_util/hamiltonian_plot.jl")
include("plot_util/ode_sol.jl")
include("plot_util/bath_plot.jl")

export solve_unitary,
    solve_schrodinger,
    solve_von_neumann,
    solve_redfield,
    solve_ame,
    solve_cgme,
    solve_ule,
    solve_lindblad,
    build_ensembles

export DEFAULT_INITIALIZER
export InstPulseCallback, PositivityCheckCallback

end # end module


text:
# __precompile__(false)
module Circuitscape
using AlgebraicMultigrid
using ArchGDAL
using SimpleWeightedGraphs
using LightGraphs
using IterativeSolvers
using GZip
using Pardiso

using LinearAlgebra
using SparseArrays
using SharedArrays
using Distributed
using DelimitedFiles
using Logging
using Dates
using SuiteSparse

include("consts.jl")
include("config.jl")
include("logging.jl")
include("utils.jl")
include("io.jl")
include("out.jl")
include("core.jl")
include("network/pairwise.jl")
include("raster/pairwise.jl")
include("raster/advanced.jl")
include("network/advanced.jl")
include("raster/onetoall.jl")
include("run.jl")
include("INIBuilder/INIBuilder.jl")

using .INIBuilder
export start

end


text:
module TestCorr
using Test,CorrDCA,DelimitedFiles

datadir = joinpath(dirname(@__FILE__),"data")
filePijtrue = joinpath(datadir,"Pij.dat")
filePitrue = joinpath(datadir,"Pi.dat")
fileal = joinpath(datadir,"pf14short.fasta.gz")
filecov = joinpath(datadir,"covariance.dat")
Z = read_fasta_alignment(fileal,0.8)[1:3,1:100]
Pijtrue = readdlm(filePijtrue)
Pitrue = readdlm(filePitrue)
Cijtrue = readdlm(filecov) # covariance matrix with pc = 0.1
W,Meff = compute_weights(Z,:auto)
Pij,Pi = Pij,Pi=compute_frequencies(Z,W)
Cij,mi = covariance_matrix(Z,W,pc=0.1)
q  = maximum(Z)
N  = size(Z,1)
s = q-1
@test sum(abs2,reshape(Pitrue,s,3)-Pi) < 1e-20
for i in 1:N
    for j in i:N
        idx =  ((i-1)*s+1):i*s
        jdx =  ((j-1)*s+1):j*s
        pmat = Pijtrue[idx, jdx]
        vp = view(Pij,:,:,i,j)
        delta = sum(abs2, pmat-vp)
        #println("i=$i j=$j delta=$delta")
        @test pmat ≈ vp
    end
end

Pijpc,Pipc=add_pseudocount(Pi,Pij,0.2)

@test permutedims(Pijpc,[2,1,4,3]) == Pijpc
for i in 1:N
    for a in 1:s
        @test Pijpc[a,a,i,i] ≈ Pipc[a,i]
    end
end

@test Cijtrue ≈ Cij

printstyled("All TestCorr passed!\n",color=:light_green,bold=true)
end # end module

max_stars_repo_name: pagnani/CorrDCA

max_issues_repo_name: pagnani/CorrDCA

max_forks_repo_name: pagnani/CorrDCA


text:
using ADI
using HCIDatasets: BetaPictoris
using LinearAlgebra
using PSFModels
using StableRNGs
using Statistics
using Test

if get(ENV, "CI", "false") == "true"
   ENV["DATADEPS_ALWAYS_ACCEPT"] = true
end

rng = StableRNG(8799)

@testset "Classic" begin include("classic.jl") end
@testset "PCA" begin include("pca.jl") end
@testset "TPCA" begin include("tpca.jl") end
@testset "NMF" begin include("nmf.jl") end
@testset "GreeDS" begin include("greeds.jl") end

@testset "Metrics" begin
    include("metrics.jl")
    include("contrast.jl")
end

@testset "SDI" begin include("sdi.jl") end
@testset "Framewise" begin include("framewise.jl") end
@testset "Geometries" begin include("geometries.jl") end


text:
using Test
using SafeTestsets

@safetestset "Polynomial basis" begin include("polynomialbasis_test.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
@testset "set operations" begin
    @testset "union" begin
        d1 = UnitDisk()
        d2 = (-.9..0.9)^2
        d3 = ProductDomain(-.5 .. -.1, ChebyshevInterval())
        d4 = (0.0..1.5)
        d5 = [1.0,3.0]

        @test convert(Domain{SVector{2,Float64}}, d3) isa Domain{SVector{2,Float64}}

        u1 = UnitDisk() ∪ (-.9..0.9)^2
        u2 = u1 ∪ d3
        @test dimension(u1) == 2
        @test dimension(u2) == 2
        @test boundingbox(u1) == ChebyshevInterval()^2

        u3 = d3 ∪ u1
        u4 = u1 ∪ u2
        @test SVector(0.,.15) ∈ u3
        @test SVector(0.,.15) ∈ u4
        @test SVector(1.1,.75) ∉ u3
        @test SVector(1.1,.75) ∉ u4

        ũ1 = UnionDomain(d1,d2)
        @test u1 == ũ1
        @test UnionDomain{SVector{2,Float64}}(d1,d2) == ũ1
        ũ1 = UnionDomain((d1,d2))
        @test u1 == ũ1
        ũ2 = UnionDomain([d1,d2])
        @test ũ2 == ũ2
        @test u1 == ũ2
        @test UnionDomain{SVector{2,Float64}}(d1) isa UnionDomain

        # Don't create a union with two identical elements
        @test UnitDisk() ∪ UnitDisk() isa UnitDisk

        # union with non-Domain type that implements domain interface
        u45 = (0.0..1.5) ∪ [1.0,3.0]
        @test u45 isa Domain{Float64}
        @test u45 isa UnionDomain
        @test eltype(component(u45,1)) == Float64
        @test eltype(component(u45,2)) == Float64
        @test 0.2 ∈ u45
        @test 1.2 ∈ u45
        @test 3 ∈ u45
        @test -1.2 ∉ u45
        @test convert(Domain{BigFloat}, u45) isa Domain{BigFloat}

        u45b = (0.0..1.5) ∪ [1,3]
        @test u45b isa Domain{Float64}
        @test component(u45b,2) isa AbstractArray{Float64}
        @test [1,3] ∪ (0.0..1.5) isa Domain{Float64}

        @test issubset([0,1], 0..1)
        @test !issubset([0,1,2], 0..1)
        @test issubset(Set([0,1]), 0..1)
        @test !issubset(Set([0,2]), 0..1)

        @test uniondomain() == EmptySpace{Any}()
        @test uniondomain(0..1) == 0..1
        @test uniondomain(0..1, [0,1]) == 0..1
        @test uniondomain([0,1], 0..1) == 0..1

        @test uniondomain([0,1], [0.0,1.0]) == [0,1]
        @test uniondomain([0,2], [0.0,1.0]) isa UnionDomain

        # larger union expressions
        @test uniondomain(0..1, 1..3, Point(0.4), 2..5, FullSpace(), Point(-0.2)) isa FullSpace
        @test uniondomain(0..1, 1..3, Point(0.4)) == 0..3
        @test uniondomain(0..1, Point(0.4), 1..3) == 0..3
        @test uniondomain(Point(0.4), 0..1, 1..3) == 0..3

        # ordering doesn't matter
        @test UnionDomain(d1,d2) == UnionDomain(d2,d1)

        @test UnionDomain((d1,d2)) == UnionDomain(d1,d2)
        @test UnionDomain(d1) isa UnionDomain
        @test UnionDomain(UnionDomain(d1,d2),d3) == UnionDomain(d3,UnionDomain(d1,d2))

        @test convert(Domain, [0..1,2..3,3..4]) isa UnionDomain{Int}
        @test convert(Domain{Float64}, [0..1,2..3,3..4]) isa UnionDomain{Float64}
        @test convert(Domain, Set([0..1,2..3,3..4])) isa UnionDomain{Int}
        @test convert(Domain{Float64}, Set([0..1,2..3,3..4])) isa UnionDomain{Float64}
        @test convert(Domain, Set([1,2,3])) isa UnionDomain{Int}
        @test 2 ∈ convert(Domain, Set([1,2,3]))
        @test 2 ∈ convert(Domain{Float64}, Set([1,2,3]))
        @test 4 ∉ convert(Domain, Set([1,2,3]))

        @test interior(uniondomain(0..1, 2..3)) == uniondomain(OpenInterval(0,1),OpenInterval(2,3))
        @test closure(uniondomain(OpenInterval(0,1),OpenInterval(2,3))) == uniondomain(0..1, 2..3)

        @test uniondomain(0..1, 2..3) \ (0.5..2.5) == uniondomain(0..0.5, 2.5..3.0)
        @test uniondomain(0..1, 2..3) \ uniondomain(0.5..2.5, -2..(-1)) == uniondomain(0..0.5, 2.5..3.0)

        @test !isempty(u1)
        show(io, textmime, u1)
        @test String(take!(io)) == "UnitDisk() ∪ ((-0.9..0.9) × (-0.9..0.9))"

        # repeated union
        @test ncomponents(uniondomain(UnitBall{Float64}(), UnitInterval(), UnitInterval())) == 2
        @test ncomponents(uniondomain(UnitInterval(), UnitBall{Float64}(), UnitInterval())) == 2
        @test ncomponents(uniondomain(UnitInterval(), UnitInterval(), UnitBall{Float64}())) == 2
    end

    @testset "intersect" begin
        @test IntersectDomain(0..1) == IntersectDomain((0..1))
        @test IntersectDomain{Float64}(0..1) == IntersectDomain(0.0..1.0)
        @test IntersectDomain{Float64}(0..1, 1..2) == IntersectDomain((0..1, 1..2))
        @test intersectdomain(0..1, 1..2) == Point(1)

        @test intersectdomain(0..1, 0.5..1.5) == (0..1) & (0.5..1.5)

        # intersection of productdomains
        i1 = intersectdomain((-.4..0.4)^2, (-.5 .. 0.5) × (-.1.. 0.1))
        @test i1 == productdomain(-0.4..0.4, -0.1..0.1)
        show(io,i1)
        @test String(take!(io)) == "(-0.4..0.4) × (-0.1..0.1)"
        @test intersectdomain(productdomain(UnitDisk(),-1..1), productdomain(-1..1, UnitDisk())) isa IntersectDomain
        i2 = UnitDisk() & (-.4..0.4)^2
        show(io, textmime, i2)
        @test String(take!(io)) == "UnitDisk() ∩ ((-0.4..0.4) × (-0.4..0.4))"
        @test dimension(i1) == 2
        @test dimension(i2) == 2

        i3 = ((-.5 .. 0.5) × (-.1.. 0.1)) & i2
        i4 = i2 & ((-.5 .. 0.5) × (-.1.. 0.1))
        i5 = i3 & i2
        @test SVector(0.,.05) ∈ i3
        @test SVector(0.,.05) ∈ i4
        @test SVector(0.,.75) ∉ i3
        @test SVector(0.,.75) ∉ i4

        d45 = IntersectDomain(0.0..1.5, [1.0,3.0])
        @test d45 isa Domain{Float64}
        @test 1.0 ∈ d45
        @test 1.1 ∉ d45
        @test convert(Domain{BigFloat}, d45) isa Domain{BigFloat}

        @test intersectdomain() == EmptySpace{Any}()
        @test intersectdomain(UnitDisk()) == UnitDisk()

        @test (0..1) ∩ [1.5] isa IntersectDomain{Float64}
        @test [0.5] ∩ (1..2) isa IntersectDomain{Float64}

        @test intersectdomain(0..1, [0,1]) == [0,1]
        @test intersectdomain([0,1], 0..1) == [0,1]
        @test intersectdomain([0,1], [0.0,1.0]) == [0,1]
        @test intersectdomain([0,2], [0.0,1.0]) isa IntersectDomain

        @test IntersectDomain(UnitDisk(), UnitSquare()) == IntersectDomain(UnitSquare(), UnitDisk())

        @test intersectdomain(uniondomain(0..1, 2..3), uniondomain(0..0.5, 3.5..4.5)) == 0..0.5
        @test intersectdomain(uniondomain(0..1, 2..3), 2..4) == 2..3
        @test intersectdomain(2..4, uniondomain(0..1, 2..3)) == 2..3

        # repeated intersection
        @test ncomponents(intersectdomain(UnitBall{Float64}(), UnitInterval(), UnitInterval())) == 2
        @test ncomponents(intersectdomain(UnitInterval(), UnitBall{Float64}(), UnitInterval())) == 2
        @test ncomponents(intersectdomain(UnitInterval(), UnitInterval(), UnitBall{Float64}())) == 2
        # larger intersection expressions
        @test intersectdomain(0..1, 1..3, Point(0.4), 2..5, FullSpace(), Point(-0.2)) isa EmptySpace
        @test intersectdomain(0..1, 1..3, Point(1.0)) == Point(1.0)
        @test intersectdomain(0..1, Point(1.0), 1..3) == Point(1.0)
        @test intersectdomain(Point(1.0), 0..1, 1..3) == Point(1.0)

        @test boundingbox(UnitSphere() ∩ 2UnitBall()) == (-1..1)^3
    end

    @testset "setdiff" begin
        d1 = UnitDisk() \ ProductDomain(-.5..0.5, -.1..0.1)
        @test dimension(d1) == 2
        @test SVector(0.,.74) ∈ d1
        @test SVector(0.,.25) ∈ d1
        @test SVector(0.,-.05) ∉ d1
        @test SVector(0.5, 0.1) ∉ d1
        @test SVector(1.01, 0.1) ∉ d1
        @test approx_in(SVector(1.01, 0.1), d1, 0.1)
        show(io, textmime, d1)
        @test String(take!(io)) == "UnitDisk() \\ ((-0.5..0.5) × (-0.1..0.1))"
        @test setdiff(d1, ProductDomain(-0.5..0.5, -.1..0.1)) == d1

        d2 = SetdiffDomain(0.0..3.0, [1.0, 2.5])
        @test d2 isa Domain{Float64}
        @test d2 isa SetdiffDomain
        @test 0.99 ∈ d2
        @test_throws MethodError approx_in(3.01, d2, 0.1)
        @test 1.0 ∉ d2
        @test convert(Domain{BigFloat}, d2) isa Domain{BigFloat}

        @test (0..1) \ [0.5] isa SetdiffDomain{Float64}
        d3 = [0,5] \ (0..3)
        @test d3 isa SetdiffDomain{Int}
        @test 0 ∉ d3
        @test 5 ∈ d3

        @test setdiff(0..1, 2..3) == setdiffdomain(0..1, 2..3)
        @test setdiff(0..1, 0.5) == setdiffdomain(0..1, 0.5)
        @test setdiff(0.5, 0..1) == setdiffdomain(0.5, 0..1)

        @test setdiff(0..1, EmptySpace()) == 0..1
        @test setdiff(0..1, 0.0..1.0) == EmptySpace()

        @test (0..1)^2 \ UnitCircle() == UnitInterval()^2 \ UnitCircle()
    end

    @testset "arithmetic" begin
        d1 = (0..1)
        d2 = (2..3)
        d = UnionDomain(d1) ∪ UnionDomain(d2)

        @test d .+ 1 == UnionDomain(d1 .+ 1) ∪ (d2 .+ 1)
        @test d .- 1 == UnionDomain(d1 .- 1) ∪ (d2 .- 1)
        @test 2 * d  == UnionDomain(2 * d1)  ∪ (2 * d2)
        @test d * 2 == UnionDomain(d1 * 2) ∪ (d2 * 2)
        @test d / 2 == UnionDomain(d1 / 2) ∪ (d2 / 2)
        @test 2 \ d == UnionDomain(2 \ d1) ∪ (2 \ d2)

        @test infimum(d) == minimum(d) == 0
        @test supremum(d) == maximum(d) == 3
    end

    @testset "different types" begin
        d̃1 = (0..1)
        d1 = (0f0.. 1f0)
        d2 = (2..3)

        @test UnionDomain(d1) ∪ d2 == UnionDomain(d̃1) ∪ d2
    end

    @testset "disk × interval" begin
        d = (0.0..1) × UnitDisk()
        @test SVector(0.1, 0.2, 0.3) ∈ d

        d = UnitDisk() × (0.0..1)
        @test SVector(0.1, 0.2, 0.3) ∈ d
    end

    @testset "Wrapped domain" begin
        d = 0..1.0
        w = DomainSets.WrappedDomain(d)
        @test w isa Domain
        @test 0.5 ∈ w

        d1 = DomainSets.WrappedDomain{Any}([0,5])
        @test 0 ∈ d1
        @test 1 ∉ d1

        @test convert(Domain{Float64}, [0.4,0.5]) isa Domain{Float64}
    end

    @testset "more set operations" begin
        D = UnitInterval()^3
        S = 2 * UnitBall()
        @testset "joint domain" begin
            DS = D ∪ S
            @test SA[0.0, 0.6, 0.0] ∈ DS
            @test SA[0.9, 0.6,-2.5] ∉ DS
        end

        @testset "domain intersection" begin
            DS = D ∩ S
            @test SA[0.1, 0.1, 0.1] ∈ DS
            @test SA[0.1, -0.1, 0.1] ∉ DS
        end
        @testset "domain difference" begin
            DS = D\S
            @test SA[0.1, 0.1, 0.1] ∉ DS

            D1 = 2 * D
            D2 = D * 2
            D3 = D / 2

            @test SA[2., 2., 2.] ∈ D1
            @test SA[0.9, 0.6,-2.5] ∉ D1
            @test SA[2., 2., 2.] ∈ D2
            @test SA[0.9, 0.6,-2.5] ∉ D2
            @test SA[.5, .4, .45] ∈ D3
            @test SA[.3, 0.6,-.2] ∉ D3
        end
    end
end


text:
module TestSamples

import Historic
using Test

using ..Samples: Samples, Default, Custom

function test_default()
    Historic.clear(Default)
    Samples.default_simple()
    Samples.default_withdata(1)
    Samples.default_withdata(2)
    events = Historic.events(Default)
    @test events.name == [:simple, :withdata, :withdata]
    defaultdata = Historic.defaultdata()
    @test events.data == [defaultdata, (; defaultdata..., k = 1), (; defaultdata..., k = 2)]
    @test events.threadid == fill(Threads.threadid(), 3)
end

function test_custom()
    Historic.clear(Custom)
    Samples.custom_simple()
    Samples.custom_withdata(1)
    Samples.custom_withdata(2)
    events = Historic.events(Custom)
    @test events.name == [:simple, :withdata, :withdata]
    defaultdata = (; foo = :bar)
    @test events.data == [defaultdata, (; defaultdata..., k = 1), (; defaultdata..., k = 2)]
    @test events.threadid == fill(Threads.threadid(), 3)
end

end  # module


text:
struct TestRoot end

#=
    Test that adding a new node inherits the lower bound from its parent
=#
@testset "Adding node after evaluating root node" begin
    function BB.get_branching_indices(tree::TestRoot)
        return [1]
    end

    function BB.evaluate_node!(tree::BnBTree{BB.DefaultNode, TestRoot}, node::BB.DefaultNode)
        return 0.0, Inf
    end

    bnb_model = BB.initialize(;
        branch_strategy = BB.MOST_INFEASIBLE(),
        root = TestRoot(),
        sense = :Min
    )
    BB.set_root!(bnb_model, NamedTuple())
    root = BB.get_next_node(bnb_model, BB.BFS())
    lb, ub = BB.evaluate_node!(bnb_model, root)
    BB.set_node_bound!(bnb_model.sense, root, lb, ub)

    # adding a new node with the root as parent
    BB.add_node!(bnb_model, root, NamedTuple())
    last_node_id = bnb_model.num_nodes
    @test bnb_model.nodes[last_node_id].lb == 0.0
end


text:
export simulate, simulate_k, oneturn

include("traps.jl")

"""
    simulate(list, dice, circular)

Returns the cost (number of dice throws) of a
a complete simulation of a Snakes and Ladders game, meaning
it starts at square 1, with a board disposition described
in `list`, a dice strategy (dice to play at each square) in `dice`
and the circularity of the board indicated in `circular`.
"""
function simulate(list::Vector{Int64}, dice::Vector{Int64}, circular::Bool)
    simulate_k(list, dice, circular, 1)
end

"""
    simulate_k(list, dice, circular, k)

Returns the cost (number of dice throws) of
a simulation of a Snakes and Ladders game, starting at square `k`,
with a board disposition described
in `list`, a dice strategy (dice to play at each square) in `dice`
and the circularity of the board indicated in `circular`.
"""
function simulate_k(list::Vector{Int64}, dice::Vector{Int64}, circular::Bool, k::Int64)
	cost = 0
	box = k
	while box != 15
		box = oneturn(list, box, dice[box], circular)
		# println(box)
		cost = cost + 1
	end
	return cost
end

"""
    oneturn(list, start, dice, circular)

This function is used in the `simulate` function and played until
the current square is not the square 15. Given the board disposition
in `list` and its circularity in `circular`, the current square `start`
and the `dice` to be played, it will return the next square (using
random variables to simulate the dice and traps behaviors).
"""
function oneturn(list::Vector{Int64}, start::Int64, dice::Int64, circular::Bool)
	startatthree = false
	if (start >= 11) && (start <= 14)
		path = 2
	else
		path = 1
	end
	if start == 3
		startatthree = true
		if rand([false, true])
			path = 2
			start = 10
		end
	end
	if dice == 1
		nextbox = start + rand(0:1)
	elseif dice == 2
		nextbox = start + rand(0:2)
	elseif dice == 3
		nextbox = start + rand(0:3)
	end
	if path == 1
		if nextbox > 10
			if !circular
				return 15
			else
				if nextbox != 11
					nextbox = nextbox - 11
				else
					return 15
				end
			end
		end
	else
		if nextbox > 14
			if !circular
				return 15
			else
				if nextbox > 15
					nextbox = nextbox - 15
				else
					return 15
				end
			end
		end
		if (startatthree) && (nextbox == 10)
			nextbox = 3
		end
	end
	if dice == 1
		return nextbox
	elseif dice == 2
		if rand([false, true])
			return nextsquare(list[nextbox], nextbox)
		else
			return nextbox
		end
	elseif dice == 3
		return nextsquare(list[nextbox], nextbox)
	end
end

function nextsquare(square_type, nextbox)
    if square_type == 0
    	return nextbox
    elseif square_type == 1
    	return applytrap1(nextbox)
    elseif square_type == 2
    	return applytrap2(nextbox)
    elseif square_type == 3
    	return applytrap3(nextbox)
    end
end

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
#=
    A module for digital communications
=#
module DigitalCommunications

using Plots
using QuadGK
using SpecialFunctions
using DocStringExtensions
using LinearAlgebra
import Base: show

include("utils.jl")
include("theoretical.jl")
include("source.jl")
include("coding.jl")
include("pulses.jl")
include("schemes.jl")
include("modulator.jl")
include("channel.jl")
include("detector.jl")

end # module


text:
module PQEcolaPoint

# include my codes -
include("Include.jl")

# export types -
export AbstractAssetModel
export AbstractDerivativeContractModel
export AbstractLatticeModel
export PutContractModel
export CallContractModel
export EquityModel
export CRRLatticeModel
export CRRContractPremiumLatticeModel
export CRRJITContractPremiumLatticeModel
export PQContractPremiumLatticePoint

# export methods -
export expiration
export build
export price
export premium
export intrinsic
export decode
export encode
export save
export load
export move
export tabu
export roll

# greeks -
export θ
export δ
export γ
export ρ
export vega

end # module


text:
type RandomPolicy <: Policy
    mdp::POMDP # contains the model
    rng::AbstractRNG
    action_space::AbstractSpace
    action::Action
end
# constructor
function RandomPolicy(mdp::POMDP, rng::AbstractRNG)
    as = actions(mdp)
    a = create_action(mdp)
    return RandomPolicy(mdp, rng, as, a)
end

function POMDPs.action(policy::RandomPolicy, state::State)
    action_space = actions(policy.mdp, state, policy.action_space)
    return rand!(policy.rng, policy.action, policy.action_space)
end

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch