text:
# Mapping from variable SymbolicRandVars to DReal Ex{}
typealias SymbToVar Dict{SymbolicRandVar, DReal.Ex}

"A Random variable represented as a expression in DReal"
type DRealRandVar{T} <: RandVar{T}
  ex::DReal.Ex{T}       # The random variable
  ctx::DReal.Context    # Contet
  sym_to_var::SymbToVar
end

dims(X::DRealRandVar) = union(map(dims, keys(X.sym_to_var))...)::Set{Int}

## Conversion between Symbolic RandVar and Dreal expression
## ========================================================

# The general rules are
# 1. If I'm a composite RandVar then expand all my children
# 2. If I'm an elementary randvar with constant arguments, do the interval trick
# 3. If I'm an elementary randvar with random parameters and closed fom quantile, then treat as 1
# 4. If I'm an elementary randvar with random parameters and no-closed form quantile then treat children as 1 and me as 1

function set_precision!(Y::DRealRandVar, precision::Float64)
  println("Setting DReal RandVar precision to $precision")
  DReal.set_precision!(Y.ctx, precision)
end

"Construct DReal.RandVar from Symbolic RandVar"
function convert{T}(::Type{DRealRandVar{T}}, X::SymbolicRandVar{T})
  ctx = Context(qf_nra)
  sym_to_var = SymbToVar()
  ex = expand(X, sym_to_var, ctx)
  DRealRandVar(ex, ctx, sym_to_var)
end

# 1. If I'm a composite RandVar then expand all my children
for (name, op) in all_functional_randvars
  eval(
  quote
  # Real
  function expand(X::$name, sym_to_var::SymbToVar, ctx::DReal.Context)
    rvargs = [expand(arg,sym_to_var,ctx) for arg in args(X)]
    ($op)(ctx, rvargs...)
  end
  end)
end

# ## Ambiguity
# function expand(X::ClosedFormQuantileRandVar, sym_to_var::SymbToVar, ctx::DReal.Context, args::ConstantRandVar...)
#   expr = quantile_expr(X)
#   expand(expr, sym_to_var, ctx::DReal.Context, args(expr)...)
# end

# 3. If I'm an elementary randvar with random parameters and closed fom quantile
function expand{T}(X::ClosedFormQuantileRandVar{T}, sym_to_var::SymbToVar, ctx::DReal.Context)
  # If all args are constant do the interval trick
  if all([isa(arg, ConstantRandVar) for arg in args(X)])
    if haskey(sym_to_var, X)
      sym_to_var[X]
    else
      # Auxilary variable unbounded
      sym_to_var[X] = Var(ctx, T)
    end
  # Otherwise expand it symbolically
  else
    expand(quantile_expr(X), sym_to_var, ctx::DReal.Context)
  end
end

# # 2. If not closedform elementary randvar , do the interval trick
# if arguments not constant then we'll get error at bound stage
function expand{T}(X::ElementaryRandVar{T}, sym_to_var::SymbToVar, ctx::DReal.Context)
  if haskey(sym_to_var, X)
    sym_to_var[X]
  else
    # Auxilary variable unbounded
    sym_to_var[X] = Var(ctx, T)
  end
end

# 5. If Constant just return value
expand(X::ConstantRandVar, sym_to_var::SymbToVar, ctx::DReal.Context) = X.val

# 6. If Omega just return value
function expand(X::OmegaRandVar, sym_to_var::SymbToVar, ctx::DReal.Context)
  if haskey(sym_to_var, X)
    sym_to_var[X]
  else
    sym_to_var[X] = Var(ctx, Float64, "omega$(X.dim)", 0.0, 1.0)
  end
end

## Calling DReal RandVars with input sets
## ======================================

# FIXME: A should be an interval or a real, need a type for this, or do i
function quantile{T<:Distribution}(::Type{T}, A, args...)
  dist = T(args...)
  quantile(dist, A)
end

"Returns lower and upper bounds for auxilary variable as function of event A"
function bounds(X::OmegaRandVar, A::AbstractOmega)
  A[X.dim]
end

"Returns lower and upper bounds for Elementary RandVar as function of event A"
function bounds{T <: ElementaryRandVar}(X::T, A::AbstractOmega)
  arg_bounds = [bounds(arg, A) for arg in args(X)]
  quantile(distribution_type(T), A[X.dim], arg_bounds...)
end

"Returns lower and upper bounds for Constant RandVar as function of event A"
bounds(X::ConstantRandVar, A::AbstractOmega) = X.val

# FIXME: maybe lb and upper bound should be T but quantiles return Float and DReal has no mk_num for Ints
function add_bound_constraints!{T}(ctx::DReal.Context, X::DReal.Ex{T}, lb::Float64, ub::Float64)
  # DReal currently has trouble with infinities, so just dont add the contraint
  # Since we know it is unbounded anyway
  if lb != -Inf
    lb_constraint = (>=)(ctx, X, lb)
    DReal.add!(ctx, lb_constraint)
    # println("lb_constraint = (>=)(ctx, X, $lb)")
    # println("DReal.add!(ctx, lb_constraint)")
  end

  if ub != Inf
    ub_constraint = (<=)(ctx, X, ub)
    DReal.add!(ctx, ub_constraint)
    # println("ub_constraint = (<=)(ctx, X, $ub)")
    # println("DReal.add!(ctx, ub_constraint)")
  end
end

function is_sat(ex::Ex{Bool}, X::DRealRandVar{Bool}, A::AbstractOmega)
  # println("Testing issat")
  ctx = X.ctx
  push_ctx!(ctx)
  # println("push_ctx!(ctx)")
  # push!(debugstring, "(push 1)")

  ## Define subset (box) of Omega using assertions
  ## For each variable - add constraints on its lower and upper bounds
  for (symb, vari) in X.sym_to_var
    interval = bounds(symb, A)
    add_bound_constraints!(ctx, vari, interval.l, interval. u)
  end

  DReal.add!(ctx, ex)
  # println("DReal.add!(ctx, ex)")
  result = is_satisfiable(ctx)
  # println("result = is_satisfiable(ctx)")

  # push!(debugstring, "(pop 1)")
  pop_ctx!(ctx) #undo from 2 to here
  # println("pop_ctx!(ctx)")
  result
end

"""For a given event ``A`` finds X(A), i.e. are there points in a \in A
  such that X(A) is true, false or both"""
function call(X::DRealRandVar{Bool}, A::AbstractOmega)
  # 1. ∃ω ∈ A ∩ X : Does A contain any point X?
  # println("pos case")
  pos_case = is_sat(X.ex, X, A)

  # ∃ω ∈ A \ X : Does A contain any point not in X?
  # println("neg case2")
  notex = (!)(X.ctx, X.ex)
  # FIXME: Store this Not Ex
  # println("ex = (!)(X.ctx, ex)")
  neg_case = is_sat(notex, X, A)
  # println("Tested both cases: $pos_case $neg_case")

  if pos_case & neg_case tf
  elseif pos_case t
  elseif neg_case f
  else
    # print(debugstring)
    error("Solver error: Query or its negation must be true")
  end
end

"Returns some subet A ⊆ init_box s.t. Y(A) == true.  Also returns probability it sampled that box"
function preimage_proposal{D <: Domain}(Y::DRealRandVar{Bool}, init_box::D; args...)
  DReal.push_ctx!(Y.ctx)

  ## Define subset (box) of Omega using assertions
  for (symb, vari) in Y.sym_to_var
    interval = bounds(symb, init_box)
    add_bound_constraints!(Y.ctx, vari, interval.l, interval. u)
  end
  DReal.add!(Y.ctx, Y.ex)

  issat = DReal.is_satisfiable()

  !issat && error("Cannot condition on unsatisfiable events")

  A = LazyBox(Float64) #FIXME Float64 too speific?
  for (symb, var) in Y.sym_to_var
    A[symb.dim] = DReal.model(Y.ctx, var)
  end

  # I believe all the boxes are going to have the same size, so assuming that's true
  # we can just return some constant for the size.

  # As for the box itself, we need to return an abstract omega.
  dummy_ex = first(Y.sym_to_var)[2]
  logq = DReal.opensmt_get_model_logprob(Y.ctx.ctx, dummy_ex.e)
  DReal.pop_ctx!(Y.ctx)
  A, logq, 1.0
end

cleanup(Y::DRealRandVar) = DReal.delete_ctx!(Y.ctx)


text:
import Base: Dates
import Base.Dates: Hour

int_expr = is(Int, Int32) ? :(Core.Int32) : :(Core.Int64)
p = @patch f(a, b::Int64, c=3, d::Integer=4; e=5, f::Int=6) = nothing
@test p.signature == :(f(a, b::Core.Int64, c=3, d::Core.Integer=4; e=5, f::$int_expr=6))
@test p.modules == Set([:Core])

p = @patch f(a::Integer...) = nothing
@test p.signature == :(f(a::Core.Integer...))
@test p.modules == Set([:Core])

patches = [
    @patch f(h::Base.Dates.Hour=Base.Dates.Hour(rand())) = nothing
    @patch f(h::Dates.Hour=Dates.Hour(rand())) = nothing
    @patch f(h::Hour=Hour(rand())) = nothing
]
for p in patches
    @test p.signature == :(f(h::Base.Dates.Hour=Base.Dates.Hour(Base.Random.rand())))
    @test p.modules == Set([:(Base.Dates), :(Base.Random)])
end

ext: jl


text:
using .DOS
using SpecialFunctions

"""
    fermi_dirac(E, kT=0.026)

Compute the value of the Fermi-Dirac distribution at energy `E` (relative to the Fermi energy) and thermal energy `kT`.
"""
fermi_dirac(E; kT = 0.026) = inv.(1 .+ exp.(E ./ kT))

abstract type KineticModel end

# dispatch for net rates
(km::KineticModel)(V_app; kT = 0.026) =
    abs.(km(V_app, true; kT = kT) - km(V_app, false; kT = kT))

# return a new one with a scaled prefactor
import Base.*
function *(c::Real, km::KineticModel)
    new_A = c*km.A
    other_params = getfield.([km], propertynames(km))[2:end]
    typeof(km)(new_A, other_params...)
end

# generic pretty printing
function Base.show(io::IO, m::KineticModel)
    s = repr(typeof(m)) * "("
    for field in propertynames(m)
        s =
            s *
            string(field) *
            "=" *
            string(round(getproperty(m, field), sigdigits = 3)) *
            ", "
    end
    s = s[1:end-2] * ")"
    print(io, s)
end

"""
    ButlerVolmer(A, α)
    ButlerVolmer(A)

Computes Butler-Volmer kinetics.

If initialized with one argument, assumes symmetric electron transfer (α=0.5) and sets this to be the prefactor A. Note that this prefactor implicitly contains information about equilibrium activation energies, as well as geometric information.
"""
struct ButlerVolmer <: KineticModel
    A::Float64
    α::Float64
end

# default to unit prefactor and symmetric response
ButlerVolmer() = ButlerVolmer(1.0, 0.5)
ButlerVolmer(A) = ButlerVolmer(A, 0.5)

function (bv::ButlerVolmer)(V_app, ox::Bool; kT::Real = 0.026)
    local exp_arg
    if ox
        exp_arg = (bv.α .* V_app) ./ kT
    else
        exp_arg = -((1 .- bv.α) .* V_app) ./ kT
    end
    bv.A .* exp.(exp_arg)
end

"""
    AsymptoticMarcusHushChidsey(A, λ)
    AsymptoticMarchusHushChidsey(λ)

Computes asymptotic solution to MHC model, as described in Zeng et al.: 10.1016/j.jelechem.2014.09.038d

If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor to 1.0.
"""
struct AsymptoticMarcusHushChidsey <: KineticModel
    A::Float64
    λ::Float64
end

# default prefactor is 1
AsymptoticMarcusHushChidsey(λ) = AsymptoticMarcusHushChidsey(1.0, λ)

function (amhc::AsymptoticMarcusHushChidsey)(V_app, ox::Bool; kT::Real = 0.026)
    η = (2 * ox - 1) .* V_app ./ kT
    λ_nondim = amhc.λ / kT
    a = 1 + sqrt(λ_nondim)
    arg = (λ_nondim .- sqrt.(a .+ η^2)) ./ (2 * sqrt(λ_nondim))
    pref = sqrt(π * λ_nondim) / (1 + exp(-η))
    return amhc.A .* pref .* erfc.(arg)
end

abstract type IntegralModel <: KineticModel end # "Marcus-like"

# check to catch missed dispatches for new types
# integrand(km::IntegralModel, V_dl, ox::Bool; kwargs...) =
#     error("An integral-based kinetic model must dispatch the `integrand` function!")

# TODO: check that this passes through both kT and V_q appropriately
# dispatch for net rates
integrand(km::IntegralModel, V_dl; kwargs...) =
    E -> abs.(
        integrand(km, V_dl, true; kwargs...)(E) - integrand(km, V_dl, false; kwargs...)(E),
    )

"""
    Marcus(A, λ)
    Marcus(λ)

Computes Marcus kinetics.

If initialized with one argument, assumes this to be the reorganization energy λ and sets the prefactor A to 1.0.
"""
struct Marcus <: KineticModel
    A::Float64
    λ::Float64
end

# default prefactor is 1
Marcus(λ) = Marcus(1.0, λ)

