text:
# definition of some distribution types

## Distribution & Make

abstract type Distribution{T} end

Base.eltype(::Type{<:Distribution{T}}) where {T} = T

abstract type Make{T} <: Distribution{T} end

struct Make0{T} <: Make{T} end

Make( ::Type{T}) where {T} = Make0{T}()
Make0(::Type{T}) where {T} = Make0{T}()
make(::Type{T}) where {T} = Make0{maketype(T)}()

# default
maketype(::Type{T}) where {T} = T

struct Make1{T,X} <: Make{T}
    x::X
end

Make{T}(x::X)      where {T,X} = Make1{T,X}(      x)
Make{T}(::Type{X}) where {T,X} = Make1{T,Type{X}}(X)

# for expliciteness (allows using Make1 instead of Make)
Make1{T}(x::X)      where {T,X} = Make1{T,X}(      x)
Make1{T}(::Type{X}) where {T,X} = Make1{T,Type{X}}(X)

make(::Type{T}, x::X)      where {T,X} = Make{maketype(T,x)}(x)
make(::Type{T}, ::Type{X}) where {T,X} = Make{maketype(T,X)}(X)

# default
maketype(::Type{T}, x) where {T} = T

find_deduced_type(::Type{T}, ::X,     ) where {T,X} = deduce_type(T, gentype(X))
find_deduced_type(::Type{T}, ::Type{X}) where {T,X} = deduce_type(T, X)

struct Make2{T,X,Y} <: Make{T}
    x::X
    y::Y
end

Make{T}(x::X,      y::Y)      where {T,X,Y} = Make2{T,X,      Y}(      x, y)
Make{T}(::Type{X}, y::Y)      where {T,X,Y} = Make2{T,Type{X},Y}(      X, y)
Make{T}(x::X,      ::Type{Y}) where {T,X,Y} = Make2{T,X,      Type{Y}}(x, Y)
Make{T}(::Type{X}, ::Type{Y}) where {T,X,Y} = Make2{T,Type{X},Type{Y}}(X, Y)

# for expliciteness (allows using Make2 instead of Make)
Make2{T}(x::X,      y::Y)      where {T,X,Y} = Make2{T,X,      Y}(      x, y)
Make2{T}(::Type{X}, y::Y)      where {T,X,Y} = Make2{T,Type{X},Y}(      X, y)
Make2{T}(x::X,      ::Type{Y}) where {T,X,Y} = Make2{T,X,      Type{Y}}(x, Y)
Make2{T}(::Type{X}, ::Type{Y}) where {T,X,Y} = Make2{T,Type{X},Type{Y}}(X, Y)

make(::Type{T}, x::X,      y::Y)      where {T,X,Y} = Make{maketype(T,x,y)}(x, y)
make(::Type{T}, ::Type{X}, y::Y)      where {T,X,Y} = Make{maketype(T,X,y)}(X, y)
make(::Type{T}, x::X,      ::Type{Y}) where {T,X,Y} = Make{maketype(T,x,Y)}(x, Y)
make(::Type{T}, ::Type{X}, ::Type{Y}) where {T,X,Y} = Make{maketype(T,X,Y)}(X, Y)

# default
maketype(::Type{T}, x, y) where {T} = T

find_deduced_type(::Type{T}, ::X,       ::Y)       where {T,X,Y} = deduce_type(T, gentype(X), gentype(Y))
find_deduced_type(::Type{T}, ::Type{X}, ::Y)       where {T,X,Y} = deduce_type(T, X,          gentype(Y))
find_deduced_type(::Type{T}, ::X,       ::Type{Y}) where {T,X,Y} = deduce_type(T, gentype(X), Y)
find_deduced_type(::Type{T}, ::Type{X}, ::Type{Y}) where {T,X,Y} = deduce_type(T, X,          Y)

struct Make3{T,X,Y,Z} <: Make{T}
    x::X
    y::Y
    z::Z
end

