text:
export EnsembleState, EnsembleStateMeas, norm, sort, sort!, member, deviation!, deviation

import Base: size, +, -, sort, sort!

import Statistics: mean, cov

import LinearAlgebra: norm

# """
#     EnsembleState
#
# An structure  for Ensemble.
#
# Fields:
# - 'S' : Array of the different members
# """
struct EnsembleState{Nx, Ne}
    " Array of the different ensemble members"
    S::Array{Float64,2}
end

function EnsembleState(Nx::Int64, Ne::Int64)
    return EnsembleState{Nx, Ne}(zeros(Nx,Ne))
end

size(ens::EnsembleState{Nx, Ne}) where {Nx, Ne} = (Nx, Ne)

function EnsembleState(States::Array{Float64,2})
    Nx, Ne = size(States)
    return EnsembleState{Nx, Ne}(States)
end

# """
#     Define addition of two EnsembleState
# """
function (+)(A::EnsembleState{Nx, Ne}, B::EnsembleState{Nx, Ne}) where {Nx, Ne}
    C = deepcopy(A)
    C.S .+= B.S
    return C
end

# """
#     Define addition of an Array and an EnsembleState
# """
function (+)(A::EnsembleState{Nx, Ne}, B::Array{Float64,1}) where {Nx, Ne}
    C = deepcopy(A)
    C .+= B
    return C
end

# """
#     Define substraction of two EnsembleState
# """
function (-)(A::EnsembleState{Nx, Ne}, B::EnsembleState{Nx, Ne}) where {Nx, Ne}
    C = deepcopy(A)
    C.S .-= B.S
    return C
end

# """
#     Define substraction of an Array from an EnsembleState
# """
function (-)(A::EnsembleState{Nx, Ne}, B::Array{Float64,1}) where {Nx, Ne}
    C = deepcopy(A)
    C .-= B
    return C
end

# """
#     Define norm of an EnsembleState variable, equivalent to norm(hcat(ens))
# """

norm(A::EnsembleState{Nx, Ne}) where {Nx, Ne} = norm(A.S)

sort(A::EnsembleState{Nx, Ne}, dims::Int64) where {Nx, Ne} = sort(A.S, dims = dims)
sort!(A::EnsembleState{Nx, Ne}, dims::Int64) where {Nx, Ne} = sort!(A.S, dims = dims)

member(ens::EnsembleState{Nx, Ne}, idx::Int64)  where{Nx, Ne}  = ens.S[:,idx]

mean(ens::EnsembleState{Nx, Ne}) where {Nx, Ne} = mean(ens.S, dims=2)[:,1]
# Rescaled by 1/(Ne-1)
cov(ens::EnsembleState{Nx, Ne}) where {Nx, Ne} = Statistics.cov(ens.S, dims=2)

function deviation!(fluc::EnsembleState{Nx, Ne}, ens::EnsembleState{Nx, Ne}) where {Nx, Ne}
S = mean(ens)
@. fluc.S = ens.S - S
return fluc
end

deviation(ens::EnsembleState{Nx, Ne}) where {Nx, Ne} = deviation!(EnsembleState(Nx, Ne), ens)

struct EnsembleStateMeas{Nx, Ny, Ne}
    " Array of state for the different ensemble members"
    state::EnsembleState{Nx, Ne}

    " Array of measurement for the different ensemble members"
    meas::EnsembleState{Ny, Ne}
end

EnsembleStateMeas(Nx, Ny, Ne) = EnsembleStateMeas{Nx, Ny, Ne}(EnsembleState(Nx, Ne), EnsembleState(Ny, Ne))

function EnsembleStateMeas(state::Array{Float64,2}, meas::Array{Float64,2})

    Nx, Nes = size(state)
    Ny, Nem = size(meas)
    @assert Nes==Nem "Error dimension of the ensemble"
    Ne = Nes

    return EnsembleStateMeas{Nx, Ny, Ne}(EnsembleState{Nx, Ne}(state), EnsembleState{Ny, Ne}(meas))
end


text:
__precompile__(true)

module ReactiveExtensions
    import Base: get, reduce

    export

    # subject
    Functor,
    Subject,
    get,
    init_task,
    _map,
    _filter,
    _reduce,

    # chaining
    map_,
    filter_,
    reduce,
    bridge,

    # sequence
    Sequence,
    from,
    to_list

    include("subject.jl")
    include("chaining.jl")
    include("sequence.jl")

end # module

ext: jl


text:
module Interpolations

export
    interpolate,
    interpolate!,
    extrapolate,
    scale,

    gradient!,
    gradient1,
    hessian!,
    hessian,
    hessian1,

    AbstractInterpolation,
    AbstractExtrapolation,

    OnCell,
    OnGrid,

    Flat,
    Line,
    Free,
    Periodic,
    Reflect,
    Natural,
    InPlace,
    InPlaceQ

    # see the following files for further exports:
    # b-splines/b-splines.jl
    # extrapolation/extrapolation.jl
    # scaling/scaling.jl

using Compat
using WoodburyMatrices, Ratios, AxisAlgorithms

import Base: convert, size, getindex, gradient, scale, promote_rule, ndims, eltype

abstract Flag
abstract InterpolationType <: Flag
immutable NoInterp <: InterpolationType end
abstract GridType <: Flag
immutable OnGrid <: GridType end
immutable OnCell <: GridType end

typealias DimSpec{T} Union{T,Tuple{Vararg{Union{T,NoInterp}}},NoInterp}

abstract AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} <: AbstractArray{T,N}
abstract AbstractInterpolationWrapper{T,N,ITPT,IT,GT} <: AbstractInterpolation{T,N,IT,GT}
abstract AbstractExtrapolation{T,N,ITPT,IT,GT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT}

immutable Throw <: Flag end
immutable Flat <: Flag end
immutable Line <: Flag end
immutable Free <: Flag end
immutable Periodic <: Flag end
immutable Reflect <: Flag end
immutable InPlace <: Flag end
# InPlaceQ is exact for an underlying quadratic. This is nice for ground-truth testing
# of in-place (unpadded) interpolation.
immutable InPlaceQ <: Flag end
typealias Natural Line

# TODO: size might have to be faster?
size{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->size(itp,i), N)::NTuple{N,Int}
size(exp::AbstractExtrapolation, d) = size(exp.itp, d)
bounds{T,N}(itp::AbstractInterpolation{T,N}) = tuple(zip(lbounds(itp), ubounds(itp))...)
bounds{T,N}(itp::AbstractInterpolation{T,N}, d) = (lbound(itp,d),ubound(itp,d))
lbounds{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->lbound(itp,i), N)::NTuple{N,T}
ubounds{T,N}(itp::AbstractInterpolation{T,N}) = ntuple(i->ubound(itp,i), N)::NTuple{N,T}
lbound{T,N}(itp::AbstractInterpolation{T,N}, d) = 1
ubound{T,N}(itp::AbstractInterpolation{T,N}, d) = size(itp, d)
itptype{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(::Type{AbstractInterpolation{T,N,IT,GT}}) = IT
itptype{ITP<:AbstractInterpolation}(::Type{ITP}) = itptype(supertype(ITP))
itptype(itp::AbstractInterpolation ) = itptype(typeof(itp))
gridtype{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(::Type{AbstractInterpolation{T,N,IT,GT}}) = GT
gridtype{ITP<:AbstractInterpolation}(::Type{ITP}) = gridtype(supertype(ITP))
gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
ndims{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(::Type{AbstractInterpolation{T,N,IT,GT}}) = N
ndims{ITP<:AbstractInterpolation}(::Type{ITP}) = ndims(supertype(ITP))
ndims(itp::AbstractInterpolation) = ndims(typeof(itp))
eltype{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(::Type{AbstractInterpolation{T,N,IT,GT}}) = T
eltype{ITP<:AbstractInterpolation}(::Type{ITP}) = eltype(supertype(ITP))
eltype(itp::AbstractInterpolation) = eltype(typeof(itp))
count_interp_dims{T<:AbstractInterpolation}(::Type{T}, N) = N

include("nointerp/nointerp.jl")
include("b-splines/b-splines.jl")
include("gridded/gridded.jl")
include("extrapolation/extrapolation.jl")
include("scaling/scaling.jl")
include("utils.jl")

end # module


text:
@recipe function f(curve::AbstractManifold{1, C}; npoints=1000) where C
    pts = [curve(t) for t in range(0.0, 0.1, length=npoints)]

    linewidth --> 2
    color     --> :orange
    legend    --> false

    xs = get.(pts, 1, 0.0)
    ys = get.(pts, 2, 0.0)
    zs = get.(pts, 3, 0.0)

    xs, ys, zs
end

@recipe function f(surface::AbstractManifold{2, C}; npoints = 100) where C
    rng = range(0.0, 1.0, length=npoints)

    xs = zeros(npoints, npoints)
    ys = zeros(npoints, npoints)
    zs = zeros(npoints, npoints)
    ws = zeros(npoints, npoints)
    for (i, θ) in enumerate(rng), (j, φ) in enumerate(rng)
        v = surface(θ, φ)
        xs[i, j] = get(v, 1, 0.0)
        ys[i, j] = get(v, 2, 0.0)
        zs[i, j] = get(v, 3, 0.0)
        ws[i, j] = get(v, 4, 0.0)
    end
    seriestype := :surface
    if !iszero(ws)
        surfacecolor := ws
    end

    xs, ys, zs
end

@userplot Points3d
@recipe function f(p::Points3d)
    if length(p.args) != 1 || !(typeof(first(p.args)) <: AbstractVector{<:AbstractVector})
        error("points3d is expecting a vector of vectors. " *
              "Got: $(typeof(p.args))")
    end
    pts = first(p.args)

    seriestype := :scatter
    markersize --> 0.5
    get.(pts, 1, 0.0), get.(pts, 2, 0.0), get.(pts, 3, 0.0)
end


text:
export Grad, NGrad, Hessian, gradient

"""
    NGrad{N,FT} <: Function

Obtain gradients `Grad(f)(Val(i), args..., kwargs...)`, where `i` is the index of loss in `args`. `Grad` object calls forward first, and then backward.

!!! note
    `Val(1)` is specially optimized, so putting the loss as the first parameter can avoid potential overhead.
```
"""
struct NGrad{N,FT} <: Function
    f::FT
end
function NGrad{N}(f::FT) where {N,FT}
    NGrad{N,FT}(f)
end

const Grad{FT} = NGrad{1,FT}
const Hessian{FT} = NGrad{2,FT}

Base.show_function(io::IO, b::NGrad{N}, compact::Bool) where {N} = print(io, "$(b.f)"*"'"^N)
Base.show_function(io::IO, ::MIME"text/plain", b::NGrad{N}, compact::Bool) where {N} = print(io, b)
Base.display(bf::NGrad) = print(bf)
(_::Type{Inv{NGrad{N}}})(f::NGrad{M}) where {M, N} = NGrad{M-N}(f.f)
(_::Type{Inv{NGrad{M}}})(f::NGrad{M}) where {M} = f.f

@i function (g::Grad)(il::Val{iloss}, args...; kwargs...) where iloss
    protectf(g).f(args...; kwargs...)
    GVar.(args)
    INC(args |> tget(iloss) |> grad)
    (~protectf(g).f)(args...; kwargs...)
end

@i function (g::Grad)(il::Val{1}, x, ys...; kwargs...)
    protectf(g).f(x, ys...; kwargs...)
    GVar(x)
    INC(x |> grad)
    GVar.(ys)
    (~protectf(g).f)(x, ys...; kwargs...)
end

@i function (g::Grad)(args...; iloss::Int, kwargs...)
    protectf(g).f(args...; kwargs...)
    GVar.(args)
    INC(args |> tget(iloss) |> grad)
    (~protectf(g).f)(args...; kwargs...)
end

@generated function gradient(::Val{iloss}, f, args::NTuple{N,Any}; kwargs...) where {iloss,N}
    newres = gensym()
    newargs = Any[:(GVar($newres[$i])) for i=1:N]
    newargs[iloss] = :(GVar($newres[$iloss], one($newres[$iloss])))
    quote
        $newres = f(args...; kwargs...)
        grad.(NiLangCore.wrap_tuple((~f)($(newargs...); kwargs...)))
    end
end

gradient(f, args; iloss::Int, kwargs...) = gradient(Val(iloss), f, args; kwargs...)


text:
# Modelo de tanque de armazenamento com dois alimentacoes
type Tanque2
	Tanque2()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in1}",
				:PosX=>0.1,
				:PosY=>0.20
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in2}",
				:PosX=>0.1,
				:PosY=>0.79
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>0.9,
				:PosY=>0.5
			)),
			[
				:(Fout.F = Fin1.F + Fin2.F),
				:(Fout.F * Fout.H = Fin1.F * Fin1.H + Fin2.F * Fin2.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
			],
			[
				"","","",
			],
			[:PP2,],
			[:Fin1,:Fin2,:Fout,]
		)
	end
	PP2::DanaPlugin
	Fin1::Corrente
	Fin2::Corrente
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Tanque2
function setEquationFlow(in::Tanque2)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function atributes(in::Tanque2,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/tanque2"
	drive!(fields,_)
	return fields
end
Tanque2(_::Dict{Symbol,Any})=begin
	newModel=Tanque2()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# trace data character codes
# /* s = 12 bit unsigned stored as short int, 0 to 4095, */
# /* q = 12 bit signed stored as short int, -2048 to 2047, */
# /* u = 16 bit unsigned stored as short int, 0 to 65535 */
# /* i = 16 bit signed stored as short int, -32767 to 32767, */
# /* 2 = 24 bit signed integer stored as long, */
# /* l = 32 bit signed integer stored as long, */
# /*  r = 12 bit data, 4 lsb time stored as short int, */
# /* f = float (32 bit IEEE real), */
# /* d = double (64 bit IEEE real), */
# /*  c = complex, */
# /* v = vector, */
# /* t = tensor */

function suds_decode(x::Array{UInt8,1}, code::UInt8)
  if code in (0x69, 0x71, 0x73, 0x75)   # 'i', 'q', 's', 'u'
    y = reinterpret(Int16, x)
    s = 2
  elseif code == 0x32 || code == 0x6c   # '2', 'l'
    y = reinterpret(Int32, x)
    s = 4
  elseif code == 0x63                   # 'c'
    y = reinterpret(Complex{Float32}, x)
    s = 8
  elseif code == 0x64                   # 'd'
    y = reinterpret(Float64, x)
    s = 8
  elseif code == 0x66                   # 'f'
    y = reinterpret(Float32, x)
    s = 4
  else
    error(string("no decoder for trace data code ", code, "!"))
  end

  return y,s
end
# No decoders for (and no idea what to do with)
# 0x72 # 'r',       would need a win32-style bits parser
# 0x74 # 't',       not a bits type, not defined in documentation
# 0x76 # 'v',       not a bits type, not defined in documentation


text:
"""
CSRTOSPARSE
-----------
    convert a matrix from compressed sparse row to a sparse matrix A.
    It returns the arrays that feed the sparse function in julia.

Functions
---------
- (nzi,nzj,nzv) = csr_to_sparse{T}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{T})
- (nzi,nzj,nzv) = csr_to_sparse{T}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{T},nrows::Int64)

Example
-------
~~~
i = [1;2;3]
j = [3;4;4]
v = [8;9;10]
(rp,ci,ai,m) = sparse_to_csr(i,j,v)
(nzi,nzj,nzv) = csr_to_sparse(rp,ci,ai)
A = sparse(nzi,nzj,nzv,length(rp)-1,maximum(ci))
B = csr_to_sparse_matrix(rp,ci,ai)
isequal(A,B)
~~~
"""

##################
#	Functions    #
##################
function csr_to_sparse{T}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{T})
    nrows = length(rp)-1
    ncols = length(ci)
    nzi = zeros(Int64,ncols)
    for i=1:nrows
        for j=rp[i]:rp[i+1]-1
            nzi[j] = i
        end
    end
    nzj = ci
    nzv = ai
    return (nzi,nzj,nzv)