function (m::Marcus)(V_app, ox::Bool; kT::Real = 0.026)
    local exp_arg
    if ox
        exp_arg = -(m.λ .+ V_app).^2 ./ (4 .* m.λ .* kT)
    else
        exp_arg = -(m.λ .- V_app).^2 ./ (4 .* m.λ .* kT)
    end
    m.A .* exp.(exp_arg)
end

"""
    MarcusHushChidsey(A, λ, average_dos)
    MarcusHushChidsey(λ, average_dos)

Computes Marcus-Hush-Chidsey kinetics: 10.1126/science.251.4996.919

Note that strictly speaking, `average_dos` and the prefactor `A` are redundant. They are both included primarily to facilitate comparisons with similarly parametrized `Marcus` models.
"""
struct MarcusHushChidsey <: IntegralModel
    A::Real
    λ::Real
    average_dos::Real
end

# default prefactor is 1
MarcusHushChidsey(λ, avg_dos) = MarcusHushChidsey(1.0, λ, avg_dos)
# convert more detailed DOS information to just pull out average
MarcusHushChidsey(A, λ, dd::DOSData) = MarcusHushChidsey(A, λ, dd.average_value)
MarcusHushChidsey(A, λ, dos_file::Union{Matrix,String}; kwargs...) =
    MarcusHushChidsey(A, λ, DOSData(dos_file; kwargs...))

# TODO: Check that both this and +DOS versions still match original paper
function integrand(mhc::MarcusHushChidsey, V_dl, ox::Bool; kT::Real = 0.026)
    function marcus_term(E)
        local exp_arg
        if ox
            exp_arg = -( ((E .- (mhc.λ .+ V_dl)) .^ 2)) ./ (4 * mhc.λ * kT)
        else
            exp_arg = -( ((E .- (mhc.λ .- V_dl)) .^2)) ./ (4 * mhc.λ * kT)
        end
        exp.(exp_arg)
    end
    E -> (mhc.A * mhc.average_dos) .* marcus_term(E) .* fermi_dirac(E; kT = kT)
end

"""
    MarcusHushChidseyDOS(A=1.0, λ, dos)
    MarcusHushChidseyDOS(A=1.0, λ, dos_file)

Computes Marcus-Hush-Chidsey + DOS kinetics as described in Kurchin and Viswanathan: 10.1063/5.0023611
"""
struct MarcusHushChidseyDOS <: IntegralModel
    A::Float64
    λ::Float64
    dos::DOSData
end

function Base.show(io::IO, mhcd::MarcusHushChidseyDOS)
    s = repr(typeof(mhcd)) * "("
    s *= "A=$(round(mhcd.A, sigdigits=3)), λ=$(round(mhcd.λ, sigdigits=3)))"
    print(io, s)
end

# default prefactor is 1
MarcusHushChidseyDOS(λ, dd::DOSData) = MarcusHushChidseyDOS(1.0, λ, dd)

MarcusHushChidseyDOS(A, λ, dos_file::Union{Matrix,String}; kwargs...) =
    MarcusHushChidseyDOS(A, λ, DOSData(dos_file; kwargs...))

MarcusHushChidseyDOS(λ, dos_file::Union{Matrix,String}; kwargs...) = MarcusHushChidseyDOS(1.0, λ, DOSData(dos_file; kwargs...))

function integrand(
    mhcd::MarcusHushChidseyDOS,
    V_dl,
    ox::Bool;
    kT::Real = 0.026,
    V_q = 0.0,
)
    function marcus_term(E)
        local exp_arg
        if ox
            exp_arg = -(( (mhcd.λ .- V_dl) .+ E) .^ 2) ./ (4 * mhcd.λ * kT)
        else
            exp_arg = -(( (mhcd.λ .+ V_dl) .- E) .^ 2) ./ (4 * mhcd.λ * kT)
        end
        exp.(exp_arg)
    end
    fd(E) = ox ? 1 .- fermi_dirac(E; kT = kT) : fermi_dirac(E; kT = kT)
    E -> mhcd.A .* mhcd.dos.interp_func.(E .+ V_q) .* fd(E) .* marcus_term(E)
end


text:
# with xr and xs
function fdtdKernel!{T}(n::Int,
		     p_out::Array{T},
		     p0::Array{T,3},
		     p1::Array{T,3},
		     p2::Array{T,3},
		     xr::Vector{Tuple{Int,Int,Int}},
		     xs::Vector{Tuple{Int,Int,Int}},
		     s::Array{T},
		     d1::T,
		     d2::T,
		     d3::T,
		     d4::T,
		     q::Vector{T},
		     G::Array{UInt8,3} )
	for i = 1:size(p0,3), m = 1:size(p0,2), l = 1:size(p0,1)
		if G[l,m,i] == 0x01
			air!(l,m,i,d1,d2,d3,d4,p0,p1,p2)
			#walls
		elseif G[l,m,i] == 0x02
			l1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x03
			le!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x04
			m1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x05
			me!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x06
			i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x07
			ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

			#edges
		elseif G[l,m,i] == 0x08
			l1m1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x09
			l1me!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0a
			l1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0b
			l1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x0c
			lem1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0d
			leme!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0e
			lei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0f
			leie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x10
			m1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x11
			m1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x12
			mei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x13
			meie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

			#corners
		elseif G[l,m,i] == 0x14
			l1m1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x15
			l1m1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x16
			l1mei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x17
			l1meie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x18
			lem1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x19
			lem1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x1a
			lemei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x1b
			lemeie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		end
		## add source
		for k in eachindex(xs)
			if xs[k][1] == l && xs[k][2] == m && xs[k][3] == i
				p2[l,m,i] += s[n,k]
			end
		end
		## get pressure
		for k in eachindex(xr)
			if xr[k][1] == l && xr[k][2] == m && xr[k][3] == i
				p_out[n,k] = p2[l,m,i]
			end
		end
	end
end

# full and xs
function fdtdKernel!{T}(n::Int,
		     p_out::Array{T,4}, #full output
		     p0::Array{T,3},
		     p1::Array{T,3},
		     p2::Array{T,3},
		     xs::Vector{Tuple{Int,Int,Int}},
		     s::Array{T},
		     d1::T,
		     d2::T,
		     d3::T,
		     d4::T,
		     q::Vector{T},
		     G::Array{UInt8,3} )
	for i = 1:size(p0,3), m = 1:size(p0,2), l = 1:size(p0,1)
		if G[l,m,i] == 0x01
			air!(l,m,i,d1,d2,d3,d4,p0,p1,p2)

			#walls
		elseif G[l,m,i] == 0x02
			l1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x03
			le!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x04
			m1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x05
			me!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x06
			i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x07
			ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

			#edges
		elseif G[l,m,i] == 0x08
			l1m1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x09
			l1me!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0a
			l1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0b
			l1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x0c
			lem1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0d
			leme!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0e
			lei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x0f
			leie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x10
			m1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x11
			m1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x12
			mei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x13
			meie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

			#corners
		elseif G[l,m,i] == 0x14
			l1m1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x15
			l1m1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x16
			l1mei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x17
			l1meie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)

		elseif G[l,m,i] == 0x18
			lem1i1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x19
			lem1ie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x1a
			lemei1!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		elseif G[l,m,i] == 0x1b
			lemeie!(l,m,i,d1,d2,d3,d4,p0,p1,p2,q)
		end
		## add source
		for k in eachindex(xs)
			if xs[k][1] == l && xs[k][2] == m && xs[k][3] == i
				p2[l,m,i] += s[n,k]
			end
		end
		p_out[l,m,i,n] = p2[l,m,i]
	end
end

# xr and full s
function fdtdKernel!{T}(n::Int,
		     p_out::Array{T},
		     p0::Array{T,3},
		     p1::Array{T,3},
		     p2::Array{T,3},
		     xr::Vector{Tuple{Int,Int,Int}},
		     s::Array{T,4},
		     d1::T,
		     d2::T,
		     d3::T,
		     d4::T,
		     q::Vector{T},
		     G::Array{UInt8,3} )
	for i = 1:size(p0,3), m = 1:size(p0,2), l = 1:size(p0,1)
		if G[l,m,i] == 0x01
			air!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,s)
			#walls
		elseif G[l,m,i] == 0x02
			l1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x03
			le!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x04
			m1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x05
			me!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x06
			i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x07
			ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

			#edges
		elseif G[l,m,i] == 0x08
			l1m1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x09
			l1me!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0a
			l1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0b
			l1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x0c
			lem1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0d
			leme!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0e
			lei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0f
			leie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x10
			m1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x11
			m1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x12
			mei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x13
			meie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

			#corners
		elseif G[l,m,i] == 0x14
			l1m1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x15
			l1m1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x16
			l1mei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x17
			l1meie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x18
			lem1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x19
			lem1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x1a
			lemei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x1b
			lemeie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		end
		## get pressure
		for k in eachindex(xr)
			if xr[k][1] == l && xr[k][2] == m && xr[k][3] == i
				p_out[n,k] = p2[l,m,i]
			end
		end
	end
end

# full p and full s
function fdtdKernel!{T}(n::Int,
		     p_out::Array{T,4}, #full output
		     p0::Array{T,3},
		     p1::Array{T,3},
		     p2::Array{T,3},
		     s::Array{T,4},
		     d1::T,
		     d2::T,
		     d3::T,
		     d4::T,
		     q::Vector{T},
		     G::Array{UInt8,3} )
	for i = 1:size(p0,3), m = 1:size(p0,2), l = 1:size(p0,1)
		if G[l,m,i] == 0x01
			air!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,s)
			#walls
		elseif G[l,m,i] == 0x02
			l1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x03
			le!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x04
			m1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x05
			me!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x06
			i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x07
			ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

			#edges
		elseif G[l,m,i] == 0x08
			l1m1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x09
			l1me!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0a
			l1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0b
			l1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x0c
			lem1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0d
			leme!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0e
			lei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x0f
			leie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x10
			m1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x11
			m1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x12
			mei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x13
			meie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

			#corners
		elseif G[l,m,i] == 0x14
			l1m1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x15
			l1m1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x16
			l1mei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x17
			l1meie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)

		elseif G[l,m,i] == 0x18
			lem1i1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x19
			lem1ie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x1a
			lemei1!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		elseif G[l,m,i] == 0x1b
			lemeie!(n,l,m,i,d1,d2,d3,d4,p0,p1,p2,q,s)
		end
		p_out[l,m,i,n] = p2[l,m,i]
	end
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
"""Check whether two `GeoArrays`s `a` and `b` are
geographically equal, although not necessarily in content."""
function equals(a::GeoArray, b::GeoArray)
    size(a) == size(b) && a.f == b.f && a.crs == b.crs
end

function Base.:-(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .- b.A, a.f, a.crs)
end

function Base.:+(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .+ b.A, a.f, a.crs)
end

function Base.:*(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .* b.A, a.f, a.crs)
end

function Base.:/(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A ./ b.A, a.f, a.crs)
end


text:
################################################################################################
#
# OpenCV_modules.jl
#
# Load the OpenCV 3.0 modules
# Currently the following are enabled:
#
# core.     The Core Functionality
#            * Basic structures
#            * Operations on Arrays
#            * Clustering
#            * Utility and System Functions and Macros
#
# imgproc.  Image Processing
#            * Image Filtering
#            * Geometric Image Transformations
#            * Miscellaneous Image Transformations
#            * Drawing Functions
#            * ColorMaps
#            * Histograms
#            * Structural Analysis and Shape Descriptors
#            * Motion Analysis and Object Tracking
#            * Feature Detection
#            * Object Detection
#
# imgcodecs. Image file reading and writing
#            * Reading and Writing Images
#
# videoio.   Media I/O
#            * Reading and Writing Video
#
# highgui.   High-level GUI and Media I/O
#            * User Interface
#
# video.     Motion Analysis and Object Tracking
#
#
# TO DO:
# calib3d.    Camera Calibration and 3D Reconstruction
# features2d. 2D Features Framework
# objdetect.  Object Detection
# ml.         Machine Learning
# flann.      Clustering and Search in Multi-Dimensional Spaces
# photo.      Computational Photography
# shape.      Shape Distance and Matching
# superres.   Super Resolution
# videostab.  Video Stabilization
# viz.        3D Visualizer
#             Binary descriptors for lines extracted from an image
# optflow.    Optical Flow Algorithms
# reg.        Image Registration
#             Saliency API
# tracking.   Tracking API
# xfeatures2d.Extra 2D Features Framework
# ximgproc.   Extended Image Processing
# xobjdetect. Extended object detection.
# xphoto.     Additional photo processing algorithms
#
#
#
#
# Sources:
# http://docs.opencv.org/trunk/modules/refman.html
# http://physics.nyu.edu/grierlab/manuals/opencv/namespacecv.html#a346f563897249351a34549137c8532a0
####################################################################################################

# Core
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_core.jl"))
# Custom module for Mat arrays (e.g., get and set methods)
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_Mat.jl"))

# Image processing
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_imgproc.jl"))

# Image file reading and writing
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_imgcodecs.jl"))

# Reading and Writing Video
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_videoio.jl"))

# High-level GUI and Media I/O
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_highgui.jl"))

# Video: Motion Analysis and Object Tracking
# Enable only if the user provides the path for opencv-contrib-master in src/OpenCV.jl lines 101-105
# Unfortunately, some of the functions in opencv-contrib-master are necessary to support the video module in OpenCV 3.0.0

