text:
#=
    The File contains `Subspacik<T>` and related funcs implementation
=#

# I guess find_basis funcs should be

#------------------------------------------------------------------------------

include("dok_sparsik.jl")
include("densik.jl")

#------------------------------------------------------------------------------

import Nemo: QQ, GF, Field

#------------------------------------------------------------------------------

# Subspacik
mutable struct Subspacik{T<:Field}
    field::T
    # the echelon form of vectors is an invariant!
    echelon_form::Dict{Int, AbstractSparseObject{T}}
end

#------------------------------------------------------------------------------

# Subspacik
mutable struct HashedSubspacik{T<:Field}
    # inherited from Subspacik
    field::T
    # pivot --> row
    # maybe we can call this "rref"
    echelon_form::Dict{Int, AbstractSparseObject{T}}

    # -----
    # new hash things

    hash_vector::AbstractSparseObject{T}
    # this vector is reduced over time
    reduced_hash_vector::AbstractSparseObject{T}

    # pivot --> hash
    hashes::Dict{Int, FieldElem}
end

#------------------------------------------------------------------------------

# convenience ctor
function Subspacik(field::T) where {T<:Field}
    return Subspacik(field, Dict{Int, AbstractSparseObject{T}}())
end

# deepcopy redefinition in order to preserve the field
function Base.deepcopy_internal(x::Subspacik{T}, stackdict::IdDict) where {T<:Field}
    y = Subspacik(base_ring(x), Base.deepcopy_internal(x.echelon_form, stackdict))
    stackdict[x] = y
    return y
end

#------------------------------------------------------------------------------

# base_ring field!
# Gleb: also suggest base_ring
# Alex: to be changed soon..
base_ring(v::Union{Subspacik, HashedSubspacik}) = v.field

#------------------------------------------------------------------------------

# returns basis vectors of V sorted by pivot indices
function basis(V::Union{Subspacik, HashedSubspacik})
    return sort(collect(deepcopy(values(V.echelon_form))), by=first_nonzero)
end

#------------------------------------------------------------------------------

# returns the dimension of the ambient space of V
function ambient_dim(V::Union{Subspacik, HashedSubspacik})
    return dim(first(values(V.echelon_form)))
end

# returns the dimension of V
function dim(V::Union{Subspacik, HashedSubspacik})
    return length(V.echelon_form)
end

#------------------------------------------------------------------------------

@enum EatCode begin
    reduced = -1
    # Gleb: I would call discarded skipped
    skipped = -2
end

# adds `new_vector` to the set of spanning vectors of V
# while modifying both
# returns -1 if `new_vector` lies in V, new pivot index otherwise

# O(kR)  if k = nnz(new_vector) and R = Σnnz(v) for v in echelon_form
function eat_sparsik!(V::Subspacik, new_vector::AbstractSparseObject; ω=1.0)
    # PRobably, we need eat_sparsik( ;verbose) which will return
    # coordinates of new_vector if V basis in case of linear dependency

    # O(kr) if k = nnz(new_vector) and r = len(echelon_form)
    for (piv, vect) in V.echelon_form
        if haskey(new_vector, piv)
            reduce!(new_vector, vect, -new_vector[piv])
        end
    end

    if iszero(new_vector)
        return reduced
    end

    # it should(must) be optimized somehow
    # O(k)
    if density(new_vector) > ω
        return skipped
    end

    # O(k)
    pivot = first_nonzero(new_vector)
    scale!(new_vector, inv(new_vector[pivot]))

    # O(kR)  if k = nnz(new_vector) and R = Σnnz(v) for v in echelon_form
    for (piv, vect) in V.echelon_form
        if haskey(vect, pivot)
            reduce!(V.echelon_form[piv], new_vector, -vect[pivot])
        end
    end

    V.echelon_form[pivot] = new_vector
    return pivot
end

#------------------------------------------------------------------------------

# returns a new Subspacik formed as a linear span
# of the given vectors
function linear_span!(vectors)
    V = Subspacik(base_ring(first(vectors)))
    for vect in vectors
        eat_sparsik!(V, vect)
    end
    return V
end

#------------------------------------------------------------------------------

 # transformes V to the smallest subspace invariant under
 # the given matrices and containing the original one
 function apply_matrices_inplace!(V::Union{Subspacik, HashedSubspacik}, matrices; ω=1.)
     global CNT
     global NONZERO

     fat_vectors = [ ]
     new_pivots = collect(keys(V.echelon_form))
     i = 0

     while !isempty(new_pivots)
         pivots_to_process = sort(deepcopy(new_pivots))
         empty!(new_pivots)

         for pivot in pivots_to_process
             for vect in matrices
                 product = apply_vector(V.echelon_form[pivot], vect)
                 CNT += 1

                 i += 1
                 # i % 500 == 0 && print(".")

                 if !iszero(product)
                     new_pivot = eat_sparsik!(V, product, ω=ω)
                     if new_pivot == skipped
                         push!(fat_vectors, product)
                     elseif new_pivot != reduced
                         push!(new_pivots, new_pivot)
                     end
                 end

             end
         end

     end

     return fat_vectors

  end

#------------------------------------------------------------------------------

# checks whether the given vector lies in V
function check_inclusion(V::Union{Subspacik, HashedSubspacik}, vector::AbstractSparseObject)
    return check_inclusion!(V, deepcopy(vector))
end

# checks whether the given vector lies in V
#
function check_inclusion!(V::Subspacik, vector::AbstractSparseObject)
    for (piv, vect) in V.echelon_form
        if ! iszero(vector[piv])
            reduce!(vector, vect, -vector[piv])
        end
    end
    return iszero(vector)
end

# checks whether the given vector lies in V
#
function check_inclusion!(V::HashedSubspacik, vector::AbstractSparseObject)
    return iszero(inner(V.reduced_hash_vector, vector))
end

#------------------------------------------------------------------------------

# checks whether the given subspace is invariant under the matrices
function check_invariance!(V::Union{Subspacik, HashedSubspacik}, matrices)
    for matr in matrices
        for vec in values(V.echelon_form)
            product = matr * vec
            if ! check_inclusion!(V, product)
                return false
            end
        end
    end
    return true
end

#------------------------------------------------------------------------------

# checks whether U ⊆ V using the provived algorithm for vector-wise check
function check_inclusion!(V::Union{Subspacik, HashedSubspacik}, U::Union{Subspacik, HashedSubspacik}; algorithm=check_inclusion!)
    return algorithm(V, basis(U))
end

#------------------------------------------------------------------------------

# check whether the given vectors lie in V
# O(K)
# where K = ∑ᵢnnz(vᵢ)  for vᵢ ∈ vectors
function check_inclusion!(V::Union{Subspacik, HashedSubspacik}, vectors::Array)
    for vec in vectors
        if ! check_inclusion!(V, vec)
            return false
        end
    end
    return true
end

#------------------------------------------------------------------------------

Base.repr(::MIME"text/plain", V::Subspacik) = "<\n$(join(map(v -> repr(v), values(V.echelon_form)), ",\n"))\n>"
Base.show(io::IO, V::Subspacik) = print(io, repr(MIME("text/plain"), V))

#------------------------------------------------------------------------------

==(U::Subspacik{T}, V::Subspacik{T}) where {T} = (
            U.field == V.field &&
            U.echelon_form == V.echelon_form;)

#------------------------------------------------------------------------------

# convenience ctor
function HashedSubspacik(hash_vector, field::T) where {T<:Field}
    return HashedSubspacik(
        field,
        Dict{Int, AbstractSparseObject{T}}(),
        deepcopy(hash_vector),
        deepcopy(hash_vector),
        Dict{Int, FieldElem}()
    )
end

# deepcopy redefinition in order to preserve the field
function Base.deepcopy_internal(x::HashedSubspacik{T}, stackdict::IdDict) where {T<:Field}
    y = HashedSubspacik(
                base_ring(x),
                Base.deepcopy_internal(x.echelon_form, stackdict),
                Base.deepcopy_internal(x.hash_vector, stackdict),
                Base.deepcopy_internal(x.reduced_hash_vector, stackdict),
                Base.deepcopy_internal(x.hashes, stackdict)
    )
    stackdict[x] = y
    return y
end

Base.repr(::MIME"text/plain", V::HashedSubspacik) = "<\n$(join(map(v -> repr(v), values(V.echelon_form)), ",\n"))\n>\nwith hash vector $(repr(V.hash_vector))"
Base.show(io::IO, V::HashedSubspacik) = print(io, repr(MIME("text/plain"), V))

#------------------------------------------------------------------------------

# adds `new_vector` to the set of spanning vectors of V
# while modifying both
# returns -1 if `new_vector` lies in V, new pivot index otherwise
function eat_sparsik!(V::HashedSubspacik, new_vector::AbstractSparseObject; ω=1.0)
    # O(k) if k = nnz(new_vector)
    reduced_hash = inner(new_vector, V.reduced_hash_vector)
    if iszero(reduced_hash)
        return reduced
    end

    for (piv, vect) in V.echelon_form
        if !iszero(new_vector[piv])
            reduce!(new_vector, vect, -new_vector[piv])
        end
    end

    if iszero(new_vector)
        return reduced
    end

    # it should(must) be optimized somehow
    # O(k)
    if density(new_vector) > ω
        return skipped
    end

    # O(k)
    pivot = first_nonzero(new_vector)
    scale!(new_vector, inv(new_vector[pivot]))
    new_hash = inner(V.hash_vector, new_vector)

    # O(kR)  if k = nnz(new_vector) and R = Σnnz(v) for v in echelon_form
    for (piv, vect) in V.echelon_form
        if !iszero(vect[pivot])
            reduce!(V.echelon_form[piv], new_vector, -vect[pivot])
            V.hashes[piv] -= new_hash * (-vect[pivot])
        end
    end

    V.hashes[pivot] = new_hash
    V.echelon_form[pivot] = new_vector

    reduce!(V.reduced_hash_vector, pivot, -new_hash)

    return pivot
end

#------------------------------------------------------------------------------

function linear_span!(vectors, hash_vector)
    V = HashedSubspacik(hash_vector, base_ring(first(vectors)))
    for vect in vectors
        eat_sparsik!(V, vect)
    end
    return V
end

#------------------------------------------------------------------------------

# returns a new Subspacik object consisting of REF elements
# each reconstructed from V.field to QQ
#
# reconstructs *not* inplace
function rational_reconstruction(V::Subspacik)
    ans = Subspacik(QQ)
    for piv in keys(V.echelon_form)
        ans.echelon_form[piv] = rational_reconstruction(V.echelon_form[piv])
    end
    return ans
end

# returns a new HashedSubspacik object consisting of REF elements
# each reconstructed from V.field to QQ
#
# hash vectors and hashes are reconstructed too
#
# reconstructs *not* inplace
function rational_reconstruction(V::HashedSubspacik)
    hash_vector = rational_reconstruction(V.hash_vector)
    reduced_hash_vector = rational_reconstruction(V.reduced_hash_vector)
    ans = HashedSubspacik(hash_vector, QQ)
    ans.reduced_hash_vector = reduced_hash_vector
    for piv in keys(V.echelon_form)
        ans.echelon_form[piv] = rational_reconstruction(V.echelon_form[piv])
    end

    tmp_type = BigInt
    if isa(base_ring(V), GaloisField)
        tmp_type = Int
    end
    ch = convert(tmp_type, characteristic(base_ring(V)))
    ans.hashes = Dict(
        i => rational_reconstruction(x, ch)
        for (i, x) in V.hashes
    )
    return ans
end


text:
# code reflection entry-points

# forward the rest to GPUCompiler with an appropriate CompilerJob
for method in (:code_typed, :code_warntype, :code_llvm, :code_native)
    # only code_typed doesn't take a io argument
    args = method == :code_typed ? (:job,) : (:io, :job)

    @eval begin
        function $method(io::IO, @nospecialize(func), @nospecialize(types);
                         kernel::Bool=false, minthreads=nothing, maxthreads=nothing,
                         blocks_per_sm=nothing, maxregs=nothing, kwargs...)
            source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
            target = MetalCompilerTarget(macos=get_macos_v();)
            params = MetalCompilerParams()
            job = CompilerJob(target, source, params)
            GPUCompiler.$method($(args...); kwargs...)
        end
        $method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $method(stdout, func, types; kwargs...)
    end
end

"""
    Metal.return_type(f, tt) -> r::Type

Return a type `r` such that `f(args...)::r` where `args::tt`.
"""
function return_type(@nospecialize(func), @nospecialize(tt))
    source = FunctionSpec(func, tt, true)
    target = MetalCompilerTarget(macos=get_macos_v();)
    params = MetalCompilerParams()
    job = CompilerJob(target, source, params)
    interp = GPUCompiler.get_interpreter(job)
    if VERSION >= v"1.8-"
        sig = Base.signature_type(job.source.f, job.source.tt)
        Core.Compiler.return_type(interp, sig)
    else
        Core.Compiler.return_type(interp, job.source.f, job.source.tt)
    end
end

#
# @device_code_* functions
#

export @device_code_lowered, @device_code_typed, @device_code_warntype,
       @device_code_llvm, @device_code_metal, @device_code