end

function csr_to_sparse_matrix{T}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{T})
    (i,j,k) = csr_to_sparse(rp,ci,ai)
    A = sparse(i,j,k)
    return A
end

function csr_to_sparse_matrix{T}(rp::Vector{Int64},ci::Vector{Int64},
                                            ai::Vector{T},nrows::Int64,ncols::Int64)
    (i,j,k) = csr_to_sparse(rp,ci,ai);
    A = sparse(i,j,k,nrows,ncols)
    return A
end


text:
# Example of running this routine:

include("run_hmc_background_student_rescale.jl")

# File containing the transit times:
fname = "timing_simon.txt"

felements = "../../data/elements_noprior_students.txt"

# Name an output file:
foutput = string("T1_run_hmc_student_ecc_lgndof_V1exp2nuinv_nstep10_eps0.1_nleap20_000.jld2")

# Call the likelihood profile function:
#run_hmc(fname,felements,foutput,2000,0.1,20,false)
run_hmc(fname,felements,foutput,10,0.1,20,false)


text:
module Artifacts

import Base: get, SHA1
using Base.BinaryPlatforms
using Artifacts
import Artifacts: artifact_names, ARTIFACTS_DIR_OVERRIDE, ARTIFACT_OVERRIDES, artifact_paths,
                  artifacts_dirs, pack_platform!, unpack_platform, load_artifacts_toml,
                  query_override, with_artifacts_directory, load_overrides
import ..get_bool_env
import ..set_readonly
import ..GitTools
import ..TOML
using ..MiniProgressBars
using ..PlatformEngines
import ..pkg_server, ..can_fancyprint, ..stderr_f, ..printpkgstyle
import ..Types: write_env_usage, parse_toml

using SHA

export create_artifact, artifact_exists, artifact_path, remove_artifact, verify_artifact,
       artifact_meta, artifact_hash, bind_artifact!, unbind_artifact!, download_artifact,
       find_artifacts_toml, ensure_artifact_installed, @artifact_str, archive_artifact,
       select_downloadable_artifacts

"""
    create_artifact(f::Function)

Creates a new artifact by running `f(artifact_path)`, hashing the result, and moving it
to the artifact store (`~/.julia/artifacts` on a typical installation).  Returns the
identifying tree hash of this artifact.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function create_artifact(f::Function)
    # Ensure the `artifacts` directory exists in our default depot
    artifacts_dir = first(artifacts_dirs())
    mkpath(artifacts_dir)

    # Temporary directory where we'll do our creation business
    temp_dir = mktempdir(artifacts_dir)

    try
        # allow the user to do their work inside the temporary directory
        f(temp_dir)

        # Calculate the tree hash for this temporary directory
        artifact_hash = SHA1(GitTools.tree_hash(temp_dir))

        # If we created a dupe, just let the temp directory get destroyed. It's got the
        # same contents as whatever already exists after all, so it doesn't matter.  Only
        # move its contents if it actually contains new contents.  Note that we explicitly
        # set `honor_overrides=false` here, as we wouldn't want to drop things into the
        # system directory by accidentally creating something with the same content-hash
        # as something that was foolishly overridden.  This should be virtually impossible
        # unless the user has been very unwise, but let's be cautious.
        new_path = artifact_path(artifact_hash; honor_overrides=false)
        if !isdir(new_path)
            # Move this generated directory to its final destination, set it to read-only
            mv(temp_dir, new_path)
            chmod(new_path, filemode(dirname(new_path)))
            set_readonly(new_path)
        end

        # Give the people what they want
        return artifact_hash
    finally
        # Always attempt to cleanup
        rm(temp_dir; recursive=true, force=true)
    end
end

"""
    remove_artifact(hash::SHA1; honor_overrides::Bool=false)

Removes the given artifact (identified by its SHA1 git tree hash) from disk.  Note that
if an artifact is installed in multiple depots, it will be removed from all of them.  If
an overridden artifact is requested for removal, it will be silently ignored; this method
will never attempt to remove an overridden artifact.

In general, we recommend that you use `Pkg.gc()` to manage artifact installations and do
not use `remove_artifact()` directly, as it can be difficult to know if an artifact is
being used by another package.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function remove_artifact(hash::SHA1)
    if query_override(hash) !== nothing
        # We never remove overridden artifacts.
        return
    end

    # Get all possible paths (rooted in all depots)
    possible_paths = artifacts_dirs(bytes2hex(hash.bytes))
    for path in possible_paths
        if isdir(path)
            rm(path; recursive=true, force=true)
        end
    end
end

"""
    verify_artifact(hash::SHA1; honor_overrides::Bool=false)

Verifies that the given artifact (identified by its SHA1 git tree hash) is installed on-
disk, and retains its integrity.  If the given artifact is overridden, skips the
verification unless `honor_overrides` is set to `true`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function verify_artifact(hash::SHA1; honor_overrides::Bool=false)
    # Silently skip overridden artifacts unless we really ask for it
    if !honor_overrides
        if query_override(hash) !== nothing
            return true
        end
    end

    # If it doesn't even exist, then skip out
    if !artifact_exists(hash)
        return false
    end

    # Otherwise actually run the verification
    return all(hash.bytes .== GitTools.tree_hash(artifact_path(hash)))
end

"""
    archive_artifact(hash::SHA1, tarball_path::String; honor_overrides::Bool=false)

Archive an artifact into a tarball stored at `tarball_path`, returns the SHA256 of the
resultant tarball as a hexidecimal string. Throws an error if the artifact does not
exist.  If the artifact is overridden, throws an error unless `honor_overrides` is set.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function archive_artifact(hash::SHA1, tarball_path::String; honor_overrides::Bool=false)
    if !honor_overrides
        if query_override(hash) !== nothing
            error("Will not archive an overridden artifact unless `honor_overrides` is set!")
        end
    end

    if !artifact_exists(hash)
        error("Unable to archive artifact $(bytes2hex(hash.bytes)): does not exist!")
    end

    # Package it up
    package(artifact_path(hash), tarball_path)

    # Calculate its sha256 and return that
    return open(tarball_path, "r") do io
        return bytes2hex(sha256(io))
    end
end

"""
    bind_artifact!(artifacts_toml::String, name::String, hash::SHA1;
                   platform::Union{AbstractPlatform,Nothing} = nothing,
                   download_info::Union{Vector{Tuple},Nothing} = nothing,
                   lazy::Bool = false,
                   force::Bool = false)

Writes a mapping of `name` -> `hash` within the given `(Julia)Artifacts.toml` file. If
`platform` is not `nothing`, this artifact is marked as platform-specific, and will be
a multi-mapping.  It is valid to bind multiple artifacts with the same name, but
different `platform`s and `hash`'es within the same `artifacts_toml`.  If `force` is set
to `true`, this will overwrite a pre-existant mapping, otherwise an error is raised.

`download_info` is an optional vector that contains tuples of URLs and a hash.  These
URLs will be listed as possible locations where this artifact can be obtained.  If `lazy`
is set to `true`, even if download information is available, this artifact will not be
downloaded until it is accessed via the `artifact"name"` syntax, or
`ensure_artifact_installed()` is called upon it.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function bind_artifact!(artifacts_toml::String, name::String, hash::SHA1;
                        platform::Union{AbstractPlatform,Nothing} = nothing,
                        download_info::Union{Vector{<:Tuple},Nothing} = nothing,
                        lazy::Bool = false,
                        force::Bool = false)
    # First, check to see if this artifact is already bound:
    if isfile(artifacts_toml)
        artifact_dict = parse_toml(artifacts_toml)

        if !force && haskey(artifact_dict, name)
            meta = artifact_dict[name]
            if !isa(meta, Vector)
                error("Mapping for '$name' within $(artifacts_toml) already exists!")
            elseif any(isequal(platform), unpack_platform(x, name, artifacts_toml) for x in meta)
                error("Mapping for '$name'/$(triplet(platform)) within $(artifacts_toml) already exists!")
            end
        end
    else
        artifact_dict = Dict{String, Any}()
    end

    # Otherwise, the new piece of data we're going to write out is this dict:
    meta = Dict{String,Any}(
        "git-tree-sha1" => bytes2hex(hash.bytes),
    )

    # If we're set to be lazy, then lazy we shall be
    if lazy
        meta["lazy"] = true
    end

    # Integrate download info, if it is given.  We represent the download info as a
    # vector of dicts, each with its own `url` and `sha256`, since different tarballs can
    # expand to the same tree hash.
    if download_info !== nothing
        meta["download"] = [
            Dict("url" => dl[1],
                 "sha256" => dl[2],
            ) for dl in download_info
        ]
    end

    if platform === nothing
        artifact_dict[name] = meta
    else
        # Add platform-specific keys to our `meta` dict
        pack_platform!(meta, platform)

        # Insert this entry into the list of artifacts
        if !haskey(artifact_dict, name)
            artifact_dict[name] = [meta]
        else
            # Delete any entries that contain identical platforms
            artifact_dict[name] = filter(
                x -> unpack_platform(x, name, artifacts_toml) != platform,
                artifact_dict[name]
            )
            push!(artifact_dict[name], meta)
        end
    end

    # Spit it out onto disk
    let artifact_dict = artifact_dict
        open(artifacts_toml, "w") do io
            TOML.print(io, artifact_dict, sorted=true)
        end
    end

    # Mark that we have used this Artifact.toml
    write_env_usage(artifacts_toml, "artifact_usage.toml")
    return
end

"""
    unbind_artifact!(artifacts_toml::String, name::String; platform = nothing)

Unbind the given `name` from an `(Julia)Artifacts.toml` file.
Silently fails if no such binding exists within the file.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function unbind_artifact!(artifacts_toml::String, name::String;
                         platform::Union{AbstractPlatform,Nothing} = nothing)
    artifact_dict = parse_toml(artifacts_toml)
    if !haskey(artifact_dict, name)
        return
    end

    if platform === nothing
        delete!(artifact_dict, name)
    else
        artifact_dict[name] = filter(
            x -> unpack_platform(x, name, artifacts_toml) != platform,
            artifact_dict[name]
        )
    end

    open(artifacts_toml, "w") do io
        TOML.print(io, artifact_dict, sorted=true)
    end
    return