# include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_video.jl"))

# OpenCL support (ocl)
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_ocl.jl"))

# OpenGL support (ogl)
# include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_opengl.jl"))

# Support for conversion and manipulation of images from external packages
include(joinpath(Pkg.dir("OpenCV"), "./src/OpenCV_ImagesSupport.jl"))


text:
function linsolve(operator, b, x₀, alg::CG, a₀::Real = 0, a₁::Real = 1)
    # Initial function operation and division defines number type
    y₀ = apply(operator, x₀)
    T = typeof(dot(b, y₀)/norm(b)*one(a₀)*one(a₁))
    α₀ = convert(T, a₀)
    α₁ = convert(T, a₁)
    # Continue computing r = b - a₀ * x₀ - a₁ * operator(x₀)
    r = copyto!(similar(b, T), b)
    r = α₀ == 0 ? r : axpy!(-α₀, x₀, r)
    r = axpy!(-α₁, y₀, r)
    x = copyto!(similar(r), x₀)
    normr = norm(r)
    S = typeof(normr)

    # Algorithm parameters
    maxiter = alg.maxiter
    tol = convert(S, alg.rtol == 0 ? alg.atol : max(alg.atol, alg.rtol*norm(b)))
    numops = 1 # operator has been applied once to determine r
    numiter = 0

    # Check for early return
    normr < tol && return (x, ConvergenceInfo(1, r, normr, numiter, numops))

    # First iteration
    ρ = normr^2
    p = copyto!(similar(r), r)
    q = apply(operator, p)
    if α₀ != zero(α₀) || α₁ != one(α₁)
        axpby!(α₀, p, α₁, q)
    end
    α = ρ / dot(p,q)
    axpy!(+α, p, x)
    axpy!(-α, q, r)
    normr = norm(r)
    ρold = ρ
    ρ = normr^2
    β = ρ/ρold
    numops += 1
    numiter += 1

    # Check for early return
    normr < tol && return (x, ConvergenceInfo(1, r, normr, numiter, numops))

    while numiter < maxiter
        axpby!(1, r, β, p)
        q = apply(operator, p)
        if α₀ != zero(α₀) || α₁ != one(α₁)
            axpby!(α₀, p, α₁, q)
        end
        α = ρ / dot(p,q)
        axpy!(+α, p, x)
        axpy!(-α, q, r)
        normr = norm(r)
        if normr < tol # recompute to account for buildup of floating point errors
            r = copyto!(similar(b, T), b)
            r = α₀ == 0 ? r : axpy!(-α₀, x, r)
            r = axpy!(-α₁, apply(operator, x), r)
            normr = norm(r)
            ρ = normr^2
            β = 0 # restart CG
        else
            ρold = ρ
            ρ = normr^2
            β = ρ/ρold
        end
        normr < tol && return (x, ConvergenceInfo(1, r, normr, numiter, numops))
        numops += 1
        numiter += 1
    end
    return (x, ConvergenceInfo(0, r, normr, numiter, numops))
end


text:
facts("Scenario generation") do
  context("Internal data structure") do
    @fact _scgen_symbol_from_function(scenarioDuplication) --> :Duplicate
    @fact _scgen_symbol_from_function(scenarioConcatenation) --> :Concatenate
    @fact _scgen_symbol_from_function(scenarioMerging) --> :Merge
    @fact _scgen_symbol_from_function(scenarioMixing) --> :Mix

    @fact ReservoirManagement._scgen_function_from_symbol(:Duplicate) --> exactly(scenarioDuplication)
    @fact ReservoirManagement._scgen_function_from_symbol(:Concatenate) --> exactly(scenarioConcatenation)
    @fact ReservoirManagement._scgen_function_from_symbol(:Merge) --> exactly(scenarioMerging)
    @fact ReservoirManagement._scgen_function_from_symbol(:Mix) --> exactly(scenarioMixing)

    @fact ReservoirManagement._scgen_string_from_symbol(:Duplicate) --> "duplicate"
    @fact ReservoirManagement._scgen_string_from_symbol(:Concatenate) --> "concatenate"
    @fact ReservoirManagement._scgen_string_from_symbol(:Merge) --> "merge"
    @fact ReservoirManagement._scgen_string_from_symbol(:Mix) --> "mix"

    @fact ReservoirManagement._scgen_string_from_function(scenarioDuplication) --> "duplicate"
    @fact ReservoirManagement._scgen_string_from_function(scenarioConcatenation) --> "concatenate"
    @fact ReservoirManagement._scgen_string_from_function(scenarioMerging) --> "merge"
    @fact ReservoirManagement._scgen_string_from_function(scenarioMixing) --> "mix"

    @fact ReservoirManagement._scgen_string_from_any(:Duplicate) --> "duplicate"
    @fact ReservoirManagement._scgen_string_from_any(:Concatenate) --> "concatenate"
    @fact ReservoirManagement._scgen_string_from_any(:Merge) --> "merge"
    @fact ReservoirManagement._scgen_string_from_any(:Mix) --> "mix"

    @fact ReservoirManagement._scgen_string_from_any(scenarioDuplication) --> "duplicate"
    @fact ReservoirManagement._scgen_string_from_any(scenarioConcatenation) --> "concatenate"
    @fact ReservoirManagement._scgen_string_from_any(scenarioMerging) --> "merge"
    @fact ReservoirManagement._scgen_string_from_any(scenarioMixing) --> "mix"
  end

  context("Duplication") do
    TS = gen_tses(52, 2000, 3)

    @fact_throws scenarioDuplication([gen_ts(length, year) for i in 1:6], 2) # Same year for each scenario.

    dTS = scenarioDuplication(TS, 2)
    @fact length(dTS) --> length(TS)
    for i in 1:length(dTS)
      @fact values(TS[i][1:52]) --> values(dTS[i][1:52])
      @fact values(TS[i][1:52]) --> values(dTS[i][53:104])
    end

    dTS = scenarioDuplication(TS, 3)
    @fact length(dTS) --> length(TS)
    for i in 1:length(dTS)
      @fact values(TS[i][1:52]) --> values(dTS[i][1:52])
      @fact values(TS[i][1:52]) --> values(dTS[i][53:104])
      @fact values(TS[i][1:52]) --> values(dTS[i][105:156])
    end
  end

  context("Duplication (river APIs)") do
    TS = gen_river("R", 52, 2000, 3)

    dTS = scenarioDuplication(TS, 2)
    @fact countScenarios(dTS) --> countScenarios(TS)
    @fact countTimeSteps(dTS) --> 2 * countTimeSteps(TS)
    for i in 1:countScenarios(dTS)
      @fact getScenario(TS, i)[1:52] --> getScenario(dTS, i)[1:52]
      @fact getScenario(TS, i)[1:52] --> getScenario(dTS, i)[53:104]
    end

    compare_rivers(scenarioGeneration(TS, :Duplicate, 2), dTS)
    compare_rivers(scenarioGeneration(TS, scenarioDuplication, 2), dTS)

    dTS = scenarioDuplication(TS, 3)
    @fact countScenarios(dTS) --> countScenarios(TS)
    @fact countTimeSteps(dTS) --> 3 * countTimeSteps(TS)
    for i in 1:countScenarios(dTS)
      @fact getScenario(TS, i)[1:52] --> getScenario(dTS, i)[1:52]
      @fact getScenario(TS, i)[1:52] --> getScenario(dTS, i)[53:104]
      @fact getScenario(TS, i)[1:52] --> getScenario(dTS, i)[105:156]
    end

    compare_rivers(scenarioGeneration(TS, :Duplicate, 3), dTS)
    compare_rivers(scenarioGeneration(TS, scenarioDuplication, 3), dTS)
  end

  context("Concatenation") do
    TS = gen_tses(52, 2000, 6)

    # The same year cannot be present twice in the same generated scenario. Hence create from scratch a series of
    # time series with the same year, and try to generate scenario from them: this input data is inconsistent.
    @fact_throws scenarioConcatenation([gen_ts(length, year) for i in 1:6], 2)

    @fact_throws scenarioConcatenation(TS, 5) # Not an integer number of input scenarios per output scenario.

    cTS = scenarioConcatenation(TS, 2)
    @fact length(cTS) --> 6 / 2
    @fact values(cTS[1][1:52]) --> values(TS[1][1:52])
    @fact values(cTS[1][53:104]) --> values(TS[2][1:52])
    @fact values(cTS[2][1:52]) --> values(TS[3][1:52])
    @fact values(cTS[2][53:104]) --> values(TS[4][1:52])
    @fact values(cTS[3][1:52]) --> values(TS[5][1:52])
    @fact values(cTS[3][53:104]) --> values(TS[6][1:52])
  end

  context("Concatenation (river APIs)") do
    TS = gen_river("R", 52, 2000, 6)

    dTS = scenarioConcatenation(TS, 2)
    @fact countScenarios(dTS) --> countScenarios(TS) / 2
    @fact countTimeSteps(dTS) --> 2 * countTimeSteps(TS)
    @fact getScenario(dTS, 1)[1:52] --> getScenario(TS, 1)[1:52]
    @fact getScenario(dTS, 1)[53:104] --> getScenario(TS, 2)[1:52]
    @fact getScenario(dTS, 2)[1:52] --> getScenario(TS, 3)[1:52]
    @fact getScenario(dTS, 2)[53:104] --> getScenario(TS, 4)[1:52]
    @fact getScenario(dTS, 3)[1:52] --> getScenario(TS, 5)[1:52]
    @fact getScenario(dTS, 3)[53:104] --> getScenario(TS, 6)[1:52]

    compare_rivers(scenarioGeneration(TS, :Concatenate, 2), dTS)
    compare_rivers(scenarioGeneration(TS, scenarioConcatenation, 2), dTS)
  end

  context("Merging") do
    TS = gen_tses(52, 2000, 4)

    # The same year cannot be present twice in the same generated scenario. Hence create from scratch a series of
    # time series with the same year, and try to generate scenario from them: this input data is inconsistent.
    @fact_throws scenarioMerging([gen_ts(length, year) for i in 1:6], 2)

    cTS = scenarioMerging(TS, 2)
    @fact length(cTS) --> 3
    @fact values(cTS[1][1:52]) --> values(TS[1][1:52])
    @fact values(cTS[1][53:104]) --> values(TS[2][1:52])
    @fact values(cTS[2][1:52]) --> values(TS[2][1:52])
    @fact values(cTS[2][53:104]) --> values(TS[3][1:52])
    @fact values(cTS[3][1:52]) --> values(TS[3][1:52])
    @fact values(cTS[3][53:104]) --> values(TS[4][1:52])
  end

  context("Merging (river APIs)") do
    TS = gen_river("R", 52, 2000, 4)

    dTS = scenarioMerging(TS, 2)
    @fact countScenarios(dTS) --> countScenarios(TS) - 1
    @fact countTimeSteps(dTS) --> 2 * countTimeSteps(TS)
    @fact getScenario(dTS, 1)[1:52] --> getScenario(TS, 1)[1:52]
    @fact getScenario(dTS, 1)[53:104] --> getScenario(TS, 2)[1:52]
    @fact getScenario(dTS, 2)[1:52] --> getScenario(TS, 2)[1:52]
    @fact getScenario(dTS, 2)[53:104] --> getScenario(TS, 3)[1:52]
    @fact getScenario(dTS, 3)[1:52] --> getScenario(TS, 3)[1:52]
    @fact getScenario(dTS, 3)[53:104] --> getScenario(TS, 4)[1:52]

    compare_rivers(scenarioGeneration(TS, :Merge, 2), dTS)
    compare_rivers(scenarioGeneration(TS, scenarioMerging, 2), dTS)
  end

  context("Mixing") do
    TS = gen_tses(52, 2000, 3)

    # The same year cannot be present twice in the same generated scenario. Hence create from scratch a series of
    # time series with the same year, and try to generate scenario from them: this input data is inconsistent.
    @fact_throws scenarioMixing([gen_ts(length, year) for i in 1:6], 2)

    cTS = scenarioMixing(TS, 2)
    @fact length(cTS) --> 9 # All pairs among three objects (can take twice the same object!).

    possible_values = Array{Float64, 1}[values(TS[i][1:52]) for i in 1:3]
    for i in 1:9
      @fact in(values(cTS[i][1:52]), possible_values) --> true
      @fact in(values(cTS[i][53:104]), possible_values) --> true
    end
  end

  context("Mixing (river APIs)") do
    TS = gen_river("R", 52, 2000, 3)

    dTS = scenarioMixing(TS, 2)
    @fact countScenarios(dTS) --> 9
    @fact countTimeSteps(dTS) --> 2 * countTimeSteps(TS)
    @fact getScenario(dTS, 1)[1:52] --> getScenario(TS, 1)[1:52]
    possible_values = Array{Float64, 1}[getScenario(TS, i)[1:52] for i in 1:3]
    for i in 1:9
      @fact in(getScenario(dTS, i)[1:52], possible_values) --> true
      @fact in(getScenario(dTS, i)[53:104], possible_values) --> true
    end

    # The different function calls must give the exact same results, as they are synonymous.
    compare_rivers(scenarioGeneration(TS, :Mix, 2), dTS)
    compare_rivers(scenarioGeneration(TS, scenarioMixing, 2), dTS)
  end
end


text:
@testset "Indices" begin

