text:
module MountainCarUtils
using GVFHordes
# import StatsBase
using StatsBase
using SparseArrays
import Random

const discount = 0.99
import ..Curiosity
import ..FeatureCreator

const SRCU = Curiosity.SRCreationUtils

import ..MountainCarConst
import ..GVFSRHordes

function load_policy(policy_name)
    Curiosity.LearnedPolicy("./src/data/MC_learned_policies/$(policy_name).bson")
end

function get_policies(parsed)
    policies = if parsed["learned_policy"]
        [load_policy(name) for name in parsed["learned_policy_names"]]
    else
        [EnergyPumpPolicy(true),EnergyPumpPolicy(true)]
    end
    return policies
end

function create_demons(parsed, fc)
    action_space = 3
    policies = get_policies(parsed)
    demons = if parsed["demon_learner"] != "SR"
        GVFHordes.Horde(
            [steps_to_wall_gvf(policies[1]),steps_to_goal_gvf(policies[2])])
    elseif parsed["demon_learner"] == "SR"
        @assert fc != nothing
        # pred_horde =  GVFHordes.Horde(
        #         [steps_to_wall_gvf(policies[1]),steps_to_goal_gvf(policies[2])])

        # pred_horde =  GVFHordes.Horde(
        #                 [GVFHordes.GVF(GVFHordes.GVFParamFuncs.FeatureCumulant(i+2),
        #                      GVFHordes.GVFParamFuncs.ConstantDiscount(0.0),
        #                      GoalPolicy(i)) for i in 1:4])

        pred_horde =  GVFHordes.Horde(
                    [GVFHordes.GVF(step_to_wall_cumulant(),
                        GVFHordes.GVFParamFuncs.ConstantDiscount(0.0),
                        policies[1]),
                    GVFHordes.GVF(step_to_goal_cumulant(),
                        GVFHordes.GVFParamFuncs.ConstantDiscount(0.0),
                        policies[2],
                        )
                     ])

        SF_policies = [g.policy for g in pred_horde.gvfs]
        SF_discounts = [GVFParamFuncs.StateTerminationDiscount(discount, is_wall, 0.0), GVFParamFuncs.StateTerminationDiscount(discount, is_goal, 0.0)]
        num_SFs = length(SF_policies)
        SF_horde = SRCU.create_SF_horde(SF_policies, SF_discounts, fc,1:action_space)
        GVFSRHordes.SRHorde(pred_horde, SF_horde, num_SFs, fc)
    else
        throw(ArgumentError("Cannot create demons"))
    end
    return demons
end

struct MountainCarStateActionCumulant <: GVFParamFuncs.AbstractCumulant
    state_num::Int
    action::Int
    state_constructor::Any
end

function Base.get(cumulant::MountainCarStateActionCumulant; kwargs...)
    state = kwargs[:constructed_state_t]
    action = kwargs[:action_t]
    if state[cumulant.state_num] == 1 && action == cumulant.action
        return 1
    else
        return 0
    end
end

# function make_behaviour_gvf(discount, state_constructor_func, learner, exploration_strategy)
function make_behaviour_gvf(behaviour_learner, γ, fc, exploration_strategy)
    function b_π(state_constructor, learner, exploration_strategy; kwargs...)
        s = state_constructor(kwargs[:state_t])
        preds = learner(s)
        return exploration_strategy(preds)[kwargs[:action_t]]
    end
    GVF_policy = GVFParamFuncs.FunctionalPolicy((;kwargs...) -> b_π(fc, behaviour_learner, exploration_strategy; kwargs...))
    BehaviourGVF = GVF(GVFParamFuncs.RewardCumulant(), goal_pseudoterm(γ), GVF_policy)
end

# function make_SF_for_policy(gvf_policy, gvf_pseudoterm, num_features, num_actions, state_constructor)
#     return GVFSRHordes.SFHorde([GVF(MountainCarStateActionCumulant(s,a,state_constructor),
#                     gvf_pseudoterm,
#                     gvf_policy) for s in 1:num_features for a in 1:num_actions])
# end
#
# function make_SF_horde(discount, num_features, num_actions, state_constructor)
#     steps_to_wall_horde = make_SF_for_policy(EnergyPumpPolicy(true),
#                 GVFParamFuncs.StateTerminationDiscount(discount, steps_to_wall_pseudoterm, 0.0),
#                 num_features, num_actions, state_constructor)
#
#     steps_to_goal_horde = make_SF_for_policy(EnergyPumpPolicy(true),
#                 GVFParamFuncs.StateTerminationDiscount(discount, steps_to_goal_pseudoterm, 0.0),
#                 num_features, num_actions, state_constructor)
#
#     SF_horde = GVFSRHordes.merge(steps_to_wall_horde,steps_to_goal_horde)
#     return SF_horde
# end

struct IdealDemonFeatures <: FeatureCreator
end

function project_features(fc::IdealDemonFeatures, state,action,state_tp1)
    new_state = sparsevec(convert(Array{Int,1},[is_wall(;state_tp1 = state_tp1), is_goal(;state_tp1 = state_tp1)] ))
    return new_state
end

(FP::IdealDemonFeatures)(state, action, next_state) = project_features(FP, state, action, next_state)
Base.size(FP::IdealDemonFeatures) = 2
Base.length(FP::IdealDemonFeatures) = 2