end

"""
    download_artifact(tree_hash::SHA1, tarball_url::String, tarball_hash::String;
                      verbose::Bool = false, io::IO=stderr)

Download/install an artifact into the artifact store.  Returns `true` on success.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function download_artifact(
    tree_hash::SHA1,
    tarball_url::String,
    tarball_hash::Union{String, Nothing} = nothing;
    verbose::Bool = false,
    quiet_download::Bool = false,
    io::IO=stderr_f(),
)
    if artifact_exists(tree_hash)
        return true
    end

    if Sys.iswindows()
        # The destination directory we're hoping to fill:
        dest_dir = artifact_path(tree_hash; honor_overrides=false)
        mkpath(dest_dir)

        # On Windows, we have some issues around stat() and chmod() that make properly
        # determining the git tree hash problematic; for this reason, we use the "unsafe"
        # artifact unpacking method, which does not properly verify unpacked git tree
        # hash.  This will be fixed in a future Julia release which will properly interrogate
        # the filesystem ACLs for executable permissions, which git tree hashes care about.
        try
            download_verify_unpack(tarball_url, tarball_hash, dest_dir, ignore_existence=true,
                                   verbose=verbose, quiet_download=quiet_download, io=io)
        catch err
            @debug "download_artifact error" tree_hash tarball_url tarball_hash err
            # Clean that destination directory out if something went wrong
            rm(dest_dir; force=true, recursive=true)

            if isa(err, InterruptException)
                rethrow(err)
            end
            return false
        end
    else
        # We download by using `create_artifact()`.  We do this because the download may
        # be corrupted or even malicious; we don't want to clobber someone else's artifact
        # by trusting the tree hash that has been given to us; we will instead download it
        # to a temporary directory, calculate the true tree hash, then move it to the proper
        # location only after knowing what it is, and if something goes wrong in the process,
        # everything should be cleaned up.  Luckily, that is precisely what our
        # `create_artifact()` wrapper does, so we use that here.
        calc_hash = try
            create_artifact() do dir
                download_verify_unpack(tarball_url, tarball_hash, dir, ignore_existence=true, verbose=verbose,
                    quiet_download=quiet_download, io=io)
            end
        catch err
            @debug "download_artifact error" tree_hash tarball_url tarball_hash err
            if isa(err, InterruptException)
                rethrow(err)
            end
            # If something went wrong during download, return false
            return false
        end

        # Did we get what we expected?  If not, freak out.
        if calc_hash.bytes != tree_hash.bytes
            msg  = "Tree Hash Mismatch!\n"
            msg *= "  Expected git-tree-sha1:   $(bytes2hex(tree_hash.bytes))\n"
            msg *= "  Calculated git-tree-sha1: $(bytes2hex(calc_hash.bytes))"
            # Since tree hash calculation is still broken on some systems, e.g. Pkg.jl#1860,
            # and Pkg.jl#2317 so we allow setting JULIA_PKG_IGNORE_HASHES=1 to ignore the
            # error and move the artifact to the expected location and return true
            ignore_hash = get_bool_env("JULIA_PKG_IGNORE_HASHES")
            if ignore_hash
                msg *= "\n\$JULIA_PKG_IGNORE_HASHES is set to 1: ignoring error and moving artifact to the expected location"
            end
            @error(msg)
            if ignore_hash
                # Move it to the location we expected
                src = artifact_path(calc_hash; honor_overrides=false)
                dst = artifact_path(tree_hash; honor_overrides=false)
                mv(src, dst; force=true)
                return true
            end
            return false
        end
    end

    return true
end

"""
    ensure_artifact_installed(name::String, artifacts_toml::String;
                              platform::AbstractPlatform = HostPlatform(),
                              pkg_uuid::Union{Base.UUID,Nothing}=nothing,
                              verbose::Bool = false,
                              quiet_download::Bool = false,
                              io::IO=stderr)

Ensures an artifact is installed, downloading it via the download information stored in
`artifacts_toml` if necessary.  Throws an error if unable to install.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function ensure_artifact_installed(name::String, artifacts_toml::String;
                                   platform::AbstractPlatform = HostPlatform(),
                                   pkg_uuid::Union{Base.UUID,Nothing}=nothing,
                                   verbose::Bool = false,
                                   quiet_download::Bool = false,
                                   io::IO=stderr_f())
    meta = artifact_meta(name, artifacts_toml; pkg_uuid=pkg_uuid, platform=platform)
    if meta === nothing
        error("Cannot locate artifact '$(name)' in '$(artifacts_toml)'")
    end

    return ensure_artifact_installed(name, meta, artifacts_toml; platform=platform,
                                     verbose=verbose, quiet_download=quiet_download, io=io)
end

function ensure_artifact_installed(name::String, meta::Dict, artifacts_toml::String;
                                   platform::AbstractPlatform = HostPlatform(),
                                   verbose::Bool = false,
                                   quiet_download::Bool = false,
                                   io::IO=stderr_f())
    hash = SHA1(meta["git-tree-sha1"])

    if !artifact_exists(hash)
        # first try downloading from Pkg server
        # TODO: only do this if Pkg server knows about this package
        if (server = pkg_server()) !== nothing
            url = "$server/artifact/$hash"
            download_success = with_show_download_info(io, name, quiet_download) do
                download_artifact(hash, url; verbose=verbose, quiet_download=quiet_download, io=io)
            end
            download_success && return artifact_path(hash)
        end

        # If this artifact does not exist on-disk already, ensure it has download
        # information, then download it!
        if !haskey(meta, "download")
            error("Cannot automatically install '$(name)'; no download section in '$(artifacts_toml)'")
        end

        # Attempt to download from all sources
        for entry in meta["download"]
            url = entry["url"]
            tarball_hash = entry["sha256"]
            download_success = with_show_download_info(io, name, quiet_download) do
                download_artifact(hash, url, tarball_hash; verbose=verbose, quiet_download=quiet_download, io=io)
            end
            download_success && return artifact_path(hash)
        end
        error("Unable to automatically install '$(name)' from '$(artifacts_toml)'")
    else
        return artifact_path(hash)
    end
end

function with_show_download_info(f, io, name, quiet_download)
    fancyprint = can_fancyprint(io)
    if !quiet_download
        fancyprint && print_progress_bottom(io)
        printpkgstyle(io, :Downloading, "artifact: $name")
    end
    try
        return f()
    finally
        if !quiet_download
            fancyprint && print(io, "\033[1A") # move cursor up one line
            fancyprint && print(io, "\033[2K") # clear line
            fancyprint && printpkgstyle(io, :Downloaded, "artifact: $name")
        end
    end
end

"""
    ensure_all_artifacts_installed(artifacts_toml::String;
                                   platform = HostPlatform(),
                                   pkg_uuid = nothing,
                                   include_lazy = false,
                                   verbose = false,
                                   quiet_download = false,
                                   io::IO=stderr)

Installs all non-lazy artifacts from a given `(Julia)Artifacts.toml` file. `package_uuid` must
be provided to properly support overrides from `Overrides.toml` entries in depots.

If `include_lazy` is set to `true`, then lazy packages will be installed as well.

This function is deprecated and should be replaced with the following snippet:

    artifacts = select_downloadable_artifacts(artifacts_toml; platform, include_lazy)
    for name in keys(artifacts)
        ensure_artifact_installed(name, artifacts[name], artifacts_toml; platform=platform)
    end

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.

!!! warning
    This function is deprecated in Julia 1.6 and will be removed in a future version.
    Use `select_downloadable_artifacts()` and `ensure_artifact_installed()` instead.
"""
function ensure_all_artifacts_installed(artifacts_toml::String;
                                        platform::AbstractPlatform = HostPlatform(),
                                        pkg_uuid::Union{Nothing,Base.UUID} = nothing,
                                        include_lazy::Bool = false,
                                        verbose::Bool = false,
                                        quiet_download::Bool = false,
                                        io::IO=stderr_f())
    # This function should not be called anymore; use `select_downloadable_artifacts()` directly.
    Base.depwarn("`ensure_all_artifacts_installed()` is deprecated; iterate over `select_downloadable_artifacts()` output with `ensure_artifact_installed()`.", :ensure_all_artifacts_installed)
    # Collect all artifacts we're supposed to install
    artifacts = select_downloadable_artifacts(artifacts_toml; platform, include_lazy, pkg_uuid)
    for name in keys(artifacts)
        # Otherwise, let's try and install it!
        ensure_artifact_installed(name, artifacts[name], artifacts_toml; platform=platform,
                                  verbose=verbose, quiet_download=quiet_download, io=io)
    end
end

"""
    extract_all_hashes(artifacts_toml::String;
                       platform = HostPlatform(),
                       pkg_uuid = nothing,
                       include_lazy = false)

Extract all hashes from a given `(Julia)Artifacts.toml` file. `package_uuid` must
be provided to properly support overrides from `Overrides.toml` entries in depots.

If `include_lazy` is set to `true`, then lazy packages will be installed as well.
"""
function extract_all_hashes(artifacts_toml::String;
                            platform::AbstractPlatform = HostPlatform(),
                            pkg_uuid::Union{Nothing,Base.UUID} = nothing,
                            include_lazy::Bool = false)
    hashes = Base.SHA1[]
    if !isfile(artifacts_toml)
        return hashes
    end

    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)

    for name in keys(artifact_dict)
        # Get the metadata about this name for the requested platform
        meta = artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)

        # If there are no instances of this name for the desired platform, skip it
        meta === nothing && continue

        # If it's a lazy one and we aren't including lazy ones, skip
        if get(meta, "lazy", false) && !include_lazy
            continue
        end

        # Otherwise, add it to the list!
        push!(hashes, Base.SHA1(meta["git-tree-sha1"]))
    end

    return hashes
end

# Support `AbstractString`s, but avoid compilers needing to track backedges for callers
# of these functions in case a user defines a new type that is `<: AbstractString`
archive_artifact(hash::SHA1, tarball_path::AbstractString; kwargs...) =
    archive_artifact(hash, string(tarball_path)::String; kwargs...)
bind_artifact!(artifacts_toml::AbstractString, name::AbstractString, hash::SHA1; kwargs...) =
    bind_artifact!(string(artifacts_toml)::String, string(name)::String, hash; kwargs...)
unbind_artifact!(artifacts_toml::AbstractString, name::AbstractString) =
    unbind_artifact!(string(artifacts_toml)::String, string(name)::String)
download_artifact(tree_hash::SHA1, tarball_url::AbstractString, args...; kwargs...) =
    download_artifact(tree_hash, string(tarball_url)::String, args...; kwargs...)
ensure_artifact_installed(name::AbstractString, artifacts_toml::AbstractString; kwargs...) =
    ensure_artifact_installed(string(name)::String, string(artifacts_toml)::String; kwargs...)
ensure_artifact_installed(name::AbstractString, meta::Dict, artifacts_toml::AbstractString; kwargs...) =
    ensure_artifact_installed(string(name)::String, meta, string(artifacts_toml)::String; kwargs...)
ensure_all_artifacts_installed(artifacts_toml::AbstractString; kwargs...) =
    ensure_all_artifacts_installed(string(name)::String; kwargs...)
extract_all_hashes(artifacts_toml::AbstractString; kwargs...) =
    extract_all_hashes(string(artifacts_toml)::String; kwargs...)

end # module Artifacts


text:
using FNVHash
using Base.Test

@test fnv32([0x01,0x47,0x6c,0x10,0xf3]) == 0
@test fnv32([0xfd, 0x45, 0x41, 0x08, 0xa0]) == 0

# Taken from http://www.isthe.com/chongo/src/fnv/test_fnv.c ...

TEST(s) = s

function TEST0(s)
    a = Vector{UInt8}(s)
    push!(a, 0x00)
end

R10(s::AbstractString) = s^10
R500(s::AbstractString) = s^500

R10(s) = repeat(s, outer=[10])
R500(s) = repeat(s, outer=[500])