# forward the rest to GPUCompiler
@eval $(Symbol("@device_code_lowered")) = $(getfield(GPUCompiler, Symbol("@device_code_lowered")))
@eval $(Symbol("@device_code_typed")) = $(getfield(GPUCompiler, Symbol("@device_code_typed")))
@eval $(Symbol("@device_code_warntype")) = $(getfield(GPUCompiler, Symbol("@device_code_warntype")))
@eval $(Symbol("@device_code_llvm")) = $(getfield(GPUCompiler, Symbol("@device_code_llvm")))
@eval $(Symbol("@device_code_metal")) = $(getfield(GPUCompiler, Symbol("@device_code_native")))
@eval $(Symbol("@device_code")) = $(getfield(GPUCompiler, Symbol("@device_code")))

ext: jl


text:
#TODO: write tests for mincubic and maxcubic
seed = rand(Int)
print_with_color(:green, "Testing basic construction")
print_with_color(:blue, "    using seed $seed\n")
function randfloat()
	return 2*(1-rand(Float64))::Float64
end

max_len = 10.0
len = rand(1:max_len)
knots = sort(unique([randfloat()*len::Float64 for i in 1:len]))
len = length(knots)

function randindex()
	return rand(1:len)
end

values = [randfloat() for i in knots]

@assert isa(knots, Array{Float64, 1})
@assert isa(values, Array{Float64, 1})

wrongLength = rand(union(1:len-1, len+1:max_len))
values_wrongLength = [randfloat() for i in 1:wrongLength]
@assert length(values_wrongLength) != length(values)

@test_throws ArgumentError Spline(knots, values_wrongLength)
s = Spline(knots, values)

#Test for refusal to work with NaN values
knotsNaN = copy(knots)
valuesNAN = copy(values)
knotsNaN[rand(1:len)] = NaN
valuesNAN[rand(1:len)] = NaN
@test_throws ArgumentError Spline(knotsNaN, values)
@test_throws ArgumentError Spline(knots, valuesNAN)
@test_throws ArgumentError Spline(knotsNaN, valuesNAN)

#Test for refusal to work with duplicated knots
knotsDup = copy(knots)
dupInd = rand(1:len)
dupInd2 = dupInd
while dupInd2 == dupInd
	dupInd2 = rand(1:len)
end
knotsDup[dupInd] = knotsDup[dupInd2]
@test_throws ArgumentError Spline(knotsDup, values)

domainMin = minimum(knots)
domainMax = maximum(knots)
@test domain(s) == (domainMin, domainMax)

value_indomain = rand(domainMin:domainMax)
value_notindomain = value_indomain + domainMax + rand(Float64)
@assert value_notindomain < domainMin || value_notindomain > domainMax

@test_throws DomainError s(value_notindomain)

for (i, knot) in enumerate(knots)
	@test_approx_eq_eps s(knot) values[i] values[i]/10000
end

function test_values()
	tol = rand(Float64)/10
	disc = uniform_discretize(s, tol)
	mesh = Splines.uniform_discretize_mesh(s, tol)
	for i in eachindex(disc)
		@test s(mesh[i]) ==  disc[i]
	end
end
test_values()

@test_approx_eq s(xmax(s))  maximum(s)

###### Now test the functions that apply to splines #######

##### recalculate!
seed = rand(Int)
print_with_color(:green, "Testing recalculate!")
print_with_color(:blue, "    using seed $seed\n")
values2 = [randfloat() for i in knots]
@test_throws ArgumentError Splines.recalculate!(s, sort(valuesNAN))
@test_throws ArgumentError Splines.recalculate!(s, values_wrongLength)
Splines.recalculate!(s, values2)
for (i, knot) in enumerate(knots)
	@test_approx_eq s(knot) values2[i]
end
test_values()
##### End recalculate! tests

##### insert!
seed = rand(Int)
print_with_color(:green, "Testing insert!")
print_with_color(:blue, "    using seed $seed\n")
# just change value if knot is already there
insert!(s, Float64(randindex()),  rand())
test_values()
#inserting a value actually adds a knot
newknot = Float64(randindex()+randfloat())
newvalue = rand()
insert!(s, newknot, newvalue)
@test newknot in s.knots
@test_approx_eq s(newknot) newvalue

#test that the array version of insert also works
knots_somesize = [Float64(i) for i in rand(3:10)]
values_someothersize = [rand() for i in 11:20]
@test_throws ArgumentError insert!(s, knots_somesize, values_someothersize) #test that arrays of different sizes are rejected
oldknots = copy(s.knots)
oldvalues = copy(s.values)
newknots = [Float64(rand(domain(s)[1]:rand():domain(s)[2])) for i in 1:rand(1:10)]

#make sure that an error is thrown if the new knots array contains duplicates
insert!(newknots, rand(1:length(newknots)), newknots[rand(1:length(newknots))])
@test_throws ArgumentError insert!(s, newknots, [rand() for i in newknots])

#now ensure that some knots are duplicated between the knots already in the spline and the array of new knots to make sure this is handled correctly
for i in 1:rand(2:len)
	insert!(newknots, rand(1:length(newknots)), oldknots[rand(1:len)])
end
newknots = unique(newknots)

newvalues = [rand() for i in newknots]
insert!(s, copy(newknots), copy(newvalues))
#test that the new knots were properly added
@test s.knots == unique(sort(union(oldknots, newknots)))
for (i, knot) in enumerate(newknots)
	@test_approx_eq s(knot) newvalues[i]
end

#test that the previous values are unchanged
for (i, knot) in enumerate(oldknots)
	if !(knot in newknots)
		@test_approx_eq s(knot) oldvalues[i]
	end
end
##### End insert! tests

##### moveknot! tests
seed = rand(Int)
print_with_color(:green, "Testing moveknot!")
print_with_color(:blue, "    using seed $seed\n")
srand(1234)

knots = sort(unique([rand(1.0:1.0:100) for i in 1:10]))
nKnots = length(knots)
values = [rand() for i in knots]
s = Spline(knots, values)
sdomain = domain(s)[1]:rand():domain(s)[2]
@test_throws DomainError moveknot!(s, 1, maximum(knots)+rand()) #can't move outside the domain
@test_throws BoundsError moveknot!(s, -1, rand(sdomain)) #can't move a knot at an out-of-bounds index
@test_throws BoundsError moveknot!(s, nKnots+rand(1:10), rand(sdomain)) #can't move a knot at an out-of-bounds index
@test_throws ArgumentError moveknot!(s, 1, knots[rand(1:nKnots)]) #can't move a knot to where there already is one
knottomove = rand(1:nKnots)
newlocation = knots[knottomove]+rand()
moveknot!(s, knottomove, newlocation)
@test_approx_eq s(newlocation) values[knottomove]

#for the array version of moveknot!
knots = sort(unique([rand(1.0:1.0:100)::Float64 for i in 1:10]))
values = [rand() for i in knots]
s = Spline(knots, values)
newlocations = [(i+rand())::Float64 for i in knots]
newlocations[end] = knots[end]-rand()
moveknot!(s, 1:length(knots), newlocations)
for i in eachindex(newlocations)
	@test_approx_eq s(newlocations[i]) values[i]
end

##### end moveknot! tests


text:
using NeRCA
using Test

const DETX = joinpath(@__DIR__, "data", "detx_v3.detx")

@testset "Calibration" begin
    hits = [Hit(1, 8, 100, 20, false),
            Hit(1, 9, 101, 21, true),
            Hit(2, 8, 112, 22, true),
            Hit(2, 9, 114, 23, false),
            Hit(1, 8, 134, 24, true),
            Hit(2, 9, 156, 25, false),
            Hit(2, 8, 133, 26, true),
            Hit(0, 8, 145, 26, false)]

    calib = Calibration(
        # detector ID
        1,
        # PMT positions
        Dict{Int32, Vector{NeRCA.Position}}(
            8 => [(1,2,3), (4,5,6), (7,8,9)],
            9 => [(10,11,12), (13,14,15),(16,17,18)]
        ),
        # PMT directions
        Dict{Int32, Vector{NeRCA.Direction}}(
            8 => [(1,2,3), (4,5,6), (7,8,9)],
            9 => [(10,11,12), (13,14,15),(16,17,18)]
        ),
        # t0
        Dict{Int32, Vector{Int32}}(
            8 => [1000,2000,3000],
            9 => [4000,5000,6000]
        ),
        # DU
        Dict{Int32, UInt8}(8 => 42, 9 => 23),
        # floor
        Dict{Int32, UInt8}(8 => 18, 9 => 5),
        # OMKeys
        Dict{NeRCA.OMKey, Int32}(OMKey(2, 5) => 123456),  # just dummy data
        # PMTs
        Dict(1 => (2, 5), 2 => (3, 7)),
        # max_z
        10,
        # n_dus
        1
    )

    chits = calibrate(calib, hits)
    @test 42 == chits[1].du
    @test NeRCA.Position(4,5,6) == chits[1].pos
    @test NeRCA.Direction(16,17,18) == chits[6].dir
    @test 5000 == chits[2].t0
    @test 18 == chits[3].floor
    @test 123456 == calib.omkeys[OMKey(2, 5)]

end

@testset "floordist()" begin
    calib = Calibration(DETX)
    @test 1 ≈ floordist(calib)
end


text:
#####################################
## Simple function fitting demo

using PyPlot, PyCall
using Distributions

# true param
W = Array([1.0, 0.0, 1.0])

# generate data
sigma = 0.5
N = 20
X = linspace(-0.4,2.4,N)
Y = [W[1] + W[2]*x + W[3]*x^2 + sigma*randn() for x in X]
X_min = minimum(X)
X_max = maximum(X)

# regression1
X_all = linspace(X_min, X_max, 100)
W1 = sum(Y.*X) / sum(X.^2)
Y1 = [W1*x for x in X_all]