function MCNorm(obs)4
    pos_limit = MountainCarConst.pos_limit
    vel_limit = MountainCarConst.vel_limit
    return Float32[(obs[1] - pos_limit[1])/(pos_limit[2] - pos_limit[1]),
                   (obs[2] - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end

function task_gvf()
    GVF(GVFParamFuncs.FunctionalCumulant(task_cumulant),
        GVFParamFuncs.StateTerminationDiscount(discount, task_pseudoterm, 0.0),
        EnergyPumpPolicy(true))
end

function steps_to_wall_gvf(policy)
    GVF(step_to_wall_cumulant(),
        wall_pseudoterm(discount),
        policy)
end

function steps_to_goal_gvf(policy)
    GVF(step_to_goal_cumulant(),
        goal_pseudoterm(discount),
        policy)
end

function leave_cumulant()
    GVFParamFuncs.FunctionalCumulant(task_cumulant)
end
function step_to_goal_cumulant()
    GVFParamFuncs.FunctionalCumulant(steps_to_goal)
end
function step_to_wall_cumulant()
    GVFParamFuncs.FunctionalCumulant(steps_to_wall)
end
function negative_step_to_goal_cumulant()
    GVFParamFuncs.FunctionalCumulant(negative_steps_to_goal)
end
function negative_step_to_wall_cumulant()
    GVFParamFuncs.FunctionalCumulant(negative_steps_to_wall)
end
function wall_pseudoterm(γ)
    GVFParamFuncs.StateTerminationDiscount(γ, is_wall, 0.0)
end
function goal_pseudoterm(γ)
    GVFParamFuncs.StateTerminationDiscount(γ, is_goal, 0.0)
end

struct EnergyPumpPolicy <: GVFParamFuncs.AbstractPolicy
    normalized::Bool
    EnergyPumpPolicy(normalized=false) = new(normalized)
end

vel_0(π::EnergyPumpPolicy) = if π.normalized == true
    0.5
else
    0.0
end

# NOTE IS THIS TP1 or just T
Base.get(π::EnergyPumpPolicy; kwargs...) =
    if kwargs[:state_t][2] >= vel_0(π)
        MountainCarConst.Accelerate == kwargs[:action_t] ? 1.0 : 0.0 # Go Right/Accelerate
    else
        MountainCarConst.Reverse == kwargs[:action_t] ? 1.0 : 0.0 # Go left/Reverse
    end

StatsBase.sample(rng::Random.AbstractRNG, π::EnergyPumpPolicy; kwargs...) =
    if kwargs[:state_t] state_t[2] >= vel_0(π)
        MountainCarConst.Accelerate # Go Right/Accelerate
    else
        MountainCarConst.Reverse # Go left/Reverse
    end
StatsBase.sample(rng::Random.AbstractRNG, π::EnergyPumpPolicy, state_t, actions) =
    sample(rng, Weights([get(π; state_t = state_t, action_t=a) for a in actions]))
StatsBase.sample(π::EnergyPumpPolicy, args...) =
    StatsBase.sample(Random.GLOBAL_RNG, π::EnergyPumpPolicy, args...)

# StatsBase.sample(rng::Random.AbstractRNG, π::EnergyPumpPolicy, state_t) =
#     sample(rng, Weights([get(π; state_t = state_t, action_t=a) for a in actions]))

function is_wall(;kwargs...)
    normed_wall_pos = 0.0
    # return obs[1] <= MountainCarConst.pos_limit[1]
    return kwargs[:state_tp1][1] <= normed_wall_pos
end
function steps_to_wall(;kwargs...)
    normed_wall_pos = 0.0
    # return obs[1] <= MountainCarConst.pos_limit[1] ? 1.0 : 0.0
    return kwargs[:state_tp1][1] <= normed_wall_pos ? 1.0 : 0.0
end
function negative_steps_to_wall(;kwargs...)
    return -1
end

function is_goal(;kwargs...)
    normed_goal_pos = 1.0
    # return obs[1] >= MountainCarConst.pos_limit[2]
    return kwargs[:state_tp1][1] >= normed_goal_pos
end

function steps_to_goal(;kwargs...)
    normed_goal_pos = 1.0
    # return obs[1] >= MountainCarConst.pos_limit[2] ? 1.0 : 0.0
    return kwargs[:state_tp1][1] >= normed_goal_pos ? 1.0 : 0.0
end
function negative_steps_to_goal(;kwargs...)
    return -1
end

function task_pseudoterm(;kwargs...)
    normed_goal_pos = 1.0
    return kwargs[:state_tp1][1] >= normed_goal_pos
end

function task_cumulant(;kwargs...)
    normed_goal_pos = 1.0
    # return obs[1] >= MountainCarConst.pos_limit[2] ? 1.0 : 0.0
    return kwargs[:state_tp1][1] >= normed_goal_pos ? 0.0 : -1.0
end

end #end MountainCarUtils

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# PROBABILISTIC CIRCUITS LIBRARY ROOT

module ProbabilisticCircuits

# USE EXTERNAL MODULES

using Reexport

using LogicCircuits
# only reexport selectively from LogicCircuits
export pos_literals, neg_literals
# circuit queries
export issmooth, isdecomposable, isstruct_decomposable,
       isdeterministic, iscanonical
# circuit status
export num_edges, num_parameters
# datasets
export twenty_datasets

include("Utils/Utils.jl")
@reexport using .Utils

include("FactorGraph/factor_graph.jl")
include("FactorGraph/fg_compile.jl")

include("abstract_prob_nodes.jl")
include("shared_prob_nodes.jl")
include("plain_prob_nodes.jl")
include("structured_prob_nodes.jl")
include("logistic_nodes.jl")
include("param_bit_circuit.jl")
include("param_bit_circuit_pair.jl")
include("parameters.jl")
include("gradient_based_learning.jl")

include("queries/likelihood.jl")
include("queries/marginal_flow.jl")
include("queries/map.jl")
include("queries/sample.jl")
include("queries/pr_constraint.jl")
include("queries/information.jl")
include("queries/expectation_rec.jl")
include("queries/expectation_graph.jl")
include("queries/expectation_bit.jl")

include("Logistic/Logistic.jl")
@reexport using .Logistic

include("mixtures/em.jl")

include("structurelearner/chow_liu_tree.jl")
include("structurelearner/init.jl")
include("structurelearner/heuristics.jl")
include("structurelearner/learner.jl")
include("structurelearner/vtree_learner.jl")

include("LoadSave/LoadSave.jl")
@reexport using .LoadSave

using Requires

function __init__()
    # optional dependency
    @require BlossomV = "6c721016-9dae-5d90-abf6-67daaccb2332" include("structurelearner/vtree_learner_blossomv.jl")
end

end

ext: jl


text:
@testset "Test base checks" begin
    unordered = [1, 4, 2, 6]

    @test_throws(PowerSystems.DataFormatError,
                 PowerSystems.check_ascending_order(unordered, "test"))

    ordered = sort(unordered)
    PowerSystems.check_ascending_order(ordered, "test")

end


text:
using DiffEqOperators

s = x, y, z = (-5:0.2:5, -5:0.2:5, -5:0.2:5)
dx = dy = dz = x[2] - x[1]

ricker(x::T, y::T, z::T) where T = (4*(x^2+y^2+z^2) - 6)*exp(-(x^2+y^2+z^2))

u0 = [ricker(X, Y, Z) for Z in z, Y in y, X in x]

Dxx = CenteredDifference{1}(2, 4, dx, length(x))
Dyy = CenteredDifference{2}(2, 4, dy, length(y))
Dzz = CenteredDifference{3}(2, 4, dz, length(z))

A = Dxx+Dyy+Dzz
Q = compose(Dirichlet0BC(Float64, length.(s))...)

dt = dx/(sqrt(3)*3e8)
t = 0.0:dt:10/3e8

f(u,p,t) = (3e8)^2 .*(A*Q*u)

function steptime(u,uold,uolder)
    return ((dt^2 .*f(u,0,0) .+ 5u .- 4uold .+ uolder)./2, u, uold)
end
let uolder = deepcopy(u0), uold = deepcopy(u0), u = deepcopy(u0)
    u,uold,uolder = steptime(u0,uold,uolder)
    # @gif for ti in t #4th order time stepper
    #         u, uold, uolder = steptime(u,uold,uolder)
    #         heatmap(u)
    # end
    for ti in t #4th order time stepper
        u, uold, uolder = steptime(u,uold,uolder)
    end
end


text:
using YaoBase

export ShiftGate, shift

"""
    ShiftGate <: PrimitiveBlock

Phase shift gate.

# Definition

```math
\\begin{pmatrix}
1 & 0
0 & e^(im θ)
\\end{pmatrix}
```
"""
mutable struct ShiftGate{T <: Real} <: PrimitiveBlock{1}
    theta::T
end

"""
    shift(θ)

Create a [`ShiftGate`](@ref) with phase `θ`.

# Example

```jldoctest
julia> shift(0.1)
shift(0.1)
```
"""
shift(θ::AbstractFloat) = ShiftGate(θ)
shift(θ::Real) = shift(Float64(θ))
mat(::Type{T}, gate::ShiftGate) where {T <: Complex} = Diagonal(T[1.0, exp(im * gate.theta)])

cache_key(gate::ShiftGate) = gate.theta

# parametric interface
niparams(::Type{<:ShiftGate}) = 1
getiparams(x::ShiftGate) = x.theta
setiparams!(r::ShiftGate, param::Real) = (r.theta = param; r)

Base.adjoint(blk::ShiftGate) = ShiftGate(-blk.theta)
Base.copy(block::ShiftGate{T}) where T = ShiftGate{T}(block.theta)
Base.:(==)(lhs::ShiftGate, rhs::ShiftGate) = lhs.theta == rhs.theta
YaoBase.isunitary(r::ShiftGate) = true


text:
using Reinforce
using Test

@testset "Reinforce" begin
  include("foo.jl")
  include("interface.jl")

  @testset "env" begin
    include("env/mountain_car.jl")
    include("env/multi-armed-bandit.jl")
  end
end


text:
using EmailScraper
using Cascadia
using Gumbo
using HTTP
using Test

const BODY = parsehtml("""

                       Jose Storopoli

Associate Professor and Researcher of the
                       Department of Computer Science at
                       Universidade Nove de Julho - UNINOVE located in São Paulo - Brazil.

Teaches undergraduate and graduate courses in Data Science, Statistics, Bayesian Statistics, Machine Learning and Deep Learning using Julia
                       , R, Python, and Stan
                       . Contributor to Julia, R and
                       Stan ecosystems. Proficient in C/C++ and
                       Rust. Has published Julia, Rust, R
                       , and Python packages in official repositories/registries.

                       Researches, publishes and advises PhD candidates on topics about Bayesian Statistical Modeling and Machine Learning applied to Decision Making. Principal Investigator of
                       LabCidades - Smart City Research Lab at UNINOVE.

                       Coauthor of Julia Data Science upcoming book
                       . Certified
                       RStudio Tidyverse Instructor.

 Julia for Data Science Book

 Bayesian Statistics using Julia and Turing

 Ciência de Dados e Computação Científica com Julia

 EmailScraper.jl

 Workshop of Bayesian Inference with Julia and Turing

 Estatística Bayesiana com R e Stan

 Ciência de Dados pandas, Scikit-Learn, TensorFlow e PyTorch

 bibjoin

 Estatística com R

 Ciência de Dados com R e tidyverse

 Rcpp

 FactorAssumptions

 bibexcel

 Website built with Franklin.jl and the Julia programming language.

                       """).root[2]

@testset "EmailScraper" begin
    @testset "Domain" begin
        @testset "request_body" begin
            @test EmailScraper.request_body("storopoli.io") isa HTMLElement
            @test isempty(text(EmailScraper.request_body("string")))
        end

        @testset "valid_links" begin
            @test EmailScraper.valid_link("https://julialang.org") ==
                  "https://julialang.org"
            @test EmailScraper.valid_link("https://julialang.org/temp.doc") isa Missing
            @test EmailScraper.valid_link("https://julialang.org/temp.docx") isa Missing
            @test EmailScraper.valid_link("https://julialang.org/temp.@@search") isa Missing
            @test EmailScraper.valid_link("https://facebook.com") isa Missing
            @test EmailScraper.valid_link("https://facebook.com") isa Missing
        end

        @testset "get_links" begin
            @test "https://github.com/storopoli" in EmailScraper.get_links(BODY)
            # if it returns true
        end
    end # Links

    @testset "Emails" begin
        @testset "extract_email" begin
            @test EmailScraper.extract_emails("LOREMIPUSM hello@hello.com LoReM Ipsum") ==
                  ["hello@hello.com"]
            @test EmailScraper.extract_emails(
                "LOREMIPUSM hello@hello.com ### !! LoReM hello@hello.com Ipsum"
            ) == ["hello@hello.com"]
            @test EmailScraper.extract_emails(
                "LOREMIPUSM hello@hello.com LoReM## !! @ hello2@hello.com Ipsum"
            ) == ["hello@hello.com", "hello2@hello.com"]
        end

        @testset "valid_email" begin
            @test EmailScraper.valid_email("hello@hello.com") == true
            @test EmailScraper.valid_email("hello@hello.com.br") == true
            @test EmailScraper.valid_email("hello@hello.co.uk") == true
            @test EmailScraper.valid_email("hello@hello.png") == false
            @test EmailScraper.valid_email("hello@hello.pdf") == false
            @test EmailScraper.valid_email("hello@hello.xpto") == false
        end
    end # Emails

    @testset "scrape_url" begin
        @test scrape_url("http://www.webscraper.io/test-sites/e-commerce/allinone") ==
              String[]
        @test scrape_url("storopoli.io") == ["thestoropoli@gmail.com"]
        @test "contact@julialang.org" in scrape_url("julialang.org/about/help/")
    end
end # EmailScraper


text:
"""
`HBlkDiag` - a homogeneous block diagonal matrix, i.e. all the diagonal blocks are the same size

A matrix consisting of k diagonal blocks of size `r×s` is stored as an `r×s×k` array.
"""
immutable HBlkDiag{T} <: AbstractMatrix{T}
    arr::Array{T,3}
end

Base.eltype{T}(A::HBlkDiag{T}) = T

Base.size(A::HBlkDiag) = ((r,s,k) = size(A.arr); (r*k,s*k))

function Base.size(A::HBlkDiag,i::Integer)
    i < 1 && throw(BoundsError())
    i > 2 && return 1
    r,s,k = size(A.arr)
    (i == 1 ? r : s)*k
end

Base.copy!{T}(d::HBlkDiag{T},s::HBlkDiag{T}) = (copy!(d.arr,s.arr); d)

Base.copy{T}(s::HBlkDiag{T}) = HBlkDiag(copy(s.arr))

function Base.LinAlg.A_ldiv_B!(R::DenseVecOrMat,A::HBlkDiag,B::DenseVecOrMat)
    Aa = A.arr
    r,s,k = size(Aa)
    r == s || throw(ArgumentError("A must be square"))
    (m = size(B,1)) == size(R,1) || throw(DimensionMismatch())
    (n = size(B,2)) == size(R,2) || throw(DimensionMismatch())
    r*k == m || throw(DimensionMismatch())
    if r == 1
        for j in 1:n, i in 1:m
            R[i,j] = B[i,j]/Aa[i]
        end
    else
        db = similar(A.arr,(r,r))       # will hold the diagonal blocks
        for b in 1:k
            copy!(db,sub(A.arr,:,:,b))
            rows = (1:r)+(b-1)*r
            rr = copy!(sub(R,rows,:),sub(B,rows,:))
            Base.LinAlg.A_ldiv_B!(ishermitian(db) ? cholfact!(db) : lufact!(db), rr)
        end
    end
    R
end

function Base.getindex{T}(A::HBlkDiag{T},i::Integer,j::Integer)
    Aa = A.arr
    r,s,k = size(Aa)
    bi,ri = divrem(i-1,r)
    bj,rj = divrem(j-1,s)
    bi == bj || return zero(T)
    Aa[ri+1,rj+1,bi+1]
end

function Base.cholfact!(A::HBlkDiag,uplo::Symbol=:U)
    Aa = A.arr
    r,s,k = size(Aa)
    r == s || throw(ArgumentError("A must be square"))
    if r == 1
        for j in 1:k
            Aa[1,1,j] = sqrt(Aa[1,1,j])
        end
    else
        for j in 1:k
            cholfact!(sub(Aa,:,:,j),uplo)
        end
    end
    A
end

function Base.full(A::HBlkDiag)
    aa = A.arr
    res = zeros(eltype(aa),size(A))
    p,q,l = size(aa)
    for b in 1:l
        bm1 = b - 1
        for j in 1:q
            for i in 1:p
                res[bm1*p+i,bm1*q+j] = aa[i,j,b]
            end
        end
    end
    res
end

"""
Equivalent to `A = A + I` without making a copy of A
"""
function inflate!(A::HBlkDiag)
    Aa = A.arr
    r,s,k = size(Aa)
    for j in 1:k, i in 1:min(r,s)
        Aa[i,i,j] += 1
    end
    A
end


text:
# Laplacian

"""
    laplacian!(v,w)

Evaluate the discrete Laplacian of `w` and return it as `v`. The data `w` can be
of type dual/primary nodes or edges; `v` must be of the same type.

# Example

```jldoctest
julia> w = Nodes(Dual,(8,6));

julia> v = deepcopy(w);

julia> w[4,3] = 1.0;

julia> laplacian!(v,w)
Nodes{Dual,8,6,Float64} data
Printing in grid orientation (lower left is (1,1))
6×8 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -4.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function laplacian!(out::Nodes{Dual,NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    view(out,2:NX-1,2:NY-1) .= view(w,2:NX-1,1:NY-2) .+ view(w,1:NX-2,2:NY-1) .+
                               view(w,3:NX,2:NY-1) .+ view(w,2:NX-1,3:NY) .-
                               4.0.*view(w,2:NX-1,2:NY-1)
    #@inbounds for y in 2:NY-1, x in 2:NX-1
    #    out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    #end
    out
end

function laplacian!(out::Nodes{Primal,NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    view(out,2:NX-2,2:NY-2) .= view(w,2:NX-2,1:NY-3) .+ view(w,1:NX-3,2:NY-2) .+
                               view(w,3:NX-1,2:NY-2) .+ view(w,2:NX-2,3:NY-1) .-
                               4.0.*view(w,2:NX-2,2:NY-2)
    #@inbounds for y in 2:NY-2, x in 2:NX-2
    #    out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    #end
    out
end

function laplacian!(out::Edges{Dual,NX, NY}, w::Edges{Dual,NX, NY}) where {NX, NY}
  view(out.u,2:NX-2,2:NY-1) .= view(w.u,2:NX-2,1:NY-2) .+ view(w.u,1:NX-3,2:NY-1) .+
                          view(w.u,3:NX-1,2:NY-1) .+ view(w.u,2:NX-2,3:NY) .-
                          4.0.*view(w.u,2:NX-2,2:NY-1)
  view(out.v,2:NX-1,2:NY-2) .= view(w.v,2:NX-1,1:NY-3) .+ view(w.v,1:NX-2,2:NY-2) .+
                               view(w.v,3:NX,2:NY-2) .+ view(w.v,2:NX-1,3:NY-1) .-
                               4.0.*view(w.v,2:NX-1,2:NY-2)
  #@inbounds for y in 2:NY-1, x in 2:NX-2
  #    out.u[x,y] = w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1]
  #end
  #@inbounds for y in 2:NY-2, x in 2:NX-1
  #    out.v[x,y] = w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1]
  #end
  out
end

function laplacian!(out::Edges{Primal,NX, NY}, w::Edges{Primal,NX, NY}) where {NX, NY}
  view(out.u,2:NX-1,2:NY-2) .= view(w.u,2:NX-1,1:NY-3) .+ view(w.u,1:NX-2,2:NY-2) .+
                               view(w.u,3:NX,2:NY-2) .+ view(w.u,2:NX-1,3:NY-1) .-
                              4.0.*view(w.u,2:NX-1,2:NY-2)
  view(out.v,2:NX-2,2:NY-1) .= view(w.v,2:NX-2,1:NY-2) .+ view(w.v,1:NX-3,2:NY-1) .+
                                view(w.v,3:NX-1,2:NY-1) .+ view(w.v,2:NX-2,3:NY) .-
                               4.0.*view(w.v,2:NX-2,2:NY-1)
  #@inbounds for y in 2:NY-2, x in 2:NX-1
  #    out.u[x,y] = w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1]
  #end
  #@inbounds for y in 2:NY-1, x in 2:NX-2
  #    out.v[x,y] = w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1]
  #end
  out
end

"""
    laplacian(w)

Evaluate the discrete Laplacian of `w`. The data `w` can be of type dual/primary
nodes or edges. The returned result is of the same type as `w`.

# Example

```jldoctest
julia> q = Edges(Primal,(8,6));

julia> q.u[2,2] = 1.0;

julia> laplacian(q)
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -4.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function laplacian(w::Nodes{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(Nodes(C,w), w)
end

function laplacian(w::Edges{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(Edges(C,w), w)
end

"""
    laplacian_symm!(v,w)

Evaluate the symmetric 5-point discrete Laplacian of `w` and return it as `v`. The data `w` can be
of type dual nodes only for now. This symmetric Laplacian also evaluates the
partial Laplacians (using only available stencil data) on the ghost nodes.
"""
function laplacian_symm!(out::Nodes{Dual,NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    @inbounds for y in 2:NY-1, x in 2:NX-1
        out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    end
    @inbounds for y in 2:NY-1
        out[1,y]  = w[1,y-1]            - 4w[1,y] + w[2,y] + w[1,y+1]
        out[NX,y] = w[NX,y-1] + w[NX-1,y]- 4w[NX,y]        + w[NX,y+1]
    end
    @inbounds for x in 2:NX-1
        out[x,1]  = w[x-1,1] + w[x+1,1] - 4w[x,1] + w[x,2]
        out[x,NY] = w[x-1,NY]+ w[x+1,NY]- 4w[x,NY] + w[x,NY-1]
    end
    out[1,1] = -4w[1,1] + w[1,2] + w[2,1]
    out[NX,1] = -4w[NX,1] + w[NX-1,1] + w[NX,2]
    out[1,NY] = -4w[1,NY] + w[1,NY-1] + w[2,NY]
    out[NX,NY] = -4w[NX,NY] + w[NX,NY-1] + w[NX-1,NY]
    out
end

"""
    plan_laplacian(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [dx=1.0],[dtype=Float64])

Constructor to set up an operator for evaluating the discrete Laplacian on
dual or primal nodal data of dimension `dims`. If the optional keyword
`with_inverse` is set to `true`, then it also sets up the inverse Laplacian
(the lattice Green's function, LGF). These can then be applied, respectively, with
`*` and `\\` operations on data of the appropriate size. The optional parameter
`dx` is used in adjusting the uniform value of the LGF to match the behavior
of the continuous analog at large distances; this is set to 1.0 by default. The
type of data on which to act is floating point by default, but can also be ComplexF64.
This is specified with the optional parameter `dtype`

Instead of the first argument, one can also supply `w::Nodes` to specify the
size of the domain.

# Example

```jldoctest
julia> w = Nodes(Dual,(5,5));

julia> w[3,3] = 1.0;

julia> L = plan_laplacian(5,5;with_inverse=true)
Discrete Laplacian (and inverse) on a (nx = 5, ny = 5) grid acting on Float64 data with spacing 1.0

julia> s = L\\w
Nodes{Dual,5,5,Float64} data
Printing in grid orientation (lower left is (1,1))
5×5 Array{Float64,2}:
 0.16707    0.129276     0.106037     0.129276    0.16707
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.106037  -0.00734343  -0.257343    -0.00734343  0.106037
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.16707    0.129276     0.106037     0.129276    0.16707

julia> L*s ≈ w
true
```
"""
function plan_laplacian end

"""
    plan_laplacian!(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [dx=1.0])

Same as [`plan_laplacian`](@ref), but operates in-place on data.
"""
function plan_laplacian! end

struct Laplacian{NX, NY, T, R, DX, inplace}
    conv::Union{CircularConvolution{NX, NY, T},Nothing}
end

for (lf,inplace) in ((:plan_laplacian,false),
                     (:plan_laplacian!,true))
    @eval function $lf(dims::Tuple{Int,Int};
                   with_inverse = false, fftw_flags = FFTW.ESTIMATE, dx = 1.0, dtype = Float64)
        NX, NY = dims
        if !with_inverse
            return Laplacian{NX, NY, dtype, false, dx, $inplace}(nothing)
        end

        G = view(LGF_TABLE, 1:NX, 1:NY)
        Laplacian{NX, NY, dtype, true, dx, $inplace}(CircularConvolution(G, fftw_flags,dtype=dtype))
    end

    @eval function $lf(nx::Int, ny::Int;
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, dx = 1.0, dtype = Float64)
        $lf((nx, ny), with_inverse = with_inverse, fftw_flags = fftw_flags, dx = dx, dtype = dtype)
    end

    @eval function $lf(nodes::Nodes{T,NX,NY};
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, dx = 1.0, dtype = Float64) where {T<:CellType,NX,NY}
        $lf(node_inds(T,(NX,NY)), with_inverse = with_inverse, fftw_flags = fftw_flags, dx = dx, dtype = dtype)
    end
end

function Base.show(io::IO, L::Laplacian{NX, NY, T, R, DX, inplace}) where {NX, NY, T, R, DX, inplace}
    nodedims = "(nx = $NX, ny = $NY)"
    inverse = R ? " (and inverse)" : ""
    isinplace = inplace ? " in-place" : ""
    print(io, "Discrete$isinplace Laplacian$inverse on a $nodedims grid acting on $T data with spacing $DX")
end

mul!(out::Nodes{C,NX,NY}, L::Laplacian, s::Nodes{C,NX,NY}) where {C<:CellType,NX,NY} = laplacian!(out, s)
*(L::Laplacian{MX,MY,T,R,DX,false}, s::Nodes{C,NX,NY}) where {MX,MY,T,R,DX,C <: CellType,NX,NY} =
      laplacian(s)
function (*)(L::Laplacian{MX,MY,T,R,DX,true}, s::Nodes{C,NX,NY}) where {MX,MY,T,R,DX,C <: CellType,NX,NY}
    laplacian!(s,deepcopy(s))
end

function ldiv!(out::Nodes{C,NX, NY,T},
                   L::Laplacian{MX, MY, T, true, DX, inplace},
                   s::Nodes{C, NX, NY,T}) where {C <: CellType, NX, NY, MX, MY, T, DX, inplace}

    mul!(out.data, L.conv, s.data)

    # Adjust the behavior at large distance to match continuous kernel
    out.data .-= (sum(s.data)/2π)*(GAMMA+log(8)/2-log(DX))
    out
end

\(L::Laplacian{MX,MY,T,R,DX,false},s::Nodes{C,NX,NY}) where {MX,MY,T,R,DX,C <: CellType,NX,NY} =
  ldiv!(Nodes(C,s), L, s)

\(L::Laplacian{MX,MY,T,R,DX,true},s::Nodes{C,NX,NY}) where {MX,MY,T,R,DX,C <: CellType,NX,NY} =
  ldiv!(s, L, deepcopy(s))


text:
using TransformVariables

export ParameterizedMeasure
abstract type ParameterizedMeasure{N} <: AbstractMeasure end

function Base.getproperty(μ::ParameterizedMeasure{N}, prop::Symbol) where {N}
    return getproperty(getfield(μ, :par), prop)
end

function Base.propertynames(μ::ParameterizedMeasure{N}) where {N}
    return N
end

function Base.show(io::IO, ::MIME"text/plain", μ::ParameterizedMeasure{()})
    io = IOContext(io, :compact => true)
    print(io, nameof(typeof(μ)), "()")
end

function Base.show(io::IO, ::MIME"text/plain", μ::ParameterizedMeasure{N}) where {N}
    io = IOContext(io, :compact => true)
    print(io, nameof(typeof(μ)))
    print(io, getfield(μ,:par))
end

export asparams

# Allow things like
#
# julia> Normal{(:μ,)}(2)
# Normal(μ = 2,)
#
function (::Type{P})(args...) where {N, P <: ParameterizedMeasure{N}}
    C = constructorof(P)
    return C(NamedTuple{N}(args...))
end

(::Type{P})(;kwargs...) where {P <: ParameterizedMeasure} = P(NamedTuple(kwargs))

"""
`asparams` build on `TransformVariables.as` to construct bijections to the
parameter space of a given parameterized measure. Because this is only possible
for continuous parameter spaces, we allow constraints to assign values to any
subset of the parameters.

--------

    asparams(::Type{<:ParameterizedMeasure}, ::Val{::Symbol})

Return a transformation for a particular parameter of a given parameterized
measure. For example,

```
julia> asparams(Normal, Val(:σ))
asℝ₊
```

-----

    asparams(::Type{<: ParameterizedMeasure{N}}, constraints::NamedTuple) where {N}

Return a transformation for a given parameterized measure subject to the named tuple
`constraints`. For example,

```
julia> asparams(Binomial{(:p,)}, (n=10,))
TransformVariables.TransformTuple{NamedTuple{(:p,), Tuple{TransformVariables.ScaledShiftedLogistic{Float64}}}}((p = as𝕀,), 1)
```

------------

    aspararams(::ParameterizedMeasure)

Return a transformation with no constraints. For example,

```
julia> asparams(Normal{(:μ,:σ)})
TransformVariables.TransformTuple{NamedTuple{(:μ, :σ), Tuple{TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}}((μ = asℝ, σ = asℝ₊), 2)
```
"""
function asparams end

asparams(μ::ParameterizedMeasure, v::Val) = asparams(constructor(μ), v)
asparams(μ, s::Symbol) = asparams(μ, Val(s))

asparams(M::Type{PM}) where {PM<:ParameterizedMeasure} = asparams(M, NamedTuple())

function asparams(::Type{M}, constraints::NamedTuple{N2}) where {N1, N2, M<: ParameterizedMeasure{N1}}
    # @show M
    thekeys = params(M, constraints)
    t1 = NamedTuple{thekeys}(asparams(M, Val(k)) for k in thekeys)
    t2 = NamedTuple{N2}(map(asConst, values(constraints)))
    C = constructorof(M)
    # @show C
    # @show constraints
    # @show transforms
    # Make sure we end up with a consistent ordering
    ordered_transforms = params(C(merge(t1, t2)))
    return TV.as(ordered_transforms)
end

asparams(μ::ParameterizedMeasure, nt::NamedTuple=NamedTuple()) = asparams(constructor(μ), nt)

export params

params(μ::ParameterizedMeasure) = getfield(μ, :par)

function params(::Type{M}, constraints::NamedTuple{N2}=NamedTuple()) where {N1, N2, M<: ParameterizedMeasure{N1}}
thekeys = tuple((k for k in N1 if k ∉ N2)...)
end

params(μ) = ()

function ConstructionBase.setproperties(d::P, nt::NamedTuple) where {P<:ParameterizedMeasure}
    return constructorof(P)(merge(params(d), nt))
end

function Accessors.set(d::ParameterizedMeasure, ::typeof(params), nt::NamedTuple)
    setproperties(d, nt)
end

function Accessors.set(d::ParameterizedMeasure{N}, ::typeof(params), p) where {N}
    setproperties(d, NamedTuple{N}(p...))
end


text:
export ExpSweep, expsweep, analyze

# ALMOST COPY OF: https://github.com/ssfrr/MeasureIR.jl/blob/master/src/expsweep.jl

struct ExpSweep{AT, SR}
    sig::AT
    w1::Float64
    w2::Float64
    gain::Float64
    samplerate::SR
end

# the default constructor doesn't auto-convert the non-parametric arguments
ExpSweep(sig::T, w1, w2, g, sr::SR) where {T, SR} = ExpSweep{T, SR}(sig, w1, w2, g, sr)

"""
    expsweep(length; options...)
Create an impulse response measurement using an exponential sinusoid sweep from
`minfreq` to `maxfreq`, following the procedure from [1]. Durations and
frequencies can be given with or without Unitful units. With the default
`samplerate` of `1` durations can be interpreted as samples and frequencies are
in cycles/sample. The default minimum frequency of 0.0004 cycles/sample
corresponds to about 19Hz at 48kHz.
## Options
- `samplerate=1`
- `minfreq=0.0004`
- `maxfreq=0.5`
- `fadein`
- `fadeout`
- `fade`
- `optimize`
- `func`
## Examples
```julia
using MeasureIR
# generate a 10000-sample sweep
expsweep(10000)
# generate a 10s sweep at 48kHz samplerate
using Unitful: s, kHz
expsweep(10s; samplerate=48kHz)
```
[1]: Farina, Angelo, Simultaneous measurement of impulse response and distortion
with a swept-sine technique, Audio Engineering Society Convention 108 (2000).
"""
function expsweep(length;
        # TODO: consider making 2π the default samplerate, which makes the
        # frequencies interpretable as radians/sample
        samplerate = 1,
        minfreq=0.0004*samplerate,
        # TODO: this probably isn't the right default. for example for
        # meas = expsweep(30s; samplerate=48000Hz), it doesn't roll off on top
        maxfreq=0.5*samplerate,
        # default fades are computed below
        fadein=nothing,
        fadeout=nothing,
        fade=nothing,
        gain=1,
        optimize=true,
        func=sin)

    # TODO: review Novak: https://ant-novak.com/posts/research/2015-10-30_JAES_Swept/
    # it was recommended by Gordon

    w1 = uconvert(NoUnits, minfreq * 2π / samplerate)
    w2 = uconvert(NoUnits, maxfreq * 2π / samplerate)
    L = round(Int, uconvert(NoUnits, length*samplerate))

    if fade !== nothing
        if fadein !== nothing || fadeout !== nothing
            throw(ArgumentError("`fade` argument cannot be used with `fadein` or `fadeout`"))
        end
        fadein = fade
        fadeout = fade
    end

    # default fade times chosen heuristically to be pretty reasonable
    # we could definitely be smarter here, especially for small L
    if fadein === nothing
        f1 = round(Int, min(L/2, 200/w1))
    else
        f1 = round(Int, uconvert(NoUnits, fadein*samplerate))
    end

    if fadeout === nothing
        f2 = round(Int, min(L/4, 800/w2))
    else
        f2 = round(Int, uconvert(NoUnits, fadeout*samplerate))
    end

    if optimize
        w1 = _optimizew1(w1, w2, L)
    end
    sig = _expsweep(func, L, w1, w2)
    winin = 0.5 .- 0.5*cos.(range(0, π; length=f1))
    winout = 0.5 .- 0.5*cos.(range(π, 0; length=f2))
    sig[1:f1] .*= winin
    sig[(end-f2+1):end] .*= winout
    ExpSweep(sig, w1, w2, gain, samplerate)
end

function _expsweep(func, L, minfreq, maxfreq)
    # switch to notation from [1]
    w1 = minfreq
    w2 = maxfreq
    T = L
    t = 0:(T-1)
    @. func(w1*T /
           log(w2/w1) * (exp(t/T * log(w2/w1))-1))
end

"""
Compute a new starting freq w < w1 that gives us an integer number
of total cycles.
"""
function _optimizew1(w1, w2, L)
    # scale to [0-1]
    f1 = w1 * L
    f2 = w2 * L
    # compute the total (integer) number of cycles we want
    M = floor((f2 - f1) /
              (2π*log(f2/f1)))

    a(f) = (f2-f)/log(f2/f) - 2π*M
    fnew = find_zero(a, f1)

    fnew/L
end

"""
    function analyze(m::ExpSweep, response::AbstractArray; noncausal=false)
When the measurement is a sweep, you can pass the `noncausal=true` keyword
argument to include the noncausal part of the impulse response, which carries
energy from nonlinear distortion in the system. In this case the zero-lag part
of the impulse response starts at index `L÷2+1`, where `L` is the length of the
IR.
"""
function analyze(m::ExpSweep, response::AbstractArray; noncausal=false)
    # we don't truncate the response here because typically sweep measurements
    # are taken without repeating, so there's no real errors we'd catch by
    # truncating.
    T = length(m.sig)
    w1 = m.w1
    w2 = m.w2
    t = 0:(T-1)
    # compute the normalizing factor from the derivative of the frequency -
    # the faster the frequency is changing the less time the signal spends
    # at that frequency so we need to boost the amplitude to compensate
    amp = @. w1/T*log(w2/w1)*exp(t/T*log(w2/w1))
    invfilt = amp .* m.sig
    # calculate an extra scaling factor to get the center of the passband to
    # unity. At some point it might make sense to look at a range instead of
    # just the center frequency, in case we catch it in the middle of some
    # ripple. Generally things should be pretty flat if there's a fade in/out
    # TODO: seems kinda silly to compute the whole FFT just to look at the
    # center frequency. We could just do the dot product of that freq
    roundtrip = xcorr(m.sig * m.gain, invfilt; padmode=:longest)
    centeridx = round(Int, (w1 + (w2-w1)/2)/π*(length(roundtrip)/2))
    scale = abs(rfft(roundtrip)[centeridx])

    ir = mapslices(response, dims=1) do v
        xcorr(v, invfilt, padmode=:longest) ./ scale
    end

    # keep the full IR (including non-causal parts representing nonlinearities)
    # if m.nonlinear is true
    zeroidx = size(ir,1)÷2+1
    if noncausal
        SampleArray(reshape(ir[zeroidx-T:zeroidx+T], :, 1), m.samplerate)
    else
        # nonlinearity detection is not working right now, but it seems like a
        # useful feature to add in the future. One issue is that we need a way
        # to estimate the expected "noise floor" of the noncausal part of the
        # IR, which is _not_ the same as the noise floor of the response. It
        # definitely helps to deconvolve just the silence part and use that to
        # estimate the IR noise floor. The other issue is that the overall
        # energy in the noncausal part just doesn't change that much, because
        # all the extra energy is pretty concentrated and it gets averaged out.
        # probably to get this to work well we either need to explicitly look
        # for peaks, and possibly look at fixed delays from the main peak that
        # correspond to harmonics.
        # deconvsilence = mapslices(resp, 1) do v
        #     xcorr(v[1:prepadding(m)], invfilt)
        # end
        # nf = noisefloor(m, deconvsilence)

        # we don't check all the way up to the center because if the impulse is
        # bandlimited it'll have energy directly preceeding it. In reality
        # how far back we should go depends on the bandwidth of the response,
        # so this number may need to be tuned.
        # fortunately this is only to heuristically decide whether we print a
        # warning, so we're never going to get to 100% accuracy anyways...
        # preimpulse = 100
        # noncausalpower = sum(ir[zeroidx-T:zeroidx-preimpulse].^2) /
        #                  (2T+1-preimpulse)
        # if any(noncausalpower .> nf .* 1.1 .+ sqrt(eps()))
        #     @warn "Energy in noncausal IR is above noise floor. Check for nonlinearity"
        # end
        SampleArray(reshape(ir[zeroidx:zeroidx+T-1], :, 1), m.samplerate)
    end
end

function analyze(m::ExpSweep, response::SampleArray; kwargs...)
    @assert nchannels(response) == 1
    analyze(m, data(response)[:]; kwargs...)
end


text:
using .CUDAKernels, .CUDAKernels.CUDA
import ModelParameters.Flatten

export CuGPU

"""
    CuGPU <: GPU

    CuGPU()
    CuGPU{threads_per_block}()

```julia
ruleset = Ruleset(rule; proc=CuGPU())
# or
output = sim!(output, rule; proc=CuGPU())
```
"""
struct CuGPU{X} <: GPU end
CuGPU() = CuGPU{32}()

# CUDA setup

kernel_setup(::CuGPU{N}) where N = CUDAKernels.CUDADevice(), (N, N)

function Adapt.adapt_structure(to, x::Union{AbstractSimData,GridData,Rule})
    Flatten.modify(A -> adapt(to, A), x, Union{CuArray,Array,AbstractDimArray}, Union{SArray,Function})
end

function Adapt.adapt_structure(to, o::Output)
    frames = map(o.frames) do f
        if f isa NamedTuple
            a = map(g -> adapt(to, g), f)
        else
            adapt(to, f)
        end
    end
    @set o.extent = adapt(to, o.extent)
    @set o.frames = frames
end

@noinline function _proc_setup(::CuGPU, obj)
    Flatten.modify(CuArray, obj, Union{Array,BitArray}, Union{CuArray,SArray,Dict,Function})
end

_copyto_output!(outgrid, grid::GridData, proc::GPU) = copyto!(outgrid, gridview(grid))

# Thread-safe atomic ops

for (f, op) in atomic_ops
    atomic_f = Symbol(:atomic_, f)
    @eval begin
        @propagate_inbounds function ($f)(d::WritableGridData{<:Any,R}, ::CuGPU, x, I...) where R
            A = parent(dest(d))
            i = Base._to_linear_index(A, (I .+ R)...)
            (CUDA.$atomic_f)(pointer(A, i), x)
        end
    end
end


text:
include("constructors.jl")
include("helpers.jl")
include("show.jl")
include("auxiliary.jl")
include("functionality.jl")
include("recycle.jl")
# include("performance.jl")


text:
using ProgressMeter, BlackBoxOptim
include("util.jl")
const datafile = "inputs/23.txt"
struct Nanobot
    x::Int64
    y::Int64
    z::Int64
    r::Int64
end

function parse_input()::Array{Nanobot,1}
    lines = getlines(datafile)
    re = r"^pos=<(-?\d+),(-?\d+),(-?\d+)>, r=(\d+)$"
    nanobots = []
    for line in lines
        matches = match(re, line)
        try
            parts = [parse(Int64, matches[i]) for i in 1:4]
            nanobot = Nanobot(parts...)
            push!(nanobots, nanobot)
        catch
            println(line)
        end
    end
    return nanobots
end

function mdist(x,y,z,x2,y2,z2)::Int64
    return abs(x-x2)+abs(y-y2)+abs(z-z2)
end

function part_1()
    nanobots = parse_input()
    radii = sort(collect(nanobots), by=x->-x.r)
    biggest_radius = radii[1]
    in_range = []
    for nanobot in nanobots
        if mdist(biggest_radius.x, biggest_radius.y, biggest_radius.z,
            nanobot.x, nanobot.y, nanobot.z) <= biggest_radius.r
            push!(in_range, nanobot)
        end
    end
    println(length(in_range))
end

function minmaxxyz(nanobots::Array{Nanobot,1},r::Int64)
    x = []
    y = []
    z = []
    for nanobot in nanobots
        push!(x, nanobot.x)
        push!(y, nanobot.y)
        push!(z, nanobot.z)
    end
    return minimum(x)-r,maximum(x)+r,minimum(y)-r,maximum(y)+r,minimum(z)-r,maximum(z)+r
end

function corners(nb::Nanobot)::Array{Tuple{Int64,Int64,Int64}}
    retval = []
    push!(retval, (nb.x+nb.r,nb.y,nb.z))
    push!(retval, (nb.x-nb.r,nb.y,nb.z))
    push!(retval, (nb.x,nb.y+nb.r,nb.z))
    push!(retval, (nb.x,nb.y-nb.r,nb.z))
    push!(retval, (nb.x,nb.y,nb.z+nb.r))
    push!(retval, (nb.x,nb.y,nb.z-nb.r))
    return retval
end

function get_neighborhood(point::Tuple{Int64,Int64,Int64})::Array{Tuple{Int64,Int64,Int64}}
    retval = []
    for z in -15:15
        for y in -15:15
            for x in -15:15
                push!(retval, (point[1]+x,point[2]+y,point[3]+z))
            end
        end
    end
    return retval
end

function score(x::Int64,y::Int64,z::Int64,nb::Array{Nanobot})::Int64
    s = 0
    for n in nb
        if mdist(x,y,z,n.x,n.y,n.z) <= n.r
            s += 1
        end
    end
    return s
end

function score()
    nanobots = parse_input()
    return (w) -> begin
        x = round(w[1])
        y = round(w[2])
        z = round(w[3])
        score = 0
        for n in nanobots
            if mdist(x,y,z,n.x,n.y,n.z) <= n.r
                score += 1
            end
        end
        return float(1000-score)+(float(mdist(x,y,z,0,0,0))/848312125.0)
    end
end

nanobots = parse_input()
radii = sort(collect(nanobots), by=x->-x.r)
biggest_radius = radii[1]
minx, maxx, miny, maxy, minz, maxz = minmaxxyz(nanobots, biggest_radius.r)
fitness_func(f) = f[1]+(f[2]/mdist(0,0,0,maxx,maxy,maxz))

function part_2()
    nanobots = parse_input()
    radii = sort(collect(nanobots), by=x->-x.r)
    biggest_radius = radii[1]
    minx, maxx, miny, maxy, minz, maxz = minmaxxyz(nanobots, biggest_radius.r)
    mi = minimum([minx,miny,minz])
    ma = maximum([maxx,maxy,maxz])
    score_func = score()
    res = bboptimize(score_func; Method=:adaptive_de_rand_1_bin_radiuslimited, SearchRange = (float(mi),float(ma)), NumDimensions = 3, MaxTime = 1800.0, PopulationSize=1000)
    return res
end

part_1()
results = part_2()
results


text:
module ParametricFunctions

using ArgCheck
using ContinuousTransformations
import ContinuousTransformations: domain
using Lazy

export                          # FIXME reorder nicely when interface stabilizes
    points, degf, basis, basis!, Partial, ValuePartial, evaluate,
    basis_matrix, fit, fit!, ParametricFamily, family, parameters, fitfun,
    ParametricFunction

abstract FunctionFamily

abstract ParametricFamily

"Domain of the parametric function or a function family."
function domain end

"Degrees of freedom (number of coefficients)."
function degf end

"Recommended points for collocation and function fitting."
function points end

"""
Basis functions evaluated at a given point, written into the third
argument which should be a vector.
"""
function basis! end

"""
Return the basis functions of `p` evaluated at `x`, as a vector.
"""
basis{T}(p::ParametricFamily, x::T) = basis!(p, x, Vector{T}(degf(p)))

"""
Partial derivative in the given coordinate.
"""
immutable Partial{T} x::T end

"""
Value and partial derivative in the given coordinate.
"""
immutable ValuePartial{T} x::T end

"""
Evaluate a parametric function with the given parameters.
"""
function evaluate end

"""
Return the basis matrix of a parametric function family evaluated at
`xs`. Note that the returned value may or may not be a dense matrix,
but will always be a conformable AbstractMatrix.
"""
function basis_matrix{T}(family, xs::AbstractVector{T})
    B = Array{T}(degf(family), length(xs))
    for (i, x) in enumerate(xs)
        basis!(family, x, @view B[i, :])
    end
    B
end

"""
Fit a parametric function to a function family, placing the
coefficients in the third argument.
"""
function fit! end

"Fit a function to a parametric family."
fit!{T <: ParametricFamily}(p::T, f::Function, θ) =  fit!(p, f.(points(p)), θ)

"""
Fit a parametric function to a function family, returning the
coefficients.
"""
fit(p::ParametricFamily, ys::AbstractVector) = fit!(p, ys, Vector{eltype(ys)}(degf(p)))

fit(p::ParametricFamily, f::Function) = fit(p, f.(points(p)))

immutable ParametricFunction{TF <: ParametricFamily, Tθ <: AbstractVector}
    family::TF
    θ::Tθ
    function ParametricFunction(fam, θ)
        @argcheck degf(fam) == length(θ) "Incompatible coefficients."
        new(fam, θ)
    end
end

ParametricFunction{TF, Tθ}(fam::TF, θ::Tθ) = ParametricFunction{TF, Tθ}(fam, θ)

"Parametric family of a parametric function."
family(f::ParametricFunction) = f.family

"Parameters of a parametric function."
parameters(f::ParametricFunction) = f.θ

@forward ParametricFunction.family domain, points, degf

"""
Fit a function from a parametric family at values `ys`, which were
evaluated at `points(family)`.
"""
fitfun(family::ParametricFamily, ys) = ParametricFunction(family, fit(family, ys))

(f::ParametricFunction)(x) = evaluate(f.family, f.θ, x)

include("Chebyshev.jl")
include("domaintrans.jl")
include("valuetrans.jl")
include("collocation.jl")
include("plots.jl")

end # module


text:
using NearestNeighbors: TreeData, get_leaf_range, isleaf, getleft, getright

const BARNES_HUT_LEAFSIZE = 16
const BARNES_HUT_THETA = 1/4

# could generalize to something like HierarchicalFactorization
# have field that determines algorithm to execute Taylor(p), BarnesHut(), FKT(p)
struct BarnesHutFactorization{T, KT, XT<:AbstractVector, YT<:AbstractVector, TT, DT, RT} <: Factorization{T}
    k::KT # kernel function, for taylor branch, needs to be evaluable as k(r²)
    x::XT
    y::YT
    Tree::TT
    D::DT # diagonal correction (i.e. noise variance for GPs)
    θ::RT # distance threshold for compression
    # expansion_type::ET # IDEA to generalize tree-based algorithms: problem: would be nice to pre-allocate moments
    # w::WT # temporary storage for weight vectors, useful to split positive and negative values
    # i::BT
    # centers_of_mass::CT # pre-allocated, has to be able to work with type of w
end

# NOTE on ball tree indices
# if reorder = true, and T = BallTree(X), we have
# X[:, T.indices[i]] = T.data[i]
# X[:, i] == T.data[invpermute!(collect(1:length(T.data)), T.indices)][i]
function BarnesHutFactorization(k, x, y = x, D = nothing; θ::Real = BARNES_HUT_THETA, leafsize::Int = BARNES_HUT_LEAFSIZE)
    xs = vector_of_static_vectors(x)
    ys = x === y ? xs : vector_of_static_vectors(y)
    Tree = BallTree(ys, leafsize = leafsize)
    m = length(y)
    XT, YT, KT, TT, DT, RT = typeof.((xs, ys, k, Tree, D, θ))
    # w = zeros(length(m))
    # i = zeros(Bool, m)
    # WT, BT = typeof(w), typeof(i)
    T = gramian_eltype(k, xs[1], ys[1])
    BarnesHutFactorization{T, KT, XT, YT, TT, DT, RT}(k, xs, ys, Tree, D, θ) #, w, i)
end
function BarnesHutFactorization(G::Gramian; θ::Real = BARNES_HUT_THETA, leafsize::Int = BARNES_HUT_LEAFSIZE)
    BarnesHutFactorization(G.k, G.x, G.y, θ = θ, leafsize = leafsize)
end
Base.size(F::BarnesHutFactorization) = (length(F.x), length(F.y))
Base.size(F::BarnesHutFactorization, i::Int) = i ≤ 2 ? size(F)[i] : 1
Base.getindex(F::BarnesHutFactorization, i::Int, j::Int) = F.k(F.x[i], F.y[j])
Base.eltype(F::BarnesHutFactorization{T}) where T = T

function LinearAlgebra.mul!(y::AbstractVector, F::BarnesHutFactorization, x::AbstractVector, α::Real = 1, β::Real = 0)
    # barneshut!(y, F, x, α, β)
    if all(≥(0), x) # if all weights are positive, proceed with regular BarnesHut
        barneshut!(y, F, x, α, β, split = false)
    else
        taylor!(y, F, x, α, β)
    end
end
function Base.:*(F::BarnesHutFactorization{<:Number}, x::AbstractVector{<:Number})
    T = promote_type(eltype(F), eltype(x))
    y = zeros(T, size(F, 1))
    mul!(y, F, x)
end

# # use cg! if it's positive definite
# keywords include maxiter = size(F, 2)
# verbose: convergence info
# log: keep track of residual norm
# Pl, Pr: left and right preconditioners (not implemented?)
function LinearAlgebra.ldiv!(x::AbstractVector, F::BarnesHutFactorization, b::AbstractVector; kwargs...)
    # cg!(x, F, b; kwargs...)
    minres!(x, F, b; kwargs...)
end
function Base.:\(F::BarnesHutFactorization, b::AbstractVector; kwargs...)
    T = promote_type(eltype(F), eltype(b))
    x = zeros(T, size(F, 2))
    ldiv!(x, F, b; initially_zero = true, kwargs...)
end

############################ Core Barnes Hut algorithm ##########################
# uses the Barnes-Hut algorithm for an approximate n*log(n) multiply with the kernel matrix G
function barneshut!(b::AbstractVector, F::BarnesHutFactorization{<:Number}, w::AbstractVector,
                    α::Number = 1, β::Number = 0, θ::Real = F.θ; split::Bool = true)
    size(F, 2) == length(w) || throw(DimensionMismatch("length of w does not match second dimension of F: $(length(w)) ≠ $(size(F, 2))"))
    eltype(b) == promote_type(eltype(F), eltype(w)) ||
            throw(TypeError("eltype of target vector b not equal to eltype of matrix-vector product: $(eltype(b)) and $(promote_type(eltype(F), eltype(w)))"))

    @. b = iszero(β) ? 0 : β * b # this avoids trouble if b is initialized with NaN's, e.g. thorugh "similar"
    if split && any(<(0), w) # if split is on, we multiply with positive and negative component of w separately, tends to increase accuracy because of center of mass calculation
        splitting_barneshut!(b, F, w, α, β, θ)
    else
        sums_w = node_sums(w, F.Tree) # IDEA: could pre-allocate
        centers_of_mass = compute_centers_of_mass(F, w)
        @threads for i in eachindex(F.x) # exactly 4 * length(y) allocations?
            Fw_i = bh_recursion(1, F.k, F.x[i], F.y, w, sums_w, θ, F.Tree, centers_of_mass)::eltype(b) # x[i] creates an allocation
            b[i] += α * Fw_i
        end
        if !isnothing(F.D) # if there is a diagonal correction, need to add it
            mul!(b, F.D, w, α, 1)
        end
    end
    return b
end

# splits weight vector a into positive and negative components, tends to have
# much better accuracy for negative a, for only ~2x time penalty
function splitting_barneshut!(b::AbstractVector, F::BarnesHutFactorization, a::AbstractVector,
                    α::Number = 1, β::Number = 0, θ::Real = F.θ)
    c = copy(a) # IDEA: pre-allocate
    i = a .< 0 # first, blend out negative indices (multiply with positive part)
    @. c[i] = 0
    barneshut!(b, F, c, α, β, θ, split = false)
    @. c = -a # multiply with -a[i] to make entries positive, and use negative α to make result correct
    @. i = !i # blend out non-negative indices
    @. c[i] = 0
    barneshut!(b, F, c, -α, 1, θ, split = false) # use β = 1 because β is already taken care of above
    return b
end

# supposing y is a single value now
# k sould be isotropic kernel with euclidean norm
# x is single target data point
# y are all source data points
# sums_a[index] contains sum(a[get_node_indices(T.tree_data, index)])
# θ < 1 is the distance threshold beyond which compression is used
# e.g. when θ = 1/10, then we compress when the current target node is 10 times further
# from the center of the hypersphere than its radius
# when θ = 1, everything is compressed, when θ = 0, nothing is
function bh_recursion(index, k, xi, y::AbstractVector, w::AbstractVector,
                      sums_w::AbstractVector, θ::Real, T::BallTree, centers_of_mass)
    h = T.hyper_spheres[index]
    if isleaf(T.tree_data.n_internal_nodes, index) # do direct computation
        elty = promote_type(eltype(xi), eltype(eltype(y)), eltype(w))
        val = zero(elty)
        @inbounds @simd for i in get_leaf_range(T.tree_data, index) # loop vectorization?
            j = T.indices[i]
            val += k(xi, y[j]) * w[j]
        end
        return val

    elseif h.r < θ * euclidean(xi, centers_of_mass[index]) # compress
        return k(xi, centers_of_mass[index]) * sums_w[index]

    else # recurse NOTE: parallelizing here is not as efficient as parallelizing over target points
        l = bh_recursion(getleft(index), k, xi, y, w, sums_w, θ, T, centers_of_mass)
        r = bh_recursion(getright(index), k, xi, y, w, sums_w, θ, T, centers_of_mass)
        return l + r
    end
end

############################# centers of mass ##################################
# this is a weighted sum, could be generalized to incorporate node_sums
function compute_centers_of_mass(w::AbstractVector, x::AbstractVector, T::BallTree)
    D = eltype(x) <: StaticVector ? length(eltype(x)) : length(x[1]) # if x is static vector
    com = [zero(MVector{D, Float64}) for _ in 1:length(T.hyper_spheres)]
    compute_centers_of_mass!(com, w, x, T)
end

function compute_centers_of_mass(F::BarnesHutFactorization, w::AbstractVector)
    compute_centers_of_mass(w, F.y, F.Tree)
end

function compute_centers_of_mass!(com::AbstractVector, w::AbstractVector, x::AbstractVector, T::BallTree)
    abs_w = abs.(w)
    weighted_node_sums!(com, abs_w, x, T)
    sum_w = node_sums(abs_w, T)
    ε = eps(eltype(w)) # ensuring division by zero it not a problem
    @. com ./= sum_w + ε
end

node_sums(x::AbstractVector, T::BallTree) = weighted_node_sums(Ones(length(x)), x, T)
function node_sums!(sums, x::AbstractVector, T::BallTree)
    weighted_node_sums!(sums, Ones(length(x)), x, T)
end

function weighted_node_sums(w::AbstractVector, x::AbstractVector, T::BallTree, index::Int = 1)
    length(x) == 0 && return zero(eltype(x))
    sums = fill(zero(w[1]'x[1]), length(T.hyper_spheres))
    weighted_node_sums!(sums, w, x, T)
end

# NOTE: x should either be vector of numbers or vector of static arrays
function weighted_node_sums!(sums::AbstractVector, w::AbstractVector,
                            x::AbstractVector, T::BallTree, index::Int = 1)
    if isleaf(T.tree_data.n_internal_nodes, index)
        i = get_leaf_range(T.tree_data, index)
        wi, xi = @views w[T.indices[i]], x[T.indices[i]]
        sums[index] = wi'xi
    else
        task = @spawn weighted_node_sums!(sums, w, x, T, getleft(index))
        weighted_node_sums!(sums, w, x, T, getright(index))
        wait(task)
        sums[index] = sums[getleft(index)] + sums[getright(index)]
    end
    return sums
end

# get centers used for compression in taylor factorization
function get_hyper_centers(F::BarnesHutFactorization)
    [hyp.center for hyp in F.Tree.hyper_spheres] # for uniform weights, doesn't loose accuracy in doing this
end


text:
module CleanEOSTable

export cleanEOSTable

using HDF5

function clip(range::Range{Int}, mask::Range{Int})::Range{Int}
    imin = max(start(range), start(mask))
    imax = min(last(range), last(mask))
    imin:imax
end

function linterp{T<:AbstractFloat}(x0::T, y0::T, x1::T, y1::T, x::T,
                                   extrap::Bool=false)::T
    if !extrap
        if x0 < x1
            if !(x0 <= x && x <= x1)
                @show x0, x1, x
            end
            @assert x0 <= x && x <= x1
        elseif x0 > x1
            @assert x1 <= x && x <= x0
        else
            # We expect the input to have a non-empty range
            @assert false
            @assert x == x0
        end
    end
    y = (x1 - x) / (x1 - x0) * y0 + (x - x0) / (x1 - x0) * y1
    if !extrap
        if y0 < y1
            @assert y0 <= y && y <= y1
        elseif y0 > y1
            @assert y1 <= y && y <= y0
        else
            @assert y == y0
        end
    end
    y
end

function linterp(x0::Number, y0::Number, x1::Number, y1::Number, x::Number,
                 extrap::Bool=false)::Number
    linterp(promote(x0, y0, x1, y1, x)..., extrap)
end

function linterp0{T<:AbstractFloat}(arr::Array{T,3}, i::Int, j::Int, k::Int)::T
    arr[i,j,k]
end
function linterp1{T<:AbstractFloat}(arr::Array{T,3}, i::Int, j::Int, k::Int,
                                    x::T)::T
    (1 - x) * linterp0(arr, i, j, k) + x * linterp0(arr, i+1, j, k)
end
function linterp2{T<:AbstractFloat}(arr::Array{T,3}, i::Int, j::Int, k::Int,
                                    x::T, y::T)::T
    (1 - y) * linterp1(arr, i, j, k, x) + y * linterp1(arr, i, j+1, k, x)
end
function linterp3{T<:AbstractFloat}(arr::Array{T,3}, i::Int, j::Int, k::Int,
                                    x::T, y::T, z::T)::T
    @assert x>=0 && x<=1
    @assert y>=0 && y<=1
    @assert z>=0 && z<=1
    (1 - z) * linterp2(arr, i, j, k, x, y) + z * linterp2(arr, i, j, k+1, x, y)
end

"""
Flag inconsistent table entries:
The internal energy needs to be a strictly monotonically increasing function
of the temperature
"""
function flag_inconsistencies{T<:AbstractFloat}(table_logenergy::Array{T, 3})::
    BitArray{3}
    println("Flagging inconsistendies in table...")
    (pointsrho, pointstemp, pointsye) = size(table_logenergy)
    outliers = falses(pointsrho, pointstemp, pointsye)
    for k in 1:pointsye, i in 1:pointsrho
        prev_j = 1
        j = 2
        while j <= pointstemp
            @assert !outliers[i,prev_j,k]
            if table_logenergy[i,j,k] <= table_logenergy[i,prev_j,k]
                # We found an inconsistency. Either this point j or the
                # previous point prev_j is an outlier that needs to be
                # removed.
                irange = clip(i-1:i+1, 1:pointsrho)
                jrange = clip(prev_j-1:j+1, 1:pointstemp)
                krange = clip(k-1:k+1, 1:pointsye)
                avg = mean(table_logenergy[irange, jrange, krange])
                # Assume that the point further from the average is the
                # outlier
                if (abs(table_logenergy[i,prev_j,k] - avg)
                    > abs(table_logenergy[i,j,k] - avg))
                    # Disable the previous point
                    outliers[i,prev_j,k] = true
                    # Next, compare the current point to the one before
                    # the previous point
                    prev_j = prev_j - 1
                    while prev_j >= 1 && outliers[i,prev_j,k]
                        prev_j = prev_j - 1
                    end
                    if prev_j >= 1
                        # keep j
                    else
                        # All points before this are disabled: move on to
                        # the next point
                        prev_j = j
                        j = j + 1
                    end
                else
                    # Disable the current point
                    outliers[i,j,k] = true
                    # keep prev_j
                    j = j + 1
                end
            else
                # Move on to the next point pair
                prev_j = j
                j = j + 1
            end
        end
        # Find the first non-disabled point
        prev_j = 1
        while prev_j <= pointstemp && outliers[i,prev_j,k]
            prev_j = prev_j + 1
        end
        # Check consistency for all following points
        for j = prev_j+1:pointstemp
            if !outliers[i,j,k]
                @assert table_logenergy[i,j,k] > table_logenergy[i,prev_j,k]
                prev_j = j
            end
        end
    end
    noutliers = count(outliers)
    npoints = length(table_logenergy)
    noutliers_percent = round(100 * noutliers / npoints, 2)
    println("    Detected $noutliers outliers ($noutliers_percent%)")
    # println("    Outliers:")
    # n = 0
    # for k in 1:pointsye, j in 2:pointstemp, i in 1:pointsrho
    #     if outlier[i,j,k]
    #         n = n + 1
    #         println("        #$n: [$i,$j,$k]")
    #     end
    # end
    return outliers
end

# Correct outliers
function correct_table{T<:AbstractFloat}(old_table_logenergy::Array{T, 3},
                                         outliers::BitArray{3})::Array{T, 3}
    println("Correcting outliers...")
    @assert size(old_table_logenergy) == size(outliers)
    (pointsrho, pointstemp, pointsye) = size(old_table_logenergy)
    table_logenergy = copy(old_table_logenergy)
    for k in 1:pointsye, j in 1:pointstemp, i in 1:pointsrho
        if outliers[i,j,k]
            # Find two valid neighbouring points
            extrap = false
            prev_j = j - 1
            while prev_j >= 1 && outliers[i,prev_j,k]
                prev_j = prev_j - 1
            end
            next_j = j + 1
            while next_j <= pointstemp && outliers[i,next_j,k]
                next_j = next_j + 1
            end
            if prev_j < 1
                extrap = true
                prev_j = next_j
                next_j = next_j + 1
                while next_j <= pointstemp && outliers[i,next_j,k]
                    next_j = next_j + 1
                end
            end
            if next_j > pointstemp
                extrap = true
                next_j = prev_j
                prev_j = prev_j - 1
                while prev_j >= 1 && outliers[i,prev_j,k]
                    prev_j = prev_j - 1
                end
            end
            @assert prev_j >= 1 && prev_j <= pointstemp
            @assert next_j >= 1 && next_j <= pointstemp
            # Linear interpolation
            avg = linterp(prev_j, table_logenergy[i,prev_j,k],
                          next_j, table_logenergy[i,next_j,k], j, extrap)
            table_logenergy[i,j,k] = avg
        end
    end
    return table_logenergy
end

# Check consistency
function check_consistency{T<:AbstractFloat}(table_logenergy::Array{T, 3})::Void
    println("Checking consistency again...")
    (pointsrho, pointstemp, pointsye) = size(table_logenergy)
    for k in 1:pointsye, j in 2:pointstemp, i in 1:pointsrho
        @assert table_logenergy[i,j,k] > table_logenergy[i,j-1,k]
    end
end

function cleanEOSTable(tablename::AbstractString,
                       new_tablename::AbstractString)
    println("CleanEOSTable")

    println("Reading EOS table \"$tablename\"...")
    table = h5open(tablename, "r")

    # The 3D tables have three indices (rho, temp, ye), in this order
    # (in Fortran, Cactus, and Julia)
    pointsrho = read(table, "pointsrho")[1]
    pointstemp = read(table, "pointstemp")[1]
    pointsye = read(table, "pointsye")[1]

    table_logrho = read(table, "logrho")
    @assert size(table_logrho) == (pointsrho,)
    table_logtemp = read(table, "logtemp")
    @assert size(table_logtemp) == (pointstemp,)
    table_ye = read(table, "ye")
    @assert size(table_ye) == (pointsye,)

    table_logenergy = read(table, "logenergy")
    @assert size(table_logenergy) == (pointsrho, pointstemp, pointsye)
    npoints = length(table_logenergy)
    println("    Read table with ($pointsrho,$pointstemp,$pointsye) = $npoints data points")

    outliers = flag_inconsistencies(table_logenergy)
    new_table_logenergy = correct_table(table_logenergy, outliers)
    check_consistency(new_table_logenergy)

    println("Writing corrected table...")
    new_table = h5open(new_tablename, "w")
    write(new_table, "logenergy", new_table_logenergy)
    close(new_table)

    println("Done.")
end

end # module


text:
using ..Atmos
using ..Atmos: MoistureModel

# Helpers to gather the thermodynamic variables across the DG grid
#

function vars_thermo(atmos::AtmosModel, FT)
    @vars begin
        temp::FT
        pres::FT
        θ_dry::FT
        e_int::FT
        h_tot::FT
        h_int::FT

        moisture::vars_thermo(atmos.moisture, FT)
    end
end
vars_thermo(::MoistureModel, FT) = @vars()
function vars_thermo(m::EquilMoist, FT)
    @vars begin
        q_liq::FT
        q_ice::FT
        q_vap::FT
        θ_vir::FT
        θ_liq_ice::FT
    end
end
num_thermo(bl, FT) = varsize(vars_thermo(bl, FT))
thermo_vars(bl, array) = Vars{vars_thermo(bl, eltype(array))}(array)

# compute thermodynamic variables visitor function, to use with `@visitQ`
function compute_thermo!(atmos::AtmosModel, state, aux, thermo)
    e_tot = state.ρe / state.ρ
    ts = thermo_state(atmos, state, aux)
    e_int = internal_energy(ts)

    thermo.temp = air_temperature(ts)
    thermo.pres = air_pressure(ts)
    thermo.θ_dry = dry_pottemp(ts)
    thermo.e_int = e_int

    R_m = gas_constant_air(ts)
    thermo.h_tot = e_tot + R_m * thermo.temp
    thermo.h_int = e_int + R_m * thermo.temp

    compute_thermo!(atmos.moisture, state, aux, ts, thermo)

    return nothing
end
function compute_thermo!(::MoistureModel, state, aux, ts, thermo)
    return nothing
end
function compute_thermo!(moist::EquilMoist, state, aux, ts, thermo)
    Phpart = PhasePartition(ts)

    thermo.moisture.q_liq = Phpart.liq
    thermo.moisture.q_ice = Phpart.ice
    thermo.moisture.q_vap = vapor_specific_humidity(ts)
    thermo.moisture.θ_vir = virtual_pottemp(ts)
    thermo.moisture.θ_liq_ice = liquid_ice_pottemp(ts)

    return nothing
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
taskSourcePath = "./sourceData"

function log_file(thisMessage::String, returnValue::Int = 0)
## when the log file is first created use append = false . . .
    #    CSV.write("logFile.csv", DataFrame(time = [get_time()], event = [thisMessage], retcode = [returnValue]), append=false)
  if true == taskLoggingOn
    CSV.write("logFile.csv", DataFrame(time = [get_time()], event = [thisMessage], retcode = [returnValue]), append=true)
  end
end

function read_log() logDF = DataFrame()
    logDF = CSV.File("logFile.csv") |> DataFrame
end

function count_log_events()
    logDF = DataFrame()
    logDF = CSV.File("logFile.csv")
    return size(logDF)[1]
end

function process_sourceFolder()
  sourceFiles = read_sourceFolder()
  for sf ∈ sourceFiles

    println("file: $sf")
  end
end

function count_CSV_files()

  sourceFiles = read_sourceFolder()

  occursin.(r".csv", sourceFiles) |>
  collect |>
  sum

end

function count_TXT_files()

  sourceFiles = read_sourceFolder()

  occursin.(r".txt", sourceFiles) |>
  collect |>
  sum

end

function show_CSV_files()

  sourceFiles = read_sourceFolder()

  return occursin.(r".csv", sourceFiles)

end

function show_TXT_files()

  sourceFiles = read_sourceFolder()

  return occursin.(r".txt", sourceFiles)

end

function read_sourceFolder()
    readdir(taskSourcePath)
end

function set_sourcePath(thisPath::String)
    global taskSourcePath = thisPath

    #    taskSourcePath = thisPath
end

function show_sourcePath()
    return taskSourcePath
end


text:
using DataDeps

register(DataDep(
  "nursery",
  "http://archive.ics.uci.edu/ml/datasets/Nursery",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/nursery/nursery.data",
  "8e0389c3dd37590248a921c2726d869ee96b817761a35eb8416afa24f31f931d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=9,
        feature_indices=1:8,
        categoric_indices=1:8,
        footerskip=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    VAE([zlength::Int,p::P] ,e::E ,d::D)

Variational Auto-Encoder.

# Arguments
* `p`: MvNormal prior p(z)
* `zlength`: length of prior MvNormal distribution
* `e`: Encoder p(z|x)
* `d`: Decoder p(x|z)

# Example
Create a vanilla VAE with mappings `f` (R⁵->(R²,R²)) and `g` (R²->(R⁵,R)) that
return mean and variance of appropriate sizes. The output of `f` (and `g`) is
passed to `MvNormal` via the `condition` function of `ConditionalDists`.
```julia-repl
julia> f   = SplitLayer(5, [2,2])
julia> enc = ConditionalMvNormal(f)
julia> g   = SplitLayer(2, [5,1])
julia> dec = ConditionalMvNormal(g)
julia> vae = VAE(2, enc, dec)
VAE:

julia> mean(vae.decoder, mean(vae.encoder, rand(5)))
5×1 Array{Float32,2}:
 -0.26742023
 -0.7905855
 -0.29494995
  0.1694059
  1.123661
```
"""
struct VAE{P<:ContinuousMultivariateDistribution,E<:ConditionalMvNormal,D<:ConditionalMvNormal} <: AbstractVAE
    prior::P
    encoder::E
    decoder::D
end

Flux.@functor VAE

function Flux.trainable(m::VAE)
    (encoder=m.encoder, decoder=m.decoder)
end

function VAE(zlength::Int, enc::ConditionalMvNormal, dec::ConditionalMvNormal)
    W = first(Flux.params(enc))
    μ = fill!(similar(W, zlength), 0)
    σ = fill!(similar(W, zlength), 1)
    prior = DistributionsAD.TuringMvNormal(μ, σ)
    VAE(prior, enc, dec)
end

"""
    elbo(m::AbstractVAE, x::AbstractArray; β=1)

Evidence lower boundary of the VAE model. `β` scales the KLD term.
"""
function elbo(m::AbstractVAE, x::AbstractArray; β=1)
    # sample latent
    z = rand(m.encoder, x)

    # reconstruction error
    llh = mean(logpdf(m.decoder, x, z))

    # KLD with `condition`ed encoder
    kld = mean(kl_divergence(condition(m.encoder, x), m.prior))

    llh - β*kld
end

"""
    mmd_mean(m::AbstractVAE, x::AbstractArray, k::Kernel)

Maximum mean discrepancy of a VAE given data `x` and kernel `k`.
Uses mean of encoded data.
"""
function mmd_mean(m::AbstractVAE, x::AbstractArray, k::AbstractKernel = GaussianKernel())
    mmd(k, mean(m.encoder,x), rand(m.prior,size(x,2)))
end

"""
    mmd_rand(m::AbstractVAE, x::AbstractArray, k::Kernel)

Maximum mean discrepancy of a VAE given data `x` and kernel `k`.
Samples from the encoder.
"""
function mmd_rand(m::AbstractVAE, x::AbstractArray, k::AbstractKernel = GaussianKernel())
    mmd(k, rand(m.encoder,x), rand(m.prior,size(x,2)))
end

function Base.show(io::IO, m::AbstractVAE)
    p = repr(m.prior)
    p = sizeof(p)>70 ? "($(p[1:70-3])...)" : p
    e = repr(m.encoder)
    e = sizeof(e)>70 ? "($(e[1:70-3])...)" : e
    d = repr(m.decoder)
    d = sizeof(d)>70 ? "($(d[1:70-3])...)" : d
    msg = """$(nameof(typeof(m))):
     prior   = $(p)
     encoder = $(e)
     decoder = $(d)
    """
    print(io, msg)
end


text:
using VLMinorityGameKit
using BSON
using DataFrames
using Plots

# helper functions -
function load_simulation_data(path_to_data::String, keyname::String)::Dict
    d = BSON.load(path_to_data)
    return d[Symbol(keyname)]
end

function compute_return_array(path_to_sim_data_file::String)::Array{Float64,2}

    keyname = "data"
    simulation_dictionary = load_simulation_data(path_to_sim_data_file, keyname)
    number_of_samples = simulation_dictionary["number_of_samples"]
    number_of_timesteps = simulation_dictionary["number_of_timesteps"]
    simulation_result_array = simulation_dictionary["simulation_result_array"]

    # initialize some space, collect the data and make the plot
    sim_return_array = Array{Float64,2}(undef, (number_of_timesteps - 1), number_of_samples)
    for sample_index = 1:number_of_samples

        # grab the sim -
        M = simulation_result_array[sample_index]

        # grab the market data, and get the price -
        P = M[!,:price]

        # the price data into the price_array -
        for time_step_index = 1:(number_of_timesteps - 1)

            # get prices -
            today_price = P[time_step_index]
            tomorrow_price = P[time_step_index + 1]

            # compute r -
            r = log(tomorrow_price / today_price)

            # package -
            sim_return_array[time_step_index, sample_index] = r
        end
    end

    return sim_return_array
end

# path to simulation file -
path_to_sim_data_file = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/results/SIM-T-M10-S5-NA1001-L3.bson"
L3_return_array = compute_return_array(path_to_sim_data_file)

# test the autocorrelation -
(number_of_timesteps, number_of_samples) = size(L3_return_array)
lag_array = range(0,stop=(number_of_timesteps - 1),step=1) |> collect
acr = compute_autocorrelation_array(L3_return_array[:,2], lag_array)


text:
import ...timeevolution: dschroedinger, dschroedinger_dynamic, check_schroedinger

"""
    stochastic.schroedinger(tspan, state0, H, Hs[; fout, ...])

Integrate stochastic Schrödinger equation.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state as Ket.
* `H`: Deterministic part of the Hamiltonian.
* `Hs`: Stochastic part(s) of the Hamiltonian (either an operator or a vector
        of operators).
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `normalize_state=false`: Specify whether or not to normalize the state after
        each time step taken by the solver.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector;
                fout::Union{Function,Nothing}=nothing,
                normalize_state::Bool=false,
                calback=nothing,
                kwargs...) where {B<:Basis,T<:Ket{B}}
    tspan_ = convert(Vector{Float64}, tspan)

    n = length(Hs)
    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    # TODO: replace checks by dispatch
    for h=Hs
        @assert isa(h, AbstractOperator{B,B})
    end

    dschroedinger_determ(t::Float64, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi)
    dschroedinger_stoch(dx::DiffArray,
            t::Float64, psi::T, dpsi::T, n::Int) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n)

    if normalize_state
        norm_func(u::Vector{ComplexF64}, t::Float64, integrator) = normalize!(u)
        ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                 func_everystep=true,
                 func_start=false)
    else
        ncb = nothing
    end

    integrate_stoch(tspan_, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n;
                    ncb=ncb,
                    kwargs...)
end
schroedinger(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, Hs::AbstractOperator{B,B}; kwargs...) where B<:Basis = schroedinger(tspan, psi0, H, [Hs]; kwargs...)

"""
    stochastic.schroedinger_dynamic(tspan, state0, fdeterm, fstoch[; fout, ...])

Integrate stochastic Schrödinger equation with dynamic Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state.
* `fdeterm`: Function `f(t, psi, u) -> H` returning the deterministic
        (time- or state-dependent) part of the Hamiltonian.
* `fstoch`: Function `f(t, psi, u, du) -> Hs` returning a vector that
        contains the stochastic terms of the Hamiltonian.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch`. If unset, the number is calculated automatically
        from the function output.
        NOTE: Set this number if you want to avoid an initial calculation of
        the function output!
* `normalize_state=false`: Specify whether or not to normalize the state after
        each time step taken by the solver.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Function;
                fout::Union{Function,Nothing}=nothing, noise_processes::Int=0,
                normalize_state::Bool=false,
                kwargs...) where T<:Ket
    tspan_ = convert(Vector{Float64}, tspan)

    if noise_processes == 0
        fs_out = fstoch(0.0, psi0)
        n = length(fs_out)
    else
        n = noise_processes
    end

    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    dschroedinger_determ(t::Float64, psi::T, dpsi::T) = dschroedinger_dynamic(t, psi, fdeterm, dpsi)
    dschroedinger_stoch(dx::DiffArray,
            t::Float64, psi::T, dpsi::T, n::Int) =
        dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n)

    if normalize_state
        norm_func(u::Vector{ComplexF64}, t::Float64, integrator) = normalize!(u)
        ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                 func_everystep=true,
                 func_start=false)
    else
        ncb = nothing
    end

    integrate_stoch(tspan, dschroedinger_determ, dschroedinger_stoch, x0, state,
            dstate, fout, n;
            ncb=ncb,
            kwargs...)