input = Any[
  TEST(""),
  TEST("a"),
  TEST("b"),
  TEST("c"),
  TEST("d"),
  TEST("e"),
  TEST("f"),
  TEST("fo"),
  TEST("foo"),
  TEST("foob"),
  TEST("fooba"),
  TEST("foobar"),
  TEST0(""),
  TEST0("a"),
  TEST0("b"),
  TEST0("c"),
  TEST0("d"),
  TEST0("e"),
  TEST0("f"),
  TEST0("fo"),
  TEST0("foo"),
  TEST0("foob"),
  TEST0("fooba"),
  TEST0("foobar"),
  TEST("ch"),
  TEST("cho"),
  TEST("chon"),
  TEST("chong"),
  TEST("chongo"),
  TEST("chongo "),
  TEST("chongo w"),
  TEST("chongo wa"),
  TEST("chongo was"),
  TEST("chongo was "),
  TEST("chongo was h"),
  TEST("chongo was he"),
  TEST("chongo was her"),
  TEST("chongo was here"),
  TEST("chongo was here!"),
  TEST("chongo was here!\n"),
  TEST0("ch"),
  TEST0("cho"),
  TEST0("chon"),
  TEST0("chong"),
  TEST0("chongo"),
  TEST0("chongo "),
  TEST0("chongo w"),
  TEST0("chongo wa"),
  TEST0("chongo was"),
  TEST0("chongo was "),
  TEST0("chongo was h"),
  TEST0("chongo was he"),
  TEST0("chongo was her"),
  TEST0("chongo was here"),
  TEST0("chongo was here!"),
  TEST0("chongo was here!\n"),
  TEST("cu"),
  TEST("cur"),
  TEST("curd"),
  TEST("curds"),
  TEST("curds "),
  TEST("curds a"),
  TEST("curds an"),
  TEST("curds and"),
  TEST("curds and "),
  TEST("curds and w"),
  TEST("curds and wh"),
  TEST("curds and whe"),
  TEST("curds and whey"),
  TEST("curds and whey\n"),
  TEST0("cu"),
  TEST0("cur"),
  TEST0("curd"),
  TEST0("curds"),
  TEST0("curds "),
  TEST0("curds a"),
  TEST0("curds an"),
  TEST0("curds and"),
  TEST0("curds and "),
  TEST0("curds and w"),
  TEST0("curds and wh"),
  TEST0("curds and whe"),
  TEST0("curds and whey"),
  TEST0("curds and whey\n"),
  TEST("hi"), TEST0("hi"),
  TEST("hello"), TEST0("hello"),
  TEST(b"\xff\x00\x00\x01"), TEST(b"\x01\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x02"), TEST(b"\x02\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x03"), TEST(b"\x03\x00\x00\xff"),
  TEST(b"\xff\x00\x00\x04"), TEST(b"\x04\x00\x00\xff"),
  TEST(b"\x40\x51\x4e\x44"), TEST(b"\x44\x4e\x51\x40"),
  TEST(b"\x40\x51\x4e\x4a"), TEST(b"\x4a\x4e\x51\x40"),
  TEST(b"\x40\x51\x4e\x54"), TEST(b"\x54\x4e\x51\x40"),
  TEST("127.0.0.1"), TEST0("127.0.0.1"),
  TEST("127.0.0.2"), TEST0("127.0.0.2"),
  TEST("127.0.0.3"), TEST0("127.0.0.3"),
  TEST("64.81.78.68"), TEST0("64.81.78.68"),
  TEST("64.81.78.74"), TEST0("64.81.78.74"),
  TEST("64.81.78.84"), TEST0("64.81.78.84"),
  TEST("feedface"), TEST0("feedface"),
  TEST("feedfacedaffdeed"), TEST0("feedfacedaffdeed"),
  TEST("feedfacedeadbeef"), TEST0("feedfacedeadbeef"),
  TEST("line 1\nline 2\nline 3"),
  TEST("chongo <Landon Curt Noll> /\\../\\"),
  TEST0("chongo <Landon Curt Noll> /\\../\\"),
  TEST("chongo (Landon Curt Noll) /\\../\\"),
  TEST0("chongo (Landon Curt Noll) /\\../\\"),
  TEST("http://antwrp.gsfc.nasa.gov/apod/astropix.html"),
  TEST("http://en.wikipedia.org/wiki/Fowler_Noll_Vo_hash"),
  TEST("http://epod.usra.edu/"),
  TEST("http://exoplanet.eu/"),
  TEST("http://hvo.wr.usgs.gov/cam3/"),
  TEST("http://hvo.wr.usgs.gov/cams/HMcam/"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/deformation.html"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/images.html"),
  TEST("http://hvo.wr.usgs.gov/kilauea/update/maps.html"),
  TEST("http://hvo.wr.usgs.gov/volcanowatch/current_issue.html"),
  TEST("http://neo.jpl.nasa.gov/risk/"),
  TEST("http://norvig.com/21-days.html"),
  TEST("http://primes.utm.edu/curios/home.php"),
  TEST("http://slashdot.org/"),
  TEST("http://tux.wr.usgs.gov/Maps/155.25-19.5.html"),
  TEST("http://volcano.wr.usgs.gov/kilaueastatus.php"),
  TEST("http://www.avo.alaska.edu/activity/Redoubt.php"),
  TEST("http://www.dilbert.com/fast/"),
  TEST("http://www.fourmilab.ch/gravitation/orbits/"),
  TEST("http://www.fpoa.net/"),
  TEST("http://www.ioccc.org/index.html"),
  TEST("http://www.isthe.com/cgi-bin/number.cgi"),
  TEST("http://www.isthe.com/chongo/bio.html"),
  TEST("http://www.isthe.com/chongo/index.html"),
  TEST("http://www.isthe.com/chongo/src/calc/lucas-calc"),
  TEST("http://www.isthe.com/chongo/tech/astro/venus2004.html"),
  TEST("http://www.isthe.com/chongo/tech/astro/vita.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/c/expert.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/calc/index.html"),
  TEST("http://www.isthe.com/chongo/tech/comp/fnv/index.html"),
  TEST("http://www.isthe.com/chongo/tech/math/number/howhigh.html"),
  TEST("http://www.isthe.com/chongo/tech/math/number/number.html"),
  TEST("http://www.isthe.com/chongo/tech/math/prime/mersenne.html"),
  TEST("http://www.isthe.com/chongo/tech/math/prime/mersenne.html#largest"),
  TEST("http://www.lavarnd.org/cgi-bin/corpspeak.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/haiku.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/rand-none.cgi"),
  TEST("http://www.lavarnd.org/cgi-bin/randdist.cgi"),
  TEST("http://www.lavarnd.org/index.html"),
  TEST("http://www.lavarnd.org/what/nist-test.html"),
  TEST("http://www.macosxhints.com/"),
  TEST("http://www.mellis.com/"),
  TEST("http://www.nature.nps.gov/air/webcams/parks/havoso2alert/havoalert.cfm"),
  TEST("http://www.nature.nps.gov/air/webcams/parks/havoso2alert/timelines_24.cfm"),
  TEST("http://www.paulnoll.com/"),
  TEST("http://www.pepysdiary.com/"),
  TEST("http://www.sciencenews.org/index/home/activity/view"),
  TEST("http://www.skyandtelescope.com/"),
  TEST("http://www.sput.nl/~rob/sirius.html"),
  TEST("http://www.systemexperts.com/"),
  TEST("http://www.tq-international.com/phpBB3/index.php"),
  TEST("http://www.travelquesttours.com/index.htm"),
  TEST("http://www.wunderground.com/global/stations/89606.html"),
  TEST(R10("21701")),
  TEST(R10("M21701")),
  TEST(R10("2^21701-1")),
  TEST(R10(b"\x54\xc5")),
  TEST(R10(b"\xc5\x54")),
  TEST(R10("23209")),
  TEST(R10("M23209")),
  TEST(R10("2^23209-1")),
  TEST(R10(b"\x5a\xa9")),
  TEST(R10(b"\xa9\x5a")),
  TEST(R10("391581216093")),
  TEST(R10("391581*2^216093-1")),
  TEST(R10(b"\x05\xf9\x9d\x03\x4c\x81")),
  TEST(R10("FEDCBA9876543210")),
  TEST(R10(b"\xfe\xdc\xba\x98\x76\x54\x32\x10")),
  TEST(R10("EFCDAB8967452301")),
  TEST(R10(b"\xef\xcd\xab\x89\x67\x45\x23\x01")),
  TEST(R10("0123456789ABCDEF")),
  TEST(R10(b"\x01\x23\x45\x67\x89\xab\xcd\xef")),
  TEST(R10("1032547698BADCFE")),
  TEST(R10(b"\x10\x32\x54\x76\x98\xba\xdc\xfe")),
  TEST(R500(b"\x00")),
  TEST(R500(b"\x07")),
  TEST(R500("~")),
  TEST(R500(b"\x7f"))
]

results = [
    0x811c9dc5,
    0x050c5d7e,
    0x050c5d7d,
    0x050c5d7c,
    0x050c5d7b,
    0x050c5d7a,
    0x050c5d79,
    0x6b772514,
    0x408f5e13,
    0xb4b1178b,
    0xfdc80fb0,
    0x31f0b262,
    0x050c5d1f,
    0x70772d5a,
    0x6f772bc7,
    0x6e772a34,
    0x6d7728a1,
    0x6c77270e,
    0x6b77257b,
    0x408f5e7c,
    0xb4b117e9,
    0xfdc80fd1,
    0x31f0b210,
    0xffe8d046,
    0x6e772a5c,
    0x4197aebb,
    0xfcc8100f,
    0xfdf147fa,
    0xbcd44ee1,
    0x23382c13,
    0x846d619e,
    0x1630abdb,
    0xc99e89b2,
    0x1692c316,
    0x9f091bca,
    0x2556be9b,
    0x628e0e73,
    0x98a0bf6c,
    0xb10d5725,
    0xdd002f35,
    0x4197aed4,
    0xfcc81061,
    0xfdf1479d,
    0xbcd44e8e,
    0x23382c33,
    0x846d61e9,
    0x1630abba,
    0xc99e89c1,
    0x1692c336,
    0x9f091ba2,
    0x2556befe,
    0x628e0e01,
    0x98a0bf09,
    0xb10d5704,
    0xdd002f3f,
    0x1c4a506f,
    0x6e772a41,
    0x26978421,
    0xe184ff97,
    0x9b5e5ac6,
    0x5b88e592,
    0xaa8164b7,
    0x20b18c7b,
    0xf28025c5,
    0x84bb753f,
    0x3219925a,
    0x384163c6,
    0x54f010d7,
    0x8cea820c,
    0xe12ab8ee,
    0x26978453,
    0xe184fff3,
    0x9b5e5ab5,
    0x5b88e5b2,
    0xaa8164d6,
    0x20b18c15,
    0xf28025a1,
    0x84bb751f,
    0x3219922d,
    0x384163ae,
    0x54f010b2,
    0x8cea8275,
    0xe12ab8e4,
    0x64411eaa,
    0x6977223c,
    0x428ae474,
    0xb6fa7167,
    0x73408525,
    0xb78320a1,
    0x0caf4135,
    0xb78320a2,
    0xcdc88e80,
    0xb78320a3,
    0x8ee1dbcb,
    0xb78320a4,
    0x4ffb2716,
    0x860632aa,
    0xcc2c5c64,
    0x860632a4,
    0x2a7ec4a6,
    0x860632ba,
    0xfefe8e14,
    0x0a3cffd8,
    0xf606c108,
    0x0a3cffdb,
    0xf906c5c1,
    0x0a3cffda,
    0xf806c42e,
    0xc07167d7,
    0xc9867775,
    0xbf716668,
    0xc78435b8,
    0xc6717155,
    0xb99568cf,
    0x7662e0d6,
    0x33a7f0e2,
    0xc2732f95,
    0xb053e78f,
    0x3a19c02a,
    0xa089821e,
    0x31ae8f83,
    0x995fa9c4,
    0x35983f8c,
    0x5036a251,
    0x97018583,
    0xb4448d60,
    0x025dfe59,
    0xc5eab3af,
    0x7d21ba1e,
    0x7704cddb,
    0xd0071bfe,
    0x0ff3774c,
    0xb0fea0ea,
    0x58177303,
    0x4f599cda,
    0x3e590a47,
    0x965595f8,
    0xc37f178d,
    0x9711dd26,
    0x23c99b7f,
    0x6e568b17,
    0x43f0245b,
    0xbcb7a001,
    0x12e6dffe,
    0x0792f2d6,
    0xb966936b,
    0x46439ac5,
    0x728d49af,
    0xd33745c9,
    0xbc382a57,
    0x4bda1d31,
    0xce35ccae,
    0x3b6eed94,
    0x445c9c58,
    0x3db8bf9d,
    0x2dee116d,
    0xc18738da,
    0x5b156176,
    0x2aa7d593,
    0xb2409658,
    0xe1489528,
    0xfe1ee07e,
    0xe8842315,
    0x3a6a63a2,
    0x06d2c18c,
    0xf8ef7225,
    0x843d3300,
    0xbb24f7ae,
    0x878c0ec9,
    0xb557810f,
    0x57423246,
    0x87f7505e,
    0xbb809f20,
    0x8932abb5,
    0x0a9b3aa0,
    0xb8682a24,
    0xa7ac1c56,
    0x11409252,
    0xa987f517,
    0xf309e7ed,
    0xc9e8f417,
    0x7f447bdd,
    0xb929adc5,
    0x57022879,
    0xdcfd2c49,
    0x6edafff5,
    0xf04fb1f1,
    0xfb7de8b9,
    0xc5f1d7e9,
    0x32c1f439,
    0x7fd3eb7d,
    0x81597da5,
    0x05eb7a25,
    0x9c0fa1b5,
    0x53ccb1c5,
    0xfabece15,
    0x4ad745a5,
    0xe5bdc495,
    0x23b3c0a5,
    0xfa823dd5,
    0x0c6c58b9,
    0xe2dbccd5,
    0xdb7f50f9
]