# regression2
X2 = zeros(3, N)
X2[1,:] = 1
X2[2,:] = X
X2[3,:] = X.^2
W2 = inv(X2*X2') * X2*Y
Y2 = [W2[1] + W2[2]*x + W2[3]*x^2 for x in X_all]

# show data
figure()
plot(X_all, Y1, "b-")
plot(X_all, Y2, "g-")
plot(X, Y, "ko")
legend(["model1","model2","data"], loc="upper left", fontsize=16)
xlabel("\$x\$", fontsize=20)
ylabel("\$y\$", fontsize=20)
show()


text:
export  PathDependent1D

mutable struct PathDependent1D
    ρ::Float64 # density
    E::Float64 # Young's modulus
    # hardening parameter, yield function = |σ - q| - (σY + Kα)
    K::Float64
    B::Float64
    σY::Float64
    α::Float64
    α_::Float64 # α to be updated in `commitHistory`
    σ0::Float64 # stress at last time step
    σ0_::Float64 # σ0 to be updated in `commitHistory`
    q::Float64 # stress at last time step
    q_::Float64 # σ0 to be updated in `commitHistory`
    ε0::Float64 # strain at last time step
    ε0_::Float64 # ε0 to be updated in `commitHistory`
end

function PathDependent1D(prop::Dict{String, Any})
    ρ = prop["rho"];
    E = prop["E"]; B = prop["B"];
    K = prop["K"]; σY = prop["sigmaY"]
    α = 0.0; α_ = 0.0
    σ0 = 0.0; σ0_ = 0.0
    q = 0.0; q_ = 0.0
    ε0 = 0.0; ε0_ = 0.0
    PathDependent1D(ρ, E, K, B, σY, α, α_, σ0, σ0_, q, q_, ε0, ε0_)
end

@doc """
    For debugging pathdependent materials, we an build arbitrary pathdependent material law

    sigma = sigma0 + (eps - eps0)**2

""" ->
function getStress(self::PathDependent1D,  strain::Float64,  Dstrain::Float64, Δt::Float64 = 0.0)

    local dΔσdΔε
    ε = strain; ε0 = Dstrain
    σ0 = self.σ0;  α0 = self.α;  q0 = self.q
    E = self.E;    K = self.K;   B = self.B;
    σY = self.σY

    σ = ε^2/20.0
    dΔσdΔε = ε/10.0

    self.α_  = self.α
    self.q_  = self.q

    self.σ0_ = σ
    self.ε0_ = ε

    return σ, dΔσdΔε
end

function getTangent(self::PathDependent1D)
    error("Not implemented")
end

function commitHistory(self::PathDependent1D)
    self.α = self.α_
    self.σ0 = self.σ0_
    self.q = self.q_
    self.ε0 = self.ε0_
end


text:
module BinningAnalysis

import Statistics: mean, var
using Reexport, Lazy, RecursiveArrayTools

# Generic functions
include("generic.jl")
export std_error, tau

# Logarithmic binning
include("log/accumulators.jl")
include("log/binning.jl")
include("log/statistics.jl")
export LogBinner, capacity
export mean, var, varN, tau, std_error, autocorrelation, autocorrelation_time
export all_vars, all_varNs, all_taus, all_std_errors, all_autocorrelations, all_autocorrelation_times
export convergence, has_converged

# "Full" binning
include("full/binning.jl")
include("full/statistics.jl")
export FullBinner

# Jackknife resampling
include("Jackknife.jl")
@reexport using .Jackknife

include("ErrorPropagation/binning.jl")
include("ErrorPropagation/statistics.jl")
export ErrorPropagator
export means, vars, varNs, taus, std_errors, covmat, all_means

include("incrementing/IncrementBinner.jl")
export IncrementBinner
export indices

include("direct_tau.jl")

end # module

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
@testset "EqualSet" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:3] <= 5, Int)
    @constraint(m, x in CS.EqualSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    # doesn't check the length
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 3])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [2, 2, 2])

    constr_indices = constraint.indices
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[2],
        5,
    )
    @test CS.fix!(com, com.search_space[constr_indices[2]], 5)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        4,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1

    # now setting it to 5 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        4,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:3] <= 5, Int)
    @constraint(m, x in CS.EqualSet())
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = com.constraints[1]
    constr_indices = constraint.indices

    # feasible and no changes
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == -5:5
    end
    @test CS.fix!(com, com.search_space[constr_indices[2]], 5)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test CS.value(com.search_space[ind]) == 5
        @test CS.isfixed(com.search_space[ind])
    end

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:3] <= 5, Int)
    @constraint(m, x in CS.EqualSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    constr_indices = constraint.indices

    # Should be synced to the other variables
    @test CS.remove_below!(com, com.search_space[constr_indices[3]], 3)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == 3:5
    end

    @test CS.rm!(com, com.search_space[constr_indices[1]], 5)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == 3:4
    end

    @test CS.rm!(com, com.search_space[constr_indices[2]], 3)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test CS.value(com.search_space[ind]) == 4
        @test CS.isfixed(com.search_space[ind])
    end
end

@testset "equalset is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:5] <= 5, Int)
    @constraint(m, x in CS.EqualSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 3; check_feasibility = false)
    @test CS.fix!(com, variables[2], 5; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:5] <= 5, Int)
    @constraint(m, x in CS.EqualSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 5; check_feasibility = false)
    @test CS.fix!(com, variables[2], 5; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end


text:
import Winston
#import Visualizer.visualize
import HDF5
import Colors

@compat function visualize{T<:AbstractString}(fnames::Array{T,1})
	p = Winston.FramedPlot()
	Winston.setattr(p.x2, "draw_axis", false)
	Winston.setattr(p.y2, "draw_axis", false)
	Winston.setattr(p.y1, "tickdir", +1)
	Winston.setattr(p.x1, "tickdir", +1)
	function func3(p,fnames,i)
		if 0 < i <= length(fnames)
			tf = TemplateFile(fnames[i])
			if tf != nothing
				plot(p, tf)
				_title = replace(fnames[i], "_", "\\_")
				Winston.setattr(p, "title", _title)
			end
		end
	end
	visualize(fnames,800,600,"Spike forms", func3, p)
end

@compat function plot(fname::AbstractString)
	p = Winston.FramedPlot()
	Winston.setattr(p.x2, "draw_axis", false)
	Winston.setattr(p.y2, "draw_axis", false)
	Winston.setattr(p.y1, "tickdir", +1)
	Winston.setattr(p.x1, "tickdir", +1)
	plot(p,fname)
	p
end

@compat function plot(p::Winston.FramedPlot,fname::AbstractString)
	tf = TemplateFile(fname)
	plot(p, tf)
	p
end

function plot(p::Winston.FramedPlot, tf::TemplateFile)
	if tf == nothing
		return p
	end
	colors = Colors.distinguishable_colors(tf.ntemplates)
	t = 1:size(tf.templates,1)
	for i=1:tf.ntemplates
		Winston.add(p, Winston.Curve(t,tf.templates[:,:,i][:];color=colors[i]))
	end
	p
end

function plot(F::Array{Features,1};kvs...)
	p = Winston.FramedPlot()
	plot(p, F;kvs...)
	p
end

function plot(p::Winston.FramedPlot, F::Array{Features,1};feature1::Symbol=:spike_width, feature2::Symbol=:spikes_in_bursts)
	f1 = zeros(length(F))
	f2 = zeros(length(F))
	if !(feature1 in names(F[1]) && feature2 in names(F[1]))
		ArgumentError("Invalid feature name")
	end
	for i in 1:length(f1)
		f1[i] = getfield(F[i],feature1)
		f2[i] = getfield(F[i],feature2)
	end
	#throw away outliers
	pf1 = percentile(f1, [1,99])
	pf2 = percentile(f2, [1,99])
	fidx = (f1.>pf1[1])&(f1.<pf1[2])&(f2.>pf2[1])&(f2.<pf2[2])
	f1 = f1[fidx]
	f2 = f2[fidx]

	Winston.add(p, Winston.Points(f1,f2))
	Winston.setattr(p.x2, "draw_axis",false)
	Winston.setattr(p.y2, "draw_axis",false)
	Winston.setattr(p.x1, "tickdir", 1)
	Winston.setattr(p.y1, "tickdir", 1)
	Winston.setattr(p, "xlabel", replace(string(feature1),"_", " "))
	Winston.setattr(p, "ylabel", replace(string(feature2),"_", " "))
	p
end


text:
#
# Get into the MGS reorthogonalization loop and see if it
# does its job.
#
function mgs_test(cond=1.e6)
(A, x0, b)=data_cook(cond)
V=zeros(3,20)
gout=kl_gmres(x0, b, matvec, V, 1.e-9; orth="mgs1",pdata=A)
gout2=kl_gmres(x0, b, matvec, V, 1.e-9; orth="mgs2",pdata=A)
del=gout.reshist-gout2.reshist
mgs2ok = (norm(del,Inf) > 1.e-12) && gout.idid && gout2.idid
mgs2ok || println("mgs_test fails")
return mgs2ok
#return(gout, gout2, mgs2ok)
end

function matvec(x,A)
return A*x
end

function data_cook(cond)
u1=[1,-2, 0]/sqrt(5.0)
u2=[0,0,1]
u3=[2,1,0]/sqrt(5.0)
U=[u1 u2 u3]
V=[u3 u1 u2]
D=diagm([1, cond, sqrt(cond)])
A=U * D * V'
xstar=ones(3,)
b=A*xstar
x0=[10.0,10.0,10.0]
return(A, x0, b)
end


text:
@testset "is_symmetric" begin
    m = Bool[
        1 0 0
        0 1 0
        0 0 1
    ]
    @test ImageMorphology.is_symmetric(m)
    @test ImageMorphology.is_symmetric(centered(m))
    mo = strel(CartesianIndex, centered(m))
    @test ImageMorphology.is_symmetric(mo)

    m = Bool[
        1 0 1
        0 1 0
        1 1 1
    ]
    mo = strel(CartesianIndex, centered(m))
    @test !ImageMorphology.is_symmetric(m)
    @test !ImageMorphology.is_symmetric(mo)

    @test ImageMorphology.is_symmetric(strel_box((3, 3)))
    @test ImageMorphology.is_symmetric(strel_diamond((3, 3)))
end

@testset "require_symmetric_strel" begin
    se = rand_se_mask(3, 2; symmetric=true)
    @test_nowarn ImageMorphology.require_symmetric_strel(se)
    se = Bool[1 0 0; 0 1 0; 0 0 0]
    msg = "structuring element must be symmetric with respect to its center"
    @test_throws ArgumentError(msg) ImageMorphology.require_symmetric_strel(se)
end

@testset "require_select_function" begin
    @test_nowarn ImageMorphology.require_select_function(max, Int)
    @test_nowarn ImageMorphology.require_select_function(min, Gray{N0f8})

    msg = "function `min` is not a well-defined select function on type `ComplexF64`: does `f(x::T, y::T)` work as expected?"
    @test_throws ArgumentError(msg) ImageMorphology.require_select_function(min, ComplexF64)
    f = (x, y) -> min(x, y)
    msg = "function `$f` is not a well-defined select function on type `RGB{Float32}`: does `f(x::T, y::T)` work as expected?"
    @test_throws ArgumentError ImageMorphology.require_select_function(f, RGB{Float32})
end


text:
import MAT
import Base.parse
using DataFrames
import Stimulus
#import Spiketrains

@compat function loadWaveformsFile(fname::AbstractString;time_conversion::Real=0.001)
	fid = open(fname,"r")
	waveforms = Array(Int16,0)
	timestamps = Array(Uint64,0)
	try
		headersize = Int(read(fid,Uint32))
		num_spikes = Int(read(fid,Uint32))
		num_chs = Int(read(fid,Uint8))
		sampling_rate = Int(read(fid,Uint32))
		timepts = Int(read(fid,Uint32))
		seek(fid,headersize)
		append!(waveforms,read(fid,Int16,num_chs*timepts*num_spikes))
		waveforms = reshape(waveforms,(num_chs,timepts,num_spikes))
		seek(fid,headersize + num_spikes*num_chs*timepts*sizeof(Int16))
		append!(timestamps,read(fid,Uint64,num_spikes))
	finally
		close(fid)
	end
	return float(waveforms),float(timestamps)*time_conversion
end

function assign_area{T,V}(sptrains::Dict{T,V},config;return_mapping=false)
	newsptrains = Dict{T,V}()
	mapping = Dict()
	for (k,v) in sptrains
            kname = assign_area(k,config)
            newsptrains[kname] = v
            mapping[k] = kname
	end
	if return_mapping
		return newsptrains,mapping
	else
		return newsptrains
	end
end

@compat function assign_area(cell::AbstractString, config::Dict)
	rr = r"g([0-9]*)c([0-9])*[s]*"
        m = match(rr,cell)
        group = int(m.captures[1])
        cellnr = int(m.captures[2])
        ad = filter((kk,vv)->group in vv, config)
        aa = collect(keys(ad))[1]
        cc = collect(values(ad))[1]
        cidx = find(cc.==group)[1]
        kname = "$(aa)_g$(cidx)c$(cellnr)"
end

@compat function getspiketrains(;session::AbstractString="",groups::Array{Int64,1}=Array(Int64,0),checkArtifact::Bool=true,verbose::Integer=0,check_trials::Bool=true,config::Dict=Dict())
	if isempty(session)
		#attempt to get session from the pl2 file
		for d in (".","..")
			files = split(chomp(readall(`find $d -name "*.pl2" -maxdepth 1`)))
			if !isempty(files)
				session = splitext(splitdir(files[1])[2])[1]
				break
			end
		end
	end
	verbose > 0 && println("Analyzing session $session ...")
	if isempty(groups)
		#get the groups from the waveforms files in the current directory
		files = split(chomp(readall(`find . -name "*waveforms.bin"`)))
	else
		files = Array(ASCIIString,0)
		for g in groups
			pp = @sprintf "%sg%04dwaveforms.bin" session g
			if isfile(pp)
				push!(files,pp)
			end
		end
	end
    if isempty(config)
        for ss in (".", "..","../..")
            if isfile("$(ss)/channel_config.csv")
                config = channel_config("$(ss)/channel_config.csv")
                break
            end
        end
    end
	SS = pmap(f->begin
					m = match(r"([A-Za-z0-9_]*)g([0-9]*)",f)
					session = m.captures[1]
					group = int(m.captures[2])
					sptrains = getspiketrains(session,group;checkArtifact=checkArtifact,verbose=verbose)
					sptrains
			end,files)
	sptrains = merge(SS...)
	#check if we have a trial structure
    if check_trials
        for ss in (".","..")
            fname = "$(ss)/event_data.mat"
            if isfile(fname)
                verbose > 0 && println("Found trial info. Checking for relevant spikes....")
                trials = Stimulus.loadTrialInfo(fname)
                rtrials = Stimulus.getTrialType(trials,:reward)
                target_time = Stimulus.gettime(rtrials,:target)
                response_time = Stimulus.gettime(rtrials,:response)
                trial_dur = 1000*maximum(response_time .- target_time) #convert from seconds to ms
                cells = Spiketrains.sortCells(sptrains)
                aligned_spikes = Spiketrains.getTrialRaster(sptrains,rtrials,:target,-200.0,trial_dur)
                #remove cells that never spike during the trial
                for cc in sort(setdiff(1:length(cells),unique(aligned_spikes.cellidx)))
                    pop!(sptrains,cells[cc])
                    verbose > 0 && println("\tRemoved cell $(cells[cc])")
                end
                break
            end
        end
    end
	if !isempty(config)
		newsptrains = assign_area(sptrains,config)
	else
		newsptrains = sptrains
	end
	return newsptrains
end

@compat function getspiketrains(session::AbstractString, group::Integer;checkArtifact::Bool=true,verbose::Integer=0)
	verbose > 0 && println("Processing group $group ...")
	wfname = @sprintf "%sg%04dwaveforms.bin" session group
	overlap_fname = @sprintf "%sg%04dwaveforms.overlap" session group
	cut_fname = @sprintf "%sg%04dwaveforms.cut" session group

	waveforms,timestamps = loadWaveformsFile(wfname)
	D = MAT.matread(overlap_fname)
	cids = map(int,readlines(open(cut_fname,"r")))
	clusters = unique(cids)
	clusters = clusters[clusters.>0]
	sptrains = Dict{ASCIIString,Array{Float64,1}}()
	c2 = 1
	for c in clusters
		spidx = (cids.==c) | ((cids.==-1)&(D["overlaps"][:,c].==1))
		cidx = cids.==c
		if checkArtifact
			if size(cids,1) == size(waveforms,3)

				mm = mean(waveforms[:,:,bool(cidx[:])],3)
				ii = indmax(-mm) #get the minimum point
				goodcluster = ii in div(size(mm,2),3):div(2*size(mm,2),3) #minimum should occur in the middle 3rd of the spike
				goodcluster = goodcluster && -minimum(mm) >= maximum(mm) #minimum should exceed maximum in absolute value
				if goodcluster
					cell = "g$(group)c$c2"
					sptrains[cell] = timestamps[spidx[:]]
					c2 += 1
				else
					verbose > 0 && println("\tCluster $c classified as artifact and removed")
				end
			else
				println("\tsize(cids,1) != size(waveforms,3); Skipping cluster $c...")
			end
		else
			cell = "g$(group)c$c"
			sptrains[cell] = timestamps[spidx[:]]
		end
	end
	return sptrains
end

@compat function parse{T<:Integer}(::Type{Range{T}},a::AbstractString)
	m = match(r"([0-9]*)[\-\:]([0-9]*)",a)
     mc = m.captures
    isempty(mc) && return nothing
	ll = map(int,mc)
    return range(ll[1],ll[2]-ll[1]+1)
end

@compat function channel_config(fname::AbstractString)
	dd = DataFrames.readtable(fname)
	#convert to ranges
	q = map(x->parse(Range{Int64},x),dd[:channels])
	dd[:channels] = q
	return Dict(dd[:area],dd[:channels])
end


text:
@testset "Wall Constraint" begin

    # Test WallConstraint
    T = Float64
    n = 4
    P = 5
    x = 4
    y = 2
    X = SVector{n,T}([13.0, 1.0, -12.0, 1.0])
    x1 = SVector{P,T}([ 0.0,  0.0,  1.0,  3.0, -2.0])
    y1 = SVector{P,T}([ 1.0, -1.0,  2.0,  2.0,  0.0])
    x2 = SVector{P,T}([ 1.0,  1.0,  2.0,  2.0,  0.0])
    y2 = SVector{P,T}([ 0.0,  0.0,  1.0,  1.0,  0.0])
    xv = SVector{P,T}([ 1.0,  1.0,  1.0,  1.0,  0.0])./sqrt(2)
    yv = SVector{P,T}([ 1.0, -1.0,  1.0, -1.0,  sqrt(2)])./sqrt(2)

    con = WallConstraint(n,x1,y1,x2,y2,xv,yv,x,y)
    @test norm(TrajectoryOptimization.evaluate(con,X) - [sqrt(2)/2, 0.0, -sqrt(2)/2, 0.0, 0.0], 1) < 1e-10
    @test (@ballocated TrajectoryOptimization.evaluate($con,$X)) == 0

    function easy_jacobian(con,X)
        function local_evaluate(X)
            return TrajectoryOptimization.evaluate(con,X)
        end
        return ForwardDiff.jacobian(local_evaluate, X)
    end
    ∇c_easy = easy_jacobian(con, X)
    ∇c = zeros(MMatrix{P,n,T,P*n})
    jacobian!(∇c,con,X)
    @test norm(∇c - ∇c_easy, 1) < 1e-10
    @test (@ballocated jacobian!($∇c,$con,$X)) == 0

    # Test Wall3DConstraint
    T = Float64
    n = 6
    P = 5
    x = 4
    y = 2
    z = 1
    X  =1
    X0 = SVector{n,T}([0.0, 0.10, -12.0, 0.10, 12.0, 11.0])
    X1 = SVector{n,T}([1.0, 0.55, -12.0, 0.55, 12.0, 11.0])
    X2 = SVector{n,T}([1.0, 1.25, -12.0, 0.75, 12.0, 11.0])
    x1 = SVector{P,T}([ 0.0,  0.0,  0.0,  0.0,  0.0])
    y1 = SVector{P,T}([ 0.0,  0.0,  0.0,  0.0,  0.0])
    z1 = SVector{P,T}([ 0.0,  0.0,  0.0,  0.0,  0.0])

    x2 = SVector{P,T}([ 1.0,  1.0,  1.0,  1.0,  1.0])
    y2 = SVector{P,T}([ 0.0,  0.0,  0.0,  0.0,  0.0])
    z2 = SVector{P,T}([ 0.0,  0.0,  1.0,  0.0,  0.0])

    x3 = SVector{P,T}([ 1.0,  1.0,  1.0,  0.0,  0.0])
    y3 = SVector{P,T}([ 1.0,  1.0,  1.0,  1.0,  1.0])
    z3 = SVector{P,T}([ 0.0,  1.0,  1.0,  0.0,  1.0])

    xv = SVector{P,T}([ 0.0,  0.0,         -1.0/sqrt(2),  0.0,  0.0])
    yv = SVector{P,T}([ 0.0, -1.0/sqrt(2),  0.0,          0.0, -1.0/sqrt(2)])
    zv = SVector{P,T}([ 1.0,  1.0/sqrt(2),  1.0/sqrt(2),  1.0,  1.0/sqrt(2)])

    con = Wall3DConstraint(n,x1,y1,z1,x2,y2,z2,x3,y3,z3,xv,yv,zv,x,y,z)
    @test norm(TrajectoryOptimization.evaluate(con,X0) - [0.0, -0.1/sqrt(2), -0.1/sqrt(2), 0.0, -0.1/sqrt(2)], 1) < 1e-10
    @test norm(TrajectoryOptimization.evaluate(con,X1) - [1.0,  0.45/sqrt(2),  0.45/sqrt(2), 1.0,  0.45/sqrt(2)], 1) < 1e-10
    @test norm(TrajectoryOptimization.evaluate(con,X2) - [0.0,  0.0,  0.0,  1.0,  -0.25/sqrt(2)], 1) < 1e-10
    @test (@ballocated TrajectoryOptimization.evaluate($con,$X0)) == 0

    ∇c_easy = easy_jacobian(con, X0)
    ∇c = zeros(MMatrix{P,n,T,P*n})
    jacobian!(∇c,con,X0)
    @test norm(∇c - ∇c_easy, 1) < 1e-10
    @test (@ballocated jacobian!($∇c,$con,$X0)) == 0
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
end


text:
using DataFrames, Gadfly

set_default_plot_size(2*3.3inch, 2*3.3inch)

D = DataFrame(x1=[76.5763, 93.611, 67.2219, 39.6499],
              x2=[31.5112, 66.0929, 58.5162, 5.66118],
              x3=[27.0953, 11.2782, 16.7029, 47.3287],
              x4=[86.1758, 61.6317, 28.7841, 46.4209])

xsc  = Scale.x_continuous(minvalue=0.0, maxvalue=100)
ysc  = Scale.y_continuous(minvalue=0.0, maxvalue=100)

p1 = plot(D, x=:x1, y=:x2, xend=:x3, yend=:x4, Geom.segment(arrow=true), xsc, ysc);
p2 = plot(D, x=:x1, y=:x2, xend=:x3, yend=:x4, color=[colorant"orange"], linestyle=[:dash], Geom.vector, xsc, ysc);

p3 = plot(z=(x,y)->x*exp(-(x^2+y^2)),
        xmin=[-2], xmax=[2], ymin=[-2], ymax=[2],
        Geom.vectorfield(scale=0.4, samples=17),
        Scale.x_continuous(minvalue=-2.0, maxvalue=2.0),
        Scale.y_continuous(minvalue=-2.0, maxvalue=2.0),
        Guide.xlabel("x"), Guide.ylabel("y"), Guide.colorkey(title="z"));

Z = fill(1.0, 5, 5)
Z[2:4, 2:4] .= 2.0
Z[3,3] = 3.0

p4 = plot(z=Z,
        Geom.vectorfield,
        Scale.x_continuous(minvalue=1.0, maxvalue=5.0),
        Scale.y_continuous(minvalue=1.0, maxvalue=5.0),
        Guide.xlabel("x"), Guide.ylabel("y"),
        Theme(key_position=:none));

gridstack([p1 p2; p3 p4])


text:
#=
You run an e-commerce website and want to record the last N order ids in a log. Implement a data structure to accomplish this, with the following API:

    record(order_id): adds the order_id to the log
    get_last(i): gets the ith last element from the log. i is guaranteed to be smaller than or equal to N.

You should be as efficient with time and space as possible.
=#

global log = []
global N = 5

function changeN(N_new::Int64)
    global N = N_new
end

function record(order_id)
    if length(log) +1 <= N
        append!(log, order_id)
    else
        popfirst!(log)
        append!(log, order_id)
    end
end

function get_last(i::Int64)
    if i <= length(log)
        return log[end - i + 1]
    else
        return "i is too large"
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
include("conv.jl")
using Functors
using UnPack

struct Grid
    cell::AbstractMatrix
    origin::Any
    r::AbstractArray
    rhat::AbstractArray
    R::AbstractArray
    dv::AbstractFloat
end
@functor Grid

Base.size(g::Grid) = size(g.r)
function Base.getproperty(g::Grid, f::Symbol)
    if f === :x
        return getindex.(g.r, 1)
    elseif f === :y
        return getindex.(g.r, 2)
    elseif f === :z
        return getindex.(g.r, 3)
    end
    getfield(g, f)
end

"""
    Grid(cell, rmax::AbstractFloat)
    Grid(
        cell::AbstractMatrix,
        sz::Union{AbstractVector,Tuple};
        origin = (sz .+ 1) ./ 2,
        )

Grid is specified by its discrete cell vectors (column-wise matrix), overall size and origin. For a uniform Cartesian 5x5 grid discretized at 0.1 with a centered origin, we get `cell = [0.1 0; 0 0.1]` & `origin = [3, 3]`. Grid cell can in general be noncartesian.
"""
function Grid(
    cell::AbstractMatrix,
    sz::Union{AbstractVector,Tuple};
    origin = (sz .+ 1) ./ 2,
)
    n = size(cell, 1)
    if n == 1
        rvecs = [cell * ([x] .- origin) for x = 1:sz[1]]
    elseif n == 2

        rvecs = [cell * ([x, y] .- origin) for x = 1:sz[1], y = 1:sz[2]]
    elseif n == 3

        rvecs = [
            cell * ([x, y, z] .- origin)
            for x = 1:sz[1], y = 1:sz[2], z = 1:sz[3]
        ]
    end
    R = norm.(rvecs)
    rhat = rvecs ./ (R .+ 1e-16)
    dv = det(cell)
    Grid(cell, origin, rvecs, rhat, R, dv)
end

function Grid(cell, rmax::AbstractFloat)
    n = size(cell, 1)
    Grid(cell, 1 .+ 2 * ceil.(rmax * (cell \ ones(n))))
end

struct Props

    nranks::Any
    ranks::Any
    parities::Any
    slices::Any
end

function Props(ranks)
    Props(; ranks)
end
function Props(; nranks = nothing, ranks = nothing)
    if ranks === nothing
        ranks = vcat([fill(i - 1, n) for (i, n) in enumerate(nranks)]...)
    elseif nranks === nothing
        nranks = zeros(Int, maximum(ranks) + 1)
        for x in ranks
            nranks[x+1] += 1
        end
    end
    ends = Int.(cumsum(2 .* ranks .+ 1))
    slices = [(i == 1 ? 1 : ends[i-1] + 1):ends[i] for i in eachindex(ends)]

    # starts = [1, (1 .+ cumsum(nranks))...]
    # grouped_slices = [slices[starts[i]:starts[i+1]-1] for i = 1:length(nranks)]
    parities = (-1) .^ ranks
    Props(nranks, ranks, parities, slices)
end

"""
    (field::AbstractArray)(rvec, grid::Grid)

Makes field objects callable. Returns interpolated tensor value of the tensor field at position rvec.
"""
function Base.get(fields::AbstractArray, props::Props, i::Int)
    fields[:, :, :, props.slices[i]]
end

"""
    Base.get(field::AbstractArray, grid::Grid, rvec::AbstractVector)
    Base.put!(
        field::AbstractArray,
        grid::Grid,
        rvec::AbstractVector,
        val::AbstractVector,
    )

With grid info we can interpolate a scalar or vector field at any location. We can also place a scalar or vector point source anywhere with automatic normalization wrt discretization. Both work via a proximity weighted average of the closest grid points (in general up to 4 in 2d and 8 in 3d).
"""
function Base.get(field::AbstractArray, grid::Grid, rvec::AbstractVector)
    sum([w * field[ix...] for (ix, w) in nearest(grid, rvec)])
end
function Base.get(
    fields::AbstractArray,
    props::Props,
    i::Int,
    rvec::AbstractVector,
)
    field = get(fields, props, i)
    get(field, props.grid, rvec)
end
function Base.get(fields::AbstractArray, props::Props, rvec::AbstractVector)
    get(fields, props, 1, rvec)
end

function nearest(grid, rvec)
    @unpack cell, origin = grid
    n = length(size(grid))
    ix = origin .+ (cell \ rvec)
    ixfloor = floor.(Int, ix)
    er = ix - ixfloor

    if n == 2
        res = [
            (ixfloor + [x, y], prod(ones(n) - abs.([x, y] - er)))
            for
            x = 0:1, y = 0:1 if [0, 0] < ixfloor + [x, y] < collect(size(grid))
        ]
    elseif n == 3
        res = [
            (ixfloor + [x, y, z], prod(ones(3) - abs.([x, y, z] - er)))
            for
            x = 0:1,
            y = 0:1,
            z = 0:1 if [0, 0, 0] < ixfloor + [x, y, z] < collect(size(grid))
        ]
    end
    res
end

function Base.put!(field::AbstractArray, grid::Grid, rvec::AbstractVector, val)
    for (ix, w) in nearest(grid, rvec)
        field[ix...] += w / grid.dv * val
    end
end

function Base.put!(f, grid, rvecs::AbstractMatrix, vals::AbstractMatrix)
    for (val, rvec) in zip(eachcol(vals), eachcol(rvecs))
        put!(f, grid, rvec, val)
    end
end


text:
# Test suite for gp.jl
@testset "Cholesky LDLT" for (J0, J, N_test) in ((1, 1, [64,1024]), (0, 16, [32,512]), (4,9,[148933]))
    # J0 = Number of real, exponential celerite kernel terms
    # J  = Total number of terms

    # Iterate until we have a positive definite kernel (defined by Sturm's theorem):
    num_pos_root = 1
    aj = rand(J)
    bj = [zeros(Float64,J0);rand(J-J0)]
    cj = rand(Float64,J)
    dj = [zeros(Float64,J0);rand(J-J0)]
    while num_pos_root > 0
        x = zeros(Float64,0)
        for j=1:J
            push!(x,aj[j],bj[j],cj[j],dj[j])
            num_pos_root = celerite.sturms_theorem(x)
        end
        if num_pos_root > 0
            aj = rand(J)
            bj = [zeros(Float64,J0);rand(J-J0)]
            cj = rand(J)
            dj = [zeros(Float64,J0);rand(J-J0)]
        end
    end

    # Run a timing test:
    ntest = length(N_test)
    time_complex = zeros(Float64,ntest)
    itest = 1
    time_prior = 0.0
    #while itest <= ntest && time_prior < 3.0
    ntrial = 2
    #while itest <= ntest
    N = N_test[itest]
    yerr = 0.1 +zeros(N)
    # Generate some random time data:
    t = sort(rand(N)).*100
#    y0 = sin.(t)
    i=1
    if J0 > 0
        kernel = celerite.RealTerm(log(aj[i]),log(cj[i]))
        while i < J0
            i +=1
            kernel = kernel + celerite.RealTerm(log(aj[i]),log(cj[i]))
        end
    end
    if i < J
        if J0 == 0
            kernel = celerite.ComplexTerm(log(aj[i]),log(bj[i]),log(cj[i]),log(dj[i]))
        end
        while i < J
            i +=1
            kernel = kernel + celerite.ComplexTerm(log(aj[i]),log(bj[i]),log(cj[i]),log(dj[i]))
        end
    end
    gp = celerite.Celerite(kernel)
    # The following was used for profiling the code:
    if itest == 9
        logdet_test = celerite.compute_ldlt!(gp, t, yerr)
        @profile logdet_test = celerite.compute_ldlt!(gp, t, yerr)
        Profile.print(format=:flat)
        #  else
    end
    # Cholesky method
    # Compute Cholesky factorization, return log of determinant:
    logdet_test = celerite.compute_ldlt!(gp, t, yerr)
    # Generate random noise, and a realization of the GP:
    noise = randn(N)
    # This simulates a Gaussian process with "gp" kernel:
    y0 = celerite.simulate_gp_ldlt(gp,noise)
    # Check that the inverse works:
    noise_test = celerite.invert_lower_ldlt(gp,y0)./sqrt.(gp.D)
#    println("noise recovered? ",maximum(abs(noise-noise_test)))
    @test isapprox(noise,noise_test)
    time_zero = tic()
    # Run a somme trials to get a timing estimate:
    for itrial = 1:ntrial
        logdet_test = celerite.compute_ldlt!(gp, t, yerr)
    end
    #  end
    time_complex[itest]=toc()/ntrial
    # Now do full solve (if N_test isn't too big):
    if N < 2000
        logdetK,K = celerite.full_solve(t,y0,aj,bj,cj,dj,yerr)
# Check the log determinant:
#        println("Log det: ",logdetK,logdet_test)
        @test isapprox(logdetK,logdet_test)
# Check that simulation works:
        K_lower = chol(K)'  # chol(K) returns upper triangular cholesky matrix, so we need to transpose.
        y0_full = *(K_lower,noise)
# Check factorization by reconstructing cholesky factor from low-rank decomposition:
        Umat = copy(gp.up)
        Wmat = copy(gp.W)
        for n=1:N
          if J0 > 0
            for j=1:J0
              Umat[j,n] *= exp(-cj[j]*t[n])
              Wmat[j,n] *= exp( cj[j]*t[n])
            end
          end
          if (J-J0) > 0
            for j=1:J-J0
              Umat[J0+2j-1,n] *= exp(-cj[j]*t[n])
              Umat[J0+2j  ,n] *= exp(-cj[j]*t[n])
              Wmat[J0+2j-1,n] *= exp( cj[j]*t[n])
              Wmat[J0+2j  ,n] *= exp( cj[j]*t[n])
            end
          end
        end
        L = tril(*(Umat', Wmat), -1)
        # Need to add in identity and then multipy by D^{1/2}:
        for i=1:N
          L[i,i] = 1.0
          for j=1:N
            L[i,j] *=sqrt(gp.D[j])
          end
        end
        # Check that lower cholesky matrix agrees with numerical cholesky decomposition of full matrix:
        @test isapprox(L, K_lower)
#        println("Cholesky error: ",maximum(abs(L - K_lower)))
        # Check that reconstructed matrix agrees with full covariance matrix:
#        println("Cholesky error: ",maximum(abs(*(L, L') - K)))
        @test isapprox(*(L, L'), K)
#        println("y0: ",maximum(abs(y0-y0_full)))
        # Check that simulated vectors agree for low-rank and full cholesky multiplication:
        @test isapprox(y0,y0_full)
# Check that the solver works:
        z = celerite.apply_inverse_ldlt(gp, y0)
        z_full = K \ y0
#        println("K \\ y0: ",maximum(abs(z_full-z)))
        @test isapprox(z_full,z)
# Check that multiplication works:
#        y_test = *(K,z)
        y_test = celerite.multiply_ldlt(gp, t, z, yerr)
#        println("Multiplication: ",maximum(abs(y0-y_test)))
        @test isapprox(y_test,y0)
# Check that the "chi-square" gives the correct value:
#        println("N: ",N," dot(y0,z): ",dot(y0,z)," dot(noise,noise): ",dot(noise,noise))
        @test isapprox(dot(y0,z),dot(noise,noise))
# Check that the log likelihood is computed correctly:
        nll = celerite.log_likelihood_ldlt(gp, y0)
        nll_full = -.5*(logdetK+N*log(2pi)+dot(y0,z))
#        println("Log likelihoods: ",nll," ",nll_full)
        @test isapprox(nll,nll_full)
        println(N_test[itest]," ",time_complex[itest])
        time_prior = time_complex[itest]
        M = N*4
        tpred = sort!(rand(M)) .* 200
        ypred = celerite.predict_ldlt!(gp, t, y0, tpred)
        ypred_full = celerite.predict_full_ldlt(gp, y0, tpred; return_cov = false)
    #end
        @test isapprox(ypred,ypred_full)
    end
    println("J0,J: ",J0," ",J," N: ",N)
end
return


text:
#jl #! format: off
# # Operations problems with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports the construction and solution of optimal power system
# scheduling problems (Operations Problems). Operations problems form the fundamental
# building blocks for [sequential simulations](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/02_sequential_simulations.ipynb).
# This example shows how to specify and customize a the mathematics that will be applied to the data with
# an `OperationsProblemTemplate`, build and execute an `OperationsProblem`, and access the results.

# ## Dependencies
using SIIPExamples

# ### Modeling Packages
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder

# ### Data management packages
using Dates
using DataFrames

# ### Optimization packages
using Cbc #solver

# ### Data
# This data depends upon the [RTS-GMLC](https://github.com/gridmod/rts-gmlc) dataset. Let's
# use [PowerSystemCaseBuilder.jl](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/10_PowerSystemCaseBuilder.ipynb) to download and build a `System`.

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

# ## Define a problem specification with an `OperationsProblemTemplate`
# You can create an empty template with:
template_uc = OperationsProblemTemplate()

# Now, you can add a `DeviceModel` for each device type to create an assignment between PowerSystems device types
# and the subtypes of `AbstractDeviceFormulation`. PowerSimulations has a variety of different
# `AbstractDeviceFormulation` subtypes that can be applied to different PowerSystems device types,
# each dispatching to different methods for populating optimization problem objectives, variables,
# and constraints.

print_tree(PSI.AbstractDeviceFormulation)

# ### Branch Formulations
# Here is an example of relatively standard branch formulations. Other formulations allow
# for selective enforcement of transmission limits and greater control on transformer settings.
set_device_model!(template_uc, Line, StaticBranch)
set_device_model!(template_uc, Transformer2W, StaticBranch)
set_device_model!(template_uc, TapTransformer, StaticBranch)

# ### Injection Device Formulations
# Here we define template entries for all devices that inject or withdraw power on the
# network. For each device type, we can define a distinct `AbstractDeviceFormulation`. In
# this case, we're defining a basic unit commitment model for thermal generators,
# curtailable renewable generators, and fixed dispatch (net-load reduction) formulations
# for `HydroDispatch` and `RenewableFix` devices.

set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, HydroDispatch, FixedOutput)
set_device_model!(template_uc, HydroEnergyReservoir, HydroDispatchRunOfRiver)
set_device_model!(template_uc, RenewableFix, FixedOutput)

# ### Service Formulations
# We have two `VariableReserve` types, parameterized by their direction. So, similar to
# creating `DeviceModel`s, we can create `ServiceModel`s. The primary difference being
# that `DeviceModel` objects define how constraints get created, while `ServiceModel` objects
# define how constraints get modified.
set_service_model!(template_uc, VariableReserve{ReserveUp}, RangeReserve)
set_service_model!(template_uc, VariableReserve{ReserveDown}, RangeReserve)

# ### Network Formulations
# Finally, we can define the transmission network specification that we'd like to model. For simplicity, we'll
# choose a copper plate formulation. But there are dozens of specifications available through
# an integration with [PowerModels.jl](https://github.com/lanl-ansi/powermodels.jl). *Note that
# many formulations will require appropriate data and may be computationally intractable*
set_transmission_model!(template_uc, CopperPlatePowerModel)

# ## `OperationsProblem`
# Now that we have a `System` and an `OperationsProblemTemplate`, we can put the two together
# to create an `OperationsProblem` that we solve.

# ### Optimizer
# It's most convenient to define an optimizer instance upfront and pass it into the
# `OperationsProblem` constructor. For this example, we can use the free Cbc solver with a
# relatively relaxed MIP gap (`ratioGap`) setting to improve speed.

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

# ### Build an `OperationsProblem`
# The construction of an `OperationsProblem` essentially applies an `OperationsProblemTemplate`
# to `System` data to create a JuMP model.

op_problem = OperationsProblem(template_uc, sys; optimizer = solver, horizon = 24)

build!(op_problem, output_dir = mktempdir())

#nb # The principal component of the `OperationsProblem` is the JuMP model. For small problems,
#nb # you can inspect it by simply printing it to the screen:
#nb # ```julia
#nb # op_problem.internal.optimization_container.JuMPmodel
#nb # ```
#nb #
#nb # For anything of reasonable size, that will be unmanageable. But you can print to a file:
#nb # ```julia
#nb # f = open("testmodel.txt","w"); print(f,op_problem.internal.optimization_container.JuMPmodel); close(f)
#nb # ```
#nb #
#nb # In addition to the JuMP model, an `OperationsProblem` keeps track of a bunch of metadata
#nb # about the problem and some references to pretty names for constraints and variables.
#nb # All of these details are contained within the `optimization_container` field.
#
#nb print_struct(typeof(op_problem.internal.optimization_container))
#
#nb # ### Solve an `OperationsProblem`
#
#nb solve!(op_problem)
#
#nb # ## Results Inspection
#nb # PowerSimulations collects the `OperationsProblem` results into a struct:
#
#nb print_struct(PSI.ProblemResults)
#
#nb res = ProblemResults(op_problem);
#
#nb # ### Optimizer Stats
#nb # The optimizer summary is included
#
#nb get_optimizer_stats(res)
#
#nb # ### Objective Function Value
#
#nb get_objective_value(res)
#
#nb # ### Variable Values
#nb # The solution value data frames for variables can be accessed by:
#
#nb variable_values = get_variables(res)
#
# ## Plotting
# Take a look at the examples in [the plotting folder.](../../notebook/3_PowerSimulations_examples/Plotting)


text:
include("./server.jl")

ENV["JULIA_MARGO_LOAD_PYPLOT"] = "no thank you"
import ClimateMARGO
using ClimateMARGO.Models
using ClimateMARGO.Optimization
using ClimateMARGO.Diagnostics

# some silly code to make our custom types work with MsgPack (JSON alternative)
import MsgPack
MsgPack.msgpack_type(::Type{ClimateModel}) = MsgPack.StructType()
MsgPack.msgpack_type(::Type{ClimateModelParameters}) = MsgPack.StructType()
MsgPack.msgpack_type(::Type{Domain}) = MsgPack.StructType()
MsgPack.msgpack_type(::Type{Economics}) = MsgPack.StructType()
MsgPack.msgpack_type(::Type{Physics}) = MsgPack.StructType()
MsgPack.msgpack_type(::Type{Controls}) = MsgPack.StructType()

function setfieldconvert!(value, name::Symbol, x)
    setfield!(value, name, convert(typeof(getfield(value, name)), x))
end
# the main margo function
@expose function opt_controls_temp(;dt=20, opt_parameters, economics=Dict(), physics=Dict())
    model_parameters = deepcopy(ClimateMARGO.IO.included_configurations["default"])::ClimateModelParameters
    model_parameters.domain = Domain(Float64(dt), 2020.0, 2200.0)
    model_parameters.economics.baseline_emissions = ramp_emissions(model_parameters.domain)
    model_parameters.economics.extra_CO₂ = zeros(size(model_parameters.economics.baseline_emissions))

    for (k, v) in economics
        setfieldconvert!(model_parameters.economics, Symbol(k), v)
    end
    for (k, v) in physics
        setfieldconvert!(model_parameters.physics, Symbol(k), v)
    end

    model = ClimateModel(model_parameters)

    parsed = Dict((Symbol(k) => v) for (k, v) in opt_parameters)
    model_optimizer = optimize_controls!(model; parsed..., print_raw_status=false)
    return Dict(
        :model_parameters => model_parameters,
        :status => ClimateMARGO.Optimization.JuMP.termination_status(model_optimizer) |> string,
        model_results(model)...
    )
end

@expose function forward_controls_temp(;dt=20, controls=Dict(), economics=Dict(), physics=Dict())
    model_parameters = deepcopy(ClimateMARGO.IO.included_configurations["default"])::ClimateModelParameters
    model_parameters.domain = Domain(Float64(dt), 2020.0, 2200.0)
    model_parameters.economics.baseline_emissions = ramp_emissions(model_parameters.domain)
    model_parameters.economics.extra_CO₂ = zeros(size(model_parameters.economics.baseline_emissions))

    for (k, v) in economics
        setfieldconvert!(model_parameters.economics, Symbol(k), v)
    end
    for (k, v) in physics
        setfieldconvert!(model_parameters.physics, Symbol(k), v)
    end

    model = ClimateModel(model_parameters)

    translations = Dict(
        :M => :mitigate,
        :R => :remove,
        :G => :geoeng,
        :A => :adapt,
    )
    for (k, v) in controls
        setfieldconvert!(model.controls, translations[Symbol(k)], v)
    end

    return Dict(
        :model_parameters => model_parameters,
        model_results(model)...
    )
end

model_results(model::ClimateModel) = Dict(
    :controls => model.controls,
    :computed => Dict(
        :temperatures => Dict(
            :baseline => T_adapt(model),
            :M => T_adapt(model; M=true),
            :MR => T_adapt(model; M=true, R=true),
            :MRG => T_adapt(model; M=true, R=true, G=true),
            :MRGA => T_adapt(model; M=true, R=true, G=true, A=true),
        ),
        :emissions => Dict(
            :baseline => effective_emissions(model),
            :M => effective_emissions(model; M=true),
            :MRGA => effective_emissions(model; M=true, R=true),
        ),
        :concentrations => Dict(
            :baseline => c(model),
            :M => c(model; M=true),
            :MRGA => c(model; M=true, R=true),
        ),
        :damages => Dict(
            :baseline => costs_dict(damage(model; discounting=true), model),
            :MRGA => costs_dict(damage(model; M=true, R=true, G=true, A=true, discounting=true), model),
        ),
        :costs => Dict(
            :M => costs_dict(cost(model; M=true, discounting=true), model),
            :R => costs_dict(cost(model; R=true, discounting=true), model),
            :G => costs_dict(cost(model; G=true, discounting=true), model),
            :A => costs_dict(cost(model; A=true, discounting=true), model),
            :MRGA => costs_dict(cost(model; M=true, R=true, G=true, A=true, discounting=true), model),
        ),
    ),
)

function costs_dict(costs, model)
    Dict(
        :discounted => costs,
        :total_discounted => sum(costs .* model.domain.dt),
    )
end

# simple function as baseline for connection speed
@expose function arithmetic(;a, b)
    Dict(
        :sum => a + b,
        :difference => a - b,
    )
end

@expose function randimage(;width, height)
    img = rand(UInt8, width * height * 4)
    Dict(
        :img => img,
    )
end


text:
import ViscousFlow: Fields, TimeMarching
import ViscousFlow.TimeMarching:RK31

using LinearAlgebra

@testset "Time Marching" begin

  @testset "Basic Unconstrained" begin

  ω = 4
  u₀ = 1.0
  uex(t) = u₀ + sin(ω*t)/ω

  Δt = 0.005
  T = 0:Δt:10
  u = [u₀]
  TimeMarching.r₁(u::Vector{Float64},t::Float64) = cos(ω*t)
  rk = RK(u,Δt,r₁,rk=TimeMarching.RK31)

  u = [u₀]
  uhist = Float64[]
  for t in T
    push!(uhist,u[1])
    t,u = rk(t,u)
  end

  @test norm(uhist-uex.(T)) ≈ 0.000494857145025

  end

  @testset "Unconstrained with integrating factor" begin

  α = 0.5
  ω = 4
  u₀ = 1.0
  uex(t) = u₀*exp(-α*t) + (α*(cos(ω*t)-exp(-α*t))+ω*sin(ω*t))/(α^2+ω^2)

  Fields.plan_intfact(t::Float64,u::Vector{Float64}) = exp(-α*t)

  Δt = 0.005
  T = 0:Δt:10
  u = [u₀]
  TimeMarching.r₁(u::Vector{Float64},t::Float64) = cos(ω*t)
  ifrk = IFRK(u,Δt,plan_intfact,r₁,rk=TimeMarching.RK31)

  u = [u₀]
  uhist = Float64[]
  for t in T
    push!(uhist,u[1])
    t,u = ifrk(t,u)
  end
  @test norm(uhist-uex.(T)) ≈ 0.0005014844449

  end

  @testset "Basic unconstrained using IFHERK" begin

  ω = 4
  u₀ = 1.0
  uex(t) = u₀ + sin(ω*t)/ω

  Δt = 0.005
  T = 0:Δt:10
  u = [u₀]
  f = Vector{Float64}()
  TimeMarching.r₁(u::Vector{Float64},t::Float64) = cos(ω*t)
  TimeMarching.r₂(u::Vector{Float64},t::Float64) = Vector{Float64}()
  TimeMarching.plan_constraints(u::Vector{Float64},t::Float64) = f -> zeros(Float64,1), u -> Vector{Float64}()
  Fields.plan_intfact(t::Float64,u::Vector{Float64}) = Matrix(1.0I,1,1)
  ifherk = IFHERK(u,f,Δt,plan_intfact,plan_constraints,(r₁,r₂),rk=TimeMarching.RK31)

  u = [u₀]
  uhist = Float64[]
  for t in T
    push!(uhist,u[1])
    t,u,_ = ifherk(t,u)
  end

  @test norm(uhist-uex.(T)) ≈ 0.0004948571450253

  end

end


text:
function make_column(array)
   if length(size(array)) == 1
      outArray = array.';
   else
      error("Must be array.");
   end
   return outArray
end


text:
"""
    u_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                        closure, U, C, K, F, pHY′, clock)

Return the tendency for the horizontal velocity in the x-direction, or the east-west
direction, ``u``, at grid point `i, j, k`.

The tendency for ``u`` is called ``G_u`` and defined via

    ``∂_t u = G_u - ∂_x ϕ_n``

where ∂_x ϕ_n is the non-hydrostatic pressure gradient in the x-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

The arguments `U`, `C`, and `K` are `NamedTuple`s with the three velocity components,
tracer fields, and precalculated diffusivities where applicable. `F` is a named tuple of
forcing functions, `pHY′` is the hydrostatic pressure anomaly.

`parameters` is a `NamedTuple` of scalar parameters for user-defined forcing functions
and `clock` is the physical clock of the model.
"""
@inline function u_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                     closure, U, C, K, F, pHY′, clock)

    return ( - div_ũu(i, j, k, grid, U)
             - x_f_cross_U(i, j, k, grid, coriolis, U)
             - ∂xᶠᵃᵃ(i, j, k, grid, pHY′)
             + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, U, K)
             + x_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, clock.time)
             + ∂t_uˢ(i, j, k, grid, surface_waves, clock.time)
             + F.u(i, j, k, grid, clock, (velocities=U, tracers=C, diffusivities=K)))
end

"""
    v_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                        closure, U, C, K, F, pHY′, clock)

Return the tendency for the horizontal velocity in the y-direction, or the north-south
direction, ``v``, at grid point `i, j, k`.

The tendency for ``v`` is called ``G_v`` and defined via

    ``∂_t v = G_v - ∂_y ϕ_n``

where ∂_y ϕ_n is the non-hydrostatic pressure gradient in the y-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

The arguments `U`, `C`, and `K` are `NamedTuple`s with the three velocity components,
tracer fields, and precalculated diffusivities where applicable. `F` is a named tuple of
forcing functions, `pHY′` is the hydrostatic pressure anomaly.

`parameters` is a `NamedTuple` of scalar parameters for user-defined forcing functions
and `clock` is the physical clock of the model.
"""
@inline function v_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                     closure, U, C, K, F, pHY′, clock)

    return ( - div_ũv(i, j, k, grid, U)
             - y_f_cross_U(i, j, k, grid, coriolis, U)
             - ∂yᵃᶠᵃ(i, j, k, grid, pHY′)
             + ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure, U, K)
             + y_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, clock.time)
             + ∂t_vˢ(i, j, k, grid, surface_waves, clock.time)
             + F.v(i, j, k, grid, clock, (velocities=U, tracers=C, diffusivities=K)))