end

function dschroedinger_stochastic(dx::D, psi::T1, Hs::Vector{T2},
            dpsi::T1, index::Int) where {D<:Vector{ComplexF64},B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}}
    recast!(dx, dpsi)
    dschroedinger(psi, Hs[index], dpsi)
end
function dschroedinger_stochastic(dx::Array{ComplexF64, 2}, psi::T1, Hs::Vector{T2},
            dpsi::T1, n::Int) where {B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}}
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(dx_i, dpsi)
        dschroedinger(psi, Hs[i], dpsi)
        recast!(dpsi, dx_i)
    end
end
function dschroedinger_stochastic(dx::DiffArray,
            t::Float64, psi::T, f::Function, dpsi::T, n::Int) where T<:Ket
    ops = f(t, psi)
    if QO_CHECKS[]
        @inbounds for h=ops
            check_schroedinger(psi, h)
        end
    end
    dschroedinger_stochastic(dx, psi, ops, dpsi, n)
end

recast!(psi::Ket, x::SubArray{ComplexF64, 1}) = (x .= psi.data)
recast!(x::SubArray{ComplexF64, 1}, psi::Ket) = (psi.data = x)


text:
using Test

function list_examples()
    res = String[]
    basedir = joinpath(@__DIR__, "..", "examples")
    for file in readdir(basedir)
        fullpath = joinpath(basedir, file)
        if isfile(fullpath) && endswith(file, ".jl")
            push!(res, fullpath)
        end
    end
    res