# Prcp1
d = collect(DateTime(2003,1,1):Day(1):DateTime(2005,12,31))

data = Array{Float64,3}(undef, 2, 2,1096)
data[1,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[1,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
Results = Array{Int64, 3}(undef, 2, 2, 3);
Results[1,1,1] = 365;
Results[1,1, 2] = 0;
Results[1,1,3] = 365;
Results[1,2,1] = 365;
Results[2,1,1] = 365;
Results[2,1,2] = 0;
Results[2,1,3] = 365;
Results[1,2,1] = 365;
Results[1,2,2] = 0;
Results[1,2,3] = 365;
Results[2,2,1] = 365;
Results[2,2,2] = 0;
Results[2,2,3] = 365;

# # Add units
# data = [data][1]u"mm"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "pr")
ind = prcp1(C)
@test ind.data.data == Results

# annualsum
data = Array{Float64,3}(undef, 2, 2, 1096)
data[1,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[1,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
Results = Array{Int64, 3}(undef, 2, 2, 3);
Results[1,1,1] = 365.;
Results[1,1,2] = 0.;
Results[1,1,3] = 365.;
Results[2,1,1] = 365.;
Results[2,1,2] = 0;
Results[2,1,3] = 365.;
Results[1,2,1] = 365.;
Results[1,2,2] = 0.;
Results[1,2,3] = 365.;
Results[2,2,1] = 365.;
Results[2,2,2] = 0.;
Results[2,2,3] = 365.;

# Add units
# data = [data][1]u"mm"
# Results = [Results][1]u"mm"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "pr")
ind = annualsum(C)
@test ind.data.data == Results

# annualmean
data = Array{Float64,3}(undef, 2, 2, 1096)
data[1,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,1,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[1,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
data[2,2,:] = vcat(ones(365,1), zeros(366,1), ones(365))
Results = Array{Int64, 3}(undef, 2, 2, 3);
Results[1,1,1] = 1.;
Results[1,1,2] = 0.;
Results[1,1,3] = 1.;
Results[2,1,1] = 1.;
Results[2,1,2] = 0;
Results[2,1,3] = 1.;
Results[1,2,1] = 1.;
Results[1,2,2] = 0.;
Results[1,2,3] = 1.;
Results[2,2,1] = 1.;
Results[2,2,2] = 0.;
Results[2,2,3] = 1.;

# Add units
# data = [data][1]u"mm"
# Results = [Results][1]u"mm"

# @test annualmean(data, d) == Results

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "pr")
ind = annualmean(C)
@test ind.data.data == Results

# Frostdays

data = Array{Float64,3}(undef, 2, 2, 1096)
data[1,1,:] = vcat(-ones(365,1), -ones(366,1), zeros(365))
data[2,1,:] = vcat(-ones(365,1), -ones(366,1), zeros(365))
data[1,2,:] = vcat(-ones(365,1), -ones(366,1), zeros(365))
data[2,2,:] = vcat(-ones(365,1), -ones(366,1), zeros(365))
Results = Array{Int64, 3}(undef, 2, 2, 3);
Results[1,1,1] = 365;
Results[1,1,2] = 366;
Results[1,1,3] = 0;
Results[2,1,1] = 365;
Results[2,1,2] = 366;
Results[2,1,3] = 0;
Results[1,2,1] = 365;
Results[1,2,2] = 366;
Results[1,2,3] = 0;
Results[2,2,1] = 365;
Results[2,2,2] = 366;
Results[2,2,3] = 0;

# Add units
# data = [data][1]u"°C"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = frostdays(C)
@test ind.data.data == Results
C = ClimateTools.ClimGrid(axisdata, variable = "tasmax")
ind = icingdays(C)
@test ind.data.data == Results

# Summer Days
d = collect(DateTime(2003,1,1):Day(1):DateTime(2007,12,31))

data = Array{Float64, 3}(undef, 2, 2, 1826)
data[1,1,:] = collect(1.0:1826.0); data[1,2,:] = collect(1.0:1826.0); data[2,1,:]=collect(1.0:1826.0); data[2,2,:] = collect(1.0:1826.0);
Results = Array{Int64, 3}(undef, 2, 2, 5);
Results[1, 1, :]=[340, 366, 365, 365, 365]'';
Results[1, 2, :] = [340, 366, 365, 365, 365]'';
Results[2, 1, :] = [340, 366, 365, 365, 365]'';
Results[2, 2, :] = [340, 366, 365, 365, 365]'';

# Add units
# data = [data][1]u"°C"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmax")
ind = summerdays(C)
@test ind.data.data == Results

data = Array{Float64, 3}(undef, 2, 2, 1826)
data[1,1,:] = collect(1.0:1826.0); data[1,2,:]=collect(1.0:1826.0); data[2,1,:]=collect(1.0:1826.0); data[2,2,:]=collect(1.0:1826.0);
Results = Array{Int64, 3}(undef, 2,2,5);
Results[1,1,:]=[345, 366, 365, 365, 365]'';
Results[1,2,:]=[345, 366, 365, 365, 365]'';
Results[2,1,:]=[345, 366, 365, 365, 365]'';
Results[2,2,:]=[345, 366, 365, 365, 365]'';

# Add units
# data = [data][1]u"°C"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = tropicalnights(C)
@test ind.data.data == Results

# Custom thresholds

data = Array{Float64, 3}(undef, 2,2,1826)
data[1,1,:] = collect(1.0:1826.0); data[1,2,:]=collect(1.0:1826.0); data[2,1,:]=collect(1.0:1826.0); data[2,2,:] = collect(1.0:1826.0);
Results = Array{Int64, 3}(undef, 2,2,5);
Results[1,1,:] = [345, 366, 365, 365, 365]'';
Results[1,2,:] = [345, 366, 365, 365, 365]'';
Results[2,1,:] = [345, 366, 365, 365, 365]'';
Results[2,2,:] = [345, 366, 365, 365, 365]'';

# Add units
# data = [data][1]u"°C"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = customthresover(C, 20)
@test ind.data.data == Results

data = Array{Float64, 3}(undef, 2,2,1826)
data[1,1,:] = collect(-800.0:1025.0); data[1,2,:] = collect(-800.0:1025.0); data[2,1,:]=collect(-800.0:1025.0); data[2,2,:] = collect(-800.0:1025.0);
Results = Array{Int64, 3}(undef, 2,2,5);
Results[1,1,:]=[365, 366, 269, 0, 0]'';
Results[1,2,:]=[365, 366, 269, 0, 0]'';
Results[2,1,:]=[365, 366, 269, 0, 0]'';
Results[2,2,:]=[365, 366, 269, 0, 0]'';

# # Add units
# data = [data][1]u"°C"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = customthresunder(C, 200)
@test ind.data.data == Results

# ANNUAL MAXIMUM
data = Array{Float64, 3}(undef, 2,2,1826)
data[1,1,:] = collect(-800.0:1025.0); data[1,2,:] = collect(-800.0:1025.0); data[2,1,:] = collect(-800.0:1025.0); data[2,2,:] = collect(-800.0:1025.0);
Results = Array{Float64, 3}(undef, 2,2,5);
Results[1,1,:]=[-436.0, -70.0, 295.0, 660.0, 1025.0]'';
Results[1,2,:]=[-436.0, -70.0, 295.0, 660.0, 1025.0]'';
Results[2,1,:]=[-436.0, -70.0, 295.0, 660.0, 1025.0]'';
Results[2,2,:]=[-436.0, -70.0, 295.0, 660.0, 1025.0]'';

# # Add units
# data = [data][1]u"K"
# Results = [Results][1]u"K"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = annualmax(C)
@test ind.data.data == Results

# ANNUAL MINIMUM
data = Array{Float64, 3}(undef, 2,2,1826)
data[1,1,:] = collect(-800.:1025.); data[1,2,:]=collect(-800.:1025.); data[2,1,:]=collect(-800.:1025.); data[2,2,:]=collect(-800.:1025.);
Results = Array{Int64, 3}(undef, 2,2,5);
Results[1,1,:]=[-800.0,-435.0,-69.0,296.0,661.0]'';
Results[1,2,:]=[-800.0,-435.0,-69.0,296.0,661.0]'';
Results[2,1,:]=[-800.0,-435.0,-69.0,296.0,661.0]'';
Results[2,2,:]=[-800.0,-435.0,-69.0,296.0,661.0]'';
 #
 # # Add units
 # data = [data][1]u"K"
 # Results = [Results][1]u"K"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, variable = "tasmin")
ind = annualmin(C)
@test ind.data.data == Results

# Period mean
data = Array{Float64, 3}(undef, 2, 2, 1826)
data[1,1,:] = collect(1.0:1826.0); data[1,2,:] = collect(-10.0:1815.0); data[2,1,:] = collect(0.0:1825.0); data[2,2,:] = collect(-1725.0:100);
Results = Array{Float64, 2}(undef, 2, 2);
Results[1, 1]= 746.0;
Results[1, 2] = 735.0;
Results[2, 1] = 745.0;
Results[2, 2] = -980.0;

# # Add units
# data = [data][1]u"K"
# Results = [Results][1]u"K"

# ClimGrid based tests
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata)
ind = periodmean(C, start_date=(2004, 06, 01), end_date=(2005, 08, 31)) #mean between June 1st 2004 and August 31st 2005
@test ind.data.data == Results

Results = Array{Float64, 2}(undef, 2, 2);
Results[1, 1]= 913.5;
Results[1, 2] = 902.5;
Results[2, 1] = 912.5;
Results[2, 2] = -812.5;

# Add units
# Results = [Results][1]u"K"

ind = periodmean(C)
@test ind.data.data == Results

# Daymean
raw = collect(0:0.25:1459.75) # 6hr data
d = NCDatasets.timedecode(raw, "days since 2000-01-01", "noleap")

data = Array{Float64}(undef, 2, 2, 5840)
data[1, 1, :] = collect(1:5840)
data[1, 2, :] = collect(1:5840)
data[2, 1, :] = collect(1:5840)
data[2, 2, :] = collect(1:5840)

axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimateTools.ClimGrid(axisdata, frequency="6h")
D = daymean(C)
@test D[1][1,1,1] == 2.5
@test D[1][1,1,end] == 5838.5

D = daysum(C)
@test D[1][1,1,1] == 10.0
@test D[1][1,1,end] == 23354.0

@test get_timevec(C)[1] == DateTimeNoLeap(2000, 01, 01)

# Test vaporpressure(specific_humidity::ClimGrid, surface_pressure::ClimGrid)
d = DateTime(2003,1,1):Day(1):DateTime(2003,1,3)
# Dummy data
data_huss = Array{Float64,3}(undef, 2,2,3)
data_huss[1,1,:] .= 0
data_huss[1,2,:] .= 0.1
data_huss[2,1,:] .= 0.5
data_huss[2,2,:] .= 1
data_ps = Array{Float64,3}(undef, 2,2,3)
data_ps[:,:,1] .= 10.0
data_ps[:,:,2] .= 100.0
data_ps[:,:,3] .= 1000.0
# Expected values
Results = Array{Float64,3}(undef, 2,2,3)
Results[1,1,1] = 0.0;
Results[1,1,2] = 0.0;
Results[1,1,3] = 0.0;
Results[2,1,1] = data_huss[2,1,1]*data_ps[2,1,1]/(0.622+data_huss[2,1,1]);
Results[2,1,2] = data_huss[2,1,2]*data_ps[2,1,2]/(0.622+data_huss[2,1,2]);
Results[2,1,3] = data_huss[2,1,3]*data_ps[2,1,3]/(0.622+data_huss[2,1,3]);
Results[1,2,1] = data_huss[1,2,1]*data_ps[1,2,1]/(0.622+data_huss[1,2,1]);
Results[1,2,2] = data_huss[1,2,2]*data_ps[1,2,2]/(0.622+data_huss[1,2,2]);
Results[1,2,3] = data_huss[1,2,3]*data_ps[1,2,3]/(0.622+data_huss[1,2,3]);
Results[2,2,1] = data_huss[2,2,1]*data_ps[2,2,1]/(0.622+data_huss[2,2,1]);
Results[2,2,2] = data_huss[2,2,2]*data_ps[2,2,2]/(0.622+data_huss[2,2,2]);
Results[2,2,3] = data_huss[2,2,3]*data_ps[2,2,3]/(0.622+data_huss[2,2,3]);
# Using the function on dummy data
axisdata_huss = AxisArray(data_huss, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_ps = AxisArray(data_ps, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_huss = ClimateTools.ClimGrid(axisdata_huss, variable = "huss")
C_ps = ClimateTools.ClimGrid(axisdata_ps, dataunits="Pa", variable = "ps")
vp = vaporpressure(C_huss, C_ps)
# Run the test
@test vp.data.data == Results

# Test  vaporpressure(specific_humidity::ClimGrid, sealevel_pressure::ClimGrid, orography::ClimGrid, daily_temperature::ClimGrid)
# Dummy data
data_huss = Array{Float64,3}(undef, 2,2,3)
data_huss[1,1,:] .= 0
data_huss[1,2,:] .= 0.1
data_huss[2,1,:] .= 0.5
data_huss[2,2,:] .= 1

data_psl = Array{Float64,3}(undef, 2,2,3)
data_psl[:,:,1] .= 10.0
data_psl[:,:,2] .= 100.0
data_psl[:,:,3] .= 1000.0
# data_psl = [data_psl][1]u"Pa"

data_orog = Array{Float64,2}(undef, 2,2)
data_orog[:,1] .= 0.0
data_orog[:,2] .= 10.0
# data_orog = [data_orog][1]u"m"

data_tas = Array{Float64,3}(undef, 2,2,3)
data_tas[1,:,:] .= 100.0
data_tas[2,:,:] .= 300.0
# data_tas = [data_tas][1]u"K"

# Expected results
Results = Array{Float64,3}(undef, 2,2,3)
Results[1,1,:] .= round(0.0, digits=10)
Results[1,2,1] = round(1.380315267090330, digits=10)
Results[1,2,2] = round(13.80315267090330, digits=10)
Results[1,2,3] = round(138.0315267090330, digits=10)
Results[2,1,1] = round(4.456327985739750, digits=10)
Results[2,1,2] = round(44.56327985739750, digits=10)
Results[2,1,3] = round(445.6327985739750, digits=10)
Results[2,2,1] = round(6.158207427095860, digits=10)
Results[2,2,2] = round(61.58207427095860, digits=10)
Results[2,2,3] = round(615.8207427095860, digits=10)
# Results = [Results][1]u"Pa"

# Create the Climgrids
axisdata_huss = AxisArray(data_huss, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_psl = AxisArray(data_psl, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_tas = AxisArray(data_tas, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_psl = AxisArray(data_psl, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_orog = AxisArray(data_orog, Axis{:lon}(1:2), Axis{:lat}(1:2))
C_huss = ClimateTools.ClimGrid(axisdata_huss, variable = "huss")
C_psl = ClimateTools.ClimGrid(axisdata_psl, dataunits = "Pa", variable = "psl")
C_tas = ClimateTools.ClimGrid(axisdata_tas, dataunits = "K", variable = "tas")
C_orog = ClimateTools.ClimGrid(axisdata_orog, dataunits = "m", variable = "orog")
# Using the function on dummy data
vp = vaporpressure(C_huss, C_psl, C_orog, C_tas)
# Run the test
@test round.(vp.data.data, digits=10) == Results

# Test for tas in Celsius
data_tas = Array{Float64,3}(undef, 2,2,3)
data_tas[1,:,:] .= 100-273.15
data_tas[2,:,:] .= 300-273.15
# data_tas = [data_tas][1]u"°C"

axisdata_tas = AxisArray(data_tas, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tas = ClimateTools.ClimGrid(axisdata_tas, dataunits = "Celsius", variable = "tas")
# Using the function on dummy data
vp = vaporpressure(C_huss, C_psl, C_orog, C_tas)
# Run the test
@test round.(vp.data.data, digits=10) == Results

# Test wbgt(diurnal_temperature::ClimGrid, vapor_pressure::ClimGrid)
d = DateTime(2003,1,1):Day(1):DateTime(2003,1,3)
# Dummy data
data_tdiu = Array{Float64,3}(undef, 2,2,3)
data_tdiu[1,1,:] .= 250.0
data_tdiu[1,2,:] .= 275.0
data_tdiu[2,1,:] .= 300.0
data_tdiu[2,2,:] .= 325.0
# data_tdiu = [data_tdiu][1]u"K"

data_vp = Array{Float64,3}(undef, 2,2,3)
data_vp[:,:,1] .= 10.0
data_vp[:,:,2] .= 100.0
data_vp[:,:,3] .= 1000.0
# data_vp = [data_vp][1]u"Pa"

# Expected results
Results = Array{Float64}(undef, 2, 2, 3)
Results[1,1,1] = -9.1467500
Results[1,2,1] = 5.02825
Results[2,1,1] = 19.20325
Results[2,2,1] = 33.37825
Results[1,1,2] = -8.7930500
Results[1,2,2] = 5.38195
Results[2,1,2] = 19.55695
Results[2,2,2] = 33.73195
Results[1,1,3] = -5.256050
Results[1,2,3] = 8.91895
Results[2,1,3] = 23.09395
Results[2,2,3] = 37.26895
# Creating climgrids
axisdata_tdiu = AxisArray(data_tdiu, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_vp = AxisArray(data_vp, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tdiu = ClimateTools.ClimGrid(axisdata_tdiu, dataunits= "K", variable = "tmean")
C_vp = ClimateTools.ClimGrid(axisdata_vp, dataunits = "Pa", variable = "vp")
# Use the function
C_wbgt = wbgt(C_tdiu, C_vp)
# Run the test
@test round.(C_wbgt.data.data, digits=10) == Results

# Test for tdiu in Celsius
data_tdiu = Array{Float64,3}(undef, 2,2,3)
data_tdiu[1,1,:] .= 250.0 - 273.15
data_tdiu[1,2,:] .= 275.0 - 273.15
data_tdiu[2,1,:] .= 300.0 - 273.15
data_tdiu[2,2,:] .= 325.0 - 273.15
# data_tdiu = [data_tdiu][1]u"°C"

axisdata_tdiu = AxisArray(data_tdiu, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tdiu = ClimateTools.ClimGrid(axisdata_tdiu, dataunits= "Celsius", variable = "tmean")
# Use the function
C_wbgt = wbgt(C_tdiu, C_vp)
# Run the test
@test round.(C_wbgt.data.data, digits=10) == Results

# Test diurnaltemperature()
d = DateTime(2003,1,1):Day(1):DateTime(2003,1,3)
# Dummy data
data_tmax = Array{Float64,3}(undef, 2,2,3)
data_tmax[1,1,:] .= 0.0
data_tmax[1,2,:] .= 10.0
data_tmax[2,1,:] .= 20.0
data_tmax[2,2,:] .= 30.0
data_tmin = Array{Float64,3}(undef, 2,2,3)
data_tmin[:,:,1] .= -10.0
data_tmin[:,:,2] .= 0.0
data_tmin[:,:,3] .= 10.0
# data_tmin = [data_tmin][1]u"°C"

α = 0.4
# Expected results
Results[1,1,1] = -4.0
Results[1,2,1] = 2.0
Results[2,1,1] = 8.0
Results[2,2,1] = 14.0
Results[1,1,2] = 0.0
Results[1,2,2] = 6.0
Results[2,1,2] = 12.0
Results[2,2,2] = 18.0
Results[1,1,3] = 4.0
Results[1,2,3] = 10.0
Results[2,1,3] = 16.0
Results[2,2,3] = 22.0
# Results = [Results][1]u"°C"

# Creating climgrids
axisdata_tmax = AxisArray(data_tmax, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_tmin = AxisArray(data_tmin, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tmax = ClimateTools.ClimGrid(axisdata_tmax, dataunits= "Celsius", variable = "tasmax")
C_tmin = ClimateTools.ClimGrid(axisdata_tmin, dataunits = "Celsius", variable = "tasmin")
# Using the function
C_tdiu = diurnaltemperature(C_tmin, C_tmax, α)
# Run the test
@test C_tdiu.data.data == Results

# Test meantemperature()
d = DateTime(2003,1,1):Day(1):DateTime(2003,1,3)
# Dummy data
data_tmax = Array{Float64,3}(undef, 2,2,3)
data_tmax[1,1,:] .= 0.0
data_tmax[1,2,:] .= 10.0
data_tmax[2,1,:] .= 20.0
data_tmax[2,2,:] .= 30.0
# data_tmax = [data_tmax][1]u"°C"

data_tmin = Array{Float64,3}(undef, 2,2,3)
data_tmin[:,:,1] .= -10.0
data_tmin[:,:,2] .= 0.0
data_tmin[:,:,3] .= 10.0
# data_tmin = [data_tmin][1]u"°C"

# Expected results
Results = Array{Float64,3}(undef, 2,2,3)
Results[1,1,1] = -5.0
Results[1,2,1] = 0.0
Results[2,1,1] = 5.0
Results[2,2,1] = 10.0
Results[1,1,2] = 0.0
Results[1,2,2] = 5.0
Results[2,1,2] = 10.0
Results[2,2,2] = 15.0
Results[1,1,3] = 5.0
Results[1,2,3] = 10.0
Results[2,1,3] = 15.0
Results[2,2,3] = 20.0
# Results = [Results][1]u"°C"

# Creating climgrids
axisdata_tmax = AxisArray(data_tmax, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
axisdata_tmin = AxisArray(data_tmin, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tmax = ClimateTools.ClimGrid(axisdata_tmax, dataunits= "Celsius", variable = "tasmax", frequency="day")
C_tmin = ClimateTools.ClimGrid(axisdata_tmin, dataunits = "Celsius", variable = "tasmin")
# Using the function
C_mean = meantemperature(C_tmin, C_tmax)
# Run the test
@test C_mean.data.data == Results

C_tmax2 = C_tmax
ens = [C_tmax2, C_tmax, C_tmax, C_tmax, C_tmax]
E = ensemble_mean(ens)
@test E[1][1, 1, 1] == C_tmax[1][1, 1, 1]

E = ensemble_std(ens)
@test E[1][1, 1, 1] == std([C_tmax[1][1,1,1], C_tmax[1][1,1,1], C_tmax[1][1,1,1], C_tmax[1][1, 1, 1], C_tmax[1][1,1,1]])

d = DateTime(2003,1,1):Day(1):DateTime(2003,1,3)
# Dummy data
data_tmax = Array{Float64,3}(undef, 2,2,3)
data_tmax[1,1,:] .= 0.0
data_tmax[1,2,:] .= 10.0
data_tmax[2,1,:] .= 20.0
data_tmax[2,2,:] .= 30.0
# data_tmax = [data_tmax][1]u"K"

axisdata_tmax = AxisArray(data_tmax, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C_tmax = ClimateTools.ClimGrid(axisdata_tmax, dataunits= "Kelvin", variable = "tasmax", frequency="day")

ens = [C_tmax2 + 2.0, C_tmax, C_tmax, C_tmax, C_tmax]
E = ensemble_max(ens)
@test E[1][1,1,1] == ens[1][1][1,1,1]

E = ensemble_min(ens)
@test E[1][1,1,1] == ens[2][1][1,1,1]

# Drought code
# core function
temp = 17.0
precip = 0.0
month = 4
dc0 = 15.0
@test drought_dc(precip, temp, dc0, month) == 19.014

dc0 = 19.014
temp = 20.0
precip = 2.4
@test round(drought_dc(precip, temp, dc0, month), digits=3) == 23.568

dc0 = 23.568
temp = 8.5
precip = 0.0
@test round(drought_dc(precip, temp, dc0, month), digits=3) == 26.052

dc0 = 40.6
temp = 7.0
precip = 9.0
@test round(drought_dc(precip, temp, dc0, month), digits=3) == 29.529

# Drought code
# ClimGrid wrapper
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
datevec = collect(DateTime(2000, 4, 13):Day(1):DateTime(2000, 4, 27))
temp = Array{Float32}(undef, 1, 1, 15)
temp[1,1,:] = [17.0, 20.0, 8.5, 6.5, 13.0, 6.0, 5.5, 8.5, 9.5, 7.0, 6.5, 6.0, 13.0, 15.5, 23.0]
precip = Array{Float32}(undef, 1, 1, 15)
precip[1,1,:] = [0.0, 2.4, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 9.0, 1.0, 0.0, 0.0, 0.0, 0.0]

tempax = AxisArray(temp, Axis{:lon}([-88.0]), Axis{:lat}([45.0]), Axis{:time}(datevec))
precax = AxisArray(precip, Axis{:lon}([-88.0]), Axis{:lat}([45.0]), Axis{:time}(datevec))

tas = ClimGrid(tempax, variable="tas", dimension_dict=dimension_dict, dataunits="Celsius")
pr = ClimGrid(precax, variable="pr", dimension_dict=dimension_dict, dataunits="mm")

dc = drought_dc(pr, tas)

res = round.([19.014, 23.567999999999998, 26.052, 28.176, 31.47, 33.504, 35.448,  37.932, 40.596000000000004, 29.52469510340648, 31.648695103406478, 33.68269510340648, NaN, NaN, NaN], digits=3)

@test sum(round.(dc[1][1,1,1:end-3], digits=3) .== res[1:end-3]) == length(res[1:end-3])

@test sum(isnan.(dc[1][1,1,:])) == 3

end


text:
function assemble_system!(problem::FDHeatProblem{TBCL, TBCR, TS}) where {TBCL<:PeriodicBC, TBCR<:PeriodicBC, TS<:BEuler}

    # build first row
    problem.A[1,1] = 1 + 2*problem.Δt/problem.Δx^2;
    problem.A[1,2] = -problem.Δt/problem.Δx^2;
    problem.A[1,end] = -problem.Δt/problem.Δx^2;
    # build rows 2-n
    for i in 2:problem.nx
        problem.A[i,i-1] = -problem.Δt/problem.Δx^2;
        problem.A[i,i] = 1 + 2*problem.Δt/problem.Δx^2;
        problem.A[i,i+1] = -problem.Δt/problem.Δx^2;
    end
    # build row nx+1
    problem.A[end,1] =  -problem.Δt/problem.Δx^2;
    problem.A[end,end-1] =  -problem.Δt/problem.Δx^2;
    problem.A[end,end] =  1 + 2*problem.Δt/problem.Δx^2;

    problem
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

using CSV, DataFrames, Luxor, Colors, Cairo

include(joinpath(dirname(dirname(@__FILE__)), "src", "messages.jl"))
include(joinpath(dirname(dirname(@__FILE__)), "src", "distance.jl"))
include(joinpath(dirname(dirname(@__FILE__)), "src", "utils.jl"))
include(joinpath(dirname(dirname(@__FILE__)), "src", "primes.jl"))

function integer_search(stop_at::Integer)::Array{Array{Number, 1}}
    A = []
    upper_bound = 10
    increment_bound = 10
	processed = []
	i = 0

    while true
        for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound
			# skip configurations that have already been processed
			# arelessthan(upper_bound - increment_bound, q, n, d) && continue
			(q, n, d) ∈ processed && continue
            hb = hamming_bound(q, n, d, no_round)
			push!(processed, (q, n, d))
			# i += 1; println("$i:\t$q, $n, $d")
            push!(A, [q, n, d, hb])
            isequal(length(A), stop_at) && return A
        end
        upper_bound += increment_bound
    end
end

function make_integer_df(stop_at::Integer)
    A = integer_search(stop_at)
    D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[])

    for i in A
        push!(D, i)
    end

	return D
end

function draw_turtle(stop_at::Integer; specific::Bool=false)
	rigidness = "non_rigid"
	if specific
		rigidness = "rigid"
	end

	df = make_integer_df(stop_at)

	golden_angle = 137.5
	angle = golden_angle
	# angle = 42

	# Drawing(600, 400, joinpath(dirname(dirname(@__FILE__)), "other", "turtle_drawing.pdf"))
	Drawing(10000, 10000, joinpath(homedir(), "projects", "CodingTheory.jl", "other", "$(rigidness)", "turtle_drawing_$(stop_at)_$(angle).pdf"))
	origin()
	# background("midnightblue")

	🐢 = Turtle()
	Pencolor(🐢, "black")
	Penwidth(🐢, 1.5)
	n = 5

	# distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances
	for row in eachrow(df)
		if specific
			if row.d < row.n && ! isone(row.d) && ! allequal(row.q, row.n, row.d) && ! iszero(mod(row.d, 2)) && ! isone(row.hamming_bound) && ! ishammingbound(BigInt(row.q), BigInt(row.n), BigInt(row.d))
				continue
			end
		else
			if row.d < row.n && ! isone(row.d)
				continue
			end
		end

		if isinteger(row.hamming_bound)
			Turn(🐢, angle)
			# HueShift(🐢)
		end
		Forward(🐢, n)
		# n += 0.75
	end

	# fontsize(20)
	# Message(🐢, "finished")
	Luxor.finish()

	return nothing
end


text:
export SemialgebraicSetsHCSolver

function ModelKit.System(V::SemialgebraicSets.AbstractAlgebraicSet; kwargs...)
    ModelKit.System(SemialgebraicSets.equalities(V); kwargs...)
end

"""
    SemialgebraicSetsHCSolver(; compile = false, options...)

Construct a `SemialgebraicSets.AbstractAlgebraicSolver` to be used in `SemialgebraicSets`.
`options` are all valid options for [`solve`](@ref).

## Example

```
julia> using HomotopyContinuation, SemialgebraicSets;

julia> solver = SemialgebraicSetsHCSolver(; show_progress = false)
SemialgebraicSetsHCSolver(; compile = :none, show_progress = false)

julia> @polyvar x y
(x, y)

julia> V = @set x^2 == 1 && y^2 == 2 solver
Algebraic Set defined by 2 equalities
 x^2 - 1.0 = 0
 y^2 - 2.0 = 0

julia> collect(V)
4-element Array{Array{Float64,1},1}:
 [1.0, 1.414213562373095]
 [1.0, -1.414213562373095]
 [-1.0, 1.414213562373095]
 [-1.0, -1.414213562373095]
```
"""
struct SemialgebraicSetsHCSolver <: SemialgebraicSets.AbstractAlgebraicSolver
    options::Any
end
SemialgebraicSetsHCSolver(; compile = :none, options...) =
    SemialgebraicSetsHCSolver((compile = compile, options...))

function Base.show(io::IO, solver::SemialgebraicSetsHCSolver)
    print(io, "SemialgebraicSetsHCSolver(; ")
    join(
        io,
        [
            "$(k) = $(v isa Symbol ? Expr(:quote, v) : v)" for
            (k, v) in pairs(solver.options)
        ],
        ", ",
    )
    print(io, ")")
end

function SemialgebraicSets.solvealgebraicequations(
    V::SemialgebraicSets.AbstractAlgebraicSet,
    solver::SemialgebraicSetsHCSolver,
)
    F = System(V)
    m, n = size(F)
    # Check that we can have isolated solutions
    m ≥ n || return nothing
    # Only return real, non-singular solutions
    return real_solutions(solve(F; solver.options...); only_nonsingular = true)
end


text:
#####
##### `sum(x)`
#####

function frule((_, ẋ), ::typeof(sum), x; dims=:)
    return sum(x; dims=dims), sum(ẋ; dims=dims)
end

function rrule(::typeof(sum), x::AbstractArray; dims=:)
    project = ProjectTo(x)
    y = sum(x; dims=dims)
    function sum_pullback(dy_raw)
        dy = unthunk(dy_raw)
        x_thunk = InplaceableThunk(
            # Protect `dy` from broadcasting, for when `x` is an array of arrays:
            dx -> dx .+= (dims isa Colon ? Ref(dy) : dy),
            @thunk project(_unsum(x, dy, dims))  # `_unsum` handles Ref internally
        )
        return (NoTangent(), x_thunk)
    end
    return y, sum_pullback
end

# This broadcasts `dy` to the shape of `x`, and should preserve e.g. CuArrays, StaticArrays.
# Ideally this would only need `typeof(x)` not `x`, but `similar` only has a suitable method
# when `eltype(x) == eltype(dy)`, which isn't guaranteed.
_unsum(x, dy, dims) = broadcast(last ∘ tuple, x, dy)
_unsum(x, dy, ::Colon) = broadcast(last ∘ tuple, x, Ref(dy))

# Allow for second derivatives of `sum`, by writing rules for `_unsum`:

function frule((_, _, dydot, _), ::typeof(_unsum), x, dy, dims)
    return _unsum(x, dy, dims), _unsum(x, dydot, dims)
end

function rrule(::typeof(_unsum), x, dy, dims)
    z = _unsum(x, dy, dims)
    function _unsum_pullback(dz)
        return (NoTangent(), NoTangent(), sum(unthunk(dz); dims=dims), NoTangent())
    end
    return z, _unsum_pullback
end

#####
##### `sum(f, x)`
#####

# Can't map over Adjoint/Transpose Vector
function rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(sum),
    f,
    xs::Union{Adjoint{<:Number,<:AbstractVector},Transpose{<:Number,<:AbstractVector}};
    kwargs...,
)
    op = xs isa Adjoint ? adjoint : transpose
    # since summing a vector we don't need to worry about dims which simplifies adjointing
    vector = parent(xs)
    y, vector_sum_pb = rrule(config, sum, f, vector; kwargs...)
    function covector_sum_pb(ȳ)
        s̄um, f̄, v̄ = vector_sum_pb(ȳ)
        return s̄um, f̄, op(v̄)
    end

    return y, covector_sum_pb
end

function rrule(
    config::RuleConfig{>:HasReverseMode}, ::typeof(sum), f, xs::AbstractArray; dims=:
)
    fx_and_pullbacks = map(x -> rrule_via_ad(config, f, x), xs)
    y = sum(first, fx_and_pullbacks; dims=dims)

    pullbacks = last.(fx_and_pullbacks)

    project = ProjectTo(xs)

    function sum_pullback(ȳ)
        call(f, x) = f(x)  # we need to broadcast this to handle dims kwarg
        f̄_and_x̄s = call.(pullbacks, ȳ)
        # no point thunking as most of work is in f̄_and_x̄s which we need to compute for both
        f̄ = if fieldcount(typeof(f)) === 0 # Then don't need to worry about derivative wrt f
            NoTangent()
        else
            sum(first, f̄_and_x̄s)
        end
        x̄s = map(unthunk ∘ last, f̄_and_x̄s) # project does not support receiving InplaceableThunks
        return NoTangent(), f̄, project(x̄s)
    end
    return y, sum_pullback
end

# https://github.com/JuliaDiff/ChainRules.jl/issues/522
# The rule above assumes `f` is callable. Arrays are not, this came up when summing
# arrays with weights in StatsBase
@opt_out ChainRulesCore.rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(sum),
    x::AbstractArray,
    y::AbstractArray;
    dims=:,
)

function frule(
    (_, _, Δx), ::typeof(sum), ::typeof(abs2), x::AbstractArray{T}; dims=:
) where {T<:Union{Real,Complex}}
    ẋ = unthunk(Δx)
    y = sum(abs2, x; dims=dims)
    ∂y = if dims isa Colon
        2 * real(dot(x, ẋ))
    elseif VERSION ≥ v"1.2" # multi-iterator mapreduce introduced in v1.2
        mapreduce(+, x, ẋ; dims=dims) do xi, dxi
            2 * _realconjtimes(xi, dxi)
        end
    else
        2 * sum(_realconjtimes.(x, ẋ); dims=dims)
    end
    return y, ∂y
end

function rrule(
    ::typeof(sum), ::typeof(abs2), x::AbstractArray{T}; dims=:
) where {T<:Union{Real,Complex}}
    y = sum(abs2, x; dims=dims)
    function sum_abs2_pullback(ȳ)
        x_thunk = InplaceableThunk(
            dx -> dx .+= 2 .* real.(ȳ) .* x, @thunk(2 .* real.(ȳ) .* x)
        )
        return (NoTangent(), NoTangent(), x_thunk)
    end
    return y, sum_abs2_pullback
end

# Fix dispatch for this pidgeon-hole optimization,
# Rules with RuleConfig dispatch with priority over without (regardless of other args).
# and if we don't specify what do do for one that HasReverseMode then it is ambigious
for Config in (RuleConfig, RuleConfig{>:HasReverseMode})
    @eval function rrule(
        ::$Config, ::typeof(sum), ::typeof(abs2), x::AbstractArray{T}; dims=:
    ) where {T<:Union{Real,Complex}}
        return rrule(sum, abs2, x; dims=dims)
    end
end

#####
##### `prod`
#####

function rrule(::typeof(prod), x::AbstractArray{T}; dims=:) where {T<:CommutativeMulNumber}
    y = prod(x; dims=dims)
    project_x = ProjectTo(x)
    # vald = dims isa Colon ? nothing : dims isa Integer ? Val(Int(dims)) : Val(Tuple(dims))
    function prod_pullback(ȳ)
        dy = unthunk(ȳ)
        x_thunk = InplaceableThunk(
            # In-place versions -- same branching
            dx -> if dims === (:)
                ∇prod!(dx, x, dy, y)
            elseif any(iszero, x)
                vald = if dims isa Colon
                    nothing
                elseif dims isa Integer
                    Val(Int(dims))
                else
                    Val(Tuple(dims))
                end
                ∇prod_dims!(dx, vald, x, dy, y)
            else
                dx .+= conj.(y ./ x) .* dy
            end,
            # Out-of-place versions
            @thunk project_x(
                if dims === (:)
                    ∇prod(x, dy, y)
                elseif any(iszero, x)  # Then, and only then, will ./x lead to NaN
                    vald = if dims isa Colon
                        nothing
                    elseif dims isa Integer
                        Val(Int(dims))
                    else
                        Val(Tuple(dims))
                    end
                    ∇prod_dims(vald, x, dy, y)  # val(Int(dims)) is about 2x faster than Val(Tuple(dims))
                else
                    conj.(y ./ x) .* dy
                end,
            )
        )
        return (NoTangent(), x_thunk)
    end
    return y, prod_pullback
end

function ∇prod_dims(vald::Val{dims}, x, dy, y=prod(x; dims=dims)) where {dims}
    T = promote_type(eltype(x), eltype(dy))
    dx = fill!(similar(x, T, axes(x)), zero(T))
    ∇prod_dims!(dx, vald, x, dy, y)
    return dx
end

function ∇prod_dims!(dx, ::Val{dims}, x, dy, y) where {dims}
    iters = ntuple(d -> d in dims ? tuple(:) : axes(x, d), ndims(x))  # Without Val(dims) this is a serious type instability
    @inbounds for ind in Iterators.product(iters...)
        jay = map(i -> i isa Colon ? 1 : i, ind)
        @views ∇prod!(dx[ind...], x[ind...], dy[jay...], y[jay...])
    end
    return dx
end

function ∇prod(x, dy::Number=1, y::Number=prod(x))
    T = promote_type(eltype(x), eltype(dy))
    dx = fill!(similar(x, T, axes(x)), zero(T)) # axes(x) makes MArray on StaticArrays, Array for structured matrices
    ∇prod!(dx, x, dy, y)
    return dx
end

function ∇prod!(dx, x, dy::Number=1, y::Number=prod(x))
    numzero = iszero(y) ? count(iszero, x) : 0
    if numzero == 0  # This can happen while y==0, if there are several small xs
        dx .+= conj.(y ./ x) .* dy
    elseif numzero == 1
        ∇prod_one_zero!(dx, x, dy)
    else
        # numzero > 1, then all first derivatives are zero
    end
    return dx
end

function ∇prod_one_zero!(dx, x, dy::Number=1)  # Assumes exactly one x is zero
    i_zero = 0
    p_rest = one(promote_type(eltype(x), typeof(dy)))
    for i in eachindex(x)
        xi = @inbounds x[i]
        p_rest *= ifelse(iszero(xi), one(xi), conj(xi))
        i_zero = ifelse(iszero(xi), i, i_zero)
    end
    dx[i_zero] += p_rest * dy
    return nothing
end

#####
##### `cumprod`
#####

function rrule(::typeof(cumprod), x::AbstractVector{<:Real}; dims::Integer=1)
    y = cumprod(x; dims=dims)  # does nothing unless dims == 1
    project_x = ProjectTo(x)
    function cumprod_pullback_1(dy_raw)
        dy = unthunk(dy_raw)
        dx_thunk = InplaceableThunk(
            dx -> if dims == 1
                ∇cumprod!(dx, x, dy, y)
            else
                dx .+= dy
            end,
            @thunk project_x(
                if dims == 1
                    ∇cumprod(x, dy, y)
                else
                    dy
                end,
            )
        )
        return (NoTangent(), dx_thunk)
    end
    return y, cumprod_pullback_1
end

function rrule(::typeof(cumprod), x::AbstractArray{<:Real}; dims::Integer)
    y = cumprod(x; dims=dims)
    project_x = ProjectTo(x)
    function cumprod_pullback_2(dy_raw)
        dy = unthunk(dy_raw)
        dx_thunk = InplaceableThunk(
            dx -> if dims <= ndims(x)
                vald = Val(Int(dims))
                ∇cumprod_dim!(dx, vald, x, dy, y)
            else
                dx .+= dy
            end,
            @thunk project_x(
                if dims <= ndims(x)
                    vald = Val(Int(dims))
                    ∇cumprod_dim(vald, x, dy, y)
                else
                    dy
                end,
            )
        )
        return (NoTangent(), dx_thunk)
    end
    return y, cumprod_pullback_2
end

function ∇cumprod_dim(
    vald::Val{dim}, x::AbstractArray, dy=fill!(zero(x), 1), y=cumprod(x; dims=dim)
) where {dim}
    T = promote_type(eltype(x), eltype(dy))
    dx = fill!(similar(x, T, axes(x)), zero(T))
    ∇cumprod_dim!(dx, vald, x, dy, y)
    return dx
end

@inline function ∇cumprod_dim!(
    dx::AbstractArray, ::Val{dim}, x::AbstractArray, dy, y
) where {dim}
    iters = ntuple(k -> k == dim ? Ref(:) : axes(x, k), ndims(x))
    for ind in Iterators.product(iters...)
        @views ∇cumprod!(dx[ind...], x[ind...], dy[ind...], y[ind...])
    end
    return dx
end

function ∇cumprod(x::AbstractVector, dy=one(x), y=cumprod(x))
    T = promote_type(eltype(x), eltype(dy))  # really needs to allow dy * y / x
    dx = fill!(similar(x, T, axes(x)), zero(T))  # axes(x) makes MArray on StaticArrays, Array for structured matrices
    ∇cumprod!(dx, x, dy, y)
    return dx
end

@inline function ∇cumprod!(dx::AbstractVector, x::AbstractVector, dy, y)
    lo, hi = firstindex(x), lastindex(x)
    z = something(findfirst(iszero, x), hi + 1)
    acc = zero(eltype(dy))
    @inbounds for k in (z - 1):-1:lo
        acc += y[k] * dy[k]
        dx[k] += acc / x[k]
    end
    @inbounds if z != hi + 1
        yk = z == 1 ? one(eltype(y)) : y[z - 1]  # will be prod(x[j] for j=1:k if j!=z)
        dx[z] += yk * dy[z]
        for k in (z + 1):hi
            yk *= x[k]
            dx[z] += yk * dy[k]
        end
    end
    return dx
end


text:
module ArrayPairsTests

using Test
using Gridap.Arrays

a = collect(1:10)
b = collect(11:20)
r = [(i,i+10) for i in 1:10]

c = pair_arrays(a,b)
test_array(c,r)

_a, _b = unpair_arrays(r)
test_array(_a,a)
test_array(_b,b)

#print_op_tree(_b)

_a, _b = unpair_arrays(c)
@test a === _a
@test b === _b

end # module


text:
#
# map
#

t = init_task(map_, x -> x + 5)
@test consume(t, Subject(1)) == Subject(6)

#
# filter
#

t = init_task(filter_, x -> x > 10)
@test consume(t, Subject(11)) == Subject(11)

#
# reduce
#

t = init_task(reduce, (x, y) -> x + y, Subject(0))
@test consume(t, Subject(1)) == Subject(1)
@test consume(t, Subject(2)) == Subject(3)
@test consume(t, Subject(3)) == Subject(6)

#
# bridge
#

t1 = init_task(map_, x -> x + 5)
t2 = init_task(map_, x -> x * 2)
t3 = init_task(bridge, t1, t2)
@test consume(t3, Subject(2)) == Subject(14)

t4 = init_task(filter_, x -> x > 10)
t5 = init_task(bridge, t3, t4)
@test consume(t5, Subject(5)) == Subject(20)

#
# Broadcast
#

parallel = [x -> x+2, x -> x*2, x -> x+2, x -> x*2, x -> x+2, x -> x*2]
answer = [Subject(13), Subject(22), Subject(13), Subject(22), Subject(13), Subject(22)]
t = init_task(bridge, x -> x+1, parallel)
@test consume(t, Subject(10)) == answer

#
# Subject chaining
#

a = Subject(10)
@test a(map_, x -> x + 2) == Subject(12)
@test a(map_, x -> x + 2)(filter_, x -> x > 10) == Subject(12)


text:
module SplitApplyCombine

using DataFrames, Query, Statistics

person_data = [
    ("Alice", "apple", "apple pie", "yes", 10),
    ("Balin", "banana", "bananas foster", "yes", 53),
    ("Cal", "carrot", "carrot cake", "no", 34),
    ("Dara", "dates", "doughnut", "yes", 83),
    ("Ezekiel", "elderberry", "elderberry tart", "no", 4),
    ("Frances", "fig", "fig bars", "no", 34)
]

person = DataFrame(name = String[], fav_fruit_veg = String[], fav_treat = String[], sweet_tooth = String[], age = Int[])
for p in person_data
    push!(person, p)
end

# DataFrames
sweet_tooth_by_age = combine(groupby(person, :sweet_tooth),
                             nrow => :age_cnt,
                             :age => mean => :age_avg)

# Query - stand alone
sweet_tooth_by_age2 = person |>
    @groupby(_.sweet_tooth) |>
    @map({sweet_tooth=key(_), age_cnt=length(_), age_avg=mean(_.age)}) |>
    DataFrame

# Query - LINQ
sweet_tooth_by_age3 = @from i in person begin
    @group i by i.sweet_tooth into g
    @select {sweet_tooth=key(g), age_cnt=length(g), age_avg=mean(g.age)}
    @collect DataFrame
end

end # module

max_stars_repo_name: curtisalexander/saspyrilia-sandboxes

max_issues_repo_name: curtisalexander/saspyrilia-sandboxes

max_forks_repo_name: curtisalexander/saspyrilia-sandboxes


text:
function Well19937c()
    return Well19937c(())
end

function Well19937c(arg0::Vector{jint})
    return Well19937c((Vector{jint},), arg0)
end

function Well19937c(arg0::jint)
    return Well19937c((jint,), arg0)
end

function Well19937c(arg0::jlong)
    return Well19937c((jlong,), arg0)
end

function next_int(obj::Well19937c)
    return jcall(obj, "nextInt", jint, ())
end


text:
#=
  Utilities for Turing.jl

  Here are utils for my Turing.jl programming.

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Turing.jl page: http://www.hakank.org/julia/turing/

=#

using Turing
using Printf, Distributions
using LinearAlgebra

"""
Return a Dict of the elements and their occurences
in the collection a.
"""
function make_hash(a)
    d = Dict()
    for e in collect(a)
        # get!(d,e,0)
        if !haskey(d,e)
            d[e] = 0
        end
        d[e] += 1
    end
    return d
end

"""
This is to simplify the conversion from WebPPL
"""
function flip(p=0.5)
    return Bernoulli(p)
end

"""
Return a normalized vector, i.e. where the sum is 1
"""
function simplex(v)
    return v./sum(v)
end

"""
Show a sorted Dict with keys, values and percentages
Example
julia> show_var_dist(chns,:d1)

Distributions of variable d1
 1 =>   13654  0.341350
 2 =>   13417  0.335425
 3 =>   12929  0.323225

Note: I'm not sure how to get more than one variables
E.g. this don't work now:
 julia> show_var_dist(chns,[:d1,:d2])
"""
function show_var_dist(chns, var)

    if var in chns.name_map.parameters
        println("Distributions of variable $var")
        len = length(vcat(chns[var]...)) # handle multiple chains
        for kv in sort(make_hash(chns[var]))
            @printf "%-3.5f => % 7d  %2.6f\n" kv[1] kv[2] kv[2]/len
        end
    else
        # println("Variable $var is not in chains")
    end
end

"""
Show distribution of a variable in a MCMCChain
Sort the dictionary in order of decreasing occurrence (percentage)
Examples:

 - show_var_dist_pct(chns, :n)       show all entries

 - show_var_dist_pct(chns, :n, 10)  show first 10 entries (e.g. for large tables)
"""
function show_var_dist_pct(chns::Chains, var, num=0)

    if var in chns.name_map.parameters
        println("Distributions of variable $var (num:$num)")
        len = length(vcat(chns[var]...)) # handle multiple chains
        c = 0
        for kv in sort(collect(make_hash(chns[var])),by=x->x[2],rev=true)
            c += 1
            if (num == 0) || (num > 0 && c <= num)
                @printf "%-3.5f => % 7d  (%2.6f)\n" kv[1] kv[2] kv[2]/len
            end
        end
    else
        # println("Variable $var is not in chains")
    end
end

"""
Show distribution of a variable in a MCMCChain

Sort the dictionary in order of decreasing occurrence (percentage)
and present the values from array `a` where the position in
`a` represents the values `1..length(a)`.

Examples:

 - show_var_dist_pct(chns, :x, ["yellow","blue", "green"])

"""
function show_var_dist_pct(chns::Chains, var, a::Array)

    if var in chns.name_map.parameters
        println("Distributions of variable $var")
        len = length(vcat(chns[var]...)) # handle multiple chains
        for kv in sort(collect(make_hash(chns[var])),by=x->x[2],rev=true)
            ix = round(Int,kv[1])
            @printf "%-10s => % 7d  (%2.6f)\n" a[ix] kv[2] kv[2]/len
        end
    else
        # println("Variable $var is not in chains")
    end
end

"""
Show distribution of the elements in a (simple) Array.
Sort the dictionary in order of decreasing occurrence (percentage)
Examples:

 - show_var_dist_pct(a)       show all entries

 - show_var_dist_pct(a, 10)  show first 10 entries (e.g. for large tables)

Full example:

julia> show_var_dist_pct(rand(Binomial(100,0.1),100000),10)

 Distributions of variable (num:10)
 10.00000 =>   13143  0.131430
 9.00000 =>   13094  0.130940
 11.00000 =>   12002  0.120020
 8.00000 =>   11403  0.114030
 12.00000 =>    9989  0.099890
 7.00000 =>    8765  0.087650
 13.00000 =>    7475  0.074750
 6.00000 =>    5968  0.059680
 14.00000 =>    5117  0.051170
 5.00000 =>    3450  0.034500
"""
function show_var_dist_pct(a::Array, num=0)
    println("Distributions of variable (num:$num)")
    c = 0
    len = length(a)
    isArray = false
    if length(a[1]) > 1 # typeof(a[1]) <: Array
        isArray = true
    end
    for kv in sort(collect(make_hash(a)),by=x->x[2],rev=true)
        c += 1
        if (num == 0) || (num > 0 && c <= num)
            if isArray
                # This is an array, vector, tuple of elements
                # len = length(kv[1])
                # format = "%-$(len)s => % 7d  %2.6f\n"
                # @printf format kv[1] kv[2] kv[2]/len
                # @printf "%-35s => % 7d  %2.6f\n" kv[1] kv[2] kv[2]/len
                println("$(kv[1])\t=>\t$(kv[2]) ($(kv[2]/len))")
            else
                @printf "%-3.5f => % 7d  (%2.6f)\n" kv[1] kv[2] kv[2]/len
            end
        end
    end
end

"""
This is ported from the WebPPL function in
https://mhtess.github.io/bdappl/chapters/03-simpleModels.html

Example in credible_interval_test.jl:

julia> credible_interval(chns, "posteriorPredictive",0.90)

  credible interval for posteriorPredictive with mass 0.9: (10.000000 .. 18.000000)
"""
function credible_interval(chns::Chains, var, credMass=0.95)
    # Using sort and vcat don't work, and this reshaping is clunky. TODO!
    sss = prod(size(chns[var]))
    chnsCat = reshape(chns[var],sss,1)
    sortedPts = sort(chnsCat,dims=1)
    len = length(sortedPts)
    ciIdxInc = round(Int,ceil(credMass*len))
    nCIs = len - ciIdxInc
    ciWidth = map(i->sortedPts[i+ciIdxInc]-sortedPts[i], 1:nCIs)
    _, i = findmin(ciWidth)
    print("Credible interval for $(var) with mass $(credMass): ")
    if typeof(sortedPts[1]) == Int64
        @printf("(%d..%d)\n", sortedPts[i],sortedPts[i+ciIdxInc])
    else
        @printf("(%f..%f)\n", sortedPts[i],sortedPts[i+ciIdxInc])
    end

end

"""
Return mean value of the variable var in the chain
"""
function mean_val(chns::Chains, var)
    mean(chns[var])
end

"""
`UniformDraw(xs)`

Draw one element from `xs` randomly with uniform probability.

For a version with handcrafted probability, see `DiscreteNonParametric`.

"""
UniformDraw(xs) = DiscreteNonParametric(xs, ones(length(xs)) ./ length(xs))

"""
`confusion_matrix(y,pred,cats,print_all=true)``

Given the true classifications `y` and the predictions `pred`,
it returns a confusion matrix as well as a bads array for the wrongly
predicted instances (in the for id=>(true_cat,pred_cat)).

Also prints info about all the instances if `print_all=true`.

`cats` is a dictionary converting the numeric category to a
number 1:number_of categories; this is for handling the
confusion matrix. (Perhaps there's a better way of doing this...)
"""
function confusion_matrix(y,pred,cats,print_all=true)
    n = length(y)
    bads = []
    num_cat = cats  |> length
    confusion = zeros(num_cat,num_cat)
    min_val, max_val = extrema(keys(cats))
    for i in 1:n
      true_cat = round(Int64,y[i])
      p = mean(pred["y[$i]"])
      diff = y[i]-p
      pred_cat = round(Int64,p)
      # println("true_cat: $true_cat pred_cat: $pred_cat")
      # Handle cases were we got predicates out of range
      # i.e. 0 when the categories is 1..
      if pred_cat < min_val
        println("pred_cat < min_val. Adjusting to min_val")
        pred_cat = min_val
      end
      if pred_cat > max_val
        println("pred_cat > max_val. Adjusting to max_val")
        pred_cat = max_val
      end
      abs_diff = abs(diff)

      if pred_cat != true_cat
        bad = "BAD!"
        # the correct cat is predicted (wrongly) as pred_cat
        push!(bads, i=>(true_cat=>pred_cat))
        confusion[cats[true_cat],cats[pred_cat]] += 1
      else
        bad = ""
        confusion[cats[true_cat],cats[true_cat]] += 1
      end
      if print_all
        println("y[$i]: true cat: $true_cat pred[$i]: ", p, " diff: ", y[i]-p, " pred cat: $pred_cat $bad")
      end
    end

    num_correct = sum(diag(confusion))
    num_total = sum(confusion)
    println("Num correct: $num_correct")
    println("Num incorrect: ", num_total - num_correct)
    println("Accuracy: ",  (num_correct / num_total))
    return confusion, bads
end

"""
`linear_regression_model(x,y=missing)`

BDA model for linear regression using MvNormal and rather strict priors.

`x` is the dataframe with the input features, `y` is the output variable (a single value).

The model can be used in prediction with the following workflow:
```
# Training: x and y is the training data
model = linear_regression_model(x,y)
chns = sample(model, NUTS(), 1_000)
display(chns)

# Prediction, note the "missing" argument instead if y
mpred = linear_regression_model(x,missing)
pred = predict(mpred,chns)
display(pred)`
```

"""
@model function linear_regression_model(x,y=missing, ::Type{T} = Float64) where {T}
    n,p = size(x)

    if y === missing
      y = Vector{T}(undef, n)
    end

    a ~ Normal(0,1)
    params ~ filldist(Normal(0,1),p)
    sigma ~ Uniform(0,1)

    mu = a .+ x*params
    y ~ MvNormal(mu, sigma)

    return y

end

"""
`split_train_test(x,y,train_pct=0.75)`

Split the datasets `x` (the inputs) and `y` (the outputs) into
training and test sets with a size of the training set of `train_pct`.

Returns `x_train, x_test, y_train, y_test`
"""
function split_train_test(x,y,train_pct=0.75)
    num_rows = length(y)
    num_train = round(Int64,num_rows*train_pct)

    train_ix = StatsBase.sample(1:num_rows,num_train,replace=false)
    test_ix = [i for i in 1:num_rows if !(i in train_ix)]

    # Split into train and test
    x_train = x[train_ix,:]
    x_test = x[test_ix,:]

    y_train = y[train_ix]
    y_test = y[test_ix]

    x_train, x_test, y_train, y_test
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#=

Log-likelihood for MuSSE type ODEs, integrating along the tree

Ignacio Quintero Mächler

t(-_-t)

Created 20 09 2017

=#

"""
    make_llf(tip_val::Dict{Int64,Array{Float64,1}},
             ed     ::Array{Int64,2},
             el     ::Array{Float64,1},
             ode_fun,
             af     ::Function,
             p      ::Array{Float64,1};
             sbrlen ::Float64 = 5.0)

Make likelihood function for a tree given an ODE function.
"""
function make_llf(tip_val::Dict{Int64,Array{Float64,1}},
                  ed     ::Array{Int64,2},
                  el     ::Array{Float64,1},
                  ode_fun,
                  af     ::Function,
                  p      ::Array{Float64,1},
                  md     ::Bool,
                  ws     ::Bool;
                  sbrlen ::Float64 = 5.0)

  k    = length(tip_val[1])::Int64
  k2   = 2k::Int64
  ke   = (k+1:2k)::UnitRange{Int64}
  ntip = length(tip_val)::Int64
  bk   = k + (k*k)
  r    = Array{Float64}(undef, k)

  # add long root branch
  ed = cat(ed, [2*ntip ntip + 1], dims = 1)
  push!(el, sbrlen)

  # get absolute times of branches as related to z(t)
  elrt = abs_time_branches(el, ed, ntip)

  ne    = size(ed,1)::Int64
  child = ed[:,2]::Array{Int64,1}
  wtp   = findall(child .<= ntip)::Array{Int64,1}

  # make trios
  trios = maketriads(ed, rev = true)::Array{Array{Int64,1},1}

  # preallocate tip likelihoods
  led = Array{Float64,1}[]
  for i in Base.OneTo(ne)
    push!(led, zeros(Float64, k2))
  end

  # assign states to terminal branches
  for wi in wtp
     led[wi][1:k] = tip_val[child[wi]]
  end

  # partial log-likelihoods
  llik = Array{Float64}(undef,k)

  # root likelihood
  rll = Array{Float64}(undef,k2)

  # log lambdas
  lλs = Array{Float64}(undef,k)

  # make ode solver
  ode_solve = make_solver(ode_fun, p)

  λevent! = make_λevent(af, k, ws, md)

  function f(p::Array{Float64})

    @inbounds begin

      ll     = 0.0::Float64
      llxtra = 0.0::Float64

      for i in Base.OneTo(k)
        lλs[i] = log(p[i])::Float64
      end

      # loop for integrating over internal branches
      for triad in trios

        pr, d1, d2 = triad::Array{Int64,1}

        ud1 = ode_solve(led[d1], p, elrt[d1,2], elrt[d1,1])::Array{Float64,1}

        if check_negs(ud1, k)
          return -Inf
        end

        ud2 = ode_solve(led[d2], p, elrt[d2,2], elrt[d2,1])::Array{Float64,1}

        if check_negs(ud2, k)
          return -Inf
        end

        # update likelihoods with speciation event
        λevent!(elrt[pr,2], llik, ud1, ud2, lλs, p)

        # loglik to sum for integration
        tosum   = minimum(llik)
        llxtra -= tosum

        # assign the remaining likelihoods
        for i in Base.OneTo(k)
          led[pr][i] = exp(llik[i] - tosum)
        end

        # assign extinction probabilities
        @views led[pr][ke] = ud1[ke]
      end

      # stem branch
      rll = ode_solve(led[ne], p, elrt[ne,2], elrt[ne,1])::Array{Float64,1}

      if check_negs(rll, k)
          return -Inf
      end

      for i in Base.OneTo(k)
        ll += rll[i]::Float64
      end

      return (log(ll) - llxtra)::Float64
    end
  end

end

"""
    make_λevent(af, k::Int64, ws::Bool, md::Bool)

Make function for speciation event likelihoods
"""
function make_λevent(af, k::Int64, ws::Bool, md::Bool)

  r  = Array{Float64}(undef, k)
  bk = k + (k*k)

  # speciation ESSE model with MULTIvariate z(t)
  function f1(t   ::Float64,
              llik::Array{Float64,1},
              ud1 ::Array{Float64,1},
              ud2 ::Array{Float64,1},
              lλs ::Array{Float64,1},
              p   ::Array{Float64,1})

    af(t, r)

    # estimate ll to sum,
    for i in Base.OneTo(k)
      llik[i] = log(ud1[i] * ud2[i]) + lλs[i] + p[bk+i] * r[i]
    end

    return nothing
  end

  # speciation ESSE model with UNIvariate z(t)
  function f2(t   ::Float64,
              llik::Array{Float64,1},
              ud1 ::Array{Float64,1},
              ud2 ::Array{Float64,1},
              lλs ::Array{Float64,1},
              p   ::Array{Float64,1})

    r = af(t)

    # estimate ll to sum,
    for i in Base.OneTo(k)
      llik[i] = log(ud1[i] * ud2[i]) + lλs[i] + p[bk+i] * r
    end

    return nothing
  end

  # other ESSE model
  function f3(t   ::Float64,
              llik::Array{Float64,1},
              ud1 ::Array{Float64,1},
              ud2 ::Array{Float64,1},
              lλs ::Array{Float64,1},
              p   ::Array{Float64,1})

    # estimate ll to sum,
    for i in Base.OneTo(k)
      llik[i] = log(ud1[i] * ud2[i]) + lλs[i]
    end

    return nothing
  end

  if ws
    if md
      return f1
    else
      return f2
    end
  else
    return f3
  end
end

"""
    make_lpf(λpriors::Float64, μpriors::Float64, qpriors::Float64, βpriors, k::Int64)

Make log-prior function.
"""
function make_lpf(λpriors::Float64,
                  μpriors::Float64,
                  qpriors::Float64,
                  βpriors::Tuple{Float64,Float64},
                  k      ::Int64,
                  T      ::Bool)

  l2::UnitRange{Int64} = (k+1):2k
  l3::UnitRange{Int64} = (2k+1):(k+k*k)
  lp::Float64          = 0.0

  if T
    l4 = ((k+k*k+1):(2k+k*k))::UnitRange{Int64}
  else
    l4 = ((k+k*k+1):(2k*k))::UnitRange{Int64}
  end

  function f(p::Array{Float64,1})
    lp = 0.0::Float64

    for i in Base.OneTo(k)
      lp += logdexp(p[i], λpriors)::Float64
    end

    for i in l2
      lp += logdexp(p[i], μpriors)::Float64
    end

    for i in l3
      lp += logdexp(p[i], qpriors)::Float64
    end

    for i in l4
      lp += logdnorm(p[i], βpriors[1], βpriors[2])::Float64
    end

    return lp::Float64
  end

  return f
end

"""
    make_lpf(λpriors::Float64, μpriors::Float64, qpriors::Float64, k::Int64)

Make log-prior function.
"""
function make_lpf(λpriors::Float64,
                  μpriors::Float64,
                  qpriors::Float64,
                  k      ::Int64)

  l2::UnitRange{Int64} = (k+1):2k
  l3::UnitRange{Int64} = (2k+1):(k+k*k)
  lp::Float64          = 0.0

  function f(p::Array{Float64,1})
    lp = 0.0::Float64

    for i in Base.OneTo(k)
      lp += logdexp(p[i], λpriors)::Float64
    end

    for i in l2
      lp += logdexp(p[i], μpriors)::Float64
    end

    for i in l3
      lp += logdexp(p[i], qpriors)::Float64
    end

    return lp::Float64
  end

  return f
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using Pkg
using SafeTestsets
using Test
using SciMLBase

@test isempty(detect_ambiguities(SciMLBase))

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

@time begin

if GROUP == "Core" || GROUP == "All"
    @time @safetestset "Display" begin include("display.jl") end
    @time @safetestset "Existence functions" begin include("existence_functions.jl") end
    @time @safetestset "Integrator interface" begin include("integrator_tests.jl") end
    @time @safetestset "Ensemble functionality" begin include("ensemble_tests.jl") end
end

if !is_APPVEYOR && GROUP == "Downstream"
    activate_downstream_env()
    @time @safetestset "Ensembles of Zero Length Solutions" begin include("downstream/ensemble_zero_length.jl") end
end
end


text:
using Documenter, Rayuela

makedocs(
    format = :html,
    sitename = "Rayuela.jl",
    modules = [Rayuela],
    pages = [
        "Home"    => "index.md",
        # "Manual"  => [
        #     "man/usage.md"
        # ],
        "Library" => [
            "PQ.md",
            "OPQ.md",
            "RVQ.md",
            "ERVQ.md",
            "ChainQ.md",
            "LSQ.md"
        ]
    ]
    # doctest = test
)

# Borrowing from CUDAnative.jl
deploydocs(
	repo = "github.com/una-dinosauria/Rayuela.jl.git",
  julia  = "",
  osname = "",
  # no need to build anything here, re-use output of `makedocs`
	target = "build",
	deps = nothing,
	make = nothing,
)