end

"""
    w_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                        closure, U, C, K, F, clock)

Return the tendency for the vertical velocity ``w`` at grid point `i, j, k`.
The tendency for ``w`` is called ``G_w`` and defined via

    ``∂_t w = G_w - ∂_z ϕ_n``

where ∂_z ϕ_n is the non-hydrostatic pressure gradient in the z-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

The arguments `U`, `C`, and `K` are `NamedTuple`s with the three velocity components,
tracer fields, and precalculated diffusivities where applicable. `F` is a named tuple of
forcing functions, `pHY′` is the hydrostatic pressure anomaly.

`parameters` is a `NamedTuple` of scalar parameters for user-defined forcing functions
and `clock` is the physical clock of the model.
"""
@inline function w_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                     closure, U, C, K, F, clock)

    return ( - div_ũw(i, j, k, grid, U)
             - z_f_cross_U(i, j, k, grid, coriolis, U)
             + ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure, U, K)
             + z_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, clock.time)
             + ∂t_wˢ(i, j, k, grid, surface_waves, clock.time)
             + F.w(i, j, k, grid, clock, (velocities=U, tracers=C, diffusivities=K)))
end

"""
    tracer_tendency(i, j, k, grid, c, tracer_index,
                    closure, buoyancy, U, C, K, Fc, clock)

Return the tendency for a tracer field `c` with index `tracer_index`
at grid point `i, j, k`.

The tendency for ``c`` is called ``G_c`` and defined via

    ``∂_t c = G_c``

`closure` and `buoyancy` are types encoding information about the prescribed
turbulence closure and buoyancy model.

The arguments `U`, `C`, and `K` are `NamedTuple`s with the three velocity components,
tracer fields, and  precalculated diffusivities where applicable.
`Fc` is the user-defined forcing function for tracer `c`.

`clock` keeps track of `clock.time` and `clock.iteration`.
"""
@inline function tracer_tendency(i, j, k, grid, c, tracer_index,
                                 closure, buoyancy, U, C, K, Fc, clock)

    return ( - div_uc(i, j, k, grid, U, c)
             + ∇_κ_∇c(i, j, k, grid, closure, c, tracer_index, K, C, buoyancy)
             + Fc(i, j, k, grid, clock, (velocities=U, tracers=C, diffusivities=K)))