end

function example_tags(fullpath)
    open(fullpath) do fp
        for l in readlines(fp)
            regex = "^#src tags:"
            isnothing(match(Regex(regex), l)) || return split(l[length(regex):end])
        end
        return Vector{String}()
    end
end

@testset "Run examples" begin
    for file in list_examples()
        "long" in example_tags(file) && continue
        @testset "$(file)" begin
            include(file)
        end
    end
end


text:
# Example from Tao et al.

using JLD2;
using Distributed;
using Printf;
using IterTools;
using Random;

@everywhere include("runit_bis.jl");
@everywhere include("../thresholds.jl");
include("../experiment_helpers.jl");

rng = MersenneTwister(1234)

# setup

dist = Gaussian();
dim = 12
num_arms = 20

# Box-Mueller
arms = Vector{Float64}[]
for k = 1:num_arms
    v = zeros(dim)
    for i = 1:dim
        v[i] = randn()
    end
    normalise = norm(v)
    v ./= norm(v)
    push!(arms, v)
end

μ = arms[1] + 0.01 * (arms[2] - arms[1]);

pep = LinearBestArm(dist, arms);

K = length(arms)
spanning_weights = ones(K)
spanning_weights[K] = 0.0
spanning_weights /= sum(spanning_weights)

srs = [
    ConvexGame(CTracking),
    LearnerK(CTracking),
    RoundRobin(),
    FixedWeights(spanning_weights),
    Static("G"),
    Static("XY"),
    XYAdaptive(),
    LinGapE(),
]
# srs = [LinGapE()]