for (i, r) in zip(input, results)
    @test fnv32(i) == r
end


text:
## tracked.jl - build ExGraph using tracked data types

import Base: +, -, *, /, log, exp, min, max, reshape, transpose, sum,
    abs, abs2, >, >=, <, <=, minimum, maximum, getindex
import Statistics: mean
# import Broadcast: broadcasted

const DEFAULT_GRAPH = Ref(ExGraph())

get_default_graph() = DEFAULT_GRAPH[]
set_default_graph(g::ExGraph) = (DEFAULT_GRAPH[] = g)
reset_default_graph!() = (DEFAULT_GRAPH[] = ExGraph())
swap_default_graph!(g::ExGraph) = (og = DEFAULT_GRAPH[]; DEFAULT_GRAPH[] = g; og)

## TRACKED REAL

struct TrackedReal <: Real
    graph::ExGraph
    var::Symbol
    val::Real
end

TrackedReal(var::Symbol, val::Real) = TrackedReal(get_default_graph(), var, val)
Base.show(io::IO, x::TrackedReal) = print(io, "Tracked($(x.var) = $(x.val))")
tracked_val(g::ExGraph, var::Symbol, val::Real) = TrackedReal(g, var, val)

struct TrackedArray{T,N} <: AbstractArray{T,N}
    graph::ExGraph
    var::Symbol
    val::AbstractArray{T,N}
end

TrackedArray(var::Symbol, val::AbstractArray) = TrackedArray(get_default_graph(), var, val)

Base.show(io::IO, x::TrackedArray{T,N}) where {T,N} =
    print(io, "Tracked{$T,$N}($(x.var) = $(x.val))")
Base.show(io::IO, ::MIME{Symbol("text/plain")}, x::TrackedArray{T,N}) where {T,N} =
    print(io, "Tracked{$T,$N}($(x.var) = $(x.val))")

tracked_val(g::ExGraph, var::Symbol, val::AbstractArray) = TrackedArray(g, var, val)

## conversion and promotion

tracked(g::ExGraph, x::TrackedReal) = x
tracked(g::ExGraph, x::TrackedArray) = x

function tracked(g::ExGraph, x::T) where {T <: Real}
    var = genname()
    push!(g, ExNode{:constant}(var, x; val=x))
    return TrackedReal(g, var, x)
end

function tracked(g::ExGraph, x::T) where {T <: AbstractArray}
    var = genname()
    push!(g, ExNode{:constant}(var, x; val=x))
    return TrackedArray(g, var, x)
end

value(x::TrackedArray) = x.val
value(x::TrackedReal) = x.val
value(x) = x

istracked(x::TrackedReal) = true
istracked(x::TrackedArray) = true
istracked(::Type{T}) where T = (T <: TrackedReal || T <: TrackedArray)
istracked(x) = false

Base.promote_rule(::Type{TrackedReal}, ::Type{<:Real}) = TrackedReal

function Base.convert(::Type{TrackedReal}, x::R) where {R <: Real}
    var = genname()
    g = get_default_graph()
    push!(g, ExNode{:constant}(var, x; val=x))
    return TrackedReal(g, var, x)
end

Base.convert(::Type{TrackedReal}, x::TrackedReal) = x
Base.convert(::Type{R}, x::TrackedReal) where {R <: Real} = x.val

## TRACKED ARRAY

# TODO: should write this to graph? presumably no
Base.size(x::TrackedArray) = size(x.val)
# TODO: and these? presumably yes, but be carefull about printing
# Base.getindex(x::TrackedArray, I...) = getindex(x.val, I...)
# Base.setindex!(x::TrackedArray, val, I...) = setindex!(x.val, val, I...)

## conversion and promotion

Base.promote_rule(::Type{TrackedArray}, ::Type{<:AbstractArray}) = TrackedArray

function Base.convert(::Type{TrackedArray}, x::A) where {A <: AbstractArray}
    var = genname()
    g = get_default_graph()
    push!(g, ExNode{:constant}(var, x; val=x))
    return TrackedArray(g, var, x)
end

Base.convert(::Type{TrackedArray}, x::TrackedArray) = x
# Base.convert(::Type{<:AbstractArray}, x::TrackedArray) = x.val

## @tracked macro

function track_call(sig)
    @assert sig.head == :call
    op = sig.args[1]
    vars_types, kw = split_params(sig)
    call_ex = nothing; ex_in_ex = nothing
    if isempty(kw)
        call_ex = Expr(:call, op,
                       [istracked(Core.eval(@__MODULE__, t)) ? :($(v).val) : v
                        for (v, t) in vars_types]...)
        ex_in_ex = Expr(:call, :Expr, QuoteNode(:call), QuoteNode(op),
                        [(istracked(Core.eval(@__MODULE__, t)) ? Expr(:., sig_var, QuoteNode(:var))
                         : sig_var) for (sig_var, t) in vars_types]...)
    else
        call_ex = Expr(:call, op, make_kw_params(kw),
                       [istracked(Core.eval(@__MODULE__, t)) ? :($(v).val) : v
                        for (v, t) in vars_types]...)
        keys = [k for (k, v) in kw]
        # we need to get this:
        # :( Expr(:call, :mult, Expr(:parameters, Expr(:kw, :pow, pow)), :(x.var)) )
        ex_in_ex = Expr(:call, :Expr, QuoteNode(:call), QuoteNode(op),
                        Expr(:call, :Expr, QuoteNode(:parameters),
                             [Expr(:call, :Expr, QuoteNode(:kw), QuoteNode(k), k)
                             for k in keys]...),
                        [(istracked(Core.eval(@__MODULE__, t)) ? Expr(:., sig_var, QuoteNode(:var))
                          : sig_var) for (sig_var, t) in vars_types]...)
    end
    defquot = quote
        function $op($(sig.args[2:end]...))
            val = $call_ex
            tv = tracked_val(x.graph, genname(), val) # TODO: x may be undefined
            ex = $ex_in_ex
            nd = ExNode{:call}(tv.var, ex; val=val)
            push!(x.graph, nd)
            return tv
        end
    end
    return defquot.args[2]
end

function track_bcast(sig)
    @assert sig.head == :.
    op = sig.args[1]
    sig_vars, types = unzip([(arg.args[1], Core.eval(@__MODULE__, arg.args[2]))
                             for arg in sig.args[2].args])
    bcast_ex = Expr(:., op, Expr(:tuple, [istracked(t) ? :($(v).val) : v
                                          for (v, t) in zip(sig_vars, types)]...))
    # we want to get this after code generation:
    # ex = :(Expr(:., :sin, Expr(:tuple, x.var, y.var)))
    ex_in_ex = Expr(:call, :Expr, QuoteNode(:.), QuoteNode(op),
                    Expr(:call, :Expr, QuoteNode(:tuple),
                         [istracked(t) ? Expr(:., sv, QuoteNode(:var)) : QuoteNode(:var)
                          for (sv, t) in zip(sig_vars, types)]...))
    defquot = quote
        function Broadcast.broadcasted(::typeof($op), $(sig.args[2].args...))
            val = $bcast_ex
            tv = tracked_val(x.graph, genname(), val) # TODO: x may be undefined
            ex = $ex_in_ex
            # println(ex)
            nd = ExNode{:bcast}(tv.var, ex; val=val)
            push!(x.graph, nd)
            return tv
        end
    end
    return defquot.args[2]
end

"""
Define a function or broadcasting rule for the specified signature which computes
the result as for ordinary (not tracked) data and writes it to the graph.

Note: this function expects at least 1 parameter of TrackedReal or TrackedArray type
with name `x`.
"""
macro tracked(sig)
    if sig.head == :call
        return track_call(sig)
    elseif sig.head == :.
        # TODO: we also need to add the op to broadcast list in `unfuse.jl`
        return track_bcast(sig)
    else
        error("Can only track calls or broadcasting")
    end
end

@tracked +(x::TrackedReal, y::TrackedReal)
@tracked -(x::TrackedReal, y::TrackedReal)
@tracked *(x::TrackedReal, y::TrackedReal)
@tracked /(x::TrackedReal, y::TrackedReal)

@tracked -(x::TrackedReal)
# @tracked (Base.:<=)(x::TrackedReal, y::TrackedReal)
# @tracked (Base.:>)(x::TrackedReal, y::TrackedReal)
# @tracked (Base.:>=)(x::TrackedReal, y::TrackedReal)
# @tracked (Base.:<)(x::TrackedReal, y::TrackedReal)

@tracked sin(x::TrackedReal)
@tracked cos(x::TrackedReal)
@tracked exp(x::TrackedReal)
@tracked log(x::TrackedReal)
@tracked abs(x::TrackedReal)
@tracked abs2(x::TrackedReal)

@tracked *(x::TrackedArray, y::TrackedArray)
@tracked maximum(x::TrackedArray)

@tracked getindex(x::TrackedArray, i::Integer)
@tracked getindex(x::TrackedArray, i::Integer, j::Integer)
@tracked getindex(x::TrackedArray, i::Integer, j::Integer, k::Integer)
@tracked getindex(x::TrackedArray, i::Integer, j::Integer, k::Integer, l::Integer)

@tracked reshape(x::TrackedArray, dims::Tuple{Int})
@tracked reshape(x::TrackedArray, dims::Tuple{Int, Int})
@tracked reshape(x::TrackedArray, dims::Tuple{Int, Int, Int})
@tracked reshape(x::TrackedArray, dims::Tuple{Int, Int, Int, Int})
@tracked reshape(x::TrackedArray, dims::Tuple{Int, Int, Int, Int, Int})
@tracked reshape(x::TrackedArray, dims::Tuple{Int, Int, Int, Int, Int, Int})

@tracked sum(x::TrackedArray)
@tracked mean(x::TrackedArray)

@tracked sin.(x::TrackedArray)
@tracked cos.(x::TrackedArray)
@tracked exp.(x::TrackedArray)
@tracked log.(x::TrackedArray)
@tracked log.(b::Integer, x::TrackedArray)

# TODO: make @nontracked macro?
# boolean operators aren't tracked
for op in [:<, :<=, :>, :>=, :(==)]
    @eval (Base.$op)(x::TrackedReal, y::TrackedReal) = $op(x.val, y.val)
end

# I couldn't find a way to overload .+ and friends as either call or broadcasting
# fortunately, the list of such unusual operations is small and fixed
function Broadcast.broadcasted(::typeof(+), x::TrackedArray, y::TrackedArray)
    val = x.val .+ y.val
    var = genname()
    nd = ExNode{:call}(var, :($(x.var) .+ $(y.var)); val=val)
    push!(x.graph, nd)
    return TrackedArray(var, val)
end
function Broadcast.broadcasted(::typeof(-), x::TrackedArray, y::TrackedArray)
    val = x.val .- y.val
    var = genname()
    nd = ExNode{:call}(var, :($(x.var) .- $(y.var)); val=val)
    push!(x.graph, nd)
    return TrackedArray(var, val)
end
function Broadcast.broadcasted(::typeof(*), x::TrackedArray, y::TrackedArray)
    val = x.val .* y.val
    var = genname()
    nd = ExNode{:call}(var, :($(x.var) .* $(y.var)); val=val)
    push!(x.graph, nd)
    return TrackedArray(var, val)
end
function Broadcast.broadcasted(::typeof(/), x::TrackedArray, y::TrackedArray)
    val = x.val ./ y.val
    var = genname()
    nd = ExNode{:call}(var, :($(x.var) ./ $(y.var)); val=val)
    push!(x.graph, nd)
    return TrackedArray(var, val)
end

# TODO: also track dot ops with scalars, e.g. x .+ 1

## utils

function tracked_exgraph(f::Function, args...)
    ctx = Dict{Any,Any}(:method => :track)
    input_vars = [genname() for i=1:length(args)]
    inputs = [iv => a for (iv, a) in zip(input_vars, args)]
    g = ExGraph(; ctx=ctx, inputs...)
    # replace default graph to capture constants to `g` as well
    og = swap_default_graph!(g)
    tr_args = [tracked_val(g, var, val) for (var, val) in inputs]
    # evaluate function with tracked args
    f(tr_args...)
    # put original graph back
    swap_default_graph!(og)
    return g
end


text:
using Documenter, PredictiveControl

makedocs(;
    modules=[PredictiveControl],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/imciner2/PredictiveControl.jl/blob/{commit}{path}#L{line}",
    sitename="PredictiveControl.jl",
    authors="Ian McInerney, Imperial College London",
    assets=String[],
)

deploydocs(;
    repo="github.com/imciner2/PredictiveControl.jl",
)


text:
module Day18
__revise_mode__ = :eval