end


text:
module Checks

import ..GitHub: @ghdef, GitHubType, namefield
export Action, Image, Annotation, Output

@ghdef mutable struct Action
    label::String
    description::String
    identifier::String
end
namefield(act::Action) = act.label

@ghdef mutable struct Image
    alt::String
    image_url::String
    caption::Union{String, Nothing}
end
namefield(img::Image) = img.alt

@ghdef mutable struct Annotation
    path::String
    start_line::Int
    end_line::Int
    start_column::Union{Int, Nothing}
    end_column::Union{Int, Nothing}
    annotation_level::String
    message::String
    title::Union{String, Nothing}
    raw_details::Union{String, Nothing}
end
namefield(ann::Annotation) = ann.filename

@ghdef mutable struct Output
    title::String
    summary::String
    text::Union{String, Nothing}
    annotations::Union{Vector{Annotation}, Nothing}
    images::Union{Vector{Image}, Nothing}
end
namefield(out::Output) = out.title

end
using .Checks


text:
using Oxygen

println("Testing Electrons")

@testset "Electrons" begin
    carbon = Oxygen.parse_config("[He]2s2 2p2")
    @test typeof(carbon.root) == String
    @test typeof(carbon.orbitals) <: Tuple
    @test carbon == Oxygen.Electrons(6)
    @test carbon == Oxygen.Electrons("C")