# δs = (0.1, 0.01, 0.0001);
δs = (0.01,);
βs = GK16.(δs);

repeats = 100;
seed = 1234;

# compute

@time data = map(  # TODO: replace by pmap (it is easier to debug with map)
    ((sr, i),) -> runit(seed + i, sr, μ, pep, βs, δs),
    Iterators.product(srs, 1:repeats),
);

dump_stats(pep, μ, δs, βs, srs, data, repeats);

# save

@save isempty(ARGS) ? "experiment_bai2.dat" : ARGS[1] dist μ pep srs data δs βs repeats seed

# visualise by loading viz_bai2.jl


text:
module IntervalAnalysis
using Plots
using Base
using Random
using Symbolics
using Documenter

export Interval,*,+,-,subdivide

"""
    Interval

Defines a single interval between two real numbers.

# Examples
```jldoctest
julia> Interval(-1,1)
-1.0 ←-→ 1.0
```
"""
struct Interval
	"Lower interval value"
	l::Real
	"Upper interval value"
	u::Real
	Interval(l,u) = l > u ? new(u,l) : new(l,u)

end

Base.show(io::IO, x::Interval) = print(io,round(x.l,digits=6)," ←-→ ",round(x.u,digits=6))

"""

	+(x₁::Interval, x₂::Interval)

Adds two intervals

# Examples
```jldoctest
julia> Interval(-1,1) + Interval(-2,4)
-3.0 ←-→ 5.0
```
"""
function Base.:+(x₁::Interval, x₂::Interval)
	Interval(x₁.l + x₂.l,x₁.u + x₂.u)
    end