Make{T}(x::X,      y::Y,      z::Z)      where {T,X,Y,Z} = Make3{T,X,      Y,      Z      }(x, y, z)
Make{T}(::Type{X}, y::Y,      z::Z)      where {T,X,Y,Z} = Make3{T,Type{X},Y,      Z      }(X, y, z)
Make{T}(x::X,      ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{T,X,      Type{Y},Z      }(x, Y, z)
Make{T}(::Type{X}, ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{T,Type{X},Type{Y},Z      }(X, Y, z)
Make{T}(x::X,      y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{T,X,      Y,      Type{Z}}(x, y, Z)
Make{T}(::Type{X}, y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{T,Type{X},Y,      Type{Z}}(X, y, Z)
Make{T}(x::X,      ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{T,X,      Type{Y},Type{Z}}(x, Y, Z)
Make{T}(::Type{X}, ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{T,Type{X},Type{Y},Type{Z}}(X, Y, Z)

# for expliciteness (allows using Make3 instead of Make)
Make3{T}(x::X,      y::Y,      z::Z)      where {T,X,Y,Z} = Make3{T,X,      Y,      Z      }(x, y, z)
Make3{T}(::Type{X}, y::Y,      z::Z)      where {T,X,Y,Z} = Make3{T,Type{X},Y,      Z      }(X, y, z)
Make3{T}(x::X,      ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{T,X,      Type{Y},Z      }(x, Y, z)
Make3{T}(::Type{X}, ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{T,Type{X},Type{Y},Z      }(X, Y, z)
Make3{T}(x::X,      y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{T,X,      Y,      Type{Z}}(x, y, Z)
Make3{T}(::Type{X}, y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{T,Type{X},Y,      Type{Z}}(X, y, Z)
Make3{T}(x::X,      ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{T,X,      Type{Y},Type{Z}}(x, Y, Z)
Make3{T}(::Type{X}, ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{T,Type{X},Type{Y},Type{Z}}(X, Y, Z)

make(::Type{T}, x::X,      y::Y,      z::Z)      where {T,X,Y,Z} = Make3{maketype(T, x, y, z)}(x, y, z)
make(::Type{T}, ::Type{X}, y::Y,      z::Z)      where {T,X,Y,Z} = Make3{maketype(T, X, y, z)}(X, y, z)
make(::Type{T}, x::X,      ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{maketype(T, x, Y, z)}(x, Y, z)
make(::Type{T}, ::Type{X}, ::Type{Y}, z::Z)      where {T,X,Y,Z} = Make3{maketype(T, X, Y, z)}(X, Y, z)
make(::Type{T}, x::X,      y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{maketype(T, x, y, Z)}(x, y, Z)
make(::Type{T}, ::Type{X}, y::Y,      ::Type{Z}) where {T,X,Y,Z} = Make3{maketype(T, X, y, Z)}(X, y, Z)
make(::Type{T}, x::X,      ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{maketype(T, x, Y, Z)}(x, Y, Z)
make(::Type{T}, ::Type{X}, ::Type{Y}, ::Type{Z}) where {T,X,Y,Z} = Make3{maketype(T, X, Y, Z)}(X, Y, Z)

# default
maketype(::Type{T}, x, y, z) where {T} = T

# deduce_type

deduce_type(::Type{T}, ::Type{X}, ::Type{Y}) where {T,X,Y} = _deduce_type(T, Val(isconcretetype(T)), X, Y)
deduce_type(::Type{T}, ::Type{X}) where {T,X} = _deduce_type(T, Val(isconcretetype(T)), X)

_deduce_type(::Type{T}, ::Val{true},  ::Type{X}, ::Type{Y}) where {T,X,Y} = T
_deduce_type(::Type{T}, ::Val{false}, ::Type{X}, ::Type{Y}) where {T,X,Y} = deduce_type(T{X}, Y)

_deduce_type(::Type{T}, ::Val{true},  ::Type{X}) where {T,X} = T
_deduce_type(::Type{T}, ::Val{false}, ::Type{X}) where {T,X} = T{X}

## Uniform

abstract type Uniform{T} <: Distribution{T} end

struct UniformType{T} <: Uniform{T} end

Uniform(::Type{T}) where {T} = UniformType{T}()

Base.getindex(::UniformType{T}) where {T} = T

struct UniformWrap{T,E} <: Uniform{E}
    val::T
end

Uniform(x::T) where {T} = UniformWrap{T,gentype(T)}(x)

Base.getindex(x::UniformWrap) = x.val

## Normal & Exponential

abstract type Normal{T} <: Distribution{T} end

struct Normal01{T} <: Normal{T} end

struct Normalμσ{T} <: Normal{T}
    μ::T
    σ::T
end

const NormalTypes = Union{AbstractFloat,Complex{<:AbstractFloat}}

Normal(::Type{T}=Float64) where {T<:NormalTypes} = Normal01{T}()
Normal(μ::T, σ::T) where {T<:NormalTypes} = Normalμσ(μ, σ)
Normal(μ::T, σ::T) where {T<:Real} = Normalμσ(AbstractFloat(μ), AbstractFloat(σ))
Normal(μ, σ) = Normal(promote(μ, σ)...)

abstract type Exponential{T} <: Distribution{T} end

struct Exponential1{T} <: Exponential{T} end

struct Exponentialθ{T} <: Exponential{T}
    θ::T
end

Exponential(::Type{T}=Float64) where {T<:AbstractFloat} = Exponential1{T}()
Exponential(θ::T) where {T<:AbstractFloat} = Exponentialθ(θ)
Exponential(θ::Real) = Exponentialθ(AbstractFloat(θ))

## floats

abstract type FloatInterval{T<:AbstractFloat} <: Uniform{T} end

abstract type CloseOpen{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenClose{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type CloseClose{T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenOpen{  T<:AbstractFloat} <: FloatInterval{T} end

struct CloseOpen12{T<:AbstractFloat} <: CloseOpen{T} end # interval [1,2)

struct CloseOpen01{ T<:AbstractFloat} <: CloseOpen{T}  end # interval [0,1)
struct OpenClose01{ T<:AbstractFloat} <: OpenClose{T}  end # interval (0,1]
struct CloseClose01{T<:AbstractFloat} <: CloseClose{T} end # interval [0,1]
struct OpenOpen01{  T<:AbstractFloat} <: OpenOpen{T}   end # interval (0,1)

struct CloseOpenAB{T<:AbstractFloat} <: CloseOpen{T} # interval [a,b)
    a::T
    b::T

    CloseOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenCloseAB{T<:AbstractFloat} <: OpenClose{T} # interval (a,b]
    a::T
    b::T

    OpenCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct CloseCloseAB{T<:AbstractFloat} <: CloseClose{T} # interval [a,b]
    a::T
    b::T

    CloseCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenOpenAB{T<:AbstractFloat} <: OpenOpen{T} # interval (a,b)
    a::T
    b::T

    OpenOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

check_interval(a, b) = a >= b && throw(ArgumentError("invalid interval specification"))

const FloatInterval_64 = FloatInterval{Float64}
const CloseOpen01_64   = CloseOpen01{Float64}
const CloseOpen12_64   = CloseOpen12{Float64}

CloseOpen01(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen12(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen12{T}()

CloseOpen(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen(a::T, b::T) where {T<:AbstractFloat} = CloseOpenAB{T}(a, b)

OpenClose(::Type{T}=Float64) where {T<:AbstractFloat} = OpenClose01{T}()
OpenClose(a::T, b::T) where {T<:AbstractFloat} = OpenCloseAB{T}(a, b)

CloseClose(::Type{T}=Float64) where {T<:AbstractFloat} = CloseClose01{T}()
CloseClose(a::T, b::T) where {T<:AbstractFloat} = CloseCloseAB{T}(a, b)

OpenOpen(::Type{T}=Float64) where {T<:AbstractFloat} = OpenOpen01{T}()
OpenOpen(a::T, b::T) where {T<:AbstractFloat} = OpenOpenAB{T}(a, b)

# convenience functions

CloseOpen(a, b) = CloseOpen(promote(a, b)...)
CloseOpen(a::T, b::T) where {T} = CloseOpen(AbstractFloat(a), AbstractFloat(b))

OpenClose(a, b) = OpenClose(promote(a, b)...)
OpenClose(a::T, b::T) where {T} = OpenClose(AbstractFloat(a), AbstractFloat(b))

CloseClose(a, b) = CloseClose(promote(a, b)...)
CloseClose(a::T, b::T) where {T} = CloseClose(AbstractFloat(a), AbstractFloat(b))

OpenOpen(a, b) = OpenOpen(promote(a, b)...)
OpenOpen(a::T, b::T) where {T} = OpenOpen(AbstractFloat(a), AbstractFloat(b))

## Bernoulli

struct Bernoulli{T<:Number} <: Distribution{T}
    p::Float64

    Bernoulli{T}(p::Real) where {T} = let pf = Float64(p)
        0.0 <= pf <= 1.0 ? new(pf) :
            throw(DomainError(p, "Bernoulli: parameter p must satisfy 0.0 <= p <= 1.0"))
    end
end

Bernoulli(p::Real=0.5) = Bernoulli(Int, p)
Bernoulli(::Type{T}, p::Real=0.5) where {T<:Number} = Bernoulli{T}(p)


text: abstract Theme <: Tile

ext: jl


text:
c0 = chain(:standard,1,10,Float64,Float64;numburnin=2)
c1 = chain(:standard,3,4,Int16,Float64)
c2 = chain(:gradient,2,3,Float32,Float32)
@test numparas(c0) == 1 && numsamples(c0) == 10
@test numparas(c1) == 3 && numsamples(c1) == 4
@test numparas(c2) == 2 && numsamples(c2) == 3

for a in [(c0,(0,1,0,1:0),(0,3,0,3:0),(0,1,0,1:0),(10,1,10,1:10)),
          (c1,(0,1,0,1:0),(0,1,0,1:0),(0,1,0,1:0),(4,1,4,1:4))]
    for (i,e)  in enumerate([:burnin,:main,:proposed,:all])
        @test GeneralizedMetropolisHastings._numsamples(Val{e},a[1]) == a[i+1][1]
        @test GeneralizedMetropolisHastings._startindex(Val{e},a[1]) == a[i+1][2]
        @test GeneralizedMetropolisHastings._endindex(Val{e},a[1]) == a[i+1][3]
        @test GeneralizedMetropolisHastings._index(a[1],e) == a[i+1][4]
    end
end

s0 = samples(:base,1,4,Float64,Float64)
s1 = samples(:base,3,2,Int16,Float64)
s2 = samples(:tensor,2,2,Float32,Float32)

#test that you can't add a base sample to a gradient chain
@test_throws MethodError store!(c2,s0,1)

#test that the number of parameters is the same for sample and chain
@test_throws AssertionError store!(c0,s1,1)
@test_throws AssertionError store!(c1,s0,1)

copy!(s0.values,1.0:4.0)
store!(c0,s0,4)
store!(c0,s0,3)
store!(c0,s0,2)
@test c0.values[:,1:3] == s0.values[:,4:-1:2]
@test isequal(c0.loglikelihood[1:3],s0.loglikelihood[4:-1:2])

for a in [(c0,(2,1,2,1:2),(1,3,3,3:3),(3,1,3,1:3),(10,1,10,1:10)),
          (c1,(0,1,0,1:0),(0,1,0,1:0),(0,1,0,1:0),(4,1,4,1:4))]
    for (i,e) in enumerate([:burnin,:main,:proposed,:all])
        @test GeneralizedMetropolisHastings._numsamples(Val{e},a[1]) == a[i+1][1]
        @test GeneralizedMetropolisHastings._startindex(Val{e},a[1]) == a[i+1][2]
        @test GeneralizedMetropolisHastings._endindex(Val{e},a[1]) == a[i+1][3]
        @test GeneralizedMetropolisHastings._index(a[1],e) == a[i+1][4]
    end
end

copy!(s1.values,[1.0 4.0;2.0 5.0;3.0 6.0])
store!(c1,s1,2)
store!(c1,s1,1)
store!(c1,s1,2)
store!(c1,s1,2)
@test c1.values == [4.0 1.0 4.0 4.0;5.0 2.0 5.0 5.0;6.0 3.0 6.0 6.0]
@test isequal(c1.loglikelihood,s1.loglikelihood[[2,1,2,2]])

for a in [(c0,(2,1,2,1:2),(1,3,3,3:3),(3,1,3,1:3),(10,1,10,1:10)),
          (c1,(0,1,0,1:0),(4,1,4,1:4),(4,1,4,1:4),(4,1,4,1:4))]
    for (i,e) in enumerate([:burnin,:main,:proposed,:all])
        @test GeneralizedMetropolisHastings._numsamples(Val{e},a[1]) == a[i+1][1]
        @test GeneralizedMetropolisHastings._startindex(Val{e},a[1]) == a[i+1][2]
        @test GeneralizedMetropolisHastings._endindex(Val{e},a[1]) == a[i+1][3]
        @test GeneralizedMetropolisHastings._index(a[1],e) == a[i+1][4]
    end
end

copy!(s2.values,[1.0 3.0;2.0 4.0])
store!(c2,s2,1)
store!(c2,s2,1)
store!(c2,s2,2)
@test c2.values == [1.0 1.0 3.0;2.0 2.0 4.0]
@test isequal(c2.loglikelihood,s2.loglikelihood[[1,1,2]])
@test isequal(c2.gradloglikelihood,s2.gradloglikelihood[:,[1,1,2]])

#these should give warnings that the chain is full
println("Two warnings should be printed below this line")
store!(c1,s1,1)
store!(c2,s2,1)

@test c0.numproposed == 3 && accepted!(c0,IndicatorStationary{Float64}(zeros(3),[1,4,3,2])) == 3
@test c1.numproposed == 4 && accepted!(c1,IndicatorStationary{Float32}(zeros(Float32,4),[1,2,1,2,2])) == 3
@test c2.numproposed == 3 && accepted!(c2,IndicatorStationary{Float64}(zeros(3),[1,1,1,2])) == 1

immutable ChainTestModel <: AbstractModel
    parameters
end

p0 = parameters([0.0],[2.0])
p1 = parameters([0.0,1.0,2.0],[10.0,5.0,3.0])
p2 = parameters([2.0,2.0],[10.0,10.0])
m0 = ChainTestModel(p0)
m1 = ChainTestModel(p1)
m2 = ChainTestModel(p2)

for a in [(c0,m0,2,1,3,10),
          (c1,m1,0,4,4,4),
          (c2,m2,0,3,3,3)]
    for (i,e) in enumerate([:burnin,:main,:proposed,:all])
        @test numsamples(a[1];store=e) == a[i+2]
        @test size(samples(a[1];store=e)) == (numparas(a[1]),a[i+2])
        @test size(samples(a[1];paras=1,store=e)) == (a[i+2],)
        @test length(loglikelihood(a[1];store=e)) == a[i+2]
        @test length(logprior(a[1],a[2];store=e)) == a[i+2]
        @test length(logposterior(a[1],a[2];store=e)) == a[i+2]
    end
end

for a in [(c0,1:2,3:3,1:3,1:10),
          (c1,1:0,1:4,1:4,1:4)]
    for (i,e) in enumerate([:burnin,:main,:proposed,:all])
        @test_approx_eq mean(a[1];store=e) mean(a[1].values[:,a[i+1]],2)
        @test_approx_eq median(a[1];store=e) length(a[i+1])>0?median(a[1].values[:,a[i+1]],2):NaN
        @test_approx_eq std(a[1];store=e) std(a[1].values[:,a[i+1]],2)
        @test_approx_eq var(a[1];store=e) var(a[1].values[:,a[i+1]],2)
    end
end

@test_approx_eq logprior(c0,m0) logprior(p0,c0.values,Float64)
@test_approx_eq logprior(c1,m1) logprior(p1,c1.values,Float64)
@test_approx_eq logprior(c2,m2) logprior(p2,c2.values,Float64)

@test_approx_eq logposterior(c0,m0) loglikelihood(c0) + logprior(c0,m0)
@test_approx_eq logposterior(c1,m1) loglikelihood(c1) + logprior(c1,m1)
@test_approx_eq logposterior(c2,m2) loglikelihood(c2) + logprior(c2,m2)

###Test show()
println()
println("====================")
println("Test show() function")
show(c0)
show(c1)
show(c2)
println("End  show() function")
println("====================")
println()
nothing


text:
function num_borders(r::Robot)
    side = Ost
    num=0
    while isborder(r,Nord)==false
        while isborder(r,side)==false
            if pass_the_borders(r,side) == true
                num += 1
                walk_by(r,side)
            end
        end
        move!(r,Nord)
        side = inverse(side)
    end
    return num
end
function pass_the_borders(r::Robot,side::HorizonSide)
    while isborder(r,Nord)==false
        if isborder(r,side) == true
            return false
        end
        move!(r,side)
    end
    return true
end
walk_by(r::Robot,side::HorizonSide) = while isborder(r,Nord)==true move!(r,side) end
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
module Geo

include("../src/Jaynes.jl")
using .Jaynes
using Distributions

geo(p::Float64) = rand(:flip, Bernoulli(p)) ? 1 : 1 + rand(:geo, geo, p)

# Define as primitive.
@primitive function logpdf(fn::typeof(geo), p, count)
    return Distributions.logpdf(Geometric(p), count)
end

# Example trace.
ret, cl = simulate(() -> rand(:geo, geo, 0.05))
display(cl.trace)

end # module


text:
mutable struct linear_solver_JULIA <: abstract_linear_system_solver
	_factor # TO DO, give type
	_factor_defined::Bool

	# options
	sym::Symbol # this determines if we use LU, cholesky or LDL.
	safe_mode::Bool
	recycle::Bool

  function linear_solver_JULIA(sym::Symbol, safe_mode::Bool, recycle::Bool)
			this = new();
			this.sym = sym
			this.safe_mode = safe_mode
			this.recycle = recycle
			this._factor_defined = false
			return this
  end
end

function ls_factor!(solver::linear_solver_JULIA, SparseMatrix::SparseMatrixCSC{Float64,Int64}, n::Int64, m::Int64, timer::class_advanced_timer)
			inertia_status_val = 1;

			start_advanced_timer(timer, "JULIA/factorize")
			if solver.sym == :unsymmetric
				# incomplete ---- how does one compute inertia????
				 solver._factor = lufact(SparseMatrix);
			elseif solver.sym == :definite
				# do a cholesky factorization
				@assert(m == 0)
				try
					if !solver.recycle || !solver._factor_defined
						#solver._factor = cholfact(Symmetric(SparseMatrix,:L));
						solver._factor = cholesky(Symmetric(SparseMatrix,:L));
						solver._factor_defined = true
					else
						cholfact!(solver._factor,Symmetric(SparseMatrix,:L));
					end
				catch(e)
					if typeof(e) == LinearAlgebra.PosDefException
						inertia_status_val = 0
					else
						println("ERROR in linear_solver_JULIA.ls_factor!")
						throw(e)
					end
				end
			elseif solver.sym == :symmetric
				# LDL factorization
				try
					if !solver.recycle || !solver._factor_defined
						#solver._factor = ldltfact(Symmetric(SparseMatrix,:L));
						solver._factor = ldlt(Symmetric(SparseMatrix,:L));
						solver._factor_defined = true
					else
						ldltfact!(solver._factor,Symmetric(SparseMatrix,:L));
					end
				catch(e)
					#PosDefException was added when updating Julia from 0.6.4 to 0.7.0 by Fadi Hamad
					#Because ArgumentError is not thrown anymore and it is being replaced by PosDefException
					#This is related to (PosDefException: matrix is not positive definite; Cholesky factorization failed.)
					if typeof(e) == ArgumentError || typeof(e) == LinearAlgebra.PosDefException || typeof(e) == ZeroPivotException
					#if typeof(e) == LinearAlgebra.PosDefException
						inertia_status_val = 0
					else
						println("ERROR in saddle_solver_JULIA.ls_factor!")
						throw(e)
					end
				end

				if inertia_status_val == 1
					# do something !!!!!!
					di = diag(solver._factor)
					tol = 1e-20
					pos_eigs =  sum(di .> tol)
					zero_eigs = sum((-tol .<= di) .& (di .<= tol))
					neg_eigs = sum(di .< -tol)
					nan_eigs = sum(isnan.(di))
					inf_eigs = sum(isinf.(di))
					if nan_eigs + inf_eigs <= 0
						inertia_status_val = inertia_status(pos_eigs, neg_eigs, zero_eigs, n, m)
					else
						if inf_eigs > 0
							warn("Inf appears in diagonal")
						end
						if nan_eigs > 0
							warn("NaN appears in diagonal")
						end
						inertia_status_val = 0
					end
				end
			else
				error("this.options.sym = " * string(solver.sym) * " not supported")
			end
			pause_advanced_timer(timer, "JULIA/factorize")

			return inertia_status_val; # inertia
end

function ls_solve!(solver::linear_solver_JULIA, my_rhs::Array{Float64,1}, my_sol::Array{Float64,1}, timer::class_advanced_timer)
	start_advanced_timer(timer, "JULIA/ls_solve")
  my_sol[1:length(my_sol)] = solver._factor \ my_rhs; #::UmfpackLU{Float64,Int64}
	pause_advanced_timer(timer, "JULIA/ls_solve")
end

function ls_solve(solver::linear_solver_JULIA, my_rhs::AbstractArray, timer::class_advanced_timer)
    start_advanced_timer(timer, "JULIA/ls_solve")
    #Vector was added when updating from Julia 0.6.4 to Julia 0.7.0 by Fadi Hamad
    #This is becasue my_rhs in somecases is being a SparseVector instead of being array.
    #Due to the code in shur.jl jac_res = eval_jac_T_prod and hess_res = hess_product are returning SparseMatrix now instead of Array
    sol = solver._factor \ Vector(my_rhs);
    pause_advanced_timer(timer, "JULIA/ls_solve")
    return sol
end


text:
module ODBC

using Printf, Dates, UUIDs, Unicode, Random
using DecFP, DBInterface, Tables
export DBInterface

include("API.jl")
include("utils.jl")
include("dbinterface.jl")
include("load.jl")

"""
    ODBC.setdebug(debug::Bool=true, tracefile::String=joinpath(tempdir(), "odbc.log"))

Turn on ODBC library call tracing. This prints debug information to `tracefile` upon every
entry and exit from calls to the underlying ODBC library (unixODBC, iODBC, or odbc32).
Debugging can be turned off by passing `false`.

Note that setting tracing on/off requires resetting the ODBC environment, which means
any open statements/connections will be closed/invalid.

Also note that due to the persistent nature of ODBC config, setting tracing will
persist acrosss Julia sessions, i.e. if you turn tracing on, then quit julia and start again
tracing will still be on, and will stay on until explicitly turned off.

The iODBC driver manager supports passing `stderr` as the `tracefile`, which will
print all tracing information into the julia session/repl.
"""
function setdebug(debug::Bool=true, tracefile::String=joinpath(tempdir(), "odbc.log"))
    if debug
        API.setdebug(true, tracefile)
        @info "Enabled tracing of odbc library calls to $tracefile"
    else
        API.setdebug(false, tracefile)
        @info "Disabled tracing of odbc library calls"
    end
    return
end

"""
    ODBC.setunixODBC(; kw...)

Set the ODBC driver manager used to unixODBC. By default, ODBC.jl sets the `ODBCINI` and
`ODBCSYSINI` environment variables to the ODBC.jl-managed "scratch" location, but users may override
these (or provide additional environment variables) via `kw...` keyword arguments
this this function. The env variables will be set before allocating the ODBC environemnt.

Note that the odbc driver shared libraries can be "sticky" with regards to changing to
system configuration files. You may need to set a `OVERRIDE_ODBCJL_CONFIG` environment
variable before starting `julia` and running `import ODBC` to ensure that no environment
variables are changed by ODBC.jl itself.

While a unixODBC driver manager shared library is available for every platform, do note
that individual ODBC driver libraries may not be compatible with unixODBC on your system;
for example, if a driver library is built against iODBC, but unixODBC is the driver manager.
"""
setunixODBC(; kw...) = API.setunixODBC(; kw...)

"""
    ODBC.setiODBC(; kw...)

Set the ODBC driver manager used to iODBC. By default, ODBC.jl sets the `ODBCINI` and
`ODBCINSTINI` environment variables to the ODBC.jl-managed "scratch" location, but users may override
these (or provide additional environment variables) via `kw...` keyword arguments
this this function. The env variables will be set before allocating the ODBC environemnt.

Note that the odbc driver shared libraries can be "sticky" with regards to changing to
system configuration files. You may need to set a `OVERRIDE_ODBCJL_CONFIG` environment
variable before starting `julia` and running `import ODBC` to ensure that no environment
variables are changed by ODBC.jl itself.

While the iODBC driver manager shared library is available for non-windows platforms, do note
that individual ODBC driver libraries may not be compatible with unixODBC on your system;
for example, if a driver library is built against unixODBC, but iODBC is the driver manager.
"""
setiODBC(; kw...) = API.setiODBC(; kw...)

"""
    ODBC.setodbc32(; kw...)

Set the ODBC driver manager used to odbc32. On windows, ODBC.jl uses the system-wide
configurations for drivers and datasources. Drivers and datasources can still be added
via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`, but you must
have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
setodbc32(; kw...) = API.setodbc32(; kw...)

# driver/dsn management
"""
    ODBC.drivers() -> Dict

List installed ODBC drivers. The primary config location for installed drivers on non-windows platforms is
a reserved "scratch" space directory, i.e. an ODBC.jl-managed
location. Other system/user locations may also be checked (and are used by default on windows)
by the underlying ODBC driver manager, but for the most consistent results, aim to allow ODBC.jl to manage
installed drivers/datasources via `ODBC.addriver`, `ODBC.removedriver`, etc.

Note that the odbc driver shared libraries can be "sticky" with regards to changing to
system configuration files. You may need to set a `OVERRIDE_ODBCJL_CONFIG` environment
variable before starting `julia` and running `import ODBC` to ensure that no environment
variables are changed by ODBC.jl itself.

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
drivers() = API.getdrivers()

"""
    ODBC.dsns() -> Dict

List installed ODBC datasources. The primary config location for installed datasources on non-windows platforms is
a reserved "scratch" space directory, i.e. an ODBC.jl-managed
location. Other system/user locations may also be checked (and are by default on windows) by the underlying ODBC
driver manager, but for the most consistent results, aim to allow ODBC.jl to manage
installed drivers/datasources via `ODBC.adddsn`, `ODBC.removedsn`, etc.

Note that the odbc driver shared libraries can be "sticky" with regards to changing to
system configuration files. You may need to set a `OVERRIDE_ODBCJL_CONFIG` environment
variable before starting `julia` and running `import ODBC` to ensure that no environment
variables are changed by ODBC.jl itself.

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
dsns() = API.getdsns()

"""
    ODBC.adddriver(name, libpath; kw...)

Install a new ODBC driver. `name` is a user-provided "friendly" name to identify
the driver. `libpath` is the absolute path to the ODBC driver shared library.
Other key-value driver properties can be provided by the `kw...` keyword arguments.

This method is provided to try and provide the simplest/easiest/most consistent setup
experience for installing a new driver. Editing configuration files by hand is error-prone
and it's easy to miss adding something that is required.

While ODBC.jl supports all 3 major ODBC driver managers (unixODBC, iODBC, and odbc32),
be aware that most DBMS ODBC driver libraries are built against only one of the 3 and
can lead to compatibility issues if a different driver manager is used. This is mainly
an issue for driver libraries built against iODBC and then tried to use with unixODBC
or vice-versa.

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
adddriver(name, libpath; kw...) = API.adddriver(name, libpath; kw...)

"""
    ODBC.removedriver(name; removedsns::Bool=true)

Remove an installed ODBC driver by `name` (as returned from `ODBC.drivers()`).
`removedsns=true` also removes any datasources that were specified to use the driver.

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
removedriver(name; removedsns::Bool=true) = API.removedriver(name, removedsns)

"""
    ODBC.adddsn(name, driver; kw...)

Install a new ODBC datasource. `name` is a user-provided "friendly" name to identify
the datasource (dsn). `driver` is the "friendly" driver name that should be used to
connect to the datasource (valid driver options can be seen from `ODBC.drivers()`).
Additional connection key-value properties can be provided by the `kw...` keyword arguments.

Datasources can be connected by calling `DBInterface.connect(ODBC.Connection, dsn, user, pwd)`,
where `dsn` is the friendly datasource name, `user` is the username, and `pwd` is the password.

An alternative approach to installing datasources is to generate a valid "connection string"
that includes all connection properties in a single string passed to `DBInterface.connect`.
[www.connectionstrings.com](https://www.connectionstrings.com/) is a convenient resource
that provides connection string templates for various database systems.

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
adddsn(name, driver; kw...) = API.adddsn(name, driver; kw...)

"""
    ODBC.removedsn(name)

Remove an installed datasource by `name` (as returned from `ODBC.dsns()`).

On windows, ODBC.jl uses the system-wide configurations for drivers and datasources. Drivers and
datasources can still be added via `ODBC.adddriver`/`ODBC.removdriver` and `ODBC.adddsn`/`ODBC.removedsn`,
but you must have administrator privileges in the Julia session. This is accomplished easiest by pressing
CTRL then right-clicking on the terminal/Julia application and choosing "Run as administrator".
"""
removedsn(name) = API.removedsn(name)

end #ODBC module


text:
using CVCore
using Base.Test

import CVCore: mat_depth, mat_channel, maketype

@testset "Scalar{T}" begin
    # Eltype should be Float64
    @test eltype(Scalar(0)) == Float64
    @test eltype(Scalar(Int(0))) == Float64

    # with explicit type paramter
    for typ in [Int32, Int, Float32, Float64]
        @test eltype(Scalar{typ}(1)) == typ
    end
end

@testset "Mat{T,N}" begin
    mat = Mat{Float64}()
    @test isempty(mat)

    mat = Mat{UInt8}(10, 20)
    @test size(mat) == (10, 20)
    @test channels(mat) == 1
    @test rows(mat) == 10
    @test cols(mat) == 20
    @test total(mat) == length(mat)
    @test isContinuous(mat)
    @test elemSize(mat) == sizeof(UInt8)*channels(mat)
    @test dims(mat) == 2
    @test !isempty(mat)

    mat = Mat{Float32}(10, 20, 2)
    @test channels(mat) == 2
    @test size(mat) == (10, 20, 2)
    @test elemSize(mat) == sizeof(Float32)*channels(mat)
    @test dims(mat) == 2

    mat = Mat{Float64}(10, 20, 3, Scalar(1,0.5,0.3))
    @test channels(mat) == 3
    @test size(mat) == (10, 20, 3)
    @test elemSize(mat) == sizeof(Float64)*channels(mat)

    matclone = clone(mat)
    @test isa(matclone, typeof(mat))
    @test matclone == mat

    @testset "Array to cv::Mat conversion" begin
        arr = rand(Float64, 10, 2)
        m = Mat(arr)
        @test eltype(m) == eltype(arr)
        @test size(arr) == reverse(size(m))
    end

    mat = Mat{UInt8}(10,10, Scalar(0))
    @test all(mat .== 0)
    mat = Mat{UInt8}(10,10, Scalar(1))
    @test all(mat .== 1)
end

@testset "UMat{T,N}" begin
    umat = UMat{UInt8}(3, 3)
    mat = Mat(umat)
    @test isa(mat, Mat)
    @test size(umat) == (3, 3)
end

macro mat_op_testdata()
    #=
    2x12 Array{Float64,2}:
     1.0  3.0  5.0  7.0   9.0  11.0  13.0  15.0  17.0  19.0  21.0  23.0
     2.0  4.0  6.0  8.0  10.0  12.0  14.0  16.0  18.0  20.0  22.0  24.0
    =#
    esc(quote
        x = map(Float64, reshape([1:24;], 2,3*4))
        xt = x'
        m = Mat(xt)
    end)
end

# Since data is stored in different order in julia arrays and cv arrays
# .== simply doesn't work as expected. The function does manually iterate
# and compare elements.
function equal2d(x1, x2)
    @assert size(x1) == size(x2)
    n, m = size(x1)
    ret = true
    for i in 1:m
        for j in 1:n
            ret = ret && x1[j,i] == x2[j,i]
        end
    end
    return ret
end

@testset "Matrix operations" begin
    @testset "+" begin
        @mat_op_testdata
        ret = m + m
        retmat = Mat(ret)
        # @test all(retmat .== 2x)
        @test equal2d(retmat, 2x)
    end

    @testset "-" begin
        @mat_op_testdata
        ret = m - m
        retmat = Mat(ret)
        @test all(retmat .== 0)
    end

    @testset ".*" begin
        @mat_op_testdata
        # TODO
        # ret = m .* 5
        ret = broadcast(*, m, 5)
        @test isa(ret, MatExpr)
        retmat = Mat(ret)
        @test equal2d(retmat, 5x)
    end

    @testset "./" begin
        @mat_op_testdata
        # TODO
        # ret = m ./ 2
        ret = broadcast(/, m, 2)
        @test isa(ret, MatExpr)
        retmat = Mat(ret)
        expected = x ./ 2
        @test equal2d(retmat, expected)
        # @test all(retmat .== expected)
    end

    @testset "*" begin
        @mat_op_testdata
        @test Mat(m * m') == x * x'
    end

    @testset "transpose" begin
        @mat_op_testdata
        retexpr = m'
        retmat = Mat(retexpr)
        @test equal2d(retmat, x')
        # @test all(retmat .== x')
    end

    @testset "inv" begin
        @mat_op_testdata
        square_mat = m * m'
        #=
        2x2 Array{Float64,2}:
        2300.0  2444.0
        2444.0  2600.0
        =#
        inv_x = inv(x * xt)
        ret = inv(square_mat)
        inv_m = Mat(ret)
        @test inv_x ≈ inv_m
    end
end

@testset "Matrix channels" begin
    @test mat_channel(CV_32FC1) == 1
    @test mat_channel(CV_32FC2) == 2
    @test mat_channel(CV_32FC3) == 3
    @test mat_channel(CV_32FC4) == 4

    @test mat_channel(CV_8UC1) == 1
    @test mat_channel(CV_8UC2) == 2
    @test mat_channel(CV_8UC3) == 3
    @test mat_channel(CV_8UC4) == 4
end

@testset "Matrix depth" begin
    for flags in [
            CV_32FC1,
            CV_32FC2,
            CV_32FC3,
            CV_32FC4
            ]
        @test mat_depth(flags) == CV_32F
    end

    for flags in [
            CV_8UC1,
            CV_8UC2,
            CV_8UC3,
            CV_8UC4
            ]
        @test mat_depth(flags) == CV_8U
    end
end

@testset "maketype" begin
    @test maketype(CV_8U, 1) == CV_8UC1
    @test maketype(CV_8U, 2) == CV_8UC2
    @test maketype(CV_8U, 3) == CV_8UC3
    @test maketype(CV_8U, 4) == CV_8UC4

    @test maketype(CV_32F, 1) == CV_32FC1
    @test maketype(CV_32F, 2) == CV_32FC2
    @test maketype(CV_32F, 3) == CV_32FC3
    @test maketype(CV_32F, 4) == CV_32FC4
end

@testset "cv::TermCriteria" begin
    criteria = TermCriteria(CVCore.TERM_CRITERIA_EPS, 30, 0.001)
    @test isa(criteria, Cxx.CppValue)
end


text:
# # import some faster exponential for fast potential assembly
# try
#     import AppleAccelerate
#     function _fast_exponential(Lc, Rc, r)
#         c2 = Rc-r
#         AppleAccelerate.rec!(c2, c2)
#         c2 = Lc * c2
#         AppleAccelerate.exp!(c2, c2)
#         return c2
#     end
# catch
#     _fast_exponential(Lc, Rc, r) = exp( Lc ./ (Rc - r) )
# end

###########################################################################
## Implementation of some basic calculators using MatSciPy.NeighborList
###########################################################################

"""`simple_binsum` : this is a placeholder for a more general function,
`binsum`, which still needs to be written! Here, it is assumed that
`size(A, 1) = 3`, and  only summation along the second dimension is allowed.
"""
function simple_binsum{TI <: Integer, TF <: AbstractFloat}(i::Vector{TI},
                                                           A::Matrix{TF})
    if size(A, 1) != 3
        error("simple_binsum: need size(A,1) = 3")
    end
    if size(A, 2) != length(i)
        error("simple_binsum: need size(A,2) = length(i)")
    end
    B = zeros(TF, 3, maximum(i))
    for m = 1:size(A,1)
        # @inbounds @simd
        for n = 1:length(i)
            B[m, i[n]] = B[m,i[n]] + A[m, n]
        end
    end
    return B
end

function simple_binsum{TI <: Integer, TF <: AbstractFloat}(i::Vector{TI},
                                                           A::Vector{TF})
    if length(A) != length(i)
        error("simple_binsum: need length(A) = length(i)")
    end
    B = zeros(TF, maximum(i))
    # this ought to be a SIMD loop. but that gives me a wrong answer! why?!
    for n = 1:length(i)
        B[i[n]] += + A[n]
    end
    return B
end

"""`PairCalculator` : basic calculator for pair potentials.
"""
type PairCalculator <: AbstractCalculator
    pp::PairPotential
end

import Potentials.cutoff
cutoff(calc::PairCalculator) = cutoff(calc.pp)

function potential_energy(at::ASEAtoms, calc::PairCalculator)
    r = neighbour_list(at, "d", cutoff(calc))
    return sum( calc.pp(r) )
end

function potential_energy_d(at::ASEAtoms, calc::PairCalculator)
    i, r, R = neighbour_list(at, "idD", cutoff(calc))
    return - 2.0 * simple_binsum(i, @GRAD calc.pp(r, R') )
end

forces(at::ASEAtoms, calc::PairCalculator) = - potential_energy_d(at, calc)

"`EAMCalculator` : basic calculator using the `EAMPotential` type"
type EAMCalculator <: AbstractCalculator
    p::EAMPotential
end

cutoff(calc::EAMCalculator) = max(cutoff(calc.p.V), cutoff(calc.p.rho))

function potential_energy(at::ASEAtoms, calc::EAMCalculator)
    i, r = neighbour_list(at, "id", cutoff(calc))
    return ( sum(calc.p.V(r))
             + sum( calc.p.embed( simple_binsum( i, calc.p.rho(r) ) ) ) )
end

function potential_energy_d(at::ASEAtoms, calc::EAMCalculator)
    i, j, r, R = neighbour_list(at, "ijdD", cutoff(calc))
    # pair potential component
    G = - 2.0 * simple_binsum(i, @GRAD calc.p.V(r, R'))
    # EAM component
    dF = @D calc.p.embed( simple_binsum(i, calc.p.rho(r)) )
    dF_drho = dF[i]' .* (@GRAD calc.p.rho(r, R'))
    G += simple_binsum(j, dF_drho) - simple_binsum(i, dF_drho)
    return G
end

# ###########################################################################
# ## Some Calculators Optimised for use with MatSciPy.NeighbourList
# ###########################################################################

# """
# `lennardjones_old(at::ASEAtoms, nlist::NeighbourList;
#                       r0=1.0, e0=1.0, quantities="EG")`

# A fast LJ assembly which exploits the `@simd` and `@inbounds` macros as well
# as the specific structure of `MatSciPy.NeighbourList`.

# **This has not yet been tested for correctness!**

# **This version does not have a cutoff radius!**
# """
# function lennardjones_old(at::ASEAtoms, nlist::NeighbourList;
#                       r0=1.0, e0=1.0, quantities="EG")
#     update!(nlist, at)
#     E = 0.0
#     r = nlist.Q['d']::Vector{Float64}
#     R = nlist.Q['D']::Array{Float64,2}
#     i = nlist.Q['i']::Vector{Int32}
#     t = Vector{Float64}(length(r))

#     @simd for n = 1:length(t) @inbounds begin
#         t[n] = r0 / r[n]
#         t[n] = t[n]*t[n]*t[n]
#         t[n] = t[n]*t[n]
#         E = E + t[n]*(t[n]-2.0)
#     end end
#     E *= e0

#     if 'G' in quantities
#         G = zeros(3, i[end])
#         @simd for n = 1:length(t) @inbounds begin
#             t[n] = e0 * 24.0 * t[n]*(t[n]-1.0) / (r[n]*r[n])
#             G[1,i[n]] = G[1,i[n]] + t[n] * R[n,1]
#             G[2,i[n]] = G[2,i[n]] + t[n] * R[n,2]
#             G[3,i[n]] = G[3,i[n]] + t[n] * R[n,3]
#         end end
#     end

#     ret = ()
#     for c in quantities
#         if c == 'E'
#             ret = tuple(ret..., E)
#         elseif c == 'G'
#             ret = tuple(ret..., G)
#         end
#     end
#     return ret
# end

# """
# `lennardjones(at::ASEAtoms, nlist::NeighbourList;
#                       r0=1.0, e0=1.0, Rc=2.7, Lc=1.0, quantities="EG")`

# A fast LJ assembly which exploits the specific structure of the
#  `MatSciPy.NeighbourList` to use  `@simd`, `@inbounds`, `@fastmath` macros
# and the `AppleAccalerate` package.

# **This has not yet been tested for correctness!**
# """
# function lennardjones(at::ASEAtoms, nlist::NeighbourList;
#                       r0=1.0, e0=1.0, Rc=2.7, Lc = 1.0,  quantities="EG")
#     update!(nlist, at)
#     E = 0.0
#     r = nlist.Q['d']::Vector{Float64}
#     R = nlist.Q['D']::Array{Float64,2}
#     i = nlist.Q['i']::Vector{Int32}
#     t = Vector{Float64}(length(r))
#     c1 = Vector{Float64}(length(r))
#     # c2 = Vector{Float64}(length(r))
#     Rc_ = Rc-0.001

#     # c2 = exp(Lc ./ (Rc-r))
#     c2 = _fast_exponential(Lc, Rc, r)

#     @fastmath @inbounds @simd for n = 1:length(t)
#         # c2[n] = exp( Lc / (Rc-r[n]) )
#         c1[n] = (1.0 / (1.0 + c2[n]))  * (r[n] >= Rc_)
#         t[n] = r0 / r[n]
#         t[n] = t[n]*t[n]*t[n]
#         t[n] = t[n]*t[n]
#         E = E + t[n] * (t[n]-2.0) * c1[n]
#     end
#     E *= e0

#     if 'G' in quantities
#         G = zeros(3, i[end])
#         @fastmath @inbounds @simd for n = 1:length(t)
#             c2[n] = -Lc * c2[n] * c1[n] / ( (Rc-r[n])*(Rc-r[n]) )
#             c2[n] = c2[n] * t[n] * (t[n]-2.0)
#             t[n] = e0 * 24.0 * t[n] * (t[n]-1.0) / (r[n]*r[n])
#             t[n] = t[n] * c1[n] + c2[n]
#             G[1,i[n]] = G[1,i[n]] + t[n] * R[n,1]
#             G[2,i[n]] = G[2,i[n]] + t[n] * R[n,2]
#             G[3,i[n]] = G[3,i[n]] + t[n] * R[n,3]
#         end
#     end

#     ret = ()
#     for c in quantities
#         if c == 'E'
#             ret = tuple(ret..., E)
#         elseif c == 'G'
#             ret = tuple(ret..., G)
#         end
#     end
#     return ret
# end


text:
include("helpers.jl")
include("eglob.jl")


text:
function spiral!(r::Robot)
    keeper = 1
    z = 0
    while z != -1
        for i in [0, 1]
            schetchik = keeper
            while schetchik != 0
                if ismarker(r)
                    break
                end
                move!(r, HorizonSide(i))
                schetchik -= 1
            end
            if ismarker(r)
                break
            end
            schetchik = keeper
        end
        if ismarker(r)
            break
        end
        keeper += 1
        for i in [2, 3]
            schetchik = keeper
            while schetchik != 0
                if ismarker(r)
                    break
                end
                move!(r, HorizonSide(i))
                schetchik -= 1
            end
            if ismarker(r)
                break
            end
            schetchik = keeper
        end
        keeper += 1
    end
end


text:
function next(side, kol)
    return HorizonSide((Int(side) + kol) % 4)
end

function moveAndKol(side, kol)
    flag = false
    count = 0
    while !isborder(robot, side)
        if (isborder(robot, next(side, kol)) && !flag)
            flag = true
            count += 1
        elseif !isborder(robot, next(side, kol))
            flag = false
        end
        move!(robot, side)
    end
    return count
end

function f()
    side = Ost
    n = 1
    count = 0
    while !isborder(robot, Nord)
        count += moveAndKol(side, n)
        side = next(side, 2)
        n = (n + 2) % 4
        move!(robot, Nord)
    end
    println(count)
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
module OneDNN

export Memory
export Conv, BatchNorm, Pooling

# stdlib
import LinearAlgebra
import Random

# temp local teps
import CachedArrays

# deps
import ChainRulesCore
import Flux
import MacroTools
import Polyester
import UnPack: @unpack
import Zygote
import ZygoteRules: ZygoteRules, _pullback, AContext, literal_getproperty, literal_getfield

function pullback_for_default_literal_getproperty(cx::AContext, x, ::Val{f}) where {f}
    return _pullback(cx, literal_getfield, x, Val{f}())
end

# Experimental - should probably not set this to `true`.
const SIMILAR_FOR_SCRATCHPAD = false

#####
##### deps
#####

# Include auto-generated wrapper for OneDNN
include("lib/lib.jl")
include("wrap.jl")
include("tiledarrays.jl")
using .TiledArrays: TiledArrays
include("utils.jl")
include("types.jl")
include("tiled.jl")
include("memory.jl")
include("bf16.jl")

# ops
include("ops/simple.jl")
include("ops/matmul.jl")
include("ops/batchnorm.jl")
include("ops/concat.jl")
include("ops/pool.jl")
include("ops/innerproduct.jl")
include("ops/convolution.jl")

# include("placeholder.jl")

# Just create a global engine and stream for everything to use for now.
const GLOBAL_ENGINE = Ref{Engine}()
const GLOBAL_STREAM = Ref{Stream}()
const GLOBAL_ATTRIBUTES = Ref{Attributes}()
const GLOBAL_THREADPOOL = Ref{Any}()

_get_in_parallel() = (Threads.threadid() != 1)
function _parallel_for(n::Cint, f::Ptr{Cvoid})
    #Threads.@threads for i in Base.OneTo(n)
    Polyester.@batch per = thread for i in Base.OneTo(n)
        Wrap.call_opaque(f, i - 1, n)
    end
end

# Initialize the default Engine and Stream.
function __init__()
    ### Engine
    engine = Engine(Lib.dnnl_cpu)
    GLOBAL_ENGINE[] = engine

    ### Thread Pool
    # The functions we create do not capture any variables, thus no need to preserve them
    # somewhere.
    threadpool = Wrap.construct_threadpool(
        @cfunction(_get_in_parallel, Bool, ()),
        @cfunction(_parallel_for, Cvoid, (Cint, Ptr{Cvoid})),
        Threads.nthreads(),
    )
    GLOBAL_THREADPOOL[] = threadpool

    ### Stream
    stream = Stream(engine)
    @apicall dnnl_threadpool_interop_stream_create(stream, engine, threadpool.cpp_object)
    GLOBAL_STREAM[] = stream

    ### Default Attributes
    GLOBAL_ATTRIBUTES[] = Attributes()

    return nothing
end

noattributes() = GLOBAL_ATTRIBUTES[]
global_engine() = GLOBAL_ENGINE[]
global_stream() = GLOBAL_STREAM[]

#####
##### Flux compat
#####

# # TODO: Turns out that updating is better parallelized across the layer dimension than the
# # within each parameter itself
function Flux.Optimise.update!(
    o::Flux.Optimise.Descent, x::Memory{T,N}, Δ::Memory{U,N}
) where {T,U,N}
    # Make sure both objects have the same memory layout.
    # TODO: Handle different types correctly.
    # mx = memorydesc(x)
    # mΔ = memorydesc(Δ)
    # if mx != mΔ
    #     error("Incompatible memory formats!")
    #     # # Enter the realm of the type unstable!
    #     # sz = logicalsize(mx, Val(N))
    #     # indexer_x = TiledArrays.TiledIndexer{layout(mx)}(sz, padded_size(mx, Val(N)))
    #     # indexer_Δ = TiledArrays.TiledIndexer{layout(mΔ)}(
    #     #     logicalsize(mΔ, Val(N)), padded_size(mΔ, Val(N))
    #     # )
    #     # update_typed!(o, parent(x), indexer_x, parent(Δ), indexer_Δ, CartesianIndices(sz))
    #     # return nothing
    # end

    eta = convert(Float32, o.eta)
    xa = vec(parent(x))
    Δa = vec(parent(Δ))
    xa .-= (eta .* Δa)
    return nothing
end

function Flux.Optimise.update!(
    o::Flux.Optimise.Descent, x::Memory{T,N}, ix, y::Memory{U,N}, iy
) where {T,U,N}
    px = parent(x)
    py = parent(y)
    eta = convert(Float32, o.eta)
    for i in eachindex(ix, iy)
        @inbounds(px[ix[i]] -= eta * py[iy[i]])
    end
    return nothing
end

#####
##### Mirrored
#####

function Flux.Optimise.update!(
    o::Flux.Optimise.Descent, x::Mirrored, Δ::Memory{T,N}
) where {T,N}
    eta = convert(Float32, o.eta)
    px = parent(x)
    pΔ = parent(Δ)
    for i in eachindex(pΔ)
        @inbounds px[i] -= eta * pΔ[i]
    end
    return nothing
end

function Flux.Optimise.update!(
    o::Flux.Optimise.Descent, x::Mirrored, ix, y::Memory{T,N}, iy
) where {T,N}
    px = parent(x)
    py = parent(y)
    eta = convert(Float32, o.eta)
    for i in eachindex(ix, iy)
        @inbounds(px[ix[i]] -= eta * py[iy[i]])
    end
    return nothing
end

#####
##### Zygote Compat
#####

Zygote.accum(x::AbstractArray, y::Memory) = +(Memory(x), y)
Zygote.accum(x::Memory, y::AbstractArray) = +(x, Memory(y))
Zygote.accum(x::Memory, y::Memory) = +(x, y)

# Special cases
function Zygote.accum(_x::SubArray, y::Memory)
    x = maybe_reorder(memorydesc(y), Memory(_x))
    binary!(+, x, x, y)
    return x
end

# function update_typed!(
#     o::Flux.Optimise.Descent,
#     x,
#     indexer_x::TiledArrays.TiledIndexer,
#     y,
#     indexer_y::TiledArrays.TiledIndexer,
#     iter,
# )
#     for i in iter
#         ix = TiledArrays.genindex(indexer_x, Tuple(i))
#         iy = TiledArrays.genindex(indexer_y, Tuple(i))
#         x[ix] -= o.eta * y[iy]
#     end
#     return nothing
# end

# function update_typed!(o::Flux.Optimise.Descent, x::Memory, Δ::Memory)
#     return x .= x .- (o.eta .* Δ)
# end

# # # Apply the negative here so we can just add together in `update!`.
# # # This is because it appears that OneDNN is lacking a binary `-`
# # Flux.Optimise.apply!(o::Flux.Optimise.Descent, x, Δ::Memory) = linear!(Δ, -o.eta)
# #
# # # Expect `Memory` objects to already be negated from the `apply!` step.
# # function Flux.Optimise.update!(o::Flux.Optimise.Descent, x::Memory, Δ::Memory)
# #     return binary!(+, x, Flux.Optimise.apply!(o, x, Δ))
# end

end # module


text:
struct WilsonClover
    cSW::Float64
    σ::Array{Array{ComplexF64,2},2}

    function WilsonClover(cSW)
        σ = make_σμν()
        return new(cSW,σ)
    end
end

function make_σμν()
    σ = Array{Array{ComplexF64,2},2}(undef,4,4)
    for μ=1:4
        γμ = γ_all[:,:,μ]
        for ν=1:4
            γν = γ_all[:,:,ν]
            σ[μ,ν] =  (γμ*γν .- γν*γμ )*(1/2)
        end
    end
    return σ
end

function Wclover!(y,U,x,A) #clover operator

end

struct WilsonClover_misc
    clover_coefficient::Float64
    internal_flags::Array{Bool,1}
    inn_table::Array{Int64,3}
    _ftmp_vectors::Array{Array{ComplexF64,3},1}
    _is1::Array{Int64,1}
    _is2::Array{Int64,1}

    function WilsonClover(U::Array{<: AbstractGaugefields{NC,Dim},1},x,clover_coefficient) where {NC,Dim}
        _,_,NN... = size(U[1])
        NV = prod(NN)
        inn_table= zeros(Int64,NV,4,2)
        internal_flags = zeros(Bool,2)
        _ftmp_vectors = Array{Array{ComplexF64,3},1}(undef,6)
        for i=1:6
            _ftmp_vectors[i] = zeros(ComplexF64,NC,NV,4)
        end

        _is1 = zeros(Int64,NV)
        _is2 = zeros(Int64,NV)

        return new(clover_coefficient,internal_flags,inn_table,_ftmp_vectors,_is1,_is2)
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
ACASX Study comparing performance of GP, MC, MCTS, GE, CE.
Single-threaded versions are used for more stable comparison.
Main entry: study_main()
"""
module ACASX_Compare

export combine_sweep_logs, combine_mc_logs, combine_mcts_logs,
    combine_ge_logs, combine_gp_logs, combine_ce_logs, combine_logs
export master_log, master_plot
export combine_and_plot

import Compat: ASCIIString, UTF8String
using GrammarExpts
using ExprSearch: GP, MC, MCTS, GE, CE
using ACASX_GP, ACASX_GE, ACASX_MC, ACASX_MCTS, ACASX_CE

using RLESUtils, Loggers, MathUtils, LatexUtils, Sweeper, LogJoiner
using DataFrames
using PGFPlots, TikzPictures

const STUDYNAME = "ACASX_Compare"
const MC_NAME = "ACASX_MC"
const MCTS_NAME = "ACASX_MCTS"
const GE_NAME = "ACASX_GE"
const GP_NAME = "ACASX_GP"
const CE_NAME = "ACASX_CE"

const CONFIGDIR = joinpath(dirname(@__FILE__), "..", "config")
const RESULTDIR = joinpath(dirname(@__FILE__), "..", "..", "..", "results")

const MASTERLOG_FILE = joinpath(RESULTDIR, STUDYNAME, "masterlog.csv.gz")
const PLOTLOG_FILE =  joinpath(RESULTDIR, STUDYNAME, "plotlog.csv.gz")
const PLOTFILEROOT = joinpath(RESULTDIR, STUDYNAME, "plots")

resultpath(dir::ASCIIString="") = joinpath(RESULTDIR, dir)
studypath(dir::ASCIIString="") = joinpath(RESULTDIR, STUDYNAME, dir)

function combine_mc_logs()
    dir = studypath(MC_NAME)
    logjoin(dir, "acasx_mc_log.txt", ["current_best", "elapsed_cpu_s"],
        joinpath(dir, "subdirjoined"))
end
function combine_mcts_logs()
    dir = studypath(MCTS_NAME)
    logjoin(dir, "acasx_mcts_log.txt", ["current_best", "elapsed_cpu_s"],
        joinpath(dir, "subdirjoined"))
end
function combine_ge_logs()
    dir = studypath(GE_NAME)
    logjoin(dir, "acasx_ge_log.txt", ["current_best", "elapsed_cpu_s"],
        joinpath(dir, "subdirjoined"))
end
function combine_gp_logs()
    dir = studypath(GP_NAME)
    logjoin(dir, "acasx_gp_log.txt", ["current_best", "elapsed_cpu_s"],
        joinpath(dir, "subdirjoined"))
end
function combine_ce_logs()
    dir = studypath(CE_NAME)
    logjoin(dir, "acasx_ce_log.txt", ["current_best", "elapsed_cpu_s"],
        joinpath(dir, "subdirjoined"))
end
function combine_sweep_logs()
    dir = studypath()
    logjoin(dir, "sweeper_log.txt", ["result"], joinpath(dir, "sweepjoined"))
end

#TODO: clean this up...
function master_log(; b_mc=true, b_mcts=true, b_ge=true, b_gp=true, b_ce=true)
    masterlog = DataFrame([Int64, Float64, Float64, UTF8String, ASCIIString, UTF8String],
        [:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name], 0)

    #MC
    if b_mc
        dir = studypath(MC_NAME)
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))
        D = join(logs["elapsed_cpu_s"], logs["current_best"], on=[:iter, :name])
        D[:algorithm] = fill("MC", nrow(D))
        rename!(D, :iter, :nevals)
        append!(masterlog, D[[:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name]])
    end

    #MCTS
    if b_mcts
        dir = studypath(MCTS_NAME)
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))
        D = join(logs["elapsed_cpu_s"], logs["current_best"], on=[:iter, :name])
        D[:algorithm] = fill("MCTS", nrow(D))
        rename!(D, :iter, :nevals)
        append!(masterlog, D[[:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name]])
    end

    #GE
    if b_ge
        dir = studypath(GE_NAME)
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))
        D = join(logs["elapsed_cpu_s"], logs["current_best"], on=[:nevals, :name])
        D[:algorithm] = fill("GE", nrow(D))
        append!(masterlog, D[[:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name]])
    end

    #GP
    if b_gp
        dir = studypath(GP_NAME)
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))
        D = join(logs["elapsed_cpu_s"], logs["current_best"], on=[:nevals, :name])
        D[:algorithm] = fill("GP", nrow(D))
        append!(masterlog, D[[:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name]])
    end

    #CE
    if b_ce
        dir = studypath(CE_NAME)
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))
        D = join(logs["elapsed_cpu_s"], logs["current_best"], on=[:nevals, :name])
        D[:algorithm] = fill("CE", nrow(D))
        append!(masterlog, D[[:nevals, :elapsed_cpu_s, :fitness, :expr, :algorithm, :name]])
    end

    writetable(MASTERLOG_FILE, masterlog)
    masterlog
end

master_plot(; kwargs...) = master_plot(readtable(MASTERLOG_FILE); kwargs...)

"""
Subsamples the collected data at 'subsample' rate to plot at a lower rate than collected
"""
function master_plot(masterlog::DataFrame; subsample::Int64=25000)
    D = masterlog

    #aggregate over seed
    D = aggregate(D[[:nevals, :elapsed_cpu_s, :fitness, :algorithm]], [:nevals, :algorithm],
        [mean, std, length, SEM])
    D = D[rem(D[:nevals], subsample) .== 0, :] #subsample data

    #workaround for naming in julia 0.4
    rename!(D, Symbol("fitness_MathUtils.SEM"), :fitness_SEM)
    rename!(D, Symbol("elapsed_cpu_s_MathUtils.SEM"), :elapsed_cpu_s_SEM)

    writetable(PLOTLOG_FILE, D)

    td = TikzDocument()
    algo_names = unique(D[:algorithm])
    n_algos = length(algo_names)

    #nevals_vs_fitness
    plotarray = Plots.Plot[]
    for i = 1:n_algos
        D1 = D[D[:algorithm].==algo_names[i], [:nevals, :fitness_mean, :fitness_SEM]]
        push!(plotarray, Plots.Linear(D1[:nevals], D1[:fitness_mean],
            errorBars=ErrorBars(; y=D1[:fitness_SEM]),
            legendentry=escape_latex(algo_names[i])))
    end
    tp = PGFPlots.plot(Axis(plotarray, xlabel="Number of Evaluations", ylabel="Fitness",
        title="Fitness vs. Number of Evaluations", legendPos="north east"))
    push!(td, tp)

    #nevals_vs_elapsed_cpu
    empty!(plotarray)
    for i = 1:n_algos
        D1 = D[D[:algorithm].==algo_names[i], [:nevals, :elapsed_cpu_s_mean,
            :elapsed_cpu_s_SEM]]
        push!(plotarray, Plots.Linear(D1[:nevals], D1[:elapsed_cpu_s_mean],
            errorBars=ErrorBars(;y=D1[:elapsed_cpu_s_SEM]),
            legendentry=escape_latex(algo_names[i])))
    end
    tp = PGFPlots.plot(Axis(plotarray, xlabel="Number of Evaluations", ylabel="Elapsed CPU Time (s)",
        title="Elapsed CPU Time vs. Number of Evaluations", legendPos="north west"))
    push!(td, tp)

    #elapsed_cpu_vs_fitness
    empty!(plotarray)
    for i = 1:n_algos
        D1 = D[D[:algorithm].==algo_names[i], [:elapsed_cpu_s_mean,
            :elapsed_cpu_s_SEM, :fitness_mean, :fitness_SEM]]
        push!(plotarray, Plots.Linear(D1[:elapsed_cpu_s_mean], D1[:fitness_mean],
            errorBars=ErrorBars(;x=D1[:elapsed_cpu_s_SEM], y=D1[:fitness_SEM]),
            legendentry=escape_latex(algo_names[i])))
    end
    tp = PGFPlots.plot(Axis(plotarray, xlabel="Elapsed CPU Time (s)", ylabel="Fitness",
        title="Fitness vs. Elapsed CPU Time", legendPos="north east"))
    push!(td, tp)

    save(PDF(PLOTFILEROOT * ".pdf"), td)
    save(TEX(PLOTFILEROOT * ".tex"), td)
end

function combine_and_plot(; subsample=25000, b_mc=true, b_mcts=true, b_ge=true, b_gp=true, b_ce=true)
    b_ge && combine_ge_logs()
    b_gp && combine_gp_logs()
    b_mc && combine_mc_logs()
    b_mcts && combine_mcts_logs()
    b_ce && combine_ce_logs()
    ml = master_log(b_mc=b_mc, b_mcts=b_mcts, b_ge=b_ge, b_gp=b_gp, b_ce=b_ce)
    master_plot(ml; subsample=subsample)
end

end #module


text:
import RecipesBase

using RecipesBase: @recipe

"""
    SplineFunctionWrapper{S<:Spline, D<:Derivative, OnlySupport} <: Function

A function type that wraps a spline of type `S`. (not exported)

Calling an object of this type evaluates the derivative `D` of the spline. The parameter
`OnlySupport` determines the value returned when evaluating the function at a point outside
of the [`support`](@ref) of the spline: If `OnlySupport === true`, a value equal to `NaN`
will be returned; if `OnlySupport === false`, a value equal to zero is returned.
"""
struct SplineFunctionWrapper{S<:Spline,D<:Derivative,OnlySupport} <: Function
    spline::S

    function SplineFunctionWrapper{S,D,OnlySupport}(spline) where {S<:Spline,D<:Derivative,OnlySupport}
        OnlySupport isa Bool || throw(ArgumentError("type parameter OnlySupport must be a Bool."))
        new(spline)
    end
end

SplineFunctionWrapper(s::Spline, d::Derivative, onlysupport::Bool=true) =
    SplineFunctionWrapper{typeof(s),typeof(d),onlysupport}(s)

SplineFunctionWrapper(s::Spline, onlysupport::Bool=true) =
    SplineFunctionWrapper(s, NoDerivative(), onlysupport)

function (f::SplineFunctionWrapper{S,D,true})(x) where {S,D}
    a,b = support(f.spline)
    val = f.spline(x, D())
    a ≤ x ≤ b ? val : convert(float(typeof(val)), NaN)
end
(f::SplineFunctionWrapper{S,D,false})(x) where {S,D} = f.spline(x, D())

support(s::SplineFunctionWrapper) = support(s.spline)

"""
    Function(spline::Spline, [deriv::Derivative], [onlysupport=true])

Wrap `spline` in an object that is a subtype of `Function`. Calling the returned function
with a single argument `x` will evaluate the spline (or one of its derivatives as specified
by `deriv`) at `x`.

If the optional argument `onlysupport` is set to `true` (the default), the returned function
will return `NaN` if evaluated outside of the [`support`](@ref) of `spline`. If
`onlysupport` is set to `false`, it will return zero there (as does calling the `spline`
directly).

Note that a `Spline` can be called without wrapping them as described here, although they
are not a subtype of `Function`. Wrapping a spline in a `Function` object is mainly intended
to aid in plotting, which is the rationale behind the `onlysupport=true` default: when using
[`Plots.jl`](https://github.com/JuliaPlots/Plots.jl), this will cause the spline to not be
drawn outside of its support.

# Examples

```jldoctest
julia> spline = approximate(sin, BSplineBasis(5, 0:5)); # create a Spline

julia> f = Function(spline, false); # f is zero outside of the interval [0,5]

julia> g = Function(spline, Derivative(1)); # g is NaN outside of the interval [0,5]

julia> f(1.5) == spline(1.5)
true

julia> g(1.5) == spline(1.5, Derivative(1))
true

julia> f(-1), g(-1)
(0.0, NaN)
```
"""
Base.Function(x::Spline, args...) = SplineFunctionWrapper(x, args...)

# Type recipe for wrapping a `Spline` in a `Function`
@recipe f(::Type{S}, s::S) where {S<:Spline} = Function(s, true)

# One argument
@recipe f(s::SplineFunctionWrapper) = (s, support(s)...)
@recipe f(s::Vector{S}) where {S<:Spline} = (Function.(s, true), joinsupport(s)...)
@recipe f(s::Vector{S}) where {S<:SplineFunctionWrapper} = (s, joinsupport(s)...)

joinsupport(splines) = mapreduce(support, _joinsupport, splines)
_joinsupport((xmin, xmax), (ymin, ymax)) = (min(xmin, ymin), max(xmax, ymax))

# Two arguments
@recipe f(x, s::Vector{S}) where {S<:Spline} = (x, Function.(s, true))
@recipe f(s::Vector{S}, x) where {S<:Spline} = (x, Function.(s, true))

# Three arguments
@recipe f(s::Vector{S}, xmin::Number, xmax::Number) where {S<:Spline} =
    (Function.(s, true), xmin, xmax)

# BSplineBasis
@recipe function f(basis::BSplineBasis, pointsperspline::Integer=100)
    label --> permutedims(eachindex(basis))
    xs = Matrix{Float64}(undef, pointsperspline, length(basis))
    ys = similar(xs)
    workspace = Vector{Float64}(undef, order(basis))
    for (j, bspline) = enumerate(basis)
        xmin, xmax = support(bspline)
        for (i, x) = enumerate(range(xmin, stop=xmax, length=pointsperspline))
            xs[i,j] = x
            ys[i,j] = bspline(x, workspace=workspace)
        end
    end
    xs, ys
end


text:
abstract type GraphGenerator end

Base.length(g::GraphGenerator) = g.count

LightGraphs.nv(g::GraphGenerator) = g.nodes

params(g::GraphGenerator) = tuple()

Base.iterate(g::GraphGenerator) = if length(g) > 0
    graph = generate(g)
    while !is_connected(graph)
        graph = generate(g)
    end
    graph, length(g)
else
    nothing
end

Base.iterate(g::GraphGenerator, i::Int) = if i > 1
    graph = generate(g)
    while !is_connected(graph)
        graph = generate(g)
    end
    graph, i - 1
else
    nothing
end

struct CycleGraphGenerator <: GraphGenerator
    count::Int
    nodes::Int
end
CycleGraphGenerator(; count::Int=1, nodes::Int=1, kwargs...) = CycleGraphGenerator(count, nodes)

generate(c::CycleGraphGenerator) = cycle_graph(c.nodes)

struct WheelGraphGenerator <: GraphGenerator
    count::Int
    nodes::Int
end
WheelGraphGenerator(; count::Int=1, nodes::Int=1, kwargs...) = WheelGraphGenerator(count, nodes)

generate(w::WheelGraphGenerator) = wheel_graph(w.nodes)

struct StarGraphGenerator <: GraphGenerator
    count::Int
    nodes::Int
end
StarGraphGenerator(; count::Int=1, nodes::Int=1, kwargs...) = StarGraphGenerator(count, nodes)

generate(s::StarGraphGenerator) = star_graph(s.nodes)

struct GridGraphGenerator <: GraphGenerator
    count::Int
    dims::Vector{Int}
    periodic::Bool
end
GridGraphGenerator(; count::Int=1, dims::Vector{Int}=[1], periodic=true, kwargs...) = GridGraphGenerator(count, dims, periodic)

LightGraphs.nv(g::GridGraphGenerator) = prod(g.dims)

generate(g::GridGraphGenerator) = grid(g.dims; periodic=g.periodic)

struct BarabasiAlbertGenerator <: GraphGenerator
    count::Int
    nodes::Int
    k::Int
end
BarabasiAlbertGenerator(; count::Int=1, nodes::Int=1, k::Int=1, kwargs...) = BarabasiAlbertGenerator(count, nodes, k)

params(ba::BarabasiAlbertGenerator) = (; k=ba.k)

generate(ba::BarabasiAlbertGenerator) = barabasi_albert(ba.nodes, ba.k)

struct ErdosRenyiGenerator <: GraphGenerator
    count::Int
    nodes::Int
    p::Float64
end
ErdosRenyiGenerator(; count::Int=1, nodes::Int=1, p::Float64=0.5, kwargs...) = ErdosRenyiGenerator(count, nodes, p)

params(er::ErdosRenyiGenerator) = (; p=er.p)

generate(er::ErdosRenyiGenerator) = erdos_renyi(er.nodes, er.p)


text:
using Fredholm
using Test
using QuadGK, Random
using LinearAlgebra: norm

include("fredholm.jl")


text:
using Test
# using JSON3

using ScottishTaxBenefitModel
using .STBParameters:
    IncomeTaxSys,weeklyise!,annualise!,load_file,load_file!,TaxBenefitSystem
using .Utils
using .GeneralTaxComponents: WEEKS_PER_YEAR
using .ExampleHelpers

@testset "IT Parameter Tests" begin
    itsysdir :: IncomeTaxSys = get_default_it_system( year=2019, scotland=true)
    itsyseng :: IncomeTaxSys = get_default_it_system( year=2019, scotland=false)
    @test itsyseng.non_savings_rates[1]≈0.20
    @test itsyseng.non_savings_thresholds[2]≈150_000/WEEKS_PER_YEAR
    it = IncomeTaxSys{Float64}()
    itweekly = deepcopy(it)
    @test itweekly.savings_rates == it.savings_rates
    weeklyise!( itweekly )
    annualise!( itweekly )
    @test isapprox( itweekly.non_savings_thresholds, it.non_savings_thresholds, rtol=0.00001 )
    @test itweekly.mca_minimum ≈ it.mca_minimum

    sys_scot = get_system( scotland=true )
    @show sys_scot
    @show sys_scot.it
    @show sys_scot.ni

    sys2 = load_file( "../params/sys_2021_22.jl" )
    weeklyise!( sys2 )
    @test sys2.it.personal_allowance ≈ 12_570/WEEKS_PER_YEAR
    @test sys2.lmt.working_tax_credit.age_50_plus_30_hrs == 2_030.00/WEEKS_PER_YEAR
    @test sys2.child_limits.max_children == 2
    @show sys2

end # example 1

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Toy model from paper
# include("../src/SFFM.jl")
# using LinearAlgebra, Plots

## Define 1st model
T = [-1.0 1.0; 1.0 -1.0]
C = [1.0; -1.0]
r = (
    r = function (x)
        [abs.((x .<= 1)-(x .> 1)) abs.(((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0))]
    end,
    R = function (x)
        [x.*abs.(((x .<= 1)-(x .> 1))) x.*abs.((((x .> 1))-2*(x .<= 1).*(x .> 0)-(x.==0)))]
    end,
)

Bounds = [0 30; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = Bounds)

## Define mesh
Δ = 5
Nodes = collect(Bounds[1,1]:Δ:Bounds[1,2])
nBases = 3
Basis = "lagrange"
mesh = SFFM.DGMesh(model, Nodes, nBases, Basis=Basis)

## Make matrices
All = SFFM.MakeAll( model, mesh, approxType = "projection")
proj = All
interp = SFFM.MakeAll(model,mesh,approxType="interpolation")

vikramD = copy(interp.D["++"]())
vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:] =
    repeat(All.Matrices.Local.V.w, SFFM.NIntervals(mesh)*SFFM.NPhases(model), 1) .* vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:]

## sims fro gound truth
x₀ = 15
NSim = 50000
IC = (φ = ones(Int, NSim), X = x₀ .* ones(NSim), Y = zeros(NSim))
y = 10
@time sims =
    SFFM.SimSFFM( model, SFFM.FirstExitY(-Inf, y), IC)

## DG stationary dist
# construct initial condition
theNodes = SFFM.CellNodes(mesh)[:,convert(Int,ceil(x₀/Δ))]
basisValues = zeros(length(theNodes))
for n in 1:length(theNodes)
    basisValues[n] = prod(x₀.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
end
initpm = [
    zeros(sum(model.C.<=0)) # LHS point mass
    zeros(sum(model.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,SFFM.NBases(mesh),SFFM.NIntervals(mesh),SFFM.NPhases(model))
initprobs[:,convert(Int,ceil(x₀/Δ)),1] = basisValues'*All.Matrices.Local.V.V*All.Matrices.Local.V.V'.*2/Δ
initdist = SFFM.SFFMDensity(
    initpm,
    initprobs,
    SFFM.CellNodes(mesh),
)
x0 = SFFM.Dist2Coeffs( model, mesh, initdist)
h = 0.0001
@time projyvals = SFFM.EulerDG( proj.D["++"](s = 0), y, x0, h = h)
@time interpyvals = SFFM.EulerDG( interp.D["++"](s = 0), y, x0, h = h)

## plotting
simdensity = SFFM.Sims2Dist(model,mesh,sims,SFFM.SFFMProbability)
projdensity = SFFM.Coeffs2Dist(model,mesh,projyvals,SFFM.SFFMProbability)
interpdensity = SFFM.Coeffs2Dist(model,mesh,interpyvals,SFFM.SFFMProbability)
println("proj error: ",SFFM.starSeminorm( projdensity, simdensity))
println("interp error: ",SFFM.starSeminorm( interpdensity, simdensity))

## plots
# plot solutions
# densities
p = SFFM.plot(model; mesh = mesh,dist = projdensity,color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p,model,mesh,interpdensity,color=:red,label="interp")#,marker=:x)
# plot sims
p = SFFM.plot!(p,model,mesh,simdensity,color=:black,label="sim")#,marker=:ltriangle)

display(p)
## probabilities
simdensity = SFFM.Sims2Dist(model,mesh,sims,SFFM.SFFMDensity)
projdensity = SFFM.Coeffs2Dist(model,mesh,projyvals,SFFM.SFFMDensity)
interpdensity = SFFM.Coeffs2Dist(model,mesh,interpyvals,SFFM.SFFMDensity)
# vikramdensity = SFFM.Coeffs2Dist(model=model,mesh=mesh,Coeffs=vikramyvals,SFFM.SFFMDensity)

## plots
# plot solutions
# densities
p = SFFM.plot(model,mesh=mesh,dist=projdensity,color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p,model,mesh,interpdensity,color=:red,label="interp")#,marker=:x)
# plot sims
p = SFFM.plot!(p,model,mesh,simdensity,color=:black,label="sim")#,marker=:ltriangle)
# savefig(p,"/Users/a1627293/Dropbox/PhD/NotesMaster/whichR/m1.png")
display(p)
## end 1st model

########################
########################
########################

## Define 2nd model
T = [-1.0 1.0; 1.0 -1.0]
C = [1.0; -1.0]
r = (
    r = function (x)
        [0.1.*x.+0.05 (x.>-1)]
    end,
    R = function (x)
        [0.1.*x.^2/2.0.+0.05.*x (x.>-1).*x]
    end,
)

Bounds = [0 30; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = Bounds)

## Define mesh
# Δ = 3
Nodes = collect(Bounds[1,1]:Δ:Bounds[1,2])
# nBases = 3
Basis = "lagrange"
mesh = SFFM.DGMesh(model, Nodes, nBases, Basis=Basis)

## Make matrices
All = SFFM.MakeAll( model, mesh, approxType = "projection")
proj = All
interp = SFFM.MakeAll( model, mesh, approxType="interpolation")

vikramD = copy(interp.D["++"]())
vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:] =
    repeat(All.Matrices.Local.V.w, SFFM.NIntervals(mesh)*SFFM.NPhases(model), 1) .* vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:]

## sims fro gound truth
x₀ = 15
NSim = 50000
IC = (φ = ones(Int, NSim), X = x₀ .* ones(NSim), Y = zeros(NSim))
y = 10
@time sims =
    SFFM.SimSFFM( model, SFFM.FirstExitY(-Inf, y), IC)

## DG stationary dist
# construct initial condition
theNodes = SFFM.CellNodes(mesh)[:,convert(Int,ceil(x₀/Δ))]
basisValues = zeros(length(theNodes))
for n in 1:length(theNodes)
    basisValues[n] = prod(x₀.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
end
initpm = [
    zeros(sum(model.C.<=0)) # LHS point mass
    zeros(sum(model.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,SFFM.NBases(mesh),SFFM.NIntervals(mesh),SFFM.NPhases(model))
initprobs[:,convert(Int,ceil(x₀/Δ)),1] = basisValues'*All.Matrices.Local.V.V*All.Matrices.Local.V.V'.*2/Δ
initdist = SFFM.SFFMDensity(
    initpm,
    initprobs,
    SFFM.CellNodes(mesh),
)
x0 = SFFM.Dist2Coeffs( model, mesh, initdist)
h = 0.0001
@time projyvals = SFFM.EulerDG(proj.D["++"](s = 0), y, x0, h = h)
@time interpyvals = SFFM.EulerDG(interp.D["++"](s = 0), y, x0, h = h)

## plotting
simdensity = SFFM.Sims2Dist( model, mesh, sims,SFFM.SFFMProbability)
projdensity = SFFM.Coeffs2Dist( model, mesh, projyvals,SFFM.SFFMProbability)
interpdensity = SFFM.Coeffs2Dist( model, mesh, interpyvals,SFFM.SFFMProbability)

println("proj error: ",SFFM.starSeminorm( projdensity, simdensity))
println("interp error: ",SFFM.starSeminorm(interpdensity, simdensity))
## plots
# plot solutions
# densities
# plot sims
p = SFFM.plot( model, mesh = mesh, dist = simdensity, color=:black,label="sim")#,marker=:ltriangle)
p = SFFM.plot!(p, model, mesh, projdensity, color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p, model, mesh, interpdensity, color=:red,label="interp")#,marker=:x)

display(p)
## probabilities
simdensity = SFFM.Sims2Dist( model, mesh, sims,SFFM.SFFMDensity)
projdensity = SFFM.Coeffs2Dist( model, mesh, projyvals,SFFM.SFFMDensity)
interpdensity = SFFM.Coeffs2Dist( model, mesh, interpyvals,SFFM.SFFMDensity)

## plots
# plot solutions
# densities
p = SFFM.plot( model, mesh = mesh, dist = projdensity, color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p, model, mesh, interpdensity, color=:red,label="interp")#,marker=:x)

# plot sims
p = SFFM.plot!(p, model, mesh, simdensity, color=:black,label="sim")#,marker=:ltriangle)
# savefig(p,"/Users/a1627293/Dropbox/PhD/NotesMaster/whichR/m2.png")
display(p)
## end 2nd model

########################
########################
########################

## Define 3rd model
T = [-1.0 1.0; 1.0 -1.0]
C = [1.0; -1.0]
r = (
    r = function (x)
        [cos.(x).+1.05 x.>-1]
    end,
    R = function (x)
        [sin.(x).+1.05.*x (x.>-1).*x]
    end,
)

Bounds = [0 10*π; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = Bounds)

## Define mesh
# Δ = 3
Nodes = collect(Bounds[1,1]:Δ:Bounds[1,2])
# nBases = 3
Basis = "lagrange"
mesh = SFFM.DGMesh(model, Nodes, nBases, Basis=Basis)

## Make matrices
All = SFFM.MakeAll( model, mesh, approxType = "projection")
proj = All
interp = SFFM.MakeAll( model, mesh, approxType="interpolation")

vikramD = copy(interp.D["++"]())
vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:] =
    repeat(All.Matrices.Local.V.w, SFFM.NIntervals(mesh)*SFFM.NPhases(model), 1) .* vikramD[sum(model.C.<=0)+1:end-sum(model.C.>=0),:]

## sims fro gound truth
x₀ = 15
NSim = 50000
IC = (φ = ones(Int, NSim), X = x₀ .* ones(NSim), Y = zeros(NSim))
y = 10
@time sims =
    SFFM.SimSFFM( model, SFFM.FirstExitY( -Inf, y), IC)

## DG stationary dist
# construct initial condition
theNodes = SFFM.CellNodes(mesh)[:,convert(Int,ceil(x₀/Δ))]
basisValues = zeros(length(theNodes))
for n in 1:length(theNodes)
    basisValues[n] = prod(x₀.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
end
initpm = [
    zeros(sum(model.C.<=0)) # LHS point mass
    zeros(sum(model.C.>=0)) # RHS point mass
]
initprobs = zeros(Float64,SFFM.NBases(mesh),SFFM.NIntervals(mesh),SFFM.NPhases(model))
initprobs[:,convert(Int,ceil(x₀/Δ)),1] = basisValues'*All.Matrices.Local.V.V*All.Matrices.Local.V.V'.*2/Δ
initdist = SFFM.SFFMDensity(
    initpm,
    initprobs,
    SFFM.CellNodes(mesh),
)
x0 = SFFM.Dist2Coeffs( model, mesh, initdist)
h = 0.0001
@time projyvals = SFFM.EulerDG( proj.D["++"](s = 0), y, x0, h = h)
@time interpyvals = SFFM.EulerDG( interp.D["++"](s = 0), y, x0, h = h)

## plotting
simdensity = SFFM.Sims2Dist( model, mesh, sims, SFFM.SFFMProbability)
projdensity = SFFM.Coeffs2Dist( model, mesh, projyvals, SFFM.SFFMProbability)
interpdensity = SFFM.Coeffs2Dist( model, mesh, interpyvals, SFFM.SFFMProbability)
println("proj error: ",SFFM.starSeminorm( projdensity, simdensity))
println("interp error: ",SFFM.starSeminorm( interpdensity, simdensity))

## plots
# plot solutions
# densities
p = SFFM.plot( model, mesh = mesh, dist = projdensity, color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p, model, mesh, interpdensity, color=:red,label="interp")#,marker=:x)
# plot sims
p = SFFM.plot!(p, model, mesh, simdensity, color=:black,label="sim")#,marker=:ltriangle)

display(p)
## probabilities
simdensity = SFFM.Sims2Dist( model, mesh, sims, SFFM.SFFMDensity)
projdensity = SFFM.Coeffs2Dist( model, mesh, projyvals, SFFM.SFFMDensity)
interpdensity = SFFM.Coeffs2Dist( model, mesh, interpyvals, SFFM.SFFMDensity)

## plots
# plot solutions
# densities
p = SFFM.plot( model, mesh = mesh, dist = projdensity, color=:blue,label="proj")#,marker=:rtriangle)
p = SFFM.plot!(p, model, mesh, interpdensity, color=:red,label="interp")#,marker=:x)
# plot sims
p = SFFM.plot!(p, model, mesh, simdensity, color=:black,label="sim")#,marker=:ltriangle)
# savefig(p,"/Users/a1627293/Dropbox/PhD/NotesMaster/whichR/m3.png")
display(p)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
function *(A::Const,B::PotLogArray)
    L=PotLogArray(B.variables,log(A.content)+B.content)
    return L
end

# Not sure I understand this but Julia automatically exports *,+,-,/
#function *(A::Const,B::Const)
#    C=Const(A.content*B.content)
#end
# must be a smarter way to do the following using macros:!
function *(A::Const,B::Const)
    return Const(A.content*B.content)
end
function /(A::Const,B::Const)
    return Const(A.content/B.content)
end
function +(A::Const,B::Const)
    return Const(A.content+B.content)
end
function -(A::Const,B::Const)
    return Const(A.content-B.content)
end
function *(P::PotArray,C::Const)
    p = deepcopy(P); p.content=C.content*p.content
    return p
end

function *(C::Const,P::PotArray)
    p = deepcopy(P); p.content=C.content*p.content
    return p
end
function /(P::PotArray,C::Const)
    p = deepcopy(P); p.content=p.content/C.content
    return p
end
function +(P::PotArray,C::Const)
    p = deepcopy(P); p.content=C.content+p.content
    return p
end
function +(C::Const,P::PotArray)
    p = deepcopy(P); p.content=C.content+p.content
    return p
end
function -(P::PotArray,C::Const)
    p = deepcopy(P); p.content=C.content-p.content
    return p
end
function -(C::Const,P::PotArray)
    p = deepcopy(P); p.content=C.content-p.content
    return p
end

function convert(::Type{PotArray},c::Const)
    return PotArray([],c.content)
end


text:
@testset "Ethernet" begin

@testset "IO" begin
    @test ethertype_string(ETHERTYPE_IPV4) == "IPv4(0800)"
    @test ethertype_string(ETHERTYPE_ARP) == "ARP(0806)"
    @test ethertype_string(0x1234) == "Unknown(1234)"

    x = EthernetHeader(
        MACAddress((0x1, 0x2, 0x3, 0x4, 0x5, 0x6)),
        MACAddress((0x6, 0x5, 0x4, 0x3, 0x2, 0x1)),
        ETHERTYPE_IPV4)
    @test string(x) == "(06:05:04:03:02:01 -> 01:02:03:04:05:06 IPv4(0800))"
end

@testset "decode" begin
    data = [
        0x1, 0x2, 0x3, 0x4, 0x5, 0x6,
        0x6, 0x5, 0x4, 0x3, 0x2, 0x1,
        0x08, 0x00,
        0xaa, 0xbb, 0xbb, 0xaa
    ]
    x = decode_ethernet(data)
    @test x.header.dst_mac == MACAddress((0x1, 0x2, 0x3, 0x4, 0x5, 0x6))
    @test x.header.src_mac == MACAddress((0x6, 0x5, 0x4, 0x3, 0x2, 0x1))
    @test x.header.ethertype == ETHERTYPE_IPV4
    @test length(x.payload) == 4
    GC.@preserve data begin
        @test x.payload[1] == 0xaa
        @test x.payload[2] == 0xbb
        @test x.payload[3] == 0xbb
        @test x.payload[4] == 0xaa
    end
end

end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols


text:
##### Generic methods #####

"""
    size(d::MatrixDistribution)

Return the size of each sample from distribution `d`.
"""
Base.size(d::MatrixDistribution)

size(d::MatrixDistribution, i) = size(d)[i]

"""
    length(d::MatrixDistribution)

The length (*i.e* number of elements) of each sample from the distribution `d`.
"""
Base.length(d::MatrixDistribution)

"""
    rank(d::MatrixDistribution)

The rank of each sample from the distribution `d`.
"""
LinearAlgebra.rank(d::MatrixDistribution)

"""
    inv(d::MatrixDistribution)

If known, returns a `MatrixDistribution` instance representing the
distribution of inv(X), where X is a random matrix with distribution `d`.
"""
Base.inv(d::MatrixDistribution)

"""
    mean(d::MatrixDistribution)

Return the mean matrix of `d`.
"""
mean(d::MatrixDistribution)

"""
    var(d::MatrixDistribution)

Compute the matrix of element-wise variances for distribution `d`.
"""
var(d::MatrixDistribution) = ((n, p) = size(d); [var(d, i, j) for i in 1:n, j in 1:p])

"""
    cov(d::MatrixDistribution)

Compute the covariance matrix for `vec(X)`, where `X` is a random matrix with distribution `d`.
"""
function cov(d::MatrixDistribution)
    M = length(d)
    V = zeros(partype(d), M, M)
    iter = CartesianIndices(size(d))
    for el1 = 1:M
        for el2 = 1:el1
            i, j = Tuple(iter[el1])
            k, l = Tuple(iter[el2])
            V[el1, el2] = cov(d, i, j, k, l)
        end
    end
    return V + tril(V, -1)'
end
cov(d::MatrixDistribution, ::Val{true}) = cov(d)

"""
    cov(d::MatrixDistribution, flattened = Val(false))

Compute the 4-dimensional array whose `(i, j, k, l)` element is `cov(X[i,j], X[k, l])`.
"""
function cov(d::MatrixDistribution, ::Val{false})
    n, p = size(d)
    [cov(d, i, j, k, l) for i in 1:n, j in 1:p, k in 1:n, l in 1:p]
end

# pdf & logpdf

# TODO: Remove or restrict - this causes many ambiguity errors...
_logpdf(d::MatrixDistribution, X::AbstractMatrix{<:Real}) = logkernel(d, X) + d.logc0

#  for testing
is_univariate(d::MatrixDistribution) = size(d) == (1, 1)
check_univariate(d::MatrixDistribution) = is_univariate(d) || throw(ArgumentError("not 1 x 1"))

##### Specific distributions #####

for fname in ["wishart.jl", "inversewishart.jl", "matrixnormal.jl",
              "matrixtdist.jl", "matrixbeta.jl", "matrixfdist.jl", "lkj.jl"]
    include(joinpath("matrix", fname))
end


text:
"""
    InterpolateAndBin{L}(f::Function, left_bin_edges, intp::L, intp_grid,
        extrapolation_bc::Union{<:Real, Interpolations.BoundaryCondition})

Indicates that a dataset consisting of both indices and values should first
be interpolated to the `intp_grid` grid using the provided `intp` scheme
(e.g. `Linear()`). After interpolating, assign the interpolated values
to the bins defined by `left_bin_edges` and summarise the values falling
in each bin using the summary function `f` (e.g. `mean`).

## Example

```julia
using UncertainData, Interpolations, StatsBase

# Assume we have the following unevenly spaced data with some `NaN` values
T = 100
time = sample(1.0:T*5, T, ordered = true, replace = false)
y1 = rand(T)
time[rand(1:T, 10)] .= NaN
y1[rand(1:T, 10)] .= NaN

# We want to first intepolate the dataset linearly to a regular time grid
# with steps of `0.1` time units.
intp = Linear()
intp_grid = 0:0.1:1000
extrapolation_bc = Flat(OnGrid())

# Then, bin the dataset in time bins `50` time units wide, collect all
# values in each bin and summarise them using `f`.
f = mean
left_bin_edges = 0:50:1000

r = InterpolateAndBin(f, left_bin_edges, intp, intp_grid, extrapolation_bc)
```
"""
struct InterpolateAndBin{L}
    f::Function
    left_bin_edges
    intp::L # Linear
    intp_grid
    extrapolation_bc::Union{<:Real, Interpolations.BoundaryCondition}
end

export InterpolateAndBin


text:
function test_thermalLoad_1()

    T = 2

    # With these values, the heater uses 0.2 power to keep a ΔT of 1
    l = DR.DER.ThermalLoad(
        index=1, T=2,
        temp_min=[20.0, 20.0],
        temp_max=[25.0, 25.0],
        temp_ext=[0.0, 0.0],
        temp_init=20.0,
        pwr_min=[0.0, 0.0],
        pwr_max=[10.0, 10.0],
        C=1.0, η=1.0, μ=0.2,
        binflag=false
    )

    # GLPK does not support -in-Interval constraints, so we bridge
    model = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

    lmin = [0.0, 0.0]
    lmax = [10.0, 10.0]
    price = [2.0, 1.0]

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [l]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 12.0

    # Temperature should be kept at minimum value, i.e., 20
    T1 = var2idx[(:thermal, l.index, :temp, 1)]
    T2 = var2idx[(:thermal, l.index, :temp, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), T1) ≈ 20.0
    @test MOI.get(model, MOI.VariablePrimal(), T2) ≈ 20.0

    # Check net load
    pnet1 = var2idx[(:thermal, l.index, :pnet, 1)]
    pnet2 = var2idx[(:thermal, l.index, :pnet, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), pnet1) ≈ 4.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet2) ≈ 4.0

    return nothing

end

function test_thermalLoad_2()

    T = 2

    # With these values, the heater uses 0.2 power to keep a ΔT of 1
    l = DR.DER.ThermalLoad(
        index=1, T=2,
        temp_min=[20.0, 20.0],
        temp_max=[25.0, 25.0],
        temp_ext=[0.0, 0.0],
        temp_init=20.0,
        pwr_min=[0.0, 0.0],
        pwr_max=[10.0, 10.0],
        C=1.0, η=1.0, μ=0.2,
        binflag=false
    )

    # GLPK does not support -in-Interval constraints, so we bridge
    model = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

    lmin = [0.0, 0.0]
    lmax = [10.0, 10.0]
    price = [1.0, 2.0]  # Now it's worth heating more in the first period

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [l]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    # T1 should be higher
    # T2 should be at minimum, i.e., 20
    T1 = var2idx[(:thermal, l.index, :temp, 1)]
    T2 = var2idx[(:thermal, l.index, :temp, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), T1) ≈ 25.0
    @test MOI.get(model, MOI.VariablePrimal(), T2) ≈ 20.0

    # Check net load
    pnet1 = var2idx[(:thermal, l.index, :pnet, 1)]
    pnet2 = var2idx[(:thermal, l.index, :pnet, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), pnet1) ≈ 9.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet2) ≈ 0.0

    # Check objective value
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 9.0
    return nothing

end

@testset "Thermal" begin
    test_thermalLoad_1()
    test_thermalLoad_2()
end


text:
using Manopt, Manifolds, Test
@testset "Stochastic Gradient Descent" begin
    M = Sphere(2)
    # 5 point mean
    p = [0.0, 0.0, 1.0]
    s = 1.0
    pts = [
        exp(M, p, X) for
        X in [zeros(3), [s, 0.0, 0.0], [-s, 0.0, 0.0], [0.0, s, 0.0], [0.0, -s, 0.0]]
    ]
    F(y) = 1 / 2 * sum([distance(M, y, x)^2 for x in pts])
    ∇F(y) = sum([-log(M, y, x) for x in pts])
    s∇F1(y) = [-log(M, y, x) for x in pts]
    s∇F2 = [(y -> -log(M, y, x)) for x in pts]
    @testset "Constructors" begin
        p1 = StochasticGradientProblem(M, s∇F1)
        p2 = StochasticGradientProblem(M, s∇F2)
        @test get_gradient(p1, 1, p) == zeros(3)
        @test get_gradient(p2, 1, p) == zeros(3)
        @test get_gradients(p1, p) == get_gradients(p2, p)
        o = StochasticGradientDescentOptions(deepcopy(p))
        o.order = collect(1:5)
        o.order_type = :Linear
        initialize_solver!(p1, o)
        o.order_type = :Linear
        step_solver!(p1, o, 1)
        @test o.x == exp(M, p, get_gradient(p1, 1, p))
    end
    @testset "Momentum and Average Processor Constructors" begin
        p1 = StochasticGradientProblem(M, s∇F1)
        p2 = GradientProblem(M, F, ∇F)
        m1 = MomentumGradient(M, p, StochasticGradient())
        m2 = MomentumGradient(p1, p)
        @test m1.direction == m2.direction #both use StochasticGradient
        m3 = MomentumGradient(M, p)
        m4 = MomentumGradient(p2, p)
        @test m3.direction == m4.direction #both use Gradient
        a1 = AverageGradient(M, p, 10, StochasticGradient())
        a2 = AverageGradient(p1, p, 10)
        @test a1.direction == a2.direction #both use StochasticGradient
        a3 = AverageGradient(M, p, 10)
        a4 = AverageGradient(p2, p, 10)
        @test a3.direction == a4.direction #both use Gradient
    end
    @testset "Comparing Stochastic Methods" begin
        x1 = stochastic_gradient_descent(M, s∇F1, p; order_type=:Linear)
        @test norm(x1) ≈ 1
        o1 = stochastic_gradient_descent(
            M, s∇F1, p; order_type=:Linear, return_options=true
        )
        x1a = get_solver_result(o1)
        @test x1 == x1a
        x2 = stochastic_gradient_descent(M, s∇F1, p; order_type=:FixedRandom)
        @test norm(x2) ≈ 1
        x3 = stochastic_gradient_descent(M, s∇F1, p; order_type=:Random)
        @test norm(x3) ≈ 1
        x4 = stochastic_gradient_descent(
            M,
            s∇F1,
            p;
            order_type=:Random,
            direction=AverageGradient(M, p, 10, StochasticGradient()),
        )
        @test norm(x4) ≈ 1
        x5 = stochastic_gradient_descent(
            M,
            s∇F1,
            p;
            order_type=:Random,
            direction=MomentumGradient(M, p, StochasticGradient()),
        )
        @test norm(x5) ≈ 1
    end
end


text:
import Example

println("Hello world!")

max_stars_repo_name: julia-actions/manifest-updater-examples

max_issues_repo_name: julia-actions/manifest-updater-examples

max_forks_repo_name: julia-actions/manifest-updater-examples


text:
using Test, Suppressor, DataFrames
using HMCResearchRandomWalks.ViewOut

const test_df = DataFrame(test_title1 = [1,2,3,4],
                        test_title2 = ["ab", "cd", "ef", "gh"])

@testset "savedata" begin
    current_df(test_df)
    current_seed(8765)
    @suppress_out begin
        global path1 = savedata()
    end
    _, ext = splitext(path1)
    @test isfile(path1)
    @test occursin("8765", path1)
    @test ext === ".csv"

    @suppress_out begin
        global path2 = savedata("test_input_name.csv")
    end
    @test isfile(path2)
    @test occursin("test_input_name.csv", path2)

    output = @capture_out begin
        global path3 = savedata("test_safe_name.csv", test_df)
    end
    correct_output = string("CSV File saved at the given path:\n", path3,"\n")
    @test isfile(path3)
    @test occursin("test_safe_name.csv", path3)
    @test output === correct_output

    correct_warning = "File name must end with the .csv extension"
    error = @capture_err begin
        savedata("test_bad_name", test_df)
    end
    @test occursin(correct_warning, error)
end

@testset "savetofolder" begin
    current_df(test_df)
    current_seed(4321)
    current_path("out/noflowdata/")
    output = @capture_out begin
        global path4 = savetofolder()
    end
    correct_output = string("CSV File saved at the given path:\n", path4,"\n")
    @test output === correct_output
    @test isfile(path4)
    @test occursin("seed_4321", path4)
end

@testset "loaddata" begin
    for p in [path1, path2, path3, path4]
        loaded_test_df = loaddata(p)
        @test loaded_test_df == test_df
    end

    easy_path = "/out/test_input_name.csv"
    easy_df = loaddata(easy_path, true)
    @test easy_df == test_df
end

@testset "clean up saved test files" begin
    for p in [path1, path2, path3, path4]
        rm(p)
        @test !isfile(p)
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
abstract type AbstractSChain{TT} end
struct SChain{TT} <: AbstractSChain{TT}
   F::TT
end

struct TypedChain{TT, IN, OUT} <: AbstractSChain{TT}
   F::TT
end

# construct a chain recursively
chain(F1, F2, args...) = chain( chain(F1, F2), args... )
# for most arguments, just form a tuple
chain(F1, F2) = SChain( (F1, F2) )
# if one of them is a chain already, then combine into a single long chain
chain(F1::SChain, F2) = SChain( tuple(F1.F..., F2) )
chain(F1, F2::SChain) = SChain( tuple(F1, F2.F...) )
chain(F1::SChain, F2::SChain) = chain( tuple(F1.F..., F2.F...) )

Base.length(c::SChain) = length(c.F)

@generated function evaluate(chain::AbstractSChain{TT}, X) where {TT}
   LEN = length(TT.types)
   code = Expr[]
   push!(code, :(X_0 = X))
   for l = 1:LEN
      push!(code, Meta.parse("F_$l = chain.F[$l]"))
      push!(code, Meta.parse("X_$l = evaluate(F_$l, X_$(l-1))"))
      push!(code, Meta.parse("release!(X_$(l-1))"))
   end
   push!(code, Meta.parse("return X_$LEN"))
   return Expr(:block, code...)
end

# TODO:
# - replace with an frule and then wrap that into an evaluate_ed
# - generated function to make this fast
# - implement the rrule

@generated function evaluate_ed(chain::AbstractSChain{TT}, X) where {TT}
   LEN = length(TT.types)
   code = Expr[]
   push!(code, :(X_0 = X))
   push!(code, Meta.parse("X_1, dF_1 = evaluate_ed(chain.F[1], X_0)"))
   for l = 2:LEN
      push!(code, Meta.parse("F_$l = chain.F[$l]"))
      push!(code, Meta.parse("X_$l, dF_$l = frule_evaluate(F_$l, X_$(l-1), dF_$(l-1))"))
      push!(code, Meta.parse("release!(X_$(l-1))"))
   end
   push!(code, Meta.parse("return X_$LEN, dF_$LEN"))
   return Expr(:block, code...)
end

# function evaluate_ed(chain::AbstractSChain{TT}, X) where {TT}
#    LEN = length(chain.F)
#    Xi = evaluate(chain.F[1], X)
#    dFi = evaluate_d(chain.F[1], X)
#    for i = 2:LEN
#       Xi, dFi = frule_evaluate(chain.F[i], Xi, dFi)
#    end
#    return Xi, dFi
# end

evaluate_d(chain::AbstractSChain, X) = evaluate_ed(chain, X)[2]

# TODO: This still needs sorting out ...
#       maybe we can no kill this?

valtype(chain::SChain) = valtype(chain.F[end])

valtype(chain::SChain, x) = valtype(chain)

gradtype(chain::SChain) = gradtype(chain.F[end])

# function gradtype(chain::SChain, x)
#    LEN = length(chain.F)
#    Ti = gradtype(chain.F[1], x)

# end

valtype(chain::TypedChain{TT, IN, OUT}, args...) where {TT, IN, OUT} = OUT

gradtype(chain::TypedChain{TT, IN, OUT}, args...) where {TT, IN, OUT} = _gradtype(OUT, IN)

_gradtype(::Type{Vector{T1}}, ::Type{T2}) where {T1 <: Number, T2 <: Number} =
      Vector{promote_type(T1, T2)}

_gradtype(::Type{Vector{T1}}, ::Type{SVector{N, T2}}) where {N, T1 <: Number, T2 <: Number} =
      Vector{SVector{N, promote_type(T1, T2)}}

##

import Base: ==

==(ch1::SChain, ch2::SChain) =
      all( F1==F2 for (F1, F2) in zip(ch1.F, ch2.F) )

write_dict(chain::SChain) = Dict(
            "__id__" => "ACE_SChain",
            "F" => write_dict.(chain.F)
         )

read_dict(::Val{:ACE_SChain}, D::Dict) =
         SChain(tuple( read_dict.(D["F"])... ))

##

#=

# Experimental AD codes

import ChainRulesCore, ChainRules
import ChainRulesCore: rrule, NoTangent

struct EVAL{TM, TX}
   m::TM
   X::TX
end

struct EVAL_D{TM, TX}
   m::TM
   X::TX
end

function (y::EVAL)(params)
   set_params!(y.m, params)
   return evaluate(y.m, y.X).val
end

function rrule(y::EVAL, params)
   set_params!(y.m, params)
   val = evaluate(y.m, y.X).val
   adj = dp -> ( NoTangent(), dp * getproperty.(grad_params(y.m, y.X), :val))
   return val, adj
end

function (y::EVAL_D)(params)
   set_params!(y.m, params)
   return grad_config(y.m, y.X)
end

function rrule(y::EVAL_D, params)
   set_params!(y.m, params)
   val = grad_config(y.m, y.X)
   adj = dp -> ( NoTangent(), adjoint_EVAL_D(y.m, y.X, dp) )
   return val, adj
end

=#


text:
# Autogenerated wrapper script for minimap2_jll for x86_64-apple-darwin
export minimap2

using Zlib_jll
JLLWrappers.@generate_wrapper_header("minimap2")
JLLWrappers.@declare_library_product(minimap2, "@rpath/libminimap2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        minimap2,
        "lib/libminimap2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()