end

@testset "Printing" begin
    carbon = Oxygen.Electrons(6)
    @test repr(carbon) == "[He] ((\"2s\", 2), (\"2p\", 2))"
end


text:
using Merlin
using Base.Threads
using LibCUDA

function bench()
    results = Array{Any}(2)
    @threads for dev = 0:1
        T = Float32
        x = rand(T,100,100)
        f = Linear(T,100,100)
        y = f(Var(x))
        results[dev+1] = Array(y.data)
    end
end

@time bench()
#@threads for i = 1:2
    # a[i] = threadid()
#    setdevice(0)
#    a[i] = threadid()
#end


text:
using Tomato
using Test
using JuMP
using GLPK

@testset "Tomato" begin
    @test 1 == 1
end


text:
mutable struct HPOptimisationCache
    theta::AbstractArray{Float64, 1}    # hyperparameters
    K::AbstractArray{Float64, 2}        # covariance matrix
    B::AbstractArray{Float64, 2}        # I + K / σ_n ^ 2
    Q::AbstractArray{Float64, 2}        # inv(B)
    L::AbstractArray{Float64, 2}        # chol(B)
    alpha::AbstractArray{Float64, 2}    # inv(B) * y
    R::AbstractArray{Float64, 2}        # α * α' - Q / σ_n ^ 2