function Base.:+(x₁::Real, x₂::Interval)
	Interval(x₁ + x₂.l,x₁ + x₂.u)
    end

function Base.:+(x₁::Interval, x₂::Real)
	Interval(x₁.l + x₂,x₁.u + x₂)
    end

 """

	-(x₁::Interval, x₂::Interval)

Subtracts one interval from another

# Examples
```jldoctest
julia> Interval(-1,1) - Interval(-2,4)
-5.0 ←-→ 3.0
```
"""
function Base.:-(x₁::Interval, x₂::Interval)
	Interval(x₁.l - x₂.u,x₁.u - x₂.l)
    end

function Base.:-(x::Interval)
	Interval(-x.l,-x.u)
    end

function Base.:-(x₁::Interval, x₂::Real)
	Interval(x₁.l - x₂,x₁.u - x₂)
    end

    function Base.:-(x₁::Real, x₂::Interval)
	Interval(x₁ - x₂.l,x₁ - x₂.u)
    end

"""

	*(x₁::Interval, x₂::Interval)

Multiplies two intervals

# Examples
```jldoctest
julia> Interval(-1,1) * Interval(-2,4)
-4.0 ←-→ 4.0
```
"""
function Base.:*(x₁::Interval, x₂::Interval)
	if 0 ≤ x₁.l
		if 0 ≤ x₂.l
			return Interval(x₁.l*x₂.l,x₁.u*x₂.u)
		elseif x₂.l < 0 < x₂.u
			return Interval(x₁.u*x₂.l,x₁.u*x₂.u)
		elseif x₂.u ≤ 0
			return Interval(x₁.u*x₂.l,x₁.l*x₂.u)
		end
	elseif x₁.l < 0 < x₁.u
		if 0 ≤ x₂.l
			return Interval(x₁.l*x₂.u,x₁.u*x₂.u)
		elseif x₂.u ≤ 0
			return Interval(x₁.u*x₂.l,x₁.l*x₂.l)
		elseif x₂.l < 0 < x₂.u
			return Interval(min(x₁.l*x₂.u,x₁.u*x₂.l),max(x₁.l*x₂.l,x₁.u*x₂.u))
		end
	elseif x₁.u ≤ 0
		if 0 ≤ x₂.l
			return Interval(x₁.l*x₂.u,x₁.u*x₂.l)
		elseif x₂.l < 0 < x₂.u
			return Interval(x₁.l*x₂.u,x₁.l*x₂.l)
		elseif x₂.u ≤ 0
			return Interval(x₁.u*x₂.u,x₁.l*x₂.l)
		end
	end