using Advent: save_input, read_charmatrix
#using Tokenize: Tokens, tokenize, untokenize
#
#function make_tokens(tokens)
#    [(Tokens.kind(x) == Tokens.INTEGER ? parse(Int, untokenize(x)) : untokenize(x)) for x in tokens if Tokens.kind(x) != Tokens.WHITESPACE]
#end
#
#function eval_expr(expr)
#    left = popfirst!(expr)
#    if left == "("
#        eval_expr(expr)
#        left = popfirst!(expr)
#        popfirst!(expr) # )
#    end
#    isempty(expr) && return left
#
#    op = popfirst!(expr)
#    op == ")" && return left
#    op == "" && return left
#    op = Symbol(op)
#
#    right = popfirst!(expr)
#    if right == "("
#        eval_expr(expr)
#        right = popfirst!(expr)
#        popfirst!(expr) # )
#    end
#
#    res = eval(:( $op($left, $right) ))
#
#    if isempty(expr)
#        res
#    else
#        pushfirst!(expr, res)
#        expr[2] == ")" && return
#        eval_expr(expr)
#    end
#end

⊗(x, y) = x + y
⊕(x, y) = x * y

function solve_a(filename)
    #sum(x -> eval_expr(make_tokens(tokenize(x))), eachline(filename))
    sum(eachline(filename)) do line
        # same precedence
        replace(line, "*" => "⊕") |> Meta.parse |> eval
    end
end

function solve_b(filename)
    sum(eachline(filename)) do line
        replace(line, "+" => "⊗") |> x -> replace(x, "*" => "⊕") |> Meta.parse |> eval
    end
end

solve(filename = save_input(2020, 18)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# Model structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

const IntTuple = Union{Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer}}
const RealTuple = Union{Tuple{Real,Real}, Tuple{Real,Real,Real}}

export Model, ModelRho

# Object for velocity/slowness models
mutable struct Model
	n::IntTuple
	d::RealTuple
	o::RealTuple
	nb::Integer	# number of absorbing boundaries points on each side
	m	# slowness squared
end

mutable struct ModelRho
	n::IntTuple
	d::RealTuple
	o::RealTuple
	nb::Integer
	m	# slowness squared
	rho
end

"""
    Model
        n::IntTuple
        d::RealTuple
        o::RealTuple
        nb::Integer
        m::Array
        rho::Array

Model structure for seismic velocity models.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`nb`: number of absorbing boundary points in each direction

`m`: velocity model in slowness squared (s^2/km^2)

`rho`: density in (g/cm^3), (optional, default is rho=1)

Constructor
===========

The parameters `n`, `d`, `o` and `m` are mandatory. All other input arguments are optional:

    Model(n, d, o, m; nb=40, rho=ones(n))

For zero Thomsen parameters, the acoustic wave equation is solved. Otherwise the TTI equation is solved.

"""
Model(n::IntTuple,d::RealTuple,o::RealTuple,m;nb=40) = Model(n,d,o,nb,m)
Model(n::IntTuple,d::RealTuple,o::RealTuple,m,rho;nb=40) = ModelRho(n,d,o,nb,m,rho)

function Model(path::String; parameter="m")
	A = load(path)
	if parameter=="m"
		return Model(A["n"], A["d"], A["o"], A["m"])
	elseif parameter=="m0"
		return Model(A["n"], A["d"], A["o"], A["m0"])
	elseif parameter=="dm"
		return Model(A["n"], A["d"], A["o"], A["dm"])
	end
end

const ModelAll = Union{Model,ModelRho}

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
struct RandomEffectsTerm <: AbstractTerm
    lhs::StatsModels.TermOrTerms
    rhs::StatsModels.TermOrTerms
    function RandomEffectsTerm(lhs,rhs)
        if isempty(intersect(StatsModels.termvars(lhs), StatsModels.termvars(rhs)))
            if !isa(rhs, Union{CategoricalTerm,InteractionTerm{<:NTuple{N,CategoricalTerm} where N}})
                throw(ArgumentError("blocking variables (those behind |) must be Categorical ($(rhs) is not)"))
            end
            new(lhs, rhs)
        else
            throw(ArgumentError("Same variable appears on both sides of |"))
        end
    end
end

function StatsModels.apply_schema(t::FunctionTerm{typeof(/)},
                                  sch::StatsModels.FullRank,
                                  Mod::Type{<:MixedModel})
    length(t.args_parsed) == 2 ||
        throw(ArgumentError("malformed nesting term: $t " *
                            "(Exactly two arguments are supported)"))

    first, second = apply_schema.(t.args_parsed, Ref(sch.schema), Mod)
    return first + first & second
end

RandomEffectsTerm(lhs, rhs::NTuple{2, AbstractTerm}) =
    (RandomEffectsTerm(lhs, rhs[1]), RandomEffectsTerm(lhs, rhs[2]))

Base.show(io::IO, t::RandomEffectsTerm) = print(io, "($(t.lhs) | $(t.rhs))")
StatsModels.is_matrix_term(::Type{RandomEffectsTerm}) = false

function StatsModels.termvars(t::RandomEffectsTerm)
    vcat(StatsModels.termvars(t.lhs), StatsModels.termvars(t.rhs))
end

function StatsModels.apply_schema(t::FunctionTerm{typeof(|)},
                                  schema::StatsModels.FullRank,
                                  Mod::Type{<:MixedModel})
    schema = StatsModels.FullRank(schema.schema)
    lhs, rhs = t.args_parsed
    if !StatsModels.hasintercept(lhs) && !StatsModels.omitsintercept(lhs)
        lhs = InterceptTerm{true}() + lhs
    end
    lhs, rhs = apply_schema.((lhs, rhs), Ref(schema), Mod)
    RandomEffectsTerm(MatrixTerm(lhs), rhs)
end

function StatsModels.modelcols(t::RandomEffectsTerm, d::NamedTuple)
    lhs = t.lhs
    z = Matrix(transpose(modelcols(lhs, d)))
    cnames = coefnames(lhs)
    T = eltype(z)
    S = size(z, 1)
    grp = t.rhs
    m = reshape(1:abs2(S), (S, S))
    inds = sizehint!(Int[], (S * (S + 1)) >> 1)
    for j in 1:S, i in j:S
        push!(inds, m[i,j])
    end
    refs, levels = _ranef_refs(grp, d)

    ReMat{T,S}(
        grp, refs, levels, isa(cnames, String) ? [cnames] : collect(cnames),
        z, z, LowerTriangular(Matrix{T}(I, S, S)), inds,
        adjA(refs, z), Matrix{T}(undef, (S, length(levels))))
end

# extract vector of refs from ranef grouping term and data
function _ranef_refs(grp::CategoricalTerm, d::NamedTuple)
    invindex = grp.contrasts.invindex
    refs = convert(Vector{Int32}, getindex.(Ref(invindex), d[grp.sym]))
    refs, grp.contrasts.levels
end

function _ranef_refs(grp::InteractionTerm{<:NTuple{N,CategoricalTerm}},
                     d::NamedTuple) where {N}
    combos = zip(getproperty.(Ref(d), [g.sym for g in grp.terms])...)
    uniques = unique(combos)
    invindex = Dict(x => i for (i,x) in enumerate(uniques))
    refs = convert(Vector{Int32}, getindex.(Ref(invindex), combos))
    refs, uniques
end

# add some syntax to manually promote to full dummy coding
fulldummy(t::AbstractTerm) =
    throw(ArgumentError("can't promote $t (of type $(typeof(t))) to full dummy " *
                        "coding (only CategoricalTerms)"))

function fulldummy(t::CategoricalTerm)
    new_contrasts = StatsModels.ContrastsMatrix(StatsModels.FullDummyCoding(),
                                                t.contrasts.levels)
    t = CategoricalTerm(t.sym, new_contrasts)
end

fulldummy(x) =
    throw(ArgumentError("fulldummy isn't supported outside of a MixedModel formula"))

function StatsModels.apply_schema(t::FunctionTerm{typeof(fulldummy)},
                                  sch::StatsModels.FullRank,
                                  Mod::Type{<:MixedModel})
    fulldummy(apply_schema.(t.args_parsed, Ref(sch), Mod)...)
end

# specify zero correlation
struct ZeroCorr <: AbstractTerm
    term::RandomEffectsTerm
end
StatsModels.is_matrix_term(::Type{ZeroCorr}) = false

"""
    zerocorr(term::RandomEffectsTerm)

Remove correlations between random effects in `term`.
"""
zerocorr(x) = ZeroCorr(x)

function StatsModels.apply_schema(t::FunctionTerm{typeof(zerocorr)},
                                  sch::StatsModels.FullRank,
                                  Mod::Type{<:MixedModel})
    ZeroCorr(apply_schema(t.args_parsed..., sch, Mod))
end

StatsModels.modelcols(t::ZeroCorr, d::NamedTuple) =
    zerocorr!(modelcols(t.term, d))


text:
@testset "Test Finite Elements" begin

using jFEMTools
import Tensors: Vec, norm, gradient
jF = jFEMTools

# Test ContinuousLagrange Element
finiteElement = ContinuousLagrange(jF.Triangle,2)
@test jF.gettopology(finiteElement) == Dict(0=>3,1=>3,2=>0)
@test jF.getgeomtopology(finiteElement) == Dict(0=>3,1=>0,2=>0)

finiteElement = ContinuousLagrange(jF.Triangle,1)
@test finiteElement.nodal_base_coefs ≈ [0.23570226039551578 0.2357022603955159 0.2357022603955158;
                                 -0.14433756729740646 0.14433756729740646 0.0; -1/12 -1/12 1/6]
@test jF.gettopology(finiteElement) == Dict(0=>3,1=>0,2=>0)

# Test Lagrange Continuous Elements
########### 2D in Triangles #######################################
function ref_value(i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 1 && return 1. - ξ_x - ξ_y
end

ξ_ref = Vec{2}((0.5,0.5))
for j = 1:3
    @test abs(ref_value(j,ξ_ref) - jF.value(finiteElement,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value(j, ξ), ξ_ref) - jF.gradient_value(finiteElement,j,ξ_ref), Inf) < eps(Float64)
end

######################## 1D in Segments #######################################3

function ref_value1(i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return 1-ξ_x
    i == 2 && return ξ_x
end
ξ_ref = Vec{1}((0.5,))
finiteElement = ContinuousLagrange(jF.Segment,1)
@test jF.gettopology(finiteElement) == Dict(0=>2,1=>0)
for j = 1:2
    @test abs(ref_value1(j,ξ_ref) - jF.value(finiteElement,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value1(j, ξ), ξ_ref) - jF.gradient_value(finiteElement,j,ξ_ref), Inf) < eps(Float64)
end

############## 3D in Tetrahedra #############################3

finiteElement = ContinuousLagrange(jF.Tetrahedron, 1)
function ref_value3(i::Int, ξ::Vec{3})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    ξ_z = ξ[3]
    i == 1 && return 1.0 - ξ_x - ξ_y - ξ_z
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 4 && return ξ_z
end

ξ_ref = Vec{3}((0.5,0.5,0.5))
for j = 1:4
    @test abs(ref_value3(j,ξ_ref) - jF.value(finiteElement,j,ξ_ref)) < eps(Float64)
end

########################### 2D in Rectangles #################

finiteElement = ContinuousLagrange(jF.Rectangle, 1)

ξ_ref = Vec{2}((0.5,0.5))
for j = 1:4
    @test abs(0.25 - jF.value(finiteElement,j,ξ_ref)) < eps(Float64)
end

########################### 3D in Hexahedra #################

finiteElement = ContinuousLagrange(jF.Hexahedron, 1)

ξ_ref = Vec{3}((0.5,0.5,0.5))
for j = 1:4
    @test abs(0.125 - jF.value(finiteElement,j,ξ_ref)) < eps(Float64)
end

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
export NonlocalLaplacian

#############
# NonlocalLaplacian implements the nonlocal integral operator:
#
#        𝓛 u(x) := \int_{B_δ(x)} ρ_δ(|x-y|) [u(y) - u(x)] dy
#
# where B_δ(x) = {y ∈ ℝᵈ : |x-y| ≤ δ} is the closed ball of radius 0 < δ < ∞,
# and the spherically symmetric kernel ρ_δ is defined by:
#
#       ρ_δ(r) = 2*Γ(d/2+1)*(d+2-α)/(π^{d/2} δ^{d+2-α}) χ_[0,δ](r)/r^α,
#
# where χ_[0,δ](⋅) is the indicator function, and 0 ≤ α < d + 2.
#############

abstract type AbstractNonlocalCalculusOperator{S, T} <: Operator{T} end

abstract type NonlocalLaplacian{S, T} <: AbstractNonlocalCalculusOperator{S, T} end

struct ConcreteNonlocalLaplacian{S, T} <: NonlocalLaplacian{S, T}
    space::S
    α::T
    δ::T
end

ConcreteNonlocalLaplacian(sp::Space, α::T, δ::T) where T = ConcreteNonlocalLaplacian{typeof(sp), T}(sp, α, δ)

ConcreteNonlocalLaplacian(sp::Space, α, δ) = ConcreteNonlocalLaplacian{typeof(sp), promote_type(typeof(α), typeof(δ))}(sp, α, δ)

NonlocalLaplacian(sp::Space, α, δ) = ConcreteNonlocalLaplacian(sp, α, δ)
NonlocalLaplacian(d::Domain, x...) = NonlocalLaplacian(Space(d), x...)

domain(L::ConcreteNonlocalLaplacian) = domain(L.space)

domainspace(L::ConcreteNonlocalLaplacian) = L.space

getindex(::ConcreteNonlocalLaplacian{UnsetSpace,T}, k::Integer, j::Integer) where T =
    error("Spaces cannot be inferred for operator")
rangespace(L::ConcreteNonlocalLaplacian{UnsetSpace,T}) where T = UnsetSpace()

rangespace(L::ConcreteNonlocalLaplacian{Hardy{s,DD,RR}}) where {s,DD,RR} = L.space
rangespace(L::ConcreteNonlocalLaplacian{Laurent{DD,RR}}) where {DD,RR} = L.space
rangespace(L::ConcreteNonlocalLaplacian{Fourier{DD,RR}}) where {DD,RR} = L.space
rangespace(L::ConcreteNonlocalLaplacian{CosSpace{DD,RR}}) where {DD,RR} = L.space
rangespace(L::ConcreteNonlocalLaplacian{SinSpace{DD,RR}}) where {DD,RR} = L.space

bandwidths(::ConcreteNonlocalLaplacian{Hardy{s,DD,RR}}) where {s,DD,RR} = 0,0
bandwidths(::ConcreteNonlocalLaplacian{Laurent{DD,RR}}) where {DD,RR} = 0,0
bandwidths(::ConcreteNonlocalLaplacian{Fourier{DD,RR}}) where {DD,RR} = 0,0
bandwidths(::ConcreteNonlocalLaplacian{CosSpace{DD,RR}}) where {DD,RR} = 0,0
bandwidths(::ConcreteNonlocalLaplacian{SinSpace{DD,RR}}) where {DD,RR} = 0,0

##TODO: Add scale for different periods.
function getindex(L::ConcreteNonlocalLaplacian{Hardy{true, DD, RR}, T},
                  k::Integer, j::Integer) where {DD <: PeriodicSegment, RR, T}
    if k == j
        fourier_lambda(k-1, L.α, L.δ, 1)
    else
        zero(T)
    end
end

function getindex(L::ConcreteNonlocalLaplacian{Hardy{false, DD, RR}, T},
                  k::Integer, j::Integer) where {DD <: PeriodicSegment, RR, T}
    if k == j
        fourier_lambda(k, L.α, L.δ, 1)
    else
        zero(T)
    end
end

for SP in (:Laurent, :Fourier)
    @eval begin
        function getindex(L::ConcreteNonlocalLaplacian{$SP{DD, RR}, T},
                          k::Integer, j::Integer) where {DD <: PeriodicSegment, RR, T}
            if k == j
                fourier_lambda(k÷2, L.α, L.δ, 1)
            else
                zero(T)
            end
        end
    end
end

function getindex(L::ConcreteNonlocalLaplacian{CosSpace{DD, RR}, T},
                  k::Integer, j::Integer) where {DD <: PeriodicSegment, RR, T}
    if k == j
        fourier_lambda(k-1, L.α, L.δ, 1)
    else
        zero(T)
    end
end

function getindex(L::ConcreteNonlocalLaplacian{SinSpace{DD, RR}, T},
                  k::Integer, j::Integer) where {DD <: PeriodicSegment, RR, T}
    if k == j
        fourier_lambda(k, L.α, L.δ, 1)
    else
        zero(T)
    end
end

function fourier_lambda(k, α, δ, d::Integer)
    d < 1 && error("The $d-dimensional Fourier spectrum is absurd.")
    0 ≤ α < d + 2 || error("The algebraic singularity strength α = $α is not permissible.")
    if abs(k*δ) < 6
        fourier_lambda_small(k, α, δ, d)
    else
        fourier_lambda_large(k, α, δ, d)
    end
end

fourier_lambda_small(k::Real, α::Real, δ::Real, d::Integer) = fourier_lambda_small(promote(k, α, δ)..., d)

function fourier_lambda_small(k::T, α::T, δ::T, d::Integer) where T <: Real
    z = -(k*δ/2)^2
    C = (d+2-α)/2
    D = 1+T(d)/2
    S₀, S₁, err, j = one(T), one(T) + C/(C+one(T))*z/(2*D), one(T), 1
    while err > 10*abs(S₀)*eps(T)
        rⱼ = ((C+j)/(C+j+1))/((D+j)*(j+2))
        S₀, S₁ = S₁, S₁ + (S₁ - S₀)*rⱼ*z
        err = abs(S₁-S₀)
        j += 1
    end
    return -k^2*S₁
end

fourier_lambda_large(k::Real, α::Real, δ::Real, d::Integer) = fourier_lambda_large(promote(k, α, δ)..., d)

function fourier_lambda_large(k::T, α::T, δ::T, d::Integer) where T <: Real
    if norm(d+2-α) < eps(T)
        return T(-k^2)
    else
        kδ = k*δ
        Td = T(d)
        scl = (d+2-α)*gamma(Td/2+1)*2/δ^2
        cst = multif((d-α)/2, 2/kδ, Td/2)/gamma(Td/2)
        tail = 2^(Td/2)*(kδ)^(α+1-d)*((d-2-α)*besselj((Td-2)/2, kδ)*lommelS2((d-2-2α)/2, (Td-4)/2, kδ) - besselj((Td-4)/2, kδ)*lommelS2((d-2α)/2, (Td-2)/2, kδ))
        return scl*(cst+tail)
    end
end

function multif(x, y, z)
    (y^(2x)*gamma(x+1)*gamma(z)/gamma(z-x)-1)/x
end

function multif(x::Float64, y::Float64, z::Float64)
    if x ≠ 0.0
        w = 2x*log(y) + HypergeometricFunctions.lanczosapprox(1.0, x) - HypergeometricFunctions.lanczosapprox(z, -x)
        return expm1(w)/x
    else
        return 2*log(y) - MathConstants.γ + digamma(z)
    end
end

function lommelS2(μ, ν, z)
    z^(μ-1)*drummond((1-μ+ν)/2, (1-μ-ν)/2, (z/2)^2)
end

function drummond(α::T, β::T, z::T) where T
    Nlo = one(T)
    Dlo = one(T)
    Tlo = Nlo/Dlo

    cst = (α+1)*(β+1)
    if norm(cst) < eps(real(T))
        return Tlo
    end
    Nmid = (z+α+β+1)/cst
    Dmid = (z+cst)/cst
    Tmid = Nmid/Dmid

    cst = (α+2)*(β+2)
    if norm(cst) < eps(real(T))
        return Tmid
    end
    Dhi = ((z+α+β+3+cst)*Dmid - (α+β+3)*Dlo)/cst
    Nhi = ((z+α+β+3+cst)*Nmid - (α+β+3)*Nlo)/cst
    Thi = Nhi/Dhi

    k = 0
    while (abs(Thi-Tmid) > 10*abs(Thi)*eps(real(T)) || abs(Tmid-Tlo) > 10*abs(Tmid)*eps(real(T))) && k < 10_000
        cst = (α+k+3)*(β+k+3)
        Nhi, Nmid, Nlo = ((z+(k+2)*(α+β+2*k+5)+cst)*Nhi - (k+2)*(α+β+3*k+6)*Nmid + (k+2)*(k+1)*Nlo)/cst, Nhi, Nmid
        Dhi, Dmid, Dlo = ((z+(k+2)*(α+β+2*k+5)+cst)*Dhi - (k+2)*(α+β+3*k+6)*Dmid + (k+2)*(k+1)*Dlo)/cst, Dhi, Dmid
        Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
        k += 1
    end
    return isnan(Thi) ? isnan(Tmid) ? Tlo : Tmid : Thi
end


text:
"""
    minor_certain(freq, n_checked, z)
Check if minor allele is certain.

- `freq`: frequency of minor or major allele
- `n_checked`: number of individuals checked so far
- `z`: cutoff of "z" value, defaults to `5.0`
"""
function minor_certain(freq::Float64, n_checked::Integer, z::Float64=5.0)
    delta = z * sqrt(freq * (1-freq) / n_checked)
    return !(freq - delta < 0.5 && freq + delta > 0.5)
end


text:
""" Matches bbls to associated variables (except for dependent variables in Regressors). """
function match_bbls_to_vars(bbls::Array,
                            vars::Array = JuMP.all_variables(bbls))
    # TODO: improve types.
    bbl_to_var = Dict()
    for bbl in bbls
        unb = false
        for var in vars
            if var in bbl.vars
                if unb
                    push!(bbl_to_var[bbl], var)
                else
                    bbl_to_var[bbl] = [var]
                    unb = true
                end
            end
        end
    end
    # Sorting so that we penalize both in the number of unbounded and total variables.
    bbl_to_var = sort(collect(bbl_to_var), by= x->length(x[2])*length(x[1].vars)^0.5)
    return bbl_to_var
end

match_bbls_to_vars(gm::GlobalModel, vars::Array = JuMP.all_variables(gm)) = match_bbls_to_vars(gm.bbls, vars)


text:
using OrdinaryDiffEq, DifferentialEquations

tspan = (0.0f0, 50.0f0)

"""
Simple model of a spiking neuron, from the paper:
    Izhikevich, Eugene M. "Simple model of spiking neurons."
    IEEE Transactions on neural networks 14, no. 6 (2003): 1569-1572.
"""
function neuron_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    Ft = U_func_time(t)
    Fs = U_func_space(u)
    a, b = p
    v, u_var = u

    du[1] = 0.04*(v^2) + 5*v + 140 - u_var
    du[2] = a*(b*v - u_var)
    # TODO: Why is this breaking?
    # du .+= Ft .+ Fs
    tmp = Ft .+ Fs
    du[1] += tmp[1]
    du[2] += tmp[2]
end

# This system also has a callback function
#   Default values define a "Regular Spiking" neuron
condition(u, t, integrator) = u[1] - 30
function affect!(integrator; c=-65, d=8)
    integrator.u[1] = c     # Reset
    integrator.u[2] += d    # Jump
end
neuron_reset = ContinuousCallback(condition, affect!,
                    save_positions = (false,false))

# Generate data, with a time component
p = [0.02, 0.2]
u0 = [-65.0, 0.0]
ts = range(tspan[1], tspan[2], length=5000)

function solve_neuron_system(;U_func_time=U_func_time_trivial,
                              U_func_space=U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->neuron_system(du, u, p, t,
                                U_func_time=U_func_time,
                                U_func_space=U_func_space),
                    u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=ts,
                callback=neuron_reset);
    return sol
end

#####
##### Turing function for Bayesian parameter estimation
#####
# Unforced version for use with Turing.jl
function neuron_system(u, p, t)
    u = convert.(eltype(p),u)

    a, b = p
    v, u = u
    du = [0.04(v^2) + 5v + 140 - u; a*(b*v - u)]
    return du
end

#####
##### True model in SINDy syntax
#####
a, b = p
#      x    y  c   xx    xy yy
 A = [[5   -1  140  0.04  0  0];
      [a*b -a  0    0     0  0]]
n = size(A, 1)
sindy_library = Dict(
    "cross_terms"=>2,
    "constant"=>nothing
);
core_dyn_true = sindyModel(ts, A,
                            convert_string2function(sindy_library),
                            ["x", "y"])

# Actually export
export ts, solve_neuron_system, neuron_system, core_dyn_true
    # ,lorenz_grad_residual


text:
println("Starting")

using MetaGen
using JSON
import YAML
using Pipe: @pipe
using Random

config_path = ARGS[1]
config = YAML.load_file(config_path)
output_dir = "results/$(config["experiment_name"])_" * ENV["SLURM_JOB_ID"]
mkdir(output_dir)

include("useful_functions.jl")
dict = []
#for i = 0:config["batches_upto"]
	#to_add =  @pipe "$(config["input_file_dir"])$(i)_data_labelled.json" |> open |> read |> String |> JSON.parse
to_add =  @pipe "$(config["input_file_dir"])data_labelled.json" |> open |> read |> String |> JSON.parse
append!(dict, to_add)
#end
#dict = @pipe "../../scratch_work_07_16_21/0_data_labelled.json" |> open |> read |> String |> JSON.parse
#dict = @pipe "../../scratch_work_07_16_21/0_data_labelled.json" |> open |> read |> String |> JSON.parse

#Random.seed!(15)
#try to make objects_observed::Array{Array{Array{Array{Detection2D}}}} of observed objects.
#outer array is for scenes, then frames, the receptive fields, then last is an array of detections

################################################################################
num_videos = config["num_videos"]
num_frames = config["num_frames"]
threshold = config["threshold"]

params = Video_Params(n_possible_objects = 5)
top_n = 5

receptive_fields = make_receptive_fields(params)
objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold, top_n)

num_particles = config["num_particles"]
mcmc_steps_outer = config["mcmc_steps_outer"]
mcmc_steps_inner = config["mcmc_steps_inner"]

################################################################################
#=
#Set up the output file
online_outfile = output_dir * "/online_output.csv"
online_file = open(online_outfile, "w")
file_header(online_file)

################################################################################
#Online MetaGen
println("start pf for online")

#@profilehtml unfold_particle_filter(false, num_particles, objects_observed, camera_trajectories, params, file)
traces, inferred_world_states, avg_v = unfold_particle_filter(nothing,
	num_particles, mcmc_steps_outer, mcmc_steps_inner, objects_observed,
	camera_trajectories, params, online_file)
close(online_file)

println("done with pf for online")

## Commented out for testing

################################################################################
#Retrospective MetaGen
println("start retrospective")

#Set up the output file
retro_outfile = output_dir * "/retrospective_output.csv"
retro_file = open(retro_outfile, "w")
file_header(retro_file)

unfold_particle_filter(avg_v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
	objects_observed, camera_trajectories, params, retro_file)
close(retro_file)

println("done with pf for retrospective")

################################################################################
#run Lesioned MetaGen
=#

shuffle_type = config["shuffle_type"]

num_videos_train = convert(Int64, num_videos/2)

if shuffle_type==0
	order = collect(1:num_videos_train)
elseif shuffle_type==1
	Random.seed!(1)
	order = shuffle(1:num_videos_train)
else shuffle_type==2
	Random.seed!(2)
	order = shuffle(1:num_videos_train)
end

training_objects_observed = objects_observed[order, :]
training_camera_trajectories = camera_trajectories[order, :]
#training set and test set
input_objects_observed = vcat(objects_observed[order, :], objects_observed[(num_videos_train+1):num_videos, :])
input_camera_trajectories = vcat(camera_trajectories[order, :], camera_trajectories[(num_videos_train+1):num_videos, :])

#Set up the output file
lesioned_V_file = open(output_dir * "/lesioned_V.csv", "w")
file_header_V(lesioned_V_file, params)
lesioned_ws_file = open(output_dir * "/lesioned_ws.csv", "w")
file_header_ws(lesioned_ws_file, params, num_particles)

v = zeros(length(params.possible_objects), 2)
v[:,1] .= 1.0
v[:,2] .= 0.5
unfold_particle_filter(v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
	input_objects_observed, input_camera_trajectories, params, lesioned_V_file,
	lesioned_ws_file)
close(lesioned_V_file)
close(lesioned_ws_file)

println("done with pf for lesioned metagen")

#=

################################################################################
#for writing an output file for a demo using MetaGen

###### add to dictionary
out = write_to_dict(dict, camera_trajectories, inferred_world_states, num_videos, num_frames)

#open("../../scratch_work_07_16_21/output_tiny_set_detections.json","w") do f
open(output_dir * "/output.json","w") do f
	JSON.print(f,out)
end

println("finished writing json")
=#

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using DrWatson
@quickactivate "EPExperiment"

using EPExperiments
using Quadrature

const LOG_DIR_NAME = "bananas_ground_truth"

function bananas_ground_truth()
    result_folder = datadir(LOG_DIR_NAME)
    exp_banana = banana()
    turing_banana = exp_banana.gamma2

    logjoint = AnnealedIS.make_log_joint_density(turing_banana)
    gamma2(x, p) = exp(logjoint((x1 = x[1], x2 = x[2])))
    prob1 = QuadratureProblem(gamma2, [-50, -50], [50, 50])
    Z2 = solve(prob1, HCubatureJL(), abstol=1e-20)
    @show Z2.u

    gamma1(x, p) = exp(logjoint((x1 = x[1], x2 = x[2]))) * banana_f(x[1], x[2])
    prob2 = QuadratureProblem(gamma1, [-30, -40], [30, 10])
    Z1 = solve(prob2, HCubatureJL(), reltol=1e-10, abstol=1e-30)
    @show Z1.u

    @show Z1.u / Z2.u

    @tagsave(joinpath(result_folder, "ground_truth.jld2"), Dict(
        "expectation_estimate" => Z1.u / Z2.u
    ))
end

bananas_ground_truth()


text:
Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
[@LorenaABarba](https://twitter.com/LorenaABarba)

12 steps to Navier–Stokes
=====
***

For a moment, recall the Navier–Stokes equations for an incompressible fluid, where $\vec{v}$ represents the velocity field:

$$
\begin{eqnarray*}
\nabla \cdot\vec{v} &=& 0 \\
\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v} &=& -\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v}
\end{eqnarray*}
$$