end

HPOptimisationCache() = HPOptimisationCache(
        zeros(0),       # theta
        zeros(0, 0),    # K
        zeros(0, 0),    # B
        zeros(0, 0),    # Q
        zeros(0, 0),    # L
        zeros(0, 0),    # alpha
        zeros(0, 0))    # R

abstract type AbstractGaussianProcess end;

"""
    GPModel

The main type that is used by most functions within the package.

All data matrices are row-major.

# Fields
- `kernel::AbstractGPKernel`: the kernel
- `gp_training_x::AbstractArray{Float64, 2}`: training `x`. Size: ``n \\times d``.
- `gp_training_y::AbstractArray{Float64, 2}`: training `y`. Size: ``n \\times 1``.
- `gp_test_x::AbstractArray{Float64, 2}`: test `x`.  Size: ``m \\times d``.
- `gp_hyperparameters::AbstractArray{Float64, 1}`: kernel hyperparameters, followed by
  standard deviation of intrinsic noise ``\\sigma_n``, which is always the last element in the array.
- `cache::HPOptimisationCache`: cache of matrices that can be re-used between calls to
  `gp_loglikelihood` and `gp_loglikelihood_grad`
"""
mutable struct GPModel <: AbstractGaussianProcess
    kernel::AbstractGPKernel

    cache::HPOptimisationCache

    gp_training_x::AbstractArray{Float64, 2}
    gp_training_y::AbstractArray{Float64, 2}
    gp_test_x::AbstractArray{Float64, 2}

    gp_hyperparameters::AbstractArray{Float64, 1}
end

"""
    GPModel(training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
            training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}
            [,test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0)])

Default constructor of [`GPModel`](@ref), that will use [`SquaredExponentialIsoKernel`](@ref).
Arguments that are passed as 1-d vectors will be reshaped into 2-d.
"""
function GPModel(training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
        training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
        test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0))
    GPModel(training_x=training_x, training_y=training_y, test_x=test_x)
end

"""
    GPModel(training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
            training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
            kernel::AbstractGPKernel
            [,test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0)])

Constructor of [`GPModel`](@ref) that allows the kernel to be specified.
Arguments that are passed as 1-d vectors will be reshaped into 2-d.
"""
function GPModel(training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
        training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}},
        kernel::AbstractGPKernel,
        test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0))
    GPModel(training_x=training_x, training_y=training_y, kernel=kernel, test_x=test_x)
end

"""
    GPModel(;training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        kernel::AbstractGPKernel=SquaredExponentialIsoKernel(),
        gp_hyperparameters::AbstractArray{Float64, 1}=Array{Float64}(0))

Constructor of [`GPModel`](@ref) with explicit arguments.
Arguments that are passed as 1-d vectors will be reshaped into 2-d.
"""
function GPModel(;training_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        training_y::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        test_x::Union{AbstractArray{Float64, 2}, AbstractArray{Float64, 1}}=zeros(0,0),
        kernel::AbstractGPKernel=SquaredExponentialIsoKernel(),
        gp_hyperparameters::AbstractArray{Float64, 1}=zeros(0))

    if ndims(training_x) == 1
        training_x = reshape(training_x, length(training_x), 1)
    end
    if size(training_x, 1) < size(training_x, 2)
        @warn "GPModel: got a $(size(training_x)) size training array. Are you using column major data instead of row major?"
    end
    if ndims(training_y) == 1
        training_y = reshape(training_y, length(training_y), 1)
    end
    if ndims(test_x) == 1
        test_x = reshape(test_x, length(test_x), 1)
    end
    if size(test_x, 1) > 0 && size(training_x, 1) > 0 && size(test_x, 2) != size(training_x, 2)
        error("Test and training x should have the same second dimension")
    end
    if size(training_y, 1) > 0 && size(training_y, 2) != 1
        error("Training y should have second dimension == 1")
    end
    expected_hypers_size = get_hyperparameters_size(kernel, training_x) + 1
    if length(gp_hyperparameters) == 0
        gp_hyperparameters = ones(expected_hypers_size)
    end
    if length(gp_hyperparameters) != expected_hypers_size
        error("Incorrect size of initial hyperparameters vector for ",
            "$(typeof(kernel)): $(length(gp_hyperparameters)). ",
            "Expected $(expected_hypers_size).")
    end
    GPModel(kernel,
        HPOptimisationCache(),
        training_x, training_y, test_x,
        gp_hyperparameters)
end;

function get_hyperparameters_size(gpem::GPModel)
    return get_hyperparameters_size(gpem.kernel, gpem.gp_training_x) + 1
end

"""
    set_hyperparameters(gpm::GPModel, hypers::AbstractArray{Float64, 1})

Set the hyperparameters of the [`GPModel`](@ref)
"""
function set_hyperparameters(gpem::GPModel, hypers::AbstractArray{Float64, 1})
    expected_hypers_size = get_hyperparameters_size(gpem)
    if length(hypers) != expected_hypers_size
        error("Incorrect size of hyperparameters vector for ",
            "$(typeof(gpem.kernel)): $(length(gpem.gp_hyperparameters)). ",
            "Expected $(expected_hypers_size).")
    end
    gpem.gp_hyperparameters = hypers
end

"""
    gp_loglikelihood_log(theta::AbstractArray{Float64, 1}, gpm::GPModel)

Log likelihood function with log hyperparameters. This is the target function of the
hyperparameters optimisation procedure. Its gradient is coputed by [`gp_loglikelihood_grad`](@ref).
"""
function gp_loglikelihood_log(theta::AbstractArray{Float64, 1}, gpem::GPModel)
    cache = gpem.cache
    update_cache!(gpem.kernel, cache, theta, gpem.gp_training_x, gpem.gp_training_y)
    sigma_n2_inv = exp(-2*theta[end])
    x = gpem.gp_training_x
    y = gpem.gp_training_y

    n = size(x, 1)
    ret = -0.5 * dot(y, cache.alpha) - sum(log.(diag(cache.L))) - 0.5 * n * log(2*π/sigma_n2_inv)
end

"""
    gp_loglikelihood(gpm::GPModel)

Compute the log likelihood function, based on the kernel and training data specified in `gpm`.
```math
log p(y \\vert X, \\theta) = - \\frac{1}{2}(y^TK^{-1}y + log \\vert K \\vert + n log 2 \\pi)
```
"""
gp_loglikelihood(gpem::GPModel) = gp_loglikelihood_log(log.(gpem.gp_hyperparameters), gpem)

"""
    gp_loglikelihood(theta::AbstractArray{Float64, 1}, gpm::GPModel)
"""
gp_loglikelihood(theta::AbstractArray{Float64, 1}, gpem::GPModel) = gp_loglikelihood_log(log.(theta), gpem)

"""
    gp_loglikelihood_grad(theta::AbstractArray{Float64, 1}, gpem::GPModel)

Gradient of the log likelihood function ([`gp_loglikelihood_log`](@ref))
with respect to logged hyperparameters.
"""
function gp_loglikelihood_grad(theta::AbstractArray{Float64, 1}, gpem::GPModel)
    cache = gpem.cache
    update_cache!(gpem.kernel, cache, theta, gpem.gp_training_x, gpem.gp_training_y)
    sigma_n2_inv = exp(-2*theta[end])
    x = gpem.gp_training_x
    y = gpem.gp_training_y
    n = size(x, 1)
    kernel_grad = covariance_grad(gpem.kernel, theta[1:end-1], x, cache.R) ./ 2
    if kernel_grad === :covariance_gradient_not_implemented
        return kernel_grad
    end
    dW = cache.Q .* cache.K
    d_sigma_n = sum(dW) * sigma_n2_inv - n + dot(cache.alpha, cache.alpha) / sigma_n2_inv
    push!(kernel_grad, d_sigma_n)
end

"""
    gp_regression(gpm::GPModel; <optional keyword arguments>)

Run the Gaussian Process Regression.

# Arguments
- `gpm`: the [`GPModel`](@ref), that contains the training data (x and y),
  the kernel, the hyperparameters and the test data for running the regression.
- `test_x`: if specified, overrides the test x in `gpm`. Size ``m \\times d``.
- `log_level::Int` (optional): log level. Default is `0`, which is no logging at all. `1`
  makes `gp_regression` print basic information to standard output.
- `full_covariance_matrix::Bool` (optional): whether we need the full
  covariance matrix, or just the variance vector. Defaults to `false` (i.e. just the variance).
- `batch_size::Int` (optional): If `full_covariance_matrix` is set to `false`,
  then the mean and variance vectors will be computed in batches of this size,
  to avoid allocating huge matrices. Defaults to 1000.
- `observation_noise::Bool` (optional): whether the observation noise (with variance
  ``\\sigma_n^2``) should be included in the output variance. Defaults to `true`.

# Return
A tuple of `(mean, var)`. `mean` is a mean vector of the output multivariate Normal
distribution, of size ``m``. `var` is either the covariance matrix of size ``m \\times m``,
or a variance vector of size ``m``, depending on `full_covariance_matrix` flag.
"""
function gp_regression(gpem::GPModel; kwargs...)
    gp_regression_log(log.(gpem.gp_hyperparameters), gpem.gp_test_x, gpem; kwargs...)
end

"""
    gp_regression(test_x::Union{AbstractArray{Float64, 1}, AbstractArray{Float64, 2}},
        gpem::GPModel; <optional keyword arguments>)
"""
function gp_regression(test_x::Union{AbstractArray{Float64, 1}, AbstractArray{Float64, 2}}, gpem::GPModel; kwargs...)
    if ndims(test_x) == 1
        test_x = reshape(test_x, length(test_x), 1)
    end
    gp_regression_log(log.(gpem.gp_hyperparameters), test_x, gpem; kwargs...)
end

function gp_regression_theta(theta::AbstractArray{Float64, 1}, gpem::GPModel; kwargs...)
    gpem.gp_hyperparameters = theta
    gp_regression_log(log.(theta), gpem.gp_test_x, gpem; kwargs...)
end