end

function Base.:*(x₁::Interval,x₂::Real)
	Interval(x₁.l*x₂,x₁.u*x₂)
end

function Base.:*(x₁::Real,x₂::Interval)
	Interval(x₁*x₂.l,x₁*x₂.u)
end

"""

	/(x₁::Interval, x₂::Interval)

Divides two intervals

# Examples
```jldoctest
julia> Interval(-1,1) / Interval(-2,4)
-0.5 ←-→ 0.5
```
"""
function Base.:/(x₁::Interval, x₂::Interval)
	x₁ * Interval(1/x₂.u,1/x₂.l)
    end

function Base.:/(x₁::Real, x₂::Interval)
	Interval(x₁/x₂.u,x₁/x₂.l)
    end

function Base.:/(x₁::Interval, x₂::Real)
	Interval(x₁.l/x₂,x₁.u/x₂)
    end

function w(x::Interval)
	x.u-x.l
end

function w(X::Vector{Interval})
	widths = w.(X)
	d,i = findmax(widths)
end

function m(x::Interval)
	(x.u+x.l)/2
end

function m(X::Vector{Interval})
	m.(X)
end

function Base.sqrt(x::Interval)
	Interval(sqrt(x.l),sqrt(x.u))
end

function Base.:^(x::Interval,p::Integer)
	if x.l > 0 || p%2 != 0
		return Interval(x.l^p,x.u^p)
	elseif x.u < 0 && p%2 == 0
		return Interval(x.u^p,x.l^p)
	else
		return Interval(0,abs(x)^p)
	end
end

function Base.exp(x::Interval)
	Interval(exp(x.l),exp(x.u))
end

function Base.log(x::Interval)
	Interval(log(x.l),log(x.u))
end

function Base.log10(x::Interval)
	Interval(log10(x.l),log10(x.u))
end

function Base.log2(x::Interval)
	Interval(log2(x.l),log2(x.u))
end

function Base.sin(x::Interval)
	if x.l ≥ -(π/2) && x.u ≤ π/2
			return Interval(sin(x.l),sin(x.u))
	end
	return Interval(-1,1)
end

function Base.cos(x::Interval)
	if x.l ≥ 0 && x.u ≤ π
			return Interval(cos(x.l),cos(x.u))
	end
	return Interval(-1,1)
end

function Base.abs(x::Interval)
	max(abs(x.l),abs(x.u))
end

function qdist(x₁::Interval,x₂::Interval)
	max(abs(x₁.l-x₂.l),abs(x₁.u-x₂.u))
end

"""

	subdivide(x::Interval,N::Integer)

Subdivides an interval into N intervals.

# Examples
```jldoctest
julia> subdivide(Interval(0,1),5)
[0.0 ←-→ 0.2, 0.2 ←-→ 0.4, 0.4 ←-→ 0.6, 0.6 ←-→ 0.8, 0.8 ←-→ 1.0]
```
"""
function subdivide(x::Interval,N::Integer)
	A = Array{Interval}(undef,N)
	w_x = w(x)/N
	for j ∈ 1:N
		A[j] = Interval(x.l+(j-1)*w_x,x.l+j*w_x)
	end
	return A
end

function subdivide_vector(X::Vector{Interval},N::Integer)
	A = Array{Vector{Interval}}(undef,length(X))
	for i ∈ 1:length(X)
		A[i] = subdivide(X[i],N)
	end
	return A
end

function bisect_box(X::Vector{Interval})
	# bisects a box over the largest axis
	d,i = w(X)
	X_s = subdivide(X[i],2)
	X₁ = copy(X)
	X₁[i] = X_s[1]
	X₂ = copy(X)
	X₂[i] = X_s[2]
	return X₁,X₂
end

function hull(X::Vector{Interval})
	l = [x.l for x ∈ X]
	u = [x.u for x ∈ X]
	Interval(minimum(l),maximum(u))
end

function refinement(x,f,N)
	X = subdivide(x,N)
	Y = f.(X)
	hull(Y)
end

function plot_refinement(x,f,∇f,N)
	plot(f,x.l,x.u,c="black",linewidth=3,label="")
	X = subdivide(x,N)
	Y = f.(X)
	Yₘᵥ = mean_value_form.(f,∇f,X)
	for i ∈ 1:length(X)
		if i == 1
			plot!([X[i].l,X[i].u],[Y[i].l,Y[i].l],fillrange=[Y[i].u,Y[i].u],linewidth=0,c="black",label="Interval Extension",fillalpha=0.15)
			plot!([X[i].l,X[i].u],[Yₘᵥ[i].l,Yₘᵥ[i].l],fillrange=[Yₘᵥ[i].u,Yₘᵥ[i].u],linewidth=0,c="red",label="Mean Value Form",fillalpha=0.25)
		else
			plot!([X[i].l,X[i].u],[Y[i].l,Y[i].l],fillrange=[Y[i].u,Y[i].u],linewidth=0,c="black",label="",fillalpha=0.15)
			plot!([X[i].l,X[i].u],[Yₘᵥ[i].l,Yₘᵥ[i].l],fillrange=[Yₘᵥ[i].u,Yₘᵥ[i].u],linewidth=0,c="red",label="",fillalpha=0.25)
		end

	end
	display(plot!(dpi=600))
end

function skelboe_moore(f,X)

	fᵤ = f(X).u # smallest upper bound (AKA the only upper bound at this stage)

	# creating a store of intervals, mean values, and interval values
	L = [Dict("X"=>X,"f_m"=>f(m(X)),"F"=>f(X))]
	Lᵣ = []
	it_count = 0

	# until tolerance...
	for it ∈ 1:30

		# initialise list of intervals to sack off
		del_index =[]
		for i in 1:length(L)
			# if the lowest interval value is above the smallest upper bound..
			if L[i]["F"].l > fᵤ
				# remove interval, we know global minimum is not here...
				push!(del_index,i)
			end
		end

		# delete all the ones we said we would
		for d_i ∈ del_index
			push!(Lᵣ,L[d_i])
		end
		deleteat!(L,del_index)

		# pick the remaining interval with smallest mean value
		s = argmin([l["f_m"] for l ∈ L])

		# divide the box into two
		X_s = bisect_box(L[s]["X"])

		# get rid of the containing interval
		deleteat!(L,s)

		for x ∈ X_s
			# evaluate interval values
			F = f(x)

			# check if we find a smallest upper bound
			if F.u < fᵤ
				fᵤ = F.u
			end
			# add new interval information to store
			push!(L,Dict("X"=>x,"f_m"=>f(m(x)),"F"=>F))

		end

	end
	return L,Lᵣ
end

function plot_boxes(L,Lᵣ)
	for l ∈ L
		x1 = l["X"][1]
		x2 = l["X"][2]
		plot!(Shape([x1.l,x1.l,x1.u,x1.u] , [x2.u,x2.l,x2.l,x2.u]),label="",linewidth=1,alpha=0,linealpha=1)
	end
	for l ∈ Lᵣ
		x1 = l["X"][1]
		x2 = l["X"][2]
		plot!(Shape([x1.l,x1.l,x1.u,x1.u] , [x2.u,x2.l,x2.l,x2.u]),c="black",alpha=0.55,label="",linewidth=1)
	end
end

function mean_value_form(f,∇f,X)
	if typeof(X) == Interval
		D = ∇f(X)[1]
	else
		D = ∇f(X)
	end
	mean = m(X)
	fₘ = f(mean)
	if typeof(X) == Interval
		return fₘ + D*(X-mean)
	else
		return fₘ + sum(D[i]*(X[i]-mean[i]) for i ∈ 1:length(X))
	end

end

function build_gradient(f,len_x)
	if len_x > 1
		@variables x[1:len_x]
		∇f = Symbolics.gradient(f(x),[x[i] for i ∈ 1:len_x])
		∇f_expr = build_function(∇f, x)
		∇f = eval(∇f_expr[1])
		return ∇f
	else
		@variables x
		∇f = Symbolics.gradient(f(x),[x])
		∇f_expr = build_function(∇f, x)
		∇f = eval(∇f_expr[1])
		return ∇f
	end

end

function g(x1,x2)
	return x1*exp(x1+x2^2)-x2^2
end

function g_vec(X)
	return g(X[1],X[2])
end

# f(x) = x^3 - x^2
# ∇f = build_gradient(f,1)
# plot_refinement(Interval(0,1),f,∇f,20)
# print(Interval(-1,1) / Interval(-2,4))

end


text:
"""
    inception_a(inplanes, pool_proj)

Create an Inception-v3 style-A module
(ref: Fig. 5 in [paper](https://arxiv.org/abs/1512.00567v3)).

# Arguments
- `inplanes`: number of input feature maps
- `pool_proj`: the number of output feature maps for the pooling projection
"""
function inception_a(inplanes, pool_proj)
  branch1x1 = Chain(conv_bn((1, 1), inplanes, 64)...)

  branch5x5 = Chain(conv_bn((1, 1), inplanes, 48)...,
                    conv_bn((5, 5), 48, 64; pad = 2)...)

  branch3x3 = Chain(conv_bn((1, 1), inplanes, 64)...,
                    conv_bn((3, 3), 64, 96; pad = 1)...,
                    conv_bn((3, 3), 96, 96; pad = 1)...)

  branch_pool = Chain(MeanPool((3, 3), pad = 1, stride = 1),
                      conv_bn((1, 1), inplanes, pool_proj)...)

  return Parallel(cat_channels,
                  branch1x1, branch5x5, branch3x3, branch_pool)
end

"""
    inception_b(inplanes)

Create an Inception-v3 style-B module
(ref: Fig. 10 in [paper](https://arxiv.org/abs/1512.00567v3)).

# Arguments
- `inplanes`: number of input feature maps
"""
function inception_b(inplanes)
  branch3x3_1 = Chain(conv_bn((3, 3), inplanes, 384; stride = 2)...)

  branch3x3_2 = Chain(conv_bn((1, 1), inplanes, 64)...,
                      conv_bn((3, 3), 64, 96; pad = 1)...,
                      conv_bn((3, 3), 96, 96; stride = 2)...)

  branch_pool = Chain(MaxPool((3, 3), stride = 2))

  return Parallel(cat_channels,
                  branch3x3_1, branch3x3_2, branch_pool)
end

"""
    inception_c(inplanes, inner_planes, n = 7)

Create an Inception-v3 style-C module
(ref: Fig. 6 in [paper](https://arxiv.org/abs/1512.00567v3)).

# Arguments
- `inplanes`: number of input feature maps
- `inner_planes`: the number of output feature maps within each branch
- `n`: the "grid size" (kernel size) for the convolution layers
"""
function inception_c(inplanes, inner_planes, n = 7)
  branch1x1 = Chain(conv_bn((1, 1), inplanes, 192)...)

  branch7x7_1 = Chain(conv_bn((1, 1), inplanes, inner_planes)...,
                      conv_bn((1, n), inner_planes, inner_planes; pad = (0, 3))...,
                      conv_bn((n, 1), inner_planes, 192; pad = (3, 0))...)

  branch7x7_2 = Chain(conv_bn((1, 1), inplanes, inner_planes)...,
                      conv_bn((n, 1), inner_planes, inner_planes; pad = (3, 0))...,
                      conv_bn((1, n), inner_planes, inner_planes; pad = (0, 3))...,
                      conv_bn((n, 1), inner_planes, inner_planes; pad = (3, 0))...,
                      conv_bn((1, n), inner_planes, 192; pad = (0, 3))...)

  branch_pool = Chain(MeanPool((3, 3), pad = 1, stride=1),
                      conv_bn((1, 1), inplanes, 192)...)

  return Parallel(cat_channels,
                  branch1x1, branch7x7_1, branch7x7_2, branch_pool)
end

"""
    inception_d(inplanes)

Create an Inception-v3 style-D module
(ref: [pytorch](https://github.com/pytorch/vision/blob/6db1569c89094cf23f3bc41f79275c45e9fcb3f3/torchvision/models/inception.py#L322)).

# Arguments
- `inplanes`: number of input feature maps
"""
function inception_d(inplanes)
  branch3x3 = Chain(conv_bn((1, 1), inplanes, 192)...,
                    conv_bn((3, 3), 192, 320; stride = 2)...)

  branch7x7x3 = Chain(conv_bn((1, 1), inplanes, 192)...,
                      conv_bn((1, 7), 192, 192; pad = (0, 3))...,
                      conv_bn((7, 1), 192, 192; pad = (3, 0))...,
                      conv_bn((3, 3), 192, 192; stride = 2)...)

  branch_pool = Chain(MaxPool((3, 3), stride=2))

  return Parallel(cat_channels,
                  branch3x3, branch7x7x3, branch_pool)
end

"""
    inception_e(inplanes)

Create an Inception-v3 style-E module
(ref: Fig. 7 in [paper](https://arxiv.org/abs/1512.00567v3)).

# Arguments
- `inplanes`: number of input feature maps
"""
function inception_e(inplanes)
  branch1x1 = Chain(conv_bn((1, 1), inplanes, 320)...)

  branch3x3_1 = Chain(conv_bn((1, 1), inplanes, 384)...)
  branch3x3_1a = Chain(conv_bn((1, 3), 384, 384; pad = (0, 1))...)
  branch3x3_1b = Chain(conv_bn((3, 1), 384, 384; pad = (1, 0))...)

  branch3x3_2 = Chain(conv_bn((1, 1), inplanes, 448)...,
                      conv_bn((3, 3), 448, 384; pad = 1)...)
  branch3x3_2a = Chain(conv_bn((1, 3), 384, 384; pad = (0, 1))...)
  branch3x3_2b = Chain(conv_bn((3, 1), 384, 384; pad = (1, 0))...)

  branch_pool = Chain(MeanPool((3, 3), pad = 1, stride = 1),
                      conv_bn((1, 1), inplanes, 192)...)

  return Parallel(cat_channels,
                  branch1x1,
                  Chain(branch3x3_1,
                        Parallel(cat_channels,
                                  branch3x3_1a, branch3x3_1b)),

                  Chain(branch3x3_2,
                        Parallel(cat_channels,
                                  branch3x3_2a, branch3x3_2b)),
                  branch_pool)
end

"""
    inception3(; nclasses = 1000)

Create an Inception-v3 model ([reference](https://arxiv.org/abs/1512.00567v3)).

# Arguments
- `nclasses`: the number of output classes

!!! warning
    `inception3` does not currently support pretrained weights.
"""
function inception3(; nclasses = 1000)
  layer = Chain(Chain(conv_bn((3, 3), 3, 32; stride = 2)...,
                      conv_bn((3, 3), 32, 32)...,
                      conv_bn((3, 3), 32, 64; pad = 1)...,
                      MaxPool((3, 3), stride = 2),
                      conv_bn((1, 1), 64, 80)...,
                      conv_bn((3, 3), 80, 192)...,
                      MaxPool((3, 3), stride = 2),
                      inception_a(192, 32),
                      inception_a(256, 64),
                      inception_a(288, 64),
                      inception_b(288),
                      inception_c(768, 128),
                      inception_c(768, 160),
                      inception_c(768, 160),
                      inception_c(768, 192),
                      inception_d(768),
                      inception_e(1280),
                      inception_e(2048)),
                Chain(AdaptiveMeanPool((1, 1)),
                      Dropout(0.2),
                      flatten,
                      Dense(2048, nclasses)))

  return layer
end

"""
    Inception3(; pretrain = false, nclasses = 1000)

Create an Inception-v3 model ([reference](https://arxiv.org/abs/1512.00567v3)).
See also [`inception3`](#).

# Arguments
- `pretrain`: set to `true` to load the pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `Inception3` does not currently support pretrained weights.
"""
struct Inception3
  layers
end

function Inception3(; pretrain = false, nclasses = 1000)
  layers = inception3(nclasses = nclasses)
  pretrain && loadpretrain!(layers, "Inception3")

  Inception3(layers)
end

@functor Inception3

(m::Inception3)(x) = m.layers(x)

backbone(m::Inception3) = m.layers[1]
classifier(m::Inception3) = m.layers[2]


text:
## Utilities

function vswap!{T <: Vector}(v1::T, v2::T, idx::Integer)
    val = v1[idx]
    v1[idx] = v2[idx]
    v2[idx] = val
end

function inmap{T}(v::T, c::Vector{T}, from::Integer = 1, to::Integer = length(c))
    exists = 0
    for j in from:to
        if exists == 0 && v == c[j]
            exists = j
        end
    end
    return exists
end

function swap!{T <: Vector}(v::T, from::Int, to::Int)
    val = v[from]
    v[from] = v[to]
    v[to] = val
end

# Inverts an objective funtion
function inverseFunc(f::Function)
  function fitnessFunc{T <: Vector}(x::T)
    return 1.0/(f(x)+eps())
  end
  return fitnessFunc
end

# Collecting interim values
function keep(interim, v, vv, col)
    if interim
        if !haskey(col, v)
            col[v] = typeof(vv)[]
        end
        push!(col[v], vv)
    end
end

# ReBalance a Permuation to start at 1
function rebalTo1{T <: Integer}(v::Vector{T})
  n = length(v)
  k = inmap(one(T),v)
  v = v[vcat(k:n,1:k-1)]
  return v
end

#Utilities
function clustmap{T}(a::AbstractArray{T})
    d = Dict{T,Vector{Int}}()
    for i = 1 : length(a)
        @inbounds k = a[i]
        if !haskey(d, k)
            d[k] = findin(a,k)
        end
    end
    return d
end

function fillZero!{T <: Integer}(v::Vector{T})
  k = 1
  for i in 1:length(v)
    if v[i] == zero(T)
      v[i] = k
      k += 1
    end
  end
  return v
end

#Graph Utilities
#check if two nodes are connected
function independent(G::Graph,dnode::Int,anode::Int)
  valid = true
  for vset in connected_components(G)
    if valid && in(dnode,vset) && in(anode,vset)
      valid = false
    end
  end
  return valid
end


text:
using SHA: sha256

"""
Double sha256 function
"""
function hash256(x::Vector{UInt8})
    return sha256(sha256(x))
end


text:
"""
    mJy_to_W( c::Cosmology.AbstractCosmology, S::Unitful.AbstractQuantity, z::Real )

Converts synchrotron flux density `S` in `[mJy]` to `[W/Hz]` for a given redshift `z` and cosmology `c`.
Conserves unit information.
"""
function mJy_to_W( c::Cosmology.AbstractCosmology, S::Unitful.AbstractQuantity, z::Real )

    # Luminsity Distance in Mpc
    dL = luminosity_dist(c, z)

    # P(ν)= S(ν) * 4π * dL^2
    # dL is the luminosity distance
    return S * 4π * dL^2 |> u"W/Hz"
end

"""
    mJy_to_W( c::Cosmology.AbstractCosmology, S::Real, z::Real )

Converts synchrotron flux density `S` in `[mJy]` to `[W/Hz]` for a given redshift `z` and cosmology `c`.
Deletes unit information.
"""
mJy_to_W( c::Cosmology.AbstractCosmology, S::Real, z::Real ) = mJy_to_W(c, S * 1.0u"mJy", z) |> ustrip

"""
    mJy_to_W(S::Union{Real, Unitful.AbstractQuantity}, h::AbstractGadgetHeader)

Converts synchrotron flux density `S` in `[mJy]` to `[W/Hz]` for a given redshift and cosmology taken from `AbstractGadgetHeader`.
"""
mJy_to_W(S::Union{Real, Unitful.AbstractQuantity}, h::AbstractGadgetHeader) = mJy_to_W( cosmology(h), S, h.z )