The first equation represents mass conservation at constant density. The second equation is the conservation of momentum. But a problem appears: the continuity equation for incompressble flow does not have a dominant variable and there is no obvious way to couple the velocity and the pressure. In the case of compressible flow, in contrast, mass continuity would provide an evolution equation for the density $\rho$, which is coupled with an equation of state relating $\rho$ and $p$.

In incompressible flow, the continuity equation $\nabla \cdot\vec{v}=0$ provides a *kinematic constraint* that requires the pressure field to evolve so that the rate of expansion $\nabla \cdot\vec{v}$ should vanish everywhere. A way out of this difficulty is to *construct* a pressure field that guarantees continuity is satisfied; such a relation can be obtained by taking the divergence of the momentum equation. In that process, a Poisson equation for the pressure shows up!

Step 10: 2D Poisson Equation
----
***

Poisson's equation is obtained from adding a source term to the right-hand-side of Laplace's equation:

$$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$

So, unlinke the Laplace equation, there is some finite value inside the field that affects the solution. Poisson's equation acts to "relax" the initial sources in the field.

In discretized form, this looks almost the same as [Step 9](./12_Step_9.ipynb), except for the source term:

$$\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}$$

As before, we rearrange this so that we obtain an equation for $p$ at point $i,j$. Thus, we obtain:

$$p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}$$

We will solve this equation by assuming an initial state of $p=0$ everywhere, and applying boundary conditions as follows:

$p=0$ at $x=0, \ 2$ and $y=0, \ 1$

and the source term consists of two initial spikes inside the domain, as follows:

$b_{i,j}=100$ at $i=\frac{1}{4}nx, j=\frac{1}{4}ny$

$b_{i,j}=-100$ at $i=\frac{3}{4}nx, j=\frac{3}{4}ny$

$b_{i,j}=0$ everywhere else.

The iterations will advance in pseudo-time to relax the initial spikes. The relaxation under Poisson's equation gets slower and slower as they progress. *Why?*

Let's look at one possible way to write the code for Poisson's equation. As always, we load our favorite Python libraries. We also want to make some lovely plots in 3D. Let's get our parameters defined and the initialization out of the way. What do you notice of the approach below?

```python
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
```

```python
# Parameters
nx = 50
ny = 50
nt  = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
p  = numpy.zeros((ny, nx))
pd = numpy.zeros((ny, nx))
b  = numpy.zeros((ny, nx))
x  = numpy.linspace(xmin, xmax, nx)
y  = numpy.linspace(xmin, xmax, ny)

# Source
b[int(ny / 4), int(nx / 4)]  = 100
b[int(3 * ny / 4), int(3 * nx / 4)] = -100
```

With that, we are ready to advance the initial guess in pseudo-time. How is the code below different from the function used in [Step 9](./12_Step_9.ipynb) to solve Laplace's equation?

```python
for it in range(nt):

    pd = p.copy()

    p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
                    (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
                    b[1:-1, 1:-1] * dx**2 * dy**2) /
                    (2 * (dx**2 + dy**2)))

    p[0, :] = 0
    p[ny-1, :] = 0
    p[:, 0] = 0
    p[:, nx-1] = 0
```

Maybe we could reuse our plotting function from [Step 9](./12_Step_9.ipynb), don't you think?

```python
def plot2D(x, y, p):
    fig = pyplot.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = numpy.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

```

```python
plot2D(x, y, p)
```

![png](output_14_0.png)

Ah! The wonders of code reuse! Now, you probably think: "Well, if I've written this neat little function that does something so useful, I want to use it over and over again. How can I do this without copying and pasting it each time? —If you are very curious about this, you'll have to learn about *packaging*. But this goes beyond the scope of our CFD lessons. You'll just have to Google it if you really want to know.

***

## Learn More

To learn more about the role of the Poisson equation in CFD, watch **Video Lesson 11** on You Tube:

```python
from IPython.display import YouTubeVideo
YouTubeVideo('ZjfxA3qq2Lg')
```

<iframe
    width="400"
    height="300"
    src="https://www.youtube.com/embed/ZjfxA3qq2Lg"
    frameborder="0"
    allowfullscreen
></iframe>

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    h2 {
        font-family: 'Fenix', serif;
    }
    h3{
		font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
       }
	h4{
		font-family: 'Fenix', serif;
       }
    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    div.text_cell_render{
        font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 135%;
        font-size: 120%;
        width:600px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro";
			font-size: 90%;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h1 {
        font-weight: 200;
        font-size: 50pt;
		line-height: 100%;
        color:#CD2305;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #CD2305;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>
<script>
    MathJax.Hub.Config({
                        TeX: {
                           extensions: ["AMSmath.js"]
                           },
                tex2jax: {
                    inlineMath: [ ['$','$'], ["\\(","\\)"] ],
                    displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
                },
                displayAlign: 'center', // Change this to 'center' to center equations.
                "HTML-CSS": {
                    styles: {'.MathJax_Display': {"margin": 4}}
                }
        });
</script>

> (The cell above executes the style for this notebook.)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
@testset "load_timeseries_data" begin
    # load one and check if all the fields are correct
    ts_input_data = load_timeseries_data(:CEP_GER1)
    @test ts_input_data.K == 366
    @test ts_input_data.T == 24
    @test all(ts_input_data.weights .== 1.)

    # load all four by name - just let them run and see if they run without error
    ts_input_data = load_timeseries_data(:DAM_GER)
    ts_input_data = load_timeseries_data(:DAM_CA)
    ts_input_data = load_timeseries_data(:CEP_GER1)
    ts_input_data = load_timeseries_data(:CEP_GER18)

    #load a folder by path
    ts_input_data = load_timeseries_data(normpath(joinpath(dirname(@__FILE__),"..","data","TS_GER_1")))

    # load single file by path
    ts_input_data = load_timeseries_data(normpath(joinpath(dirname(@__FILE__),"..","data","TS_GER_1","solar.csv")))
    @test all(["solar-germany"] .== keys(ts_input_data.data) )
end


text:
# This file is a part of DataDeps.jl. License is MIT.

function unpack_cmd(file,directory,extension,secondary_extension)
    p7zip() do exe7z
        if secondary_extension == ".tar" || extension == ".tgz" || extension == ".tbz"
            # special handling for compressed tarballs
            return pipeline(`$exe7z x $file -y -so`, `$exe7z x -si -y -ttar -o$directory`)
        else
            return `$exe7z x $file -y -o$directory`
        end
    end
end

"""
    unpack(f; keep_originals=false)

Extracts the content of an archive in the current directory;
deleting the original archive, unless the `keep_originals` flag is set.
"""
function unpack(f; keep_originals=false)
    run(unpack_cmd(f, pwd(), last(splitext(f)), last(splitext(first(splitext(f))))))
    rm("pax_global_header"; force=true)  # Non-compliant tarball extractors dump out this file. It is meaningless (google it's filename for more)
    !keep_originals && rm(f)
end


text:
macro maintain_type(expr)    # ensures that subtypes of StaticArrays.FieldVector maintain their type upon manipulation
    full_name = expr.args[2].args[1]                                     # MyType{T} or MyType
    short_name = full_name isa Symbol ? full_name : full_name.args[1]    # MyType
    @assert expr.args[2].args[2].args[1] == :FieldVector "@maintain_type is defined only for subtypes of FieldVector"
    N = expr.args[2].args[2].args[2]                                     # N in FieldVector{N,T}
    T = expr.args[2].args[2].args[3]                                     # T in FieldVector{N,T}
    st_method = if short_name == full_name
        :(StaticArrays.similar_type(::Type{$short_name}, ::Type{$T}, s::Size{($N,)}) = $short_name)
    else
        :(StaticArrays.similar_type(::Type{<:$short_name}, ::Type{$T}, s::Size{($N,)}) where {$T} = $full_name)
    end
    quote
        $(esc(expr))
        $(esc(st_method))
    end
end

function ode_heun(fn, y0, Tf, Ti=zero(Tf), N=10)
    dt = (Tf - Ti)/N
    y = y0
    t = Ti
    for i in 1:N
        k1 = dt*fn(y, t)
        k2 = dt*fn(y + k1, t + dt)
        y = y + (k1 + k2)/2
        t = t + dt
    end
    y
end

function ode_rk4(fn, y0, Tf, Ti=zero(Tf), N=10)
    dt = (Tf - Ti)/N
    y = y0
    t = Ti
    for i in 1:N
        k1 = dt*fn(y, t)
        k2 = dt*fn(y + k1/2, t + dt/2)
        k3 = dt*fn(y + k2/2, t + dt/2)
        k4 = dt*fn(y + k3, t + dt)
        y = y + (k1 + 2*k2 + 2*k3 + k4)/6
        t = t + dt
    end
    y
end

@inline mod2piF(x::T) where {T} = mod(x, T(2*pi))
@inline wrap_pi(x) = mod2piF(x + pi) - pi
@inline adiff(x, y) = wrap_pi(x - y)