function gp_regression_log(theta::AbstractArray{Float64, 1},
        x_s::AbstractArray{Float64, 2}, gpem::GPModel;
        log_level::Int=0,
        full_covariance_matrix::Bool=false,
        batch_size::Int=1000,
        observation_noise::Bool=true)
    if (log_level > 0)
        println("Regression log hyperparameters = $(theta)")
    end
    gpem.gp_test_x = x_s

    theta_kernel = theta[1:end-1]
    sigma_n_inv = exp(-theta[end])
    n_s = size(x_s, 1)

    expected_hypers_size = get_hyperparameters_size(gpem)
    if length(theta) != expected_hypers_size
        error("Incorrect size of hyperparameters vector for ",
            "$(typeof(gpem.kernel)): $(length(theta)). ",
            "Expected $(expected_hypers_size).")
    end
    cache = gpem.cache
    update_cache!(gpem.kernel, cache, theta, gpem.gp_training_x, gpem.gp_training_y)

    xs_idx_start = 1
    multiple_batches = !full_covariance_matrix && n_s > batch_size
    xs_idx_end = full_covariance_matrix ? n_s : min(n_s, batch_size)
    ret_mean = Array{Float64}(undef, n_s)
    ret_var = Array{Float64}(undef, n_s)
    batch_counter = 1
    while xs_idx_start <= n_s
        if log_level > 0
            println("Regression batch #$(batch_counter) $(xs_idx_start):$(xs_idx_end)")
            batch_counter += 1
        end
        batch = x_s[xs_idx_start:xs_idx_end, :]
        K_s = covariance(gpem.kernel, theta_kernel, gpem.gp_training_x, batch)'
        v = cache.L' \ (sigma_n_inv * K_s')
        batch_mean = K_s * cache.alpha
        batch_mean = reshape(batch_mean, length(batch_mean))
        if multiple_batches || !full_covariance_matrix
            K_ss = covariance_diagonal(gpem.kernel, theta_kernel, batch)
            if observation_noise
                K_ss += ones(size(batch, 1), 1) * exp(2 * theta[end])
            end
            batch_var = K_ss - sum(v .^ 2, dims=1)'
            batch_var = reshape(batch_var, length(batch_var))
        else
            K_ss = covariance(gpem.kernel, theta_kernel, batch, batch)
            if observation_noise
                K_ss += I * exp(2 * theta[end])
            end
            batch_var = K_ss - v' * v
            batch_var = 0.5 * (batch_var' + batch_var)
        end

        if multiple_batches
            ret_mean[xs_idx_start:xs_idx_end] = batch_mean
            ret_var[xs_idx_start:xs_idx_end] = batch_var
            xs_idx_start = xs_idx_end + 1
            xs_idx_end = min(n_s, xs_idx_start + batch_size - 1)
        else
            return batch_mean, batch_var
        end
    end
    return (ret_mean, ret_var)
end

function update_cache!(kernel::AbstractGPKernel, cache::HPOptimisationCache,
        theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, y::AbstractArray{Float64, 2})
    if cache.theta != theta
        n = size(x, 1)
        if size(cache.theta) != size(theta)
            cache.theta = copy(theta)
        else
            copyto!(cache.theta, theta)
        end
        sigma_n2_inv = exp(-2*theta[end])
        cache.K = covariance_training(kernel, theta[1:end-1], x)
        cache.B = cache.K * sigma_n2_inv + I
        cache.L = cholesky(Hermitian(cache.B)).U
        cache.Q = cache.L \ (cache.L' \ I)
        cache.alpha = (cache.L \ (cache.L' \ y)) .* sigma_n2_inv
        cache.R = cache.alpha * cache.alpha' - cache.Q .* sigma_n2_inv
    end
    0
end

"""
    gp_regression_sample(test_x::Union{AbstractArray{Float64, 1}, AbstractArray{Float64, 2}}, n_samples::Int64, gpem::GPModel)

Return `n_samples` random samples from the Gaussian process posterior, evaluated at `test_x`. [`gp_regression`](@ref).

# Arguments
- `test_x`: if specified, overrides the test x in `gpm`. Size ``m \\times d``.
- `n_samples`: integer specifying the number of posterior samples.
- `gpm`: the [`GPModel`](@ref), that contains the training data (x and y),
  the kernel, the hyperparameters and the test data for running the regression.
- `full_cov_matrix`: whether to use the full covariance matrix or just its diagonal elements (default `true`).
# Return
An array of posterior samples with shape ``m \\times`` `n_samples` if `n_samples`>1 and ``m`` otherwise.
"""
function gp_regression_sample(test_x::Union{AbstractArray{Float64, 1}, AbstractArray{Float64, 2}}, n_samples::Int64, gpem::GPModel, full_cov_matrix=true)
    mu, sigma = gp_regression(test_x, gpem, full_covariance_matrix=full_cov_matrix)
    post_samples = rand(MvNormal(mu, sigma), n_samples)
    if n_samples==1
        post_samples = dropdims(post_samples, dims=2)
    end
    return post_samples
end


text:
#!/usr/bin/env julia

# Create rectilinear grid VTK file.

using DataFrames, WriteVTK
import Compat.UTF8String
const FloatType = Float32

"""
## vtk

Plots displacements and directions

### Function
```julia
mesh(data, fm_df, sigma_df, dir, fname)
```

### Arguments
```julia
* data::Dict                 : Input dictionary
* fm_df::DataFrame           : Forces and moments DataFrame
* sigma_df::DataFrame        : Stresses DataFrame
* dir::AbstractString        : Project directory
* fname::AbstractString      : Output VTK file name
```
"""
function vtk(data::Dict, fm_df, sigma_df, dir, fname)
  vtk_filename_noext = dir*"/"*fname
  dim = data[:struc_el].ndim

  # Define grid

  Ni, Nj = length(data[:x_coords]), length(data[:y_coords])
  if dim == 2
    Nk = 1
  elseif dim == 3
    Nk = length(data[:z_coords])
  end

  x = data[:x_coords]
  y = data[:y_coords]
  z = zeros(FloatType, Nk)
  if dim == 3
    z = data[:z_coords]
  end

  # Extract node data

  disp = Vector{FloatType}(convert(Array, fm_df[:disp]))
  rotx = Vector{FloatType}(convert(Array, fm_df[:rotx]))
  roty = Vector{FloatType}(convert(Array, fm_df[:roty]))
  twistxy = Vector{FloatType}(convert(Array, fm_df[:twistxy]))

  # Extract cell data
  # Note that in structured grids, the cells are the hexahedra (3D) or quads (2D)
  # formed between grid points.

  local cdata, sigx, sigy, tauxy
  if dim == 2
      sigx = reshape(convert(Array{FloatType}, sigma_df[:sigx]), 2, 2)
      sigy = reshape(convert(Array{FloatType}, sigma_df[:sigy]), 2, 2)
      tauxy = reshape(convert(Array{FloatType}, sigma_df[:tauxy]), 2, 2)
  elseif dim == 3
      cdata = zeros(FloatType, Ni-1, Nj-1, Nk-1)
      for k = 1:Nk-1, j = 1:Nj-1, i = 1:Ni-1
          cdata[i, j, k] = 2i + 3k * sin(3*pi * (j-1) / (Nj-2))
      end
  end

  # Test extents (this is optional!!)
  ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] + 42

  # Initialise new vtr file (rectilinear grid).
  @time begin
      local vtk
      if dim == 2
          vtk = vtk_grid(vtk_filename_noext*"_$(dim)D", x, y; extent=ext)
      elseif dim == 3
          vtk = vtk_grid(vtk_filename_noext*"_$(dim)D", x, y, z;
                         extent=ext)
      end

      # Add data.
      vtk_point_data(vtk, disp, "disp")
      vtk_point_data(vtk, rotx, "rotx")
      vtk_point_data(vtk, roty, "roty")
      vtk_point_data(vtk, twistxy, "twistxy")
      vtk_cell_data(vtk, sigx, "sigx")
      vtk_cell_data(vtk, sigy, "sigy")
      vtk_cell_data(vtk, tauxy, "tauxy")

      # Save and close vtk file.
      vtk_save(vtk)
  end
end


text:
module StochasticProcesses

export cumsim, sim, initial
export SimResult, CumsimResult

export AItoProcess,
       BrownianBridge,
       BrownianMotion,
       BrownianMotionWithDrift,
       CompositeProcess,
       FirstTime,
       GeometricBrownianMotion,
       ItoIntegral,
       ItoProcess,
       PoissonProcess,
       SDE

# Basic processes information
export distribution, solution

# Internal stochastic process interface
abstract AProcess

include("brownian_bridge.jl")
include("first_time.jl")
include("poisson.jl")
include("results.jl")
include("solvers.jl")

using Distributions
import Base.convert, Base.size, Base.rand

function Base.rand{P <: AProcess}(process::P, t, k::Int=1)
  # todo: better way to code this?
  if method_exists(distribution, (typeof(process), Float64))
    return Base.rand(distribution(process, t[end]), k)
  end
  return sim(process, t, k)
end

"Ito Process(f,y0): process satisfying equation dy = f(t, dt, b, db, y)."
immutable ItoProcess{F, Y} <: AProcess
  f::F
  y0::Y
end

function init!{F,Y}(process::ItoProcess{F,Y}, g::Generator)
  g.y = repeat(vwrap(initial(process)), outer=(1, g.k))
  g.y
end

@inline function next!{F,Y}(process::ItoProcess{F,Y}, g::Generator, i)
  dy = process.f(g.t[i-1], g.dt[i-1], g.b, g.db, g.y)
  add!(g.y, dy)
  g.y
end

initial(process::ItoProcess) = process.y0

# CompositeProcess
# The resulting process is f(t, y), where y is the value of p
immutable CompositeProcess{P <: AProcess, F} <: AProcess
  p::P
  f::F
end

init!{P, F}(process::CompositeProcess{P,F}, g::Generator) =
    process.f.(0., init!(process.p, g))

@inline function next!{P, F}(process::CompositeProcess{P,F}, g::Generator, i)
  next!(process.p, g, i)
  process.f.(g.t[i], g.y)
end

initial{P, F}(process::CompositeProcess{P,F}) = initial(process.p)

size{P, F}(process::CompositeProcess{P,F}) = size(process.p)

"A process that is convertible to ItoProcess"
abstract AItoProcess <: AProcess

init!{P <: AItoProcess}(process::P, state) =
    init!(convert(ItoProcess, process), state)

next!{P <: AItoProcess}(process::P, state, i) =
    next!(convert(ItoProcess, process), state, i)

size{P <: AItoProcess}(process::P) = size(convert(ItoProcess, process))

initial{P <: AItoProcess}(process::P) = initial(convert(ItoProcess, process))

include("basic_processes.jl")

# Utilities

# Distribution with only one value
type Constant
  x
end

Base.rand(c::Constant, k) = fill(c.x, k)

function collapse_but_first(A::AbstractArray)
    s = size(A)
    dims = tuple(find([d == 1 && i != 1 for (i,d) in enumerate(s)])...)
    return squeeze(A,dims)
end

vwrap(x::AbstractArray) = x
vwrap(x::Number) = [x]

@inline function add!(y, dy)
  for i in eachindex(y)
    @inbounds y[i] += dy[i]
  end
end

@inline function cmul!(y, dy)
  for i in eachindex(y)
    @inbounds y[i] *= dy
  end
end

end


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
comm_size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

dest = mod(rank+1, comm_size)
src  = mod(rank-1, comm_size)

@testset "contiguous" begin
    X = ArrayType(rank .+ collect(reshape(1.0:16.0, 4, 4)))
    Y = ArrayType(zeros(4))
    req_send = MPI.Isend(@view(X[:,1]), comm; dest=dest)
    req_recv = MPI.Irecv!(Y, comm; source=src)

    MPI.Waitall([req_send, req_recv])

    @test Y == X[:,1] .- rank .+ src

    Y = ArrayType(zeros(2))

    req_send = MPI.Isend(Y, comm; dest=dest, tag=1)
    req_recv = MPI.Irecv!(@view(X[3:4,1]), comm; source=src, tag=1)

    MPI.Waitall([req_send, req_recv])

    @test X[3:4,1] == Y
end

@testset "strided" begin
    X = ArrayType(rank .+ collect(reshape(1.0:16.0, 4, 4)))
    Y = ArrayType(zeros(4))
    req_send = MPI.Isend(@view(X[2,:]), comm; dest=dest, tag=0)
    req_recv = MPI.Irecv!(Y, comm; source=src, tag=0)

    MPI.Waitall([req_send, req_recv])

    @test Y == X[2,:] .- rank .+ src

    Y = ArrayType(zeros(2))

    req_send = MPI.Isend(Y, comm; dest=dest, tag=1)
    req_recv = MPI.Irecv!(@view(X[3,1:2]), comm; source=src, tag=1)

    MPI.Waitall([req_send, req_recv])

    @test X[3,1:2] == Y
end

@testset "dense subarray" begin
    X = ArrayType(rank .+ collect(reshape(1.0:16.0, 4, 4)))
    Y = ArrayType(zeros(2,2))
    req_send = MPI.Isend(@view(X[2:3,3:4]), comm; dest=dest, tag=0)
    req_recv = MPI.Irecv!(Y, comm; source=src, tag=0)

    MPI.Waitall([req_send, req_recv])

    @test Y == X[2:3,3:4] .- rank .+ src

    Y = ArrayType(zeros(2,2))

    req_send = MPI.Isend(Y, comm; dest=dest, tag=1)
    req_recv = MPI.Irecv!(@view(X[3:4,1:2]), comm; source=src, tag=1)

    MPI.Waitall([req_send, req_recv])

    @test X[3:4,1:2] == Y
end

GC.gc()
MPI.Finalize()
@test MPI.Finalized()


text:
using Distributed

nprocs()

addprocs(5)

nprocs()

function count_heads(n)
    c::Int = 0
    for i = 1:n
        c += rand(Bool)
    end
    c
end

a = @spawnat :any count_heads(100000000)

fetch(a)

@everywhere function count_heads(n)
    c::Int = 0
    for i = 1:n
        c += rand(Bool)
    end
    c
end

a = @spawnat :any count_heads(100000000)

fetch(a)

b = @spawnat :any count_heads(100000000)

fetch(a) + fetch(b)

nheads = @distributed (+) for i = 1:200000000
    Int(rand(Bool))
end

a = zeros(100000)

@distributed for i = 1:100000
    a[i] = i
end

a

using SharedArrays

a = SharedArray{Float64}(10)

@distributed for i = 1:10
    a[i] = i
end

a

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP