text:
import Tracker: TrackedArray, TrackedReal

# const TrArray = Union{Array{<:Real}, TrackedArray}
# const TrArray1 = Union{Array{<:Real, 1}, TrackedArray{<:Real, 1, <:AbstractArray{<:Real, 1}}}
# const TrArray2 = Union{Array{<:Real, 2}, TrackedArray{<:Real, 1, <:AbstractArray{<:Real, 2}}}

#const TrArray{T, N} = Union{Array{T, N} where {T<:Real} where {N<:Int}, TrackedArray{T, N, AbstractArray{T, N}} where {T<:Real} where {N<:Int}}
#TrArray = TrArray{T, N} where {T, N}
const TrReal = Union{<:Real, TrackedReal}

abstract type AbstractTDiagonalNormal <: ContinuousUnivariateDistribution end
abstract type AbstractTMVDiagonalNormal <: ContinuousMultivariateDistribution end

# Distributions.rand(rng::AbstractRNG, d::AbstractTMVDiagonalNormal) = Distributions._rand!(rng, d, Vector{eltype(d)})(length(d))
# Distributions.rand(rng::AbstractRNG, d::AbstractTMVDiagonalNormal, n::Int) = Distributions._rand!(rng, d, Matrix{eltype(d)}(length(d), n))
# Distributions.rand!(rng::AbstractRNG, d::AbstractTMVDiagonalNormal, x::VecOrMat) = Distributions._rand!(rng, d, x)

struct TMVDiagonalNormal{T<:AbstractArray} <: AbstractTMVDiagonalNormal
    μ::T
    logσ::T
    function TMVDiagonalNormal{T}(μ::T, logσ::T) where T<:AbstractArray{S, 1} where {S}
        @assert size(μ) == size(logσ)
        new(μ, logσ)
    end
    function TMVDiagonalNormal{T}(μ::T, logσ::T) where T<:AbstractArray{S, 2} where {S}
        @assert size(μ,2) == 1
        @assert size(logσ,2) == 1
        TMVDiagonalNormal(dropdims(μ, dims=2), dropdims(logσ, dims=2))
    end
end

TMVDiagonalNormal(μ::T, logσ::T) where {T<:Real} = TMVDiagonalNormal{Array{T, 1}}([μ], [logσ]) # Convert to multivariate case
TMVDiagonalNormal(μ::T, logσ::T) where {T<:AbstractArray} = TMVDiagonalNormal{T}(μ, logσ)

# https://github.com/JuliaLang/julia/pull/26601
Base.broadcastable(t::AbstractTMVDiagonalNormal) = Ref(t)

Base.convert(::Type{<:MvNormal}, d::TMVDiagonalNormal{<:Array{<:Real}}) = MvNormal(d.μ, exp.(d.logσ))
#Base.convert(::Type{<:MvNormal}, d::TMVDiagonalNormal{<:TrackedArray}) = MvNormal(d.μ.data, exp.(d.logσ.data))

# TTMVDiagonalNormal: https://juliastats.github.io/Distributions.jl/latest/extends.html
Base.length(d::TMVDiagonalNormal) = length(d.μ)
sampler(d::TMVDiagonalNormal) = d
#Distributions._rand!(d::TMVDiagonalNormal, x::AbstractArray) = Distributions._rand!(convert(MvNormal, d), x)
Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{<:Array{<:Real}}) = Distributions.rand(rng, convert(MvNormal, d))
Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{<:Array{<:Real}}, n::Int) = Distributions.rand(rng, convert(MvNormal, d), n)
Distributions._logpdf(d::TMVDiagonalNormal, x::AbstractArray) = Distributions._logpdf(convert(MvNormal, d), x)
Distributions.mean(d::TMVDiagonalNormal) = d.μ
Distributions.var(d::TMVDiagonalNormal) = exp.(2 * d.logσ)
Distributions.entropy(d::TMVDiagonalNormal, args...) = Distributions.entropy(convert(MvNormal, d), args...)
Distributions.cov(d::TMVDiagonalNormal) = diagm(0 => Distributions.var(d))

Distributions._rand!(d::TMVDiagonalNormal, x::VecOrMat) = Distributions._rand!(Random.GLOBAL_RNG, convert(MvNormal, d), x)
Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal, x::VecOrMat) = Distributions._rand!(rng, convert(MvNormal, d), x)

Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal, x::AbstractMatrix) = Distributions._rand!(rng, convert(MvNormal, d), x)
Distributions._rand!(d::TMVDiagonalNormal, x::AbstractMatrix) = Distributions._rand!(Random.GLOBAL_RNG, convert(MvNormal, d), x)D
Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal, x::AbstractVector) = Distributions._rand!(rng, convert(MvNormal, d), x)
Distributions._rand!(d::TMVDiagonalNormal, x::AbstractVector) = Distributions._rand!(Random.GLOBAL_RNG, convert(MvNormal, d), x)

# For sampling with TrackedArrays, we want to use the reparmaterisation trick:
Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{<:TrackedArray}) = sample(rng, d)
Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{<:TrackedArray}, n::Int) = sample(rng, d, n)

Distributions._rand!(d::TMVDiagonalNormal{<:TrackedArray}, x::VecOrMat) = error("Not Implemented Error")
Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal{<:TrackedArray}, x::VecOrMat) = error("Not Implemented Error")

Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal{<:TrackedArray}, x::AbstractMatrix) = error("Not Implemented Error")
Distributions._rand!(d::TMVDiagonalNormal{<:TrackedArray}, x::AbstractMatrix) = error("Not Implemented Error")
Distributions._rand!(rng::AbstractRNG, d::TMVDiagonalNormal{<:TrackedArray}, x::AbstractVector) = error("Not Implemented Error")
Distributions._rand!(d::TMVDiagonalNormal{<:TrackedArray}, x::AbstractVector) = error("Not Implemented Error")

function data(
    d::TrackedDistributions.TMVDiagonalNormal{T}
) where {N, S<:Real, T<:TrackedArray{S, N, Array{S,N}}}
    TrackedDistributions.TMVDiagonalNormal(
        Array{Float64}(mean(d).data),
        Array{Float64}(log.(sqrt.(var(d)).data)),
    )
end

# Getters
logσ(d::TMVDiagonalNormal) = d.logσ

"""
    kl_q_p(q::DiagonalNormal, p::DiagonalNormal)

    Compute KL divergence between multidimenional diagonal Gaussian distributions
    with the given means and log-sigmas: KL(N(μ_q, logσ_q) || N(μ_p, logσ_p)).
    Reference: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Likelihood_function
"""
function kl_q_p(q::TMVDiagonalNormal, p::TMVDiagonalNormal)
    μ_q, logσ_q = q.μ, q.logσ
    μ_p, logσ_p = p.μ, p.logσ
    @assert(size(μ_q) == size(μ_p))
    @assert(ndims(μ_q) == 1)
    D = length(μ_q)
    return - 0.5 * D + 0.5 .* sum( 2 .* (logσ_p .- logσ_q) .+ exp.(2 .* (logσ_q .- logσ_p)) .+ ( (μ_q .- μ_p).^2 ./ exp.(2 .* logσ_p) ) )
end
#
#
"""
    sample(rng, d::DiagonalNormal)

    Used for the Reparametrization trick.  Note we pass in rng.
    TODO: Convert to sampler and rand
"""
function sample(rng::AbstractRNG, d::TMVDiagonalNormal)
    μ = d.μ
    logσ = d.logσ
    return μ .+ exp.(logσ) .* randn(rng, Float64, size(μ))
end

Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{TrackedArray{T, N, Array{T, N}}}) where {T<:Real, N} = sample(rng, d)
Distributions.rand(d::TMVDiagonalNormal{TrackedArray{T, N, Array{T, N}}}) where {T<:Real, N} = error("Seeding Enforced. use rand(rng, d)")
Distributions.rand(rng::AbstractRNG, d::TMVDiagonalNormal{TrackedArray{T, N, Array{T, N}}}, n::Int) where {T<:Real, N} = error("Not Implemented Error")

Distributions._logpdf(d::TMVDiagonalNormal{TrackedArray{T, N, Array{T, N}}}, x::AbstractArray) where {T<:Real, N} = log_pdf(d, x)

function log_pdf(d::TMVDiagonalNormal{TrackedArray{T, 1, Array{T, 1}}}, x::AbstractArray) where {T<:Real}
    μ = d.μ
    logσ = d.logσ
    n = 1
    D = length(μ)
    return -.5 * ( n * (D * log(2π) + 2 * sum(logσ)) + sum((x .- μ).^2 .* exp.(-2 .* logσ)) )
end


text:
using MPBUtils
using Test

@testset "MPBUtils.jl" begin
    # write tests here
end

ext: jl


text:
module Flight
using ..OneD, ..Pmm, ..Mpmm, ..SixDof, ...Types, ..PmmZeroDrag
using DifferentialEquations
export timeOfFlight

function rangeAtImpact()
end
function impactVelocity()
end
function impactAngle()
end
function slantRange()
end

end


text:
abstract RustModelCore

#Dense, Infinite-horizon RustModelCore
immutable DIRMC <: RustModelCore
	beta::Float64						# discount factor
	u::Array{Float64,2}					# u[k,a]: flow utilities
	v::Array{Float64,2}					# v[k,a]: choice-specific discounted value
	V::Array{Float64,1}					# V[k]: expected value
	pi::Array{Array{Float64,2},1}						# pi[a][k,k']: condition state transitions
	p::Array{Float64,2}      			# p[k,a]: conditional choice probabilities (CCP)
end

#Sparse, Infinite-horizon RustModelCore
immutable SIRMC <: RustModelCore
	beta::Float64						# discount factor
	u::Array{Float64,2}					# u[k,a]: flow utilities
	v::Array{Float64,2}					# v[k,a]: choice-specific discounted value
	V::Array{Float64,1}					# V[k]: expected value
	pi::Array{SparseMatrixCSC{Float64,Int},1}						# pi[a][k,k']: condition state transitions
	p::Array{Float64,2}      			# p[k,a]: conditional choice probabilities (CCP)
end

#Dense, Finite-horizon RustModelCore
immutable DFRMC <: RustModelCore
	beta::Float64						# discount factor
	u::Array{Float64,2}					# u[k,a]: flow utilities
	v::Array{Float64,2}					# v[k,a]: choice-specific discounted value
	V::Array{Float64,1}					# V[k]: expected value
	pi::Array{Array{Float64,2},1}						# pi[a][k,k']: condition state transitions
	p::Array{Float64,2}      			# p[k,a]: conditional choice probabilities (CCP)

	# for finite horizon models only
	horizonindex::Int 					# index of the FIRST state of the LAST day
end

#Sparse, Finite-horizon RustModelCore
immutable SFRMC <: RustModelCore
	beta::Float64						# discount factor
	u::Array{Float64,2}					# u[k,a]: flow utilities
	v::Array{Float64,2}					# v[k,a]: choice-specific discounted value
	V::Array{Float64,1}					# V[k]: expected value
	pi::Array{SparseMatrixCSC{Float64,Int},1}						# pi[a][k,k']: condition state transitions
	p::Array{Float64,2}      			# p[k,a]: conditional choice probabilities (CCP)

	# for finite horizon models only
	horizonindex::Int 					# index of the FIRST state of the LAST day

	# for :SparseFiniteHorizon only
	kindices::Array{Array{Int,1},2}		#holds the non-zero k' at a given [k,a]
end

function rustmodelcore(beta,pi::Array{Array{Float64,2},1})
	da=length(pi)
	dk=size(pi[1],1)
	DIRMC(
		beta,
		Array(Float64,(dk,da)),			#u
		zeros(dk,da),					#v  (needs to be zeros()?)
		Array(Float64,dk),				#V
		pi,
		Array(Float64,(dk,da))			#p
	)
end

function rustmodelcore(beta,pi::Array{SparseMatrixCSC{Float64,Int},1})
	da=length(pi)
	dk=size(pi[1],1)
	SIRMC(
		beta,
		Array(Float64,(dk,da)),			#u
		zeros(dk,da),					#v  (needs to be zeros()?)
		Array(Float64,dk),				#V
		pi,
		Array(Float64,(dk,da))			#p
	)
end

function rustmodelcore(beta,pi::Array{Array{Float64,2},1},horizonindex::Int)
	da=length(pi)
	dk=size(pi[1],1)
	DFRMC(
		beta,
		Array(Float64,(dk,da)),			#u
		zeros(dk,da),					#v  (needs to be zeros()?)
		Array(Float64,dk),				#V
		pi,
		Array(Float64,(dk,da)),			#p
		horizonindex
	)
end

function rustmodelcore(beta,pi::Array{SparseMatrixCSC{Float64,Int},1},horizonindex::Int)
	da=length(pi)
	dk=size(pi[1],1)
	kindices=Array(Array{Int64,1},dk,da)
	for ia=1:da
		temp=full(pi[ia]')
		for ik=1:dk
			kindices[ik,ia]=find(temp[:,ik].!=0)
		end
	end
	SFRMC(
		beta,
		Array(Float64,(dk,da)),			#u
		zeros(dk,da),					#v  (needs to be zeros()?)
		Array(Float64,dk),				#V
		pi,
		Array(Float64,(dk,da)),			#p
		horizonindex,
		kindices
	)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    OrdinaryKriging(γ)
    OrdinaryKriging(data, γ)

Ordinary Kriging with variogram model `γ`.

Optionally, pass the geospatial `data` to the [`fit`](@ref) function.
"""
struct OrdinaryKriging{G<:Variogram} <: KrigingEstimator
  γ::G
end

OrdinaryKriging(γ) = OrdinaryKriging{typeof(γ)}(γ)

OrdinaryKriging(data, γ) = GeoStatsBase.fit(OrdinaryKriging(γ), data)

nconstraints(::OrdinaryKriging) = 1

function set_constraints_lhs!(::OrdinaryKriging, LHS::AbstractMatrix, domain)
  T = eltype(LHS)
  LHS[end,:]   .= one(T)
  LHS[:,end]   .= one(T)
  LHS[end,end]  = zero(T)
  nothing
end

factorize(::OrdinaryKriging, LHS::AbstractMatrix) = bunchkaufman(Symmetric(LHS), check=false)

function set_constraints_rhs!(fitted::FittedKriging{<:OrdinaryKriging}, pₒ)
  RHS = fitted.state.RHS
  RHS[end] = one(eltype(RHS))
  nothing
end


text:
function RARL_DQN_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
         -𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* maximum(value(π, 𝒟[:sp]), dims=1)
end

function RARL_TD3_target(π, 𝒫, 𝒟, γ::Float32; i)
    ap, _ = exploration(𝒫[:π_smooth], 𝒟[:sp], π_on=π, i=i)
    y = -𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* min.(value(π, 𝒟[:sp], ap)...)
end

RARL_DQN(;kwargs...) = DQN(;c_loss=td_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=RARL_DQN_target)
RARL_TD3(;kwargs...) = DQN(;c_loss=double_Q_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=RARL_TD3_target)

function RARL(;𝒮_pro,
               𝒮_ant,
               px,
               log::NamedTuple=(;),
               train_pro_every::Int=1,
               train_ant_every::Int=1,
               buffer_size=1000, # Size of the buffer
               required_columns=Symbol[:x, :fail],
               buffer::ExperienceBuffer=ExperienceBuffer(𝒮_pro.S, 𝒮_pro.agent.space, buffer_size, required_columns), # The replay buffer
               buffer_init::Int=max(max(𝒮_pro.c_opt.batch_size, 𝒮_ant.c_opt.batch_size), 200) # Number of observations to initialize the buffer with
               )
        # Set buffer parameters to be consistent between solvers (Since buffer is shared)
        𝒮_pro.required_columns = required_columns
        𝒮_pro.buffer = buffer
        𝒮_pro.buffer_init = buffer_init
        𝒮_ant.required_columns = required_columns
        𝒮_ant.buffer = buffer
        𝒮_ant.buffer_init = buffer_init
        AdversarialOffPolicySolver(;𝒮_pro=𝒮_pro,
                                    𝒮_ant=𝒮_ant,
                                    px=px,
                                    train_pro_every=train_pro_every,
                                    train_ant_every=train_ant_every,
                                    log=LoggerParams(;dir="log/rarl", log...),)
end


text:
#===============================================================================

    state_functions.jl

    Some functions that states support regardless of the specific implementation.

===============================================================================#
function mb_state_energy(orbital::T, state::SpinState)::AbstractFloat where T<:Orbital
    """ Retrieves the energy for a many-body state (single spin species).
    """
    e = 0
    for k = 1:length(state)
        if state[k]>0
            e += orbital(k)
        end
    end
    return e
end

function mb_state_energy(orbital_up::T, orbital_down::S, state::FullState)::AbstractFloat where {T<:Orbital,S<:Orbital}
    """ Many-body energy for a full state (wraps the single-particle routine)..
    """
    state_up, state_down = f_to_s(state)
    return mb_state_energy(orbital_up, state_up) + mb_state_energy(orbital_down, state_down)
end


text:
function rotate_matrix!(mat :: Array{Int32, 2})
  N :: Int = size(mat, 1)
  for row = 1:int(N / 2)
    for col = row:(N - row)
      mat[row, col], mat[col, N - row + 1], mat[N - row + 1, N - col + 1], mat[N - col + 1, row] =
      mat[N - col + 1, row], mat[row, col], mat[col, N - row + 1], mat[N - row + 1, N - col + 1]
    end
  end
end

# Time  : O(n^2)
# Space : O(1)

using Base: Test

mat = Array{Int32, 2}([[1 2 3];[4 5 6];[7 8 9]])
rotate_matrix!(mat)
@test mat == Array{Int32, 2}([[7 4 1];[8 5 2];[9 6 3]])
rotate_matrix!(mat)
@test mat == Array{Int32, 2}([[9 8 7];[6 5 4];[3 2 1]])

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots
using QuadGK

function f(x)
    return (x^2-μ)/(2*pi)
end

D = 4
k = 1.
μ = pi^2
multiple = false
test_2fermion = false

σ⁺ = [0. 1.; 0. 0.]

if multiple
    σᶻ = [1. 0.; 0. -1.]
    Id = 1*Matrix(I,2,2)
    Id2 = 1*Matrix(I,D,D)
    KL = Constant(randn(4*D^2,4*D^2))
    KL = 0.5*(KL-KL')
    R1 = Constant(kron(kron(kron(randn(D,D),Id2),σ⁺),Id))
    R2 = Constant(kron(kron(kron(Id2,randn(D,D)),σᶻ),σ⁺))
    RLs = (R1,R2)
    QL = KL
    for R in RLs
        mul!(QL, R', R, -1/2, 1)
    end
    Ψ = InfiniteCMPS(QL, (R1,R2); gauge = :left)
    h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2])
else
    KL = Constant(randn(2*D,2*D))
    KL = 0.5*(KL-KL')
    R1 = Constant(kron(randn(D,D),σ⁺))
    RLs = (R1,)
    QL = KL
    for R in RLs
        mul!(QL, R', R, -1/2, 1)
    end
    Ψ = InfiniteCMPS(QL, (R1,); gauge = :left)
    h = k * (∂ψ[1]'*∂ψ[1]) - μ * (ψ[1]'*ψ[1])
end

H = ∫(h, (-Inf,+Inf))

alg1 = LBFGS(; verbosity = 2, maxiter = 1000000, gradtol = 1e-4);

if multiple
    Ψ, ρR, E, e, normgrad, numfg, history = groundstate5bis(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    # αs,fs, dfs1, dfs2 = groundstate5bis(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    # display(plot(αs,[dfs1,dfs2]))
    # gui()
else
    Ψ, ρR, E, e, normgrad, numfg, history = groundstate5(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    #αs,fs, dfs1, dfs2 = groundstate5(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    #display(plot(αs,[dfs1,dfs2]))
    #gui()
end

if test_2fermion
    Q = Ψ.Q
    R1 = Ψ.Rs[1]
    @show expval(ψ[1]*ψ[1],Ψ)[]
    @show expval(ψ[1]'*ψ[1],Ψ)[]
    @show E, quadgk(f,-sqrt(μ),sqrt(μ))

    σᶻ = [1. 0.; 0. -1.];
    Id = 1*Matrix(I,2,2);
    Id2 = 1*Matrix(I,D,D);
    q = Ψ.Q;
    r = Ψ.Rs[1];
    Q = Constant(kron(q[],Id2) + kron(Id2,q[]));
    R1 = Constant(kron(r[],Id2))
    R2 = Constant(kron(Id2,r[]))
    # Q = Constant(kron(q[],kron(Id2,Id)) + kron(kron(Id2,Id),q[]));
    # R1 = Constant(kron(r[],kron(Id2,Id)));
    # R2 = Constant(kron(kron(Id2,σᶻ),r[]));
    @show norm((R1*R1)[]),norm((R2*R2)[]),norm((R1*R2+R2*R1)[]),norm((R1*R2)[])
    @show norm(Q[]' + Q[] + (R1'*R1 + R2'*R2)[])
    Ψ = InfiniteCMPS(Q, (R1,R2); gauge = :left)
    h1 = k * (∂ψ[1]'*∂ψ[1]) - μ * (ψ[1]'*ψ[1])
    h2 = k * (∂ψ[2]'*∂ψ[2]) - μ * (ψ[2]'*ψ[2]);
    h = h1 + h2;
    @show expval(h1,Ψ)[]
    @show expval(h2,Ψ)[]
    @show expval(h,Ψ)[]
    @show expval(ψ[1]'*ψ[1],Ψ)[]
    @show expval(ψ[2]'*ψ[2],Ψ)[]
end

Λ = 10000.
KL = Constant(randn(2*D,2*D))
KL = 0.5*(KL-KL')
R1 = Constant(randn(2*D,2*D))
R2 = Constant(randn(2*D,2*D))
#Rls = (R1,)
RLs = (R1,R2)
QL = KL
for R in RLs
    mul!(QL, R', R, -1/2, 1)
end
Ψ = InfiniteCMPS(QL, RLs; gauge = :left)
#h = k * (∂ψ[1]'*∂ψ[1]) - μ * (ψ[1]'*ψ[1]) + Λ * (ψ[1]'*ψ[1]'*ψ[1]*ψ[1])
h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + Λ * (2*ψ[1]'*ψ[1]'*ψ[1]*ψ[1] + 2*ψ[2]'*ψ[2]'*ψ[2]*ψ[2] + 4*(ψ[1]*ψ[2] + ψ[2]*ψ[1])'*(ψ[1]*ψ[2]+ψ[2]*ψ[1]))
H = ∫(h, (-Inf,+Inf))
Ψ, ρR, E, e, normgrad, numfg, history = groundstate(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))

@show Q = Ψ.Q
@show R1 = Ψ.Rs[1]
@show expval(ψ[1]*ψ[1],Ψ)[]
@show E, quadgk(f,-sqrt(μ),sqrt(μ))


text:
function print_table(a::Array{T, 2},
                     header::Vector{String}) where T <: Union{Char, String}
    print(" | ")
    for i = 1:length(header)
        print(string(header[i]), " | ")
    end
    println()

    print(" | ")
    for i = 1:length(header)
        print("---", " | ")
    end
    println()

    for i = 1:size(a)[1]
        print(" | ")
        for j = 1:size(a)[2]
            print(string(a[i,j]), " | ")
        end
        println()
    end
end

alphabet = [char for char = 'a':'z']
offsets = Int.([0, 0, 2, 14, 18, 21, 24])
alphabet_array = Array{Char}(undef, 26, length(offsets))

for i = 1:length(offsets)
    alphabet_array[:,i] = vcat(alphabet[offsets[i]+1:26],alphabet[1:offsets[i]])
end

header = vcat(string.(offsets))
print_table(alphabet_array, header)


text:
# Logistic loss function

export LogisticLoss

"""
**Logistic loss**

    LogisticLoss(y, μ=1)

Returns the function
```math
f(x) = μ⋅∑_i log(1+exp(-y_i⋅x_i))
```
where `y` is an array and `μ` is a positive parameter.
"""
struct LogisticLoss{T <: AbstractArray, R <: Real} <: ProximableFunction
    y::T
    mu::R
    function LogisticLoss{T, R}(y::T, mu::R) where {T, R}
        if mu <= R(0)
            error("parameter mu must be positive")
        end
        new(y, mu)
    end
end

LogisticLoss(y::T, mu::R=1) where {R, T <: AbstractArray} = LogisticLoss{T, R}(y, mu)

is_separable(f::LogisticLoss) = true
is_convex(f::LogisticLoss) = true
is_smooth(f::LogisticLoss) = true
is_prox_accurate(f::LogisticLoss) = false

# f(x)  =  mu log(1 + exp(-y x))

function (f::LogisticLoss)(x::AbstractArray{R}) where R
    val = R(0)
    for k in eachindex(x)
        expyx = exp(f.y[k] * x[k])
        val += log(R(1) + R(1) / expyx)
    end
    return f.mu * val
end

# f'(x) = -mu y exp(-y x) / (1 + exp(-y x))
#       = -mu y / (1 + exp(y x))
#
# Lipschitz constant of gradient: (mu y)

function gradient!(g::AbstractArray{R}, f::LogisticLoss, x::AbstractArray{R}) where R
    val = R(0)
    for k in eachindex(x)
        expyx = exp(f.y[k] * x[k])
        g[k] = -f.mu * f.y[k] / (R(1) + expyx)
        val += log(R(1) + R(1) / expyx)
    end
    return f.mu * val
end

# Computing proximal operator:
# z = prox(f, x, gamma)
# <==> f'(z) + (z - x)/gamma = 0
# <==> (z - x)/gamma - mu y / (1 + exp(y z)) = 0
# <==> z - x - mu gamma y / (1 + exp(y z)) = 0
#
# Indicating the above condition as F(z) = 0, then
# ==> F'(z) = 1 - (mu gamma y^2 exp(y z))/(1+exp(y z))^2
#
# Newton's method (no damping) to compute z reads:
# z_{k+1} = z_k - F(z_k)/F'(z_k)
#
# To ensure convergence of Newton's method a damping is required.
# The damping coefficient could be computed by backtracking.
#
# Alternatively we can use gradient methods with constant step size.

function prox!(z::AbstractArray{R}, f::LogisticLoss, x::AbstractArray{R}, gamma::R=R(1)) where R
    c = R(1) / gamma # convexity modulus
    L = maximum(abs, f.mu .* f.y) + c # Lipschitz constants
    z .= x
    expyz = similar(z)
    Fz = similar(z)
    F1z = similar(z)
    for k = 1:20
        expyz .= exp.(f.y .* z)
        Fz .= z .- x .- f.mu * gamma * (f.y ./ (1 .+ expyz))
        z .-= Fz ./ L
    end
    expyz .= exp.(f.y .* z)
    val = R(0)
    for k in eachindex(expyz)
        val += log(R(1) + R(1)/expyz[k])
    end
    return f.mu * val
end


text:
export mwcs, mwcs_dvv

function getCoherence(dcs,ds1,ds2)
  coh = zeros(size(dcs))
  valids = intersect(findall(x -> abs.(x) > 0,ds1),findall(x -> abs.(x) > 0,ds2))
  coh[valids] = dcs[valids] ./ (ds1[valids] .* ds2[valids])
  coh[coh .> 1.0 ] .= 1.0
  return coh
end

"""

  mwcs(ref, cur, fmin, fmax, fs, tmin, window_length, window_step,
       smoothing_half_win)

Change in velocity measurement using the Moving Window Cross-Spectrum technique.

The current correlation `cur` is compared to the reference correlation `ref`.
Both time series are sliced in several overlapping windows.
Each slice is mean-adjusted and cosine-tapered (85% taper) before being Fourier-
transformed to the frequency domain.
``F_{cur}(ν)`` and ``F_{ref}(ν)`` are the first halves of the
Hermitian symmetric Fourier-transformed segments. The cross-spectrum
``X(ν)`` is defined as
``X(ν) = F_{ref}(ν) F_{cur}^*(ν)``
in which ``{}^*`` denotes the complex conjugation.
``X(ν)`` is then smoothed by convolution with a Hanning window.
The similarity of the two time-series is assessed using the cross-coherency
between energy densities in the frequency domain:
``C(ν) = \\frac{|\\overline{X(ν))}|}{\\sqrt{|\\overline{F_{ref}(ν)|^2} |\\overline{F_{cur}(ν)|^2}}}``
in which the over-line here represents the smoothing of the energy spectra for
``F_{ref}`` and ``F_{cur}`` and of the spectrum of ``X``. The mean
coherence for the segment is defined as the mean of ``C(ν)`` in the
frequency range of interest. The time-delay between the two cross correlations
is found in the unwrapped phase, ``ϕ(ν)``, of the cross spectrum and is
linearly proportional to frequency:
``ϕ_j = m. ν_j, m = 2 π δ t``
The time shift for each window between two signals is the slope ``m`` of a
weighted linear regression of the samples within the frequency band of interest.
The weights are those introduced by [Clarke2011]_,
which incorporate both the cross-spectral amplitude and cross-coherence, unlike
[Poupinet1984]_. The errors are estimated using the weights (thus the
coherence) and the squared misfit to the modelled slope:
``e_m = \\sqrt{\\sum_j{(\\frac{w_j ν_j}{\\sum_i{w_i ν_i^2}})^2}σ_{ϕ}^2}``
where ``w`` are weights, ``ν`` are cross-coherences and
``σ_{ϕ}^2`` is the squared misfit of the data to the modelled slope
and is calculated as ``σ_{ϕ}^2 = \\frac{\\sum_j(ϕ_j - m ν_j)^2}{N-1}``
The output of this process is a table containing, for each moving window: the
central time lag, the measured delay, its error and the mean coherence of the
segment.

# Arguments
- `ref::AbstractArray`: Reference correlation.
- `cur::AbstractArray`: Current correlation.
- `fmin:Float64`: minimum frequency in the correlation [Hz]
- `fmax:Float64`: maximum frequency in the correlation [Hz]
- `fs:Float64`: Sampling rate of `ref` and `cur` [Hz]
- `tmin:Float64`: The leftmost time lag [s]
- `window_len:Float64`: The moving window length [s]
- `window_step:Float64`: The step to jump for the moving window [s]
- `smoothing_half_win:Int`: Defines the half length of the smoothing hanning
                            window.

# Returns
- `time_axis:Array{Float64,1}`: Central time of each moving window [s]
- `dt:Array{Float64,1}`: dt for each moving window
- `err:Array{Float64,1}`: Errors for each moving window
- `mcoh:Array{Float64,1}`: Mean coherence for each moving window

This code is a Julia translation of the Python code from [MSNoise](https://github.com/ROBelgium/MSNoise/blob/master/msnoise/move2obspy.py).
"""
function mwcs(ref::AbstractArray,cur::AbstractArray,fmin::Float64,
                   fmax::Float64,fs::Float64,tmin::Float64,
                   window_length::Float64,window_step::Float64,
                   smoothing_half_win::Int)

    # create time axis for mwcs
    time_axis = Array(tmin + window_length / 2. : window_step : tmin +
                      length(ref) / fs - window_length / 2.)

    window_length_samples = convert(Int,window_length * fs)
    window_step_samples = convert(Int,window_step * fs)
    minind = 1:window_step_samples:length(ref) - window_length_samples
    padd = convert(Int,2 ^ (ceil(log2(abs(window_length_samples))) + 2))

    N = length(minind)
    dt = zeros(N)
    err = zeros(N)
    coh = zeros(N)
    time_axis = time_axis[1:N]

    # Find values in frequency range of interest
    freq_vec = FFTW.rfftfreq(padd,fs)
    index_range = findall(x -> x >= fmin && x <= fmax,freq_vec)
    cci = zeros(window_length_samples,N)
    cri = zeros(window_length_samples,N)

    # fill matrices
    for ii = 1:N
        cci[:,ii] = cur[minind[ii]:minind[ii]+window_length_samples-1]
        cri[:,ii] = ref[minind[ii]:minind[ii]+window_length_samples-1]
    end

    # preprocess
    SeisNoise.demean!(cci)
    SeisNoise.detrend!(cci)
    SeisNoise.taper!(cci,fs,max_percentage=0.425)
    SeisNoise.demean!(cri)
    SeisNoise.detrend!(cri)
    SeisNoise.taper!(cri,fs,max_percentage=0.425)

    # take fourier transform with padding
    fcur = rfft([cci;zeros(eltype(cci),padd-size(cci)[1],N)],1)
    fref = rfft([cri;zeros(eltype(cri),padd-size(cri)[1],N)],1)

    fcur2 = real(fcur).^2 + imag(fcur).^2
    fref2 = real(fref).^2 + imag(fref).^2

    # calculate cross-spectrum
    X = fref .* conj.(fcur)

    if smoothing_half_win != 0
        dcur = sqrt.(SeisNoise.smooth(fcur2, smoothing_half_win))
        dref = sqrt.(SeisNoise.smooth(fref2, smoothing_half_win))
        X = SeisNoise.smooth(X, smoothing_half_win)
    else
        dcur = sqrt.(fcur2)
        dref = sqrt.(fref2)
    end

    dcs = abs.(X)

    # Get coherence and its mean value
    coh = getCoherence(dcs,dref,dcur)
    mcoh = reshape(mean(coh[index_range,:],dims=1),N)

    # Get Weights
    w = 1. ./ (1. ./ coh[index_range,:].^2  .- 1.)
    w[coh[index_range,:] .>= 0.99] .= 1. / (1. / 0.9801 - 1.)
    w = sqrt.(w .* sqrt.(dcs[index_range,:]))
    w = real(w)

    # frequency array
    v = real(freq_vec[index_range]) .* 2π

    # phase
    phi = angle.(X)
    phi[1,:] .= 0.
    unwrap!(phi,dims=1)
    phi = phi[index_range,:]

    # Calculate the slope with a weighted least square linear regression
    # forced through the origin
    for ii = 1:N
        model = glm(@formula(Y ~0 + X),DataFrame(X=v,Y=phi[:,ii]),Normal(),
                    IdentityLink(),wts=w[:,ii])
        # a much faster, unweighted version of this is: m = v \ phi
        dt[ii] = coef(model)[1]

        # Errors
        e = sum((phi[:,ii] .- dt[ii] .* v).^2) / (length(v) - 1)
        s2x2 = sum(v.^2 .* w[:,ii].^2)
        sx2 = sum(w[:,ii] .* v.^2)
        err[ii] = sqrt(e * s2x2 / sx2^2)
    end

    return time_axis, dt, err, mcoh
end

"""
    mwcs_dvv(time_axis,dt,err,coh,dtt_lag,dist,dtt_v,dtt_minlag,dtt_width,
             dtt_sides,min_coh, max_err, max_dt)

Regresses dv/v from dt/t measurements.

Solves dt = a + bt, where b = dv/v, a = instrumental drift, dt = time lag
at time t. Solves with a weighted linear regression with weights equal to
1/error**2.

# Arguments
- `time_axis:Array{Float64,1}`: Central time of each moving window [s]
- `dt:Array{Float64,1}`: dt for each moving window
- `err:Array{Float64,1}`: Errors for each moving window
- `coh:Array{Float64,1}`: Mean coherence for each moving window
- `dist:Float64`: Distance between stations [km]
- `dtt_lag:String`: Type of time lag 'dynamic' or 'static'. When `dtt_lag`
    is set to "dynamic", the inter-station distance is used to determine
    the minimum time lag
- `dtt_v:Float64`: Velocity for minumum time lag. The velocity is determined by
    the user so that the minlag doesn't include the ballistic waves.
    If ballistic waves are visible with a velocity of 2 km/s, one could
    configure dtt_v=1.0. If stations are located 15 km apart, the minimum
    lag time will be set to 15 s.
- `minlag:Float64`: Statically set minimum lag [s]
- `dtt_width:Float64`: Width of the lag window used [s] i.e., if
    dtt_width = 30, and min_lag = 15, window = 15s - 45s
- `dtt_sides:String`: Sides of cross-correlation function to use. Either
    'Both' or 'left'.
- `min_coh:Float64`: Minimum allowed coherency between reference and current
                     correlation
- `max_err:Float64`: Maximum allowed error in dt/t regression
-  `max_dt:Float64`: Maximum allowed dt from MWCS [s]

# Returns
- `m::Float64`: dt/t for current correlation
- `em::Float64`: Error for calculation of `m`
- `a::Float64`: Intercept for regression calculation
- `ea::Float64`: Error on intercept
- `m0::Float64`: dt/t for current correlation with no intercept
- `em0::Float64`: Error for calculation of `m0`
"""
function mwcs_dvv(time_axis::AbstractArray, dt::AbstractArray,
                     err::AbstractArray, coh::AbstractArray, dtt_lag::String,
                     dist::Float64, dtt_v::Float64, dtt_minlag::Float64,
                     dtt_width::Float64, dtt_sides::String;
                     min_coh::Float64=0.7, max_err::Float64=1.,
                     max_dt::Float64=1.)

    # lag bins
    if dtt_lag == "static"
        lmlag = -dtt_minlag
        rmlag = dtt_minlag
    else
        lmlag = - dist / dtt_v
        rmlag = dist / dtt_v
    end

    lMlag = lmlag - dtt_width
    rMlag = rmlag + dtt_width

    # subset correlation by sides
    if dtt_sides == "both"
        tindex = findall(x -> abs.(x) .>= rmlag && abs.(x) .<= rMlag, time_axis)
    elseif dtt_sides == "left"
        tindex = findall(x -> x .>= lMlag && x .<= lmlag, time_axis)
    else
        tindex = findall(x -> x .>= rmlag && x .<= rMlag, time_axis)
    end

    # index to subset data
    dt, err, coh = dt[tindex], err[tindex], coh[tindex]
    index = intersect(findall(x -> x .<= min_coh,coh),
                      findall(x -> x .>= max_err,err),
                      findall(x -> abs.(x) .>= max_dt,dt))
    dt[index] .= 0.
    err[index] .= 1.
    coh[index] .= 1.

    # weight statistics
    w = 1. ./ err .^ 2

    # regress data using least squares
    VecXfilt = time_axis[tindex]
    model0 = glm(@formula(Y ~0 + X),DataFrame(X=VecXfilt,Y=dt),Normal(),
                IdentityLink(),wts=w)
    model = glm(@formula(Y ~ X),DataFrame(X=VecXfilt,Y=dt),Normal(),
                IdentityLink(),wts=w)

    a,m = coef(model)
    ea, em = stderror(model)

    m0 = coef(model0)[1]
    em0 = stderror(model0)[1]
    return m, em, a, ea, m0, em0
end


text:
using JuMP, Gurobi, OneClassSampling, SVDD, OneClassActiveLearning
include(joinpath(@__DIR__, "config.jl"))

experiment_name="dami_baseline_gt"
data_parameter_file = joinpath(data_input_root, "parameters_gt.jser")

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
quality_metrics = Dict(:mcc => matthews_corr, :kappa => cohens_kappa, :f1 => f1_score)

model_init_strategy = Dict{Symbol, Any}(
                            :svdd => Dict{Symbol, Any}(),
                            :sample_qe => Dict{Symbol, Any}(
                                :threshold_strategies => []
                        ))

sampling_strategies = Vector{Dict{Symbol, Any}}()

push!(sampling_strategies, Dict{Symbol, Any}(
    :method => RandomRatioSampler(1.0),
    :seed => 0
))

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
using CSV
using DataFrames
using Plots
using Statistics

pyplot()

function convert_to_array(str)
    els = split(str, ',')

    if "n/a" in els
        @warn "Benchmark contains 'n/a' values. These will be ignored."
        filter!(el -> el != "n/a", els)
    end

    return parse.(Float64, els)
end

labels = Dict(
              "gemmkernels_complex" => "Our implementation (complex)",
              "gemmkernels_dual" => "Our implementation (dual)",
              "cutlass" => "CUTLASS (complex)",
              "cudajl_complex" => "CUDA.jl (complex)",
             )

seriesnumbers = Dict(
              "gemmkernels_complex" => 1,
              "gemmkernels_dual" => 2,
              "cutlass" => 3,
              "cudajl_complex" => 4,
             )

flop_factors = Dict(
              "gemmkernels_complex" => 8,
              "gemmkernels_dual" => 6,
              "cutlass" => 8,
              "cudajl_complex" => 8
                   )

for file in readdir()
    if isfile(file) && splitext(basename(file))[2] == ".csv"
        filename = splitext(basename(file))[1]

        implementation = filename # e.g. gemmkernels_complex

        label = labels[implementation]
        seriesnr = seriesnumbers[implementation]
        markershape = :circle

        df = DataFrame(CSV.File(file))

        N = df[!, :N]
        runtime_arr = convert_to_array.(df[!, :runtime]) .* 1e3 # in ps
        flops_factor = flop_factors[implementation]
        tflops = [flops_factor * N[i] ^ 3 ./ runtime_arr[i] for i = 1 : length(N)]

        plot!(N, mean.(tflops), seriescolor=seriesnr, label=label, xscale=:log2, markershape=markershape, ribbon=std.(tflops), fillalpha=.5)
    end
end

title!("Performance of complex/dual GEMMs")
xlabel!("N")
ylabel!("TFLOPS")
savefig("plot.pdf")


text:
struct PRESS
    # Metadata
    echo_time # Echo time

    # Indices into underlying acqusition list
    ref_scans::Vector{Int}
    press_scans::Matrix{Int}

    navigator::Vector{Int}

    # Raw acquisitions
    acquisitions

    cycle_length
end

standard_metadata(p::PRESS) = standard_metadata(p.acquisitions)

"""
A standard Point Resolved Spectroscopy (PRESS) experiment with `num_averages`
acquisitions, possibly with navigator and reference scans included.

This is implemented as "single voxel spectroscopy" and is the standard product
sequence for in vivo MR spectroscopy on Siemens scanners as of early 2020.
"""
function PRESS(
    echo_time,
    ref_scans::Vector{Int},
    press_scans::Matrix{Int},
    navigator::Vector{Int},
    acquisitions
)
    PRESS(
        echo_time,
        ref_scans,
        press_scans,
        navigator,
        acquisitions,
        count_cycles(acquisitions)
    )
end

function Base.show(io::IO, press::PRESS)
    println(io, """
                PRESS experiment:
                  size(press_scans) = $(size(press.press_scans))
                  length(ref_scans) = $(length(press.ref_scans))
                Navigator experiment:
                  length(navigator) = $(length(press.navigator))""")
    show(io, press.acquisitions)
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# Tests of data checking and transformation code

@testset "test idempotent units transformations" begin
    @testset "3-bus case" begin
        data = PowerModels.parse_file("../test/data/case3.m")
        data_base = deepcopy(data)

        PowerModels.make_mixed_units(data)
        PowerModels.make_per_unit(data)

        @test compare_dict(data, data_base)
    end
    @testset "5-bus case" begin
        data = PowerModels.parse_file("../test/data/case5_asym.m")
        data_base = deepcopy(data)

        PowerModels.make_mixed_units(data)
        PowerModels.make_per_unit(data)

        @test compare_dict(data, data_base)
    end
    @testset "5-bus case with pwl costs" begin
        data = PowerModels.parse_file("../test/data/case5_pwlc.m")
        data_base = deepcopy(data)

        PowerModels.make_mixed_units(data)
        PowerModels.make_per_unit(data)

        @test compare_dict(data, data_base)
    end
    @testset "24-bus case" begin
        data = PowerModels.parse_file("../test/data/case24.m")
        data_base = deepcopy(data)

        PowerModels.make_mixed_units(data)
        PowerModels.make_per_unit(data)

        @test compare_dict(data, data_base)
    end

    @testset "3-bus case solution" begin
        result = run_ac_opf("../test/data/case3.m", ipopt_solver)
        result_base = deepcopy(result)

        PowerModels.make_mixed_units(result["solution"])
        PowerModels.make_per_unit(result["solution"])

        @test compare_dict(result, result_base)
    end
    @testset "5-bus case solution" begin
        result = run_ac_opf("../test/data/case5_asym.m", ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))
        result_base = deepcopy(result)

        PowerModels.make_mixed_units(result["solution"])
        PowerModels.make_per_unit(result["solution"])

        @test compare_dict(result, result_base)
    end
    @testset "24-bus case solution" begin
        result = run_ac_opf("../test/data/case24.m", ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))
        result_base = deepcopy(result)

        PowerModels.make_mixed_units(result["solution"])
        PowerModels.make_per_unit(result["solution"])

        @test compare_dict(result, result_base)
    end

    @testset "5-bus case solution with duals" begin
        result = run_dc_opf("../test/data/case5.m", ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true, "duals" => true)))
        result_base = deepcopy(result)

        PowerModels.make_mixed_units(result["solution"])
        PowerModels.make_per_unit(result["solution"])

        @test compare_dict(result, result_base)
    end

end

@testset "test user ext init" begin
    @testset "3-bus case" begin
        pm = build_generic_model("../test/data/case3.m", ACPPowerModel, PowerModels.post_opf, ext = Dict(:some_data => "bloop"))

        #println(pm.ext)

        @test haskey(pm.ext, :some_data)
        @test pm.ext[:some_data] == "bloop"

        result = solve_generic_model(pm, IpoptSolver(print_level=0))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 5907; atol = 1e0)
    end
end


text:
type HullWhite{AffineModelType, T <: TermStructure} <: OneFactorModel{AffineModelType}
  modT::AffineModelType
  r0::Float64
  a::ConstantParameter
  sigma::ConstantParameter
  phi::HullWhiteFittingParameter
  ts::T
  privateConstraint::PrivateConstraint
  common::ModelCommon
end

function HullWhite{T <: TermStructure}(ts::T, a::Float64 = 0.1, sigma::Float64 = 0.01)
  _rate = forward_rate(ts, 0.0, 0.0, ContinuousCompounding(), NoFrequency())
  r0 = _rate.rate
  a_const = ConstantParameter([a], PositiveConstraint())
  sigma_const = ConstantParameter([sigma], PositiveConstraint())

  privateConstraint = PrivateConstraint(ConstantParameter[a_const, sigma_const])

  phi  = HullWhiteFittingParameter(a, sigma, ts)

  return HullWhite(AffineModelType(), r0, a_const, sigma_const, phi, ts, privateConstraint, ModelCommon()) #, ShortRateModelCommon())
end

## Dynamics ##
type HullWhiteDynamics{P <: Parameter} <: ShortRateDynamics
  process::OrnsteinUhlenbeckProcess
  fitting::P
  a::Float64
  sigma::Float64
end

HullWhiteDynamics{P <: Parameter}(fitting::P, a::Float64, sigma::Float64) = HullWhiteDynamics(OrnsteinUhlenbeckProcess(a, sigma), fitting, a, sigma)

short_rate(dynamic::HullWhiteDynamics, t::Float64, x::Float64) = x + dynamic.fitting(t)

generate_arguments!(m::HullWhite) = m.phi = HullWhiteFittingParameter(get_a(m), get_sigma(m), m.ts)

get_dynamics(m::HullWhite) = HullWhiteDynamics(m.phi, get_a(m), get_sigma(m))

## Tree methods ##
function tree(model::HullWhite, grid::TimeGrid)
  phi = TermStructureFittingParameter(model.ts)
  numericDynamics = HullWhiteDynamics(phi, get_a(model), get_sigma(model))
  trinomial = TrinomialTree(numericDynamics.process, grid)
  numericTree = OneFactorShortRateTree{HullWhiteDynamics}(trinomial, numericDynamics, grid)

  reset_param_impl!(phi)

  @simd for i = 1:length(grid.times) - 1
    @inbounds discountBond = discount(model.ts, grid.times[i + 1])
    statePrices = get_state_prices!(numericTree, i)
    sz = get_size(numericTree, i)
    @inbounds dt = grid.dt[i]
    @inbounds dx = trinomial.dx[i]
    x = get_underlying(trinomial, i, 1)
    val = 0.0
    for j = 1:sz
      @inbounds val += statePrices[j] * exp(-x * dt)
      x += dx
    end
    val = log(val / discountBond) / dt
    @inbounds set_params!(phi, grid.times[i], val)
  end

  return numericTree
end

# Additional methods #
function B(model::HullWhite, t::Float64, T::Float64)
  _a = get_a(model)
  if _a < sqrt(eps())
    return T - t
  else
    return (-expm1(-_a * (T - t))) / _a
  end
end

function A(model::HullWhite, t::Float64, T::Float64)
  discount1 = discount(model.ts, t)
  discount2 = discount(model.ts, T)

  forward = forward_rate(model.ts, t, t, ContinuousCompounding(), NoFrequency())
  temp = get_sigma(model) * B(model, t, T)
  val = B(model, t, T) * forward.rate - 0.25 * temp * temp * B(model, 0.0, 2.0* t)

  return exp(val) * discount2 / discount1
end

function discount_bond_option{O <: OptionType}(model::HullWhite, optionType::O, strike::Float64, maturity::Float64, bondStart::Float64, bondMaturity::Float64)
  _a = get_a(model)
  if _a < sqrt(eps())
    v = get_sigma(model) * B(model, bondStart, bondMaturity) * sqrt(maturity)
  else
    v = get_sigma(model) / (_a * sqrt(2.0 * _a)) * sqrt(exp(-2.0 * _a * (bondStart - maturity)) - exp(-2.0 * _a * bondStart) - 2.0 * (exp(-_a * (bondStart + bondMaturity - 2.0 * maturity))
        - exp(-_a * (bondStart + bondMaturity))) + exp(-2.0 * _a * (bondMaturity - maturity)) - exp(-2.0 * _a * bondMaturity))
  end

   f = discount(model.ts, bondMaturity)
   k = discount(model.ts, bondStart) * strike

   return black_formula(optionType, k, f, v)
end


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

unsafe_write(::Val{:precision}, key, name, range, vals) =
    @fame_call_check(fame_write_precisions,
        (Cint, Cstring, Ptr{FameRange}, Ref{Float64}),
        key, name, range, vals)

unsafe_write(::Val{:numeric}, key, name, range, vals) =
    @fame_call_check(fame_write_numerics,
        (Cint, Cstring, Ptr{FameRange}, Ref{Float32}),
        key, name, range, vals)

unsafe_write(::Val{:boolean}, key, name, range, vals) =
    @fame_call_check(fame_write_booleans,
        (Cint, Cstring, Ptr{FameRange}, Ref{Int32}),
        key, name, range, vals)

@inline function unsafe_write(::Val{TY}, key, name, range, vals) where {TY}
    @assert(haskey(fame_freq, TY), "Cannot write data of type $TY.")
    @fame_call_check(fame_write_dates,
        (Cint, Cstring, Ptr{FameRange}, Cint, Ref{FameDate}),
        key, name, range, val_to_int(TY, fame_freq), vals)
end

@inline function unsafe_write(::Val{:string}, key, name, range, vals::Ref{String})
    tmp = [vals[]]
    @fame_call_check(fame_write_strings,
        (Cint, Cstring, Ptr{FameRange}, Ptr{Ptr{UInt8}}),
        key, name, range, tmp)

end

@inline function unsafe_write(::Val{:string}, key, name, range, vals::AbstractVector{String})
    @fame_call_check(fame_write_strings,
        (Cint, Cstring, Ptr{FameRange}, Ptr{Ptr{UInt8}}),
        key, name, range, vals)
end

@inline function unsafe_write(::Val{:namelist}, key, name, range, vals::Ref{String})
    tmp = uppercase(vals[])
    @cfm_call_check(cfmwtnl,
        (Cint, Cstring, Cint, Ptr{UInt8}),
        key, name, HNLALL, tmp)
end

@inline function do_delete_object(key, name)
    try
        @cfm_call_check(cfmdlob, (Cint, Cstring), key, name)
    catch e
        (e isa HLIError && e.status == HNOOBJ) || rethrow()
    end
end

@inline int_class_type_freq(v::FameObject{CL,TY,FR}) where {CL,TY,FR} =
    (
        val_to_int(CL, fame_class),
        try
            val_to_int(TY, fame_type)
        catch
            val_to_int(TY, fame_freq)
        end,
        val_to_int(FR, fame_freq)
    )

@inline do_new_object(key, fo::FameObject) = begin
    (cl, ty, fr) = int_class_type_freq(fo)
    @cfm_call_check(cfmnwob, (Cint, Cstring, Cint, Cint, Cint, Cint, Cint),
        key, fo.name, cl, fr, ty, fame_basis.daily,
        fame_observed[eltype(fo.data) <: AbstractFloat ? :summed : :undefined])
end

_get_range(::FameObject{:scalar}) = C_NULL
_get_range(fo::FameObject{:series,TY,FR}) where {TY,FR} =
    Ref(FameRange(val_to_int(FR, fame_freq), fo.first_index[], fo.last_index[]))

"""
    do_write(obj::FameObject, db::FameDatabase)

Write the given [`FameObject`](@ref) to the given [`FameDatabase`](@ref). If an
object by the same name already exists, it is deleted before the new object is
written.
"""
function do_write end
export do_write

function do_write(obj::FameObject{CL,FT}, db::FameDatabase) where {CL,FT}
    do_delete_object(db.key, obj.name)
    do_new_object(db.key, obj)
    range = _get_range(obj)
    (CL == :series) && (range[].r_start == FAME_INDEX_NC || range[].r_end == FAME_INDEX_NC) && return
    unsafe_write(Val(FT), db.key, obj.name, range, obj.data)
end


text:
using LightGraphs; const lg=LightGraphs
using SimpleWeightedGraphs
using LinearAlgebra
using GraphPlot
using Colors
using JSON
include("pmg_graphs.jl")

# Type to store atomic graphs
# TO CONSIDER: store ref to featurization rather than the thing itself? Does this matter for any performance we care about?
# TO CONSIDER: store an ID of some kind? (e.g. mp-123 or whatever)
mutable struct AtomGraph <: lg.AbstractGraph{Float32}
    graph::SimpleWeightedGraph{Int32,Float32}
    elements::Vector{String} # list of elemental symbols corresponding to each node
    lapl::Matrix{Float32} # graph laplacian (normalized)
    features::Matrix{Float32} # feature matrix (size (# features, # nodes))
    featurization::Vector{AtomFeat} # featurization scheme in the form of a list of AtomFeat objects
end

# basic constructor
function AtomGraph(gr::SimpleWeightedGraph{Int32,Float32}, el_list::Vector{String}, features::Matrix{Float32}, featurization::Vector{AtomFeat})
    # check that el_list is the right length
    num_atoms = size(gr)[1]
    @assert length(el_list)==num_atoms "Element list length doesn't match graph size!"

    # check that features is the right dimensions (# features x # nodes)
    expected_feature_length = sum([f.num_bins for f in featurization])
    @assert size(features) == (expected_feature_length, num_atoms) "Feature matrix is of wrong dimension! It should be of size (# features, # nodes)"

    # if all these are good, calculate laplacian and build the thing
    lapl = normalized_laplacian(gr)
    AtomGraph(gr, el_list, lapl, features, featurization)
end

# one without features or featurization initialized yet
function AtomGraph(gr::SimpleWeightedGraph{Int32,Float32}, el_list::Vector{String})
    # check that el_list is the right length
    num_atoms = size(gr)[1]
    @assert length(el_list)==num_atoms "Element list length doesn't match graph size!"

    lapl = Float32.(normalized_laplacian(gr))
    AtomGraph(gr, el_list, lapl, zeros(Float32, 1, num_atoms), AtomFeat[])
end

# TODO, maybe: constructor where you give adjacency matrix and it builds the graph for you also

# TODO: fix so that it will print feature length even if featurization isn't initialized
# pretty printing, short version
function Base.show(io::IO, g::AtomGraph)
    st = "AtomGraph with $(nv(g)) nodes, $(ne(g)) edges"
    if length(g.featurization)!=0
        st = string(st, ", feature vector length $(size(g.features)[1])")
    end
    print(io, st)
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", g::AtomGraph)
    st = "AtomGraph with $(nv(g)) nodes, $(ne(g)) edges\n   atoms: $(g.elements)\n   feature vector length: "
    if length(g.featurization)==0
        st = string(st, "uninitialized\n   encoded features: uninitialized")
    else
        st = string(st, "$(size(g.features)[1])\n   encoded features: ",  [string(f.name, ", ") for f in g.featurization]...)[1:end-2]
    end
    print(io, st)
end

# now define some functions...

# first, the LightGraphs required ones (easy, just run them on the graph contents)
lg.edges(g::AtomGraph) = lg.edges(g.graph)
lg.eltype(g::AtomGraph) = lg.eltype(g.graph)
lg.edgetype(g::AtomGraph) = lg.edgetype(g.graph)
lg.ne(g::AtomGraph) = lg.ne(g.graph)
lg.nv(g::AtomGraph) = lg.nv(g.graph)
lg.weights(g::AtomGraph) = lg.weights(g.graph)
lg.is_directed(g::AtomGraph) = lg.is_directed(g.graph)

lg.outneighbors(g::AtomGraph, node) = lg.outneighbors(g.graph, node)
lg.inneighbors(g::AtomGraph, node) = lg.inneighbors(g.graph, node)
lg.has_vertex(g::AtomGraph, v::Integer) = lg.has_vertex(g.graph, v)
lg.has_edge(g::AtomGraph, i, j) = lg.has_edge(g.graph, i, j)

lg.zero(AtomGraph) = AtomGraph(zero(SimpleWeightedGraph{Int32,Float32}), String[])

# TODO: maybe some subgraph stuff for cutting up graphs later, will need to make sure to account for feature matrix and element list properly as well as recompute laplacian, would be cool to be able to e.g. filter on particular features and only pull matching nodes

# this cribbed from GeometricFlux
function normalized_laplacian(g::G) where G<:lg.AbstractGraph
    a = adjacency_matrix(g)
    d = vec(sum(a, dims=1))
    inv_sqrt_d = diagm(0=>d.^(-0.5))
    Float32.(I - inv_sqrt_d * a * inv_sqrt_d)
end

normalized_laplacian(g::AtomGraph) = g.lapl

# function to add/change node features, note that featurization scheme must be specified!
function add_features!(g::AtomGraph, features::Matrix{Float32}, featurization::Vector{AtomFeat})
    num_atoms = nv(g)

    # check that features is the right dimensions (# features x # nodes)
    expected_feature_length = sum([f.num_bins for f in featurization])
    @assert size(features) == (expected_feature_length, num_atoms) "Feature matrix is of wrong dimension! It should be of size (# features, # nodes)"

    # okay now we can set the features
    g.features = features
    g.featurization = featurization
end

# alternate version where it builds the features too, you have to pass in the results of the make_feature_vectors function
function add_features!(g::AtomGraph, atom_feature_vecs::Dict{String, Vector{Float32}}, featurization::Vector{AtomFeat})
    @assert Set(String.(g.elements)) <= Set(keys(atom_feature_vecs)) "Some atoms in your graph do not have corresponding feature vectors! This could be because some features you requested had missing values for these atoms."
    feature_mat = Float32.(hcat([atom_feature_vecs[e] for e in g.elements]...))
    add_features!(g, feature_mat, featurization)
end

# and finally the ones where it makes the feature vectors too...(defined for both signatures of the make_feature_vectors function just for completeness)
function add_features!(g::AtomGraph, featurization::Vector{AtomFeat})
    feature_vecs, featurization = make_feature_vectors(featurization)
    add_features!(g, feature_vecs, featurization)
end

# is there a clever way to roll this into the previous one since the syntax is identical?
function add_features!(g::AtomGraph, feature_names::Vector{Symbol})
    feature_vecs, featurization = make_feature_vectors(feature_names)
    add_features!(g, feature_vecs, featurization)
end

function add_features!(g::AtomGraph, feature_names::Vector{Symbol}, nbins::Vector{<:Integer}, logspaced=false)
    feature_vecs, featurization = make_feature_vectors(feature_names, nbins=nbins, logspaced=logspaced)
    add_features!(g, feature_vecs, featurization)
end

# now visualization stuff...

"Get a list of colors to use for graph visualization."
function graph_colors(atno_list, seed_color=colorant"cyan4")
    atom_types = unique(atno_list)
    atom_type_inds = Dict(atom_types[i]=>i for i in 1:length(atom_types))
    color_inds = [atom_type_inds[i] for i in atno_list]
    colors = distinguishable_colors(length(atom_types), seed_color)
    return colors[color_inds]
end

"Compute edge widths (proportional to weights on graph) for graph visualization."
function graph_edgewidths(ag::AtomGraph)
    edgewidths = []
    # should be able to do this as
    for e in edges(ag)
        append!(edgewidths, ag.weights[e.src, e.dst])
    end
    return edgewidths
end

"Visualize a given graph."
function visualize_graph(ag::AtomGraph)
    # gplot doesn't work on weighted graphs
    sg = SimpleGraph(adjacency_matrix(ag))
    plt = gplot(sg, nodefillc=graph_colors(ag.elements), nodelabel=ag.elements, edgelinewidth=graph_edgewidths(sg, weights(ag)))
    display(plt)
end


text:
#Simple methods

export eigvals_power, eigvals_ii, eigvals_rqi

function powm_method!{T}(log::ConvergenceHistory, K::KrylovSubspace{T};
    tol::Real=eps(T)*K.n^3, maxiter::Int=K.n
    )
    θ = zero(T)
    v = Array(T, K.n)
    for iter=1:maxiter
        nextiter!(log,mvps=1)
        v = lastvec(K)
        y = nextvec(K)
        θ = dot(v, y)
        resnorm = real(norm(y - θ*v))
        push!(log, :resnorm, resnorm)
        resnorm <= tol*abs(θ) && (setconv(log, resnorm >= 0); break)
        appendunit!(K, y)
    end
    Eigenpair(θ, v)
end

function eigvals_power(A, x=nothing; tol::Real=size(A,1)^3*eps(), maxiter::Int=size(A,1))
    K = KrylovSubspace(A, 1)
    x==nothing ? initrand!(K) : init!(K, x/norm(x))
    eigvals_power(K; tol=tol, maxiter=maxiter)
end
function eigvals_power(K::KrylovSubspace, x=nothing;
    tol::Real=size(A,1)^3*eps(), maxiter::Int=size(A,1)
    )
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    pair = powm_method!(history, K; tol=tol, maxiter=maxiter)
    pair, history
end

#Inverse iteration/inverse power method
function invpowm_method!{T}(log::ConvergenceHistory, K::KrylovSubspace{T}, σ::Number=zero(T);
    tol::Real=eps(T)*K.n^3, maxiter::Int=K.n
    )
    θ = zero(T)
    v = Array(T, K.n)
    y = Array(T, K.n)
    σ = convert(T, σ)
    for iter=1:maxiter
        nextiter!(log,mvps=1)
        v = lastvec(K)
        y = (K.A-σ*eye(K))\v
        θ = dot(v, y)
        resnorm = norm(y - θ*v)
        push!(log, :resnorm, resnorm)
        resnorm <= tol*abs(θ) && (setconv(log, resnorm >= 0); break)
        appendunit!(K, y)
    end
    shrink!(log)
    Eigenpair(σ+1/θ, y/θ)
end

function eigvals_ii(A, σ::Number, x=nothing; tol::Real=eps(), maxiter::Int=size(A,1))
    K = KrylovSubspace(A, 1)
    x==nothing ? initrand!(K) : init!(K, x/norm(x))
    eigvals_ii(K, σ; tol=tol, maxiter=maxiter)
end
function eigvals_ii(K::KrylovSubspace, σ::Number, x=nothing;
    tol::Real=size(A,1)^3*eps(), maxiter::Int=size(A,1)
    )
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    pair = invpowm_method!(history, K, σ; tol=tol, maxiter=maxiter)
    pair, history
end

#Rayleigh quotient iteration
#XXX Doesn't work well
function rqi_method!{T}(log::ConvergenceHistory, K::KrylovSubspace{T}, σ::Number; tol::Real=eps(T), maxiter::Int=K.n)
    v = lastvec(K)
    ρ = dot(v, nextvec(K))
    for iter=1:maxiter
        nextiter!(log,mvps=1)
        y = (K.A-ρ*eye(K))\v
        θ = norm(y)
        ρ += dot(y,v)/θ^2
        v = y/θ
        resnorm=1/θ
        push!(log,:resnorm,resnorm)
        θ >= 1/tol && (setconv(log, resnorm >= 0); break)
    end
    shrink!(history)
    Eigenpair(ρ, v)
end

function eigvals_rqi(A, σ::Number, x=nothing; tol::Real=eps(), maxiter::Int=size(A,1))
    K = KrylovSubspace(A, 1)
    x==nothing ? initrand!(K) : init!(K, x/norm(x))
    eigvals_rqi(K, σ; tol=tol, maxiter=maxiter)
end
function eigvals_rqi{T}(K::KrylovSubspace{T}, σ::Number, x=nothing; tol::Real=eps(T), maxiter::Int=K.n)
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    rqi_method(history, K, σ, tol=tol, maxiter=maxiter)
    x, history
end


text:
module ThreadingUtilities

using VectorizationBase:
    pause, StaticInt, StridedPointer, stridedpointer, offsets, cache_linesize, align, __vload, __vstore!, num_threads, assume, False, register_size, NativeTypes

@enum ThreadState::UInt begin
    SPIN = 0   # 0: spinning
    WAIT = 1   # 1: waiting, check if ≤ 1 to see if task is free and waiting
    TASK = 2   # 2: task available
end
const TASKS = Task[]
const THREADBUFFERSIZE = 512
const THREADPOOL = UInt[]
const THREADPOOLPTR =  Ref{Ptr{UInt}}(C_NULL);

include("atomics.jl")
include("threadtasks.jl")
include("utils.jl")
include("warnings.jl")

function initialize_task(tid::Int)
    t = Task(ThreadTask(taskpointer(tid)));
    t.sticky = true # create and pin
    # set to tid, we have tasks 2...nthread, from 1-based ind perspective
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), t, tid % Cint)
    TASKS[tid] = t
    wake_thread!(tid) # task should immediately sleep
end
function reinitialize_tasks!(verbose::Bool = true)
    for (tid,t) ∈ enumerate(TASKS)
        if istaskfailed(t)
            verbose && dump(t)
            _atomic_store!(taskpointer(tid), reinterpret(UInt, SPIN))
            initialize_task(tid)
        end
    end
end

function __init__()
    _print_exclusivity_warning()
    nt = min(Threads.nthreads(),(Sys.CPU_THREADS)::Int) - 1
    resize!(THREADPOOL, (THREADBUFFERSIZE ÷ sizeof(UInt)) * nt + (cache_linesize() ÷ sizeof(UInt)) - 1)
    THREADPOOL .= 0
    Threads.atomic_fence() # ensure 0-initialization
    resize!(TASKS, nt)
    GC.@preserve THREADPOOL begin
        THREADPOOLPTR[] = align(pointer(THREADPOOL)) - THREADBUFFERSIZE
        foreach(initialize_task, 1:nt)
        for tid ∈ 1:nt
            # wait for it to sleep, to be sure
            while _atomic_load(taskpointer(tid)) ≠ reinterpret(UInt, WAIT)
                yield()
                # pause()
            end
        end
    end
end

end # module


text:
export bfs, bfsdistances

"""
    bfs(s::Stinger, source::Int64)

This version is slower as it allocates the maximum possible vertices in the Stinger
graph. If you know the maximum number of active vertices, call `bfs(s, source, nv)`
which is faster.
"""
function bfs(s::Stinger, source::Int64)
    bfs(s, source, s[max_nv])
end

"""
    bfs(s::Stinger, source::Int64, nv::Int64)

`bfs` returns a parents array of length `nv`. An empty array is returned on
failure.
"""
function bfs(s::Stinger, source::Int64, nv::Int64)
    if source>=nv
        return zeros(Int64, 0)
    end
    bfskernel(s, source, nv)
end

function bfskernel(s::Stinger, source::Int64, nv::Int64)
    parents = fill(-2, nv) #Initialize parents array with -2's.
    parents[source+1] = -1 #Set source to -1
    next = Vector{Int64}([source])
    sizehint!(next, nv)
    while !isempty(next)
        src = shift!(next) #Get first element
        foralledges(s, src) do edge, src, etype
            direction, neighbor = edgeparse(edge)
            if (direction != 1 && parents[neighbor+1] == -2)
                parents[neighbor+1] = src
                push!(next, neighbor)
            end
        end
    end
    parents
end

"""
    bfsdistances(s::Stinger, source::Int64)

This version is slower as it allocates the maximum possible vertices in the Stinger
graph. If you know the maximum number of active vertices, call `bfsdistances(s, source, nv)`
which is faster.
"""
function bfsdistances(s::Stinger, source::Int64)
    bfsdistances(s, source, s[max_nv])
end

"""
    bfsdistances(s::Stinger, source::Int64, nv::Int64)

Obtain both the `parents` array as well as the `distances` from source.
Empty arrays are returned on failure.
"""
function bfsdistances(s::Stinger, source::Int64, nv::Int64)
    if source>=nv
        return zeros(Int64, 0), zeros(Int64, 0)
    end
    bfsdistanceskernel(s, source, nv)
end

function bfsdistanceskernel(s::Stinger, source::Int64, nv::Int64)
    parents = fill(-2, nv) #Initialize parents array with -2's.
    distances = fill(-1, nv)
    parents[source+1]=-1 #Set source to -1
    distances[source+1]=0
    next = Vector{Int64}([source])
    sizehint!(next, nv)
    while !isempty(next)
        src = shift!(next) #Get first element
        foralledges(s, src) do edge, src, etype
            direction, neighbor = edgeparse(edge)
            #If edge direction is out and the parent is not set yet.
            if (direction != 1 && parents[neighbor+1] == -2)
                parents[neighbor+1] = src
                distances[neighbor+1] = distances[src+1] + 1
                push!(next, neighbor)
            end
        end
    end
    parents, distances
end


text:
function  generaMatCeldasConf(dc)

    M=dc.MATCELDASCONF;
    M=replace!(M, NaN=>-1);
    M=convert(Array{Int64,2},M);

    NumConfig, NumCeldasSelec =size(M);

    numCeldas=maximum(M);

    MATCELDASCONF=zeros(NumConfig,numCeldas);
    MATCONFHOR=zeros(NumConfig,numCeldas);
    MATCONFVERT=zeros(NumConfig,numCeldas);
    for i=1:NumConfig
        for j=1:NumCeldasSelec
            celdaSelec_i=dc.MATCELDASCONF[i,j];
            celdaSelec_i=floor(Int, celdaSelec_i);
            if celdaSelec_i>=1
                MATCELDASCONF[i,celdaSelec_i]=1;
                if dc.MATCONFHOR[i,j]==1
                    MATCONFHOR[i,celdaSelec_i]=1;
                else
                    MATCONFVERT[i,celdaSelec_i]=1;
                end
            end
        end
    end

    MATCELDASCONF, MATCONFHOR, MATCONFVERT

end


text:
@safetestset "ishash" begin
    using Validator

    # length = 32
    algorithms = ["md5", "md4", "ripemd128", "tiger128"]
    for algorithm in algorithms
        @test Validator.ishash("d94f3f016ae679c3008de268209132f2", algorithm) === true
        @test Validator.ishash("751adbc511ccbe8edf23d486fa4581cd", algorithm) === true
        @test Validator.ishash("88dae00e614d8f24cfd5a8b3f8002e93", algorithm) === true
        @test Validator.ishash("0bf1c35032a71a14c2f719e5a14c1e96", algorithm) === true
        @test Validator.ishash("d94f3F016Ae679C3008de268209132F2", algorithm) === true
        @test Validator.ishash("88DAE00e614d8f24cfd5a8b3f8002E93", algorithm) === true

        @test Validator.ishash("q94375dj93458w34", algorithm) === false
        @test Validator.ishash("39485729348", algorithm) === false
        @test Validator.ishash("%&FHKJFvk", algorithm) === false
        @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a1", algorithm) === false
        @test Validator.ishash("", algorithm) === false
        @test Validator.ishash("d94f3f016ae679c3008de268209132f&", algorithm) === false
    end

    algorithms = ["crc32", "crc32b"]
    for algorithm in algorithms
        @test Validator.ishash("d94f3f01", algorithm) === true
        @test Validator.ishash("751adbc5", algorithm) === true
        @test Validator.ishash("88dae00e", algorithm) === true
        @test Validator.ishash("0bf1c350", algorithm) === true
        @test Validator.ishash("88DAE00e", algorithm) === true
        @test Validator.ishash("751aDBc5", algorithm) === true

        @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", algorithm) === false
        @test Validator.ishash("q94375dj93458w34", algorithm) === false
        @test Validator.ishash("q943", algorithm) === false
        @test Validator.ishash("39485729348", algorithm) === false
        @test Validator.ishash("%&FHKJFvk", algorithm) === false
    end

    algorithms = ["sha1", "tiger160", "ripemd160"]
    for algorithm in algorithms
        @test Validator.ishash("3ca25ae354e192b26879f651a51d92aa8a34d8d3", algorithm) === true
        @test Validator.ishash("aaf4c61ddcc5e8a2dabede0f3b482cd9aea9434d", algorithm) === true
        @test Validator.ishash("beb8c3f30da46be179b8df5f5ecb5e4b10508230", algorithm) === true
        @test Validator.ishash("efd5d3b190e893ed317f38da2420d63b7ae0d5ed", algorithm) === true
        @test Validator.ishash("AAF4c61ddCC5e8a2dabede0f3b482cd9AEA9434D", algorithm) === true
        @test Validator.ishash("3ca25AE354e192b26879f651A51d92aa8a34d8D3", algorithm) === true

        @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", algorithm) === false
        @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1dsjkjkjkjkkjk", algorithm) === false
        @test Validator.ishash("q94375dj93458w34", algorithm) === false
        @test Validator.ishash("39485729348", algorithm) === false
        @test Validator.ishash("%&FHKJFvk", algorithm) === false
    end

    @test Validator.ishash("2cf24dba5fb0a30e26e83b2ac5b9e29e1b161e5c1fa7425e73043362938b9824", "sha256") === true
    @test Validator.ishash("1d996e033d612d9af2b44b70061ee0e868bfd14c2dd90b129e1edeb7953e7985", "sha256") === true
    @test Validator.ishash("80f70bfeaed5886e33536bcfa8c05c60afef5a0e48f699a7912d5e399cdcc441", "sha256") === true
    @test Validator.ishash("579282cfb65ca1f109b78536effaf621b853c9f7079664a3fbe2b519f435898c", "sha256") === true
    @test Validator.ishash("2CF24dba5FB0a30e26E83b2AC5b9E29E1b161e5C1fa7425E73043362938b9824", "sha256") === true
    @test Validator.ishash("80F70bFEAed5886e33536bcfa8c05c60aFEF5a0e48f699a7912d5e399cdCC441", "sha256") === true

    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", "sha256") === false
    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1dsjkjkjkjkkjk", "sha256") === false
    @test Validator.ishash("q94375dj93458w34", "sha256") === false
    @test Validator.ishash("39485729348", "sha256") === false
    @test Validator.ishash("%&FHKJFvk", "sha256") === false

    @test Validator.ishash("3fed1f814d28dc5d63e313f8a601ecc4836d1662a19365cbdcf6870f6b56388850b58043f7ebf2418abb8f39c3a42e31", "sha384") === true
    @test Validator.ishash("b330f4e575db6e73500bd3b805db1a84b5a034e5d21f0041d91eec85af1dfcb13e40bb1c4d36a72487e048ac6af74b58", "sha384") === true
    @test Validator.ishash("bf547c3fc5841a377eb1519c2890344dbab15c40ae4150b4b34443d2212e5b04aa9d58865bf03d8ae27840fef430b891", "sha384") === true
    @test Validator.ishash("fc09a3d11368386530f985dacddd026ae1e44e0e297c805c3429d50744e6237eb4417c20ffca8807b071823af13a3f65", "sha384") === true
    @test Validator.ishash("3fed1f814d28dc5d63e313f8A601ecc4836d1662a19365CBDCf6870f6b56388850b58043f7ebf2418abb8f39c3a42e31", "sha384") === true
    @test Validator.ishash("b330f4E575db6e73500bd3b805db1a84b5a034e5d21f0041d91EEC85af1dfcb13e40bb1c4d36a72487e048ac6af74b58", "sha384") === true

    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", "sha384") === false
    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1dsjkjkjkjkkjk", "sha384") === false
    @test Validator.ishash("q94375dj93458w34", "sha384") === false
    @test Validator.ishash("39485729348", "sha384") === false
    @test Validator.ishash("%&FHKJFvk", "sha384") === false

    @test Validator.ishash("9b71d224bd62f3785d96d46ad3ea3d73319bfbc2890caadae2dff72519673ca72323c3d99ba5c11d7c7acc6e14b8c5da0c4663475c2e5c3adef46f73bcdec043", "sha512") === true
    @test Validator.ishash("83c586381bf5ba94c8d9ba8b6b92beb0997d76c257708742a6c26d1b7cbb9269af92d527419d5b8475f2bb6686d2f92a6649b7f174c1d8306eb335e585ab5049", "sha512") === true
    @test Validator.ishash("45bc5fa8cb45ee408c04b6269e9f1e1c17090c5ce26ffeeda2af097735b29953ce547e40ff3ad0d120e5361cc5f9cee35ea91ecd4077f3f589b4d439168f91b9", "sha512") === true
    @test Validator.ishash("432ac3d29e4f18c7f604f7c3c96369a6c5c61fc09bf77880548239baffd61636d42ed374f41c261e424d20d98e320e812a6d52865be059745fdb2cb20acff0ab", "sha512") === true
    @test Validator.ishash("9B71D224bd62f3785D96d46ad3ea3d73319bFBC2890CAAdae2dff72519673CA72323C3d99ba5c11d7c7ACC6e14b8c5DA0c4663475c2E5c3adef46f73bcDEC043", "sha512") === true
    @test Validator.ishash("432AC3d29E4f18c7F604f7c3c96369A6C5c61fC09Bf77880548239baffd61636d42ed374f41c261e424d20d98e320e812a6d52865be059745fdb2cb20acff0ab", "sha512") === true

    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", "sha512") === false
    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1dsjkjkjkjkkjk", "sha512") === false
    @test Validator.ishash("q94375dj93458w34", "sha512") === false
    @test Validator.ishash("39485729348", "sha512") === false
    @test Validator.ishash("%&FHKJFvk", "sha512") === false

    @test Validator.ishash("6281a1f098c5e7290927ed09150d43ff3990a0fe1a48267c", "tiger192") === true
    @test Validator.ishash("56268f7bc269cf1bc83d3ce42e07a85632394737918f4760", "tiger192") === true
    @test Validator.ishash("46fc0125a148788a3ac1d649566fc04eb84a746f1a6e4fa7", "tiger192") === true
    @test Validator.ishash("7731ea1621ae99ea3197b94583d034fdbaa4dce31a67404a", "tiger192") === true
    @test Validator.ishash("6281A1f098c5e7290927ed09150d43ff3990a0fe1a48267C", "tiger192") === true
    @test Validator.ishash("46FC0125a148788a3AC1d649566fc04eb84A746f1a6E4fa7", "tiger192") === true

    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1", "tiger192") === false
    @test Validator.ishash("KYT0bf1c35032a71a14c2f719e5a14c1dsjkjkjkjkkjk", "tiger192") === false
    @test Validator.ishash("q94375dj93458w34", "tiger192") === false
    @test Validator.ishash("39485729348", "tiger192") === false
    @test Validator.ishash("%&FHKJFvk", "tiger192") === false
end


text:
using Test
using WaveProp
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh
using LinearAlgebra

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
"""
    PropertyBatchLens(names)

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> lens = PropertyBatchLens(:a, :b, :c);

julia> get((a=1, b=2, c=3, d=4), lens)
(a = 1, b = 2, c = 3)

julia> set((a=1, b=2, c=3, d=4), lens, (a=10, b=20, c=30))
(a = 10, b = 20, c = 30, d = 4)
```
"""
PropertyBatchLens

"""
    KeyBatchLens(names)

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> lens = KeyBatchLens(:a, :b, :c);

julia> get((a=1, b=2, c=3, d=4), lens)
(a = 1, b = 2, c = 3)

julia> set((a=1, b=2, c=3, d=4), lens, Dict(:a=>10, :b=>20, :c=>30))
(a = 10, b = 20, c = 30, d = 4)
```
"""
KeyBatchLens

"""
    IndexBatchLens(names)

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> lens = IndexBatchLens(:a, :b, :c);

julia> get((a=1, b=2, c=3, d=4), lens)
(1, 2, 3)

julia> set((a=1, b=2, c=3, d=4), lens, (10, 20, 30))
(a = 10, b = 20, c = 30, d = 4)
```
"""
IndexBatchLens

@enum Accessor INDEX KEY PROPERTY

struct BatchLens{names, objacc, valueacc} <: Lens end
const PropertyBatchLens{names} = BatchLens{names, PROPERTY, PROPERTY}
const KeyBatchLens{names} = BatchLens{names, PROPERTY, KEY}
const IndexBatchLens{names} = BatchLens{names, PROPERTY, INDEX}

lenstypenameof(::BatchLens{<:Any, PROPERTY, valueacc}) where valueacc =
    Dict(
        PROPERTY => :PropertyBatchLens,
        KEY => :KeyBatchLens,
        INDEX => :IndexBatchLens,
    )[valueacc]

setter(objacc::Accessor) = setter(Val(objacc))
setter(::Val{INDEX}) = setindices
setter(::Val{KEY}) = setkeys
setter(::Val{PROPERTY}) = setproperties

Setfield.set(obj, ::BatchLens{names, objacc, INDEX}, val) where {names, objacc} =
    setter(objacc)(obj, NamedTuple{names}(Tuple(val)))

Setfield.set(obj, ::BatchLens{names, objacc, KEY}, val) where {names, objacc} =
    setter(objacc)(obj, NamedTuple{names}(map(n -> val[n], names)))

Setfield.set(obj, ::BatchLens{names, objacc, PROPERTY}, val) where {names, objacc} =
    setter(objacc)(obj, NamedTuple{names}(map(n -> getproperty(val, n), names)))

_get(::Val{INDEX}, obj, (i, name)) = getindex(obj, i)
_get(::Val{KEY}, obj, (i, name)) = getindex(obj, name)
_get(::Val{PROPERTY}, obj, (i, name)) = getproperty(obj, name)
_get(acc::Accessor, obj) = x -> _get(Val(acc), obj, x)

getastuple(obj, ::BatchLens{names, objacc}) where {names, objacc} =
    _map(_get(objacc, obj), _enumerate(names))

Setfield.get(obj, lens::BatchLens{<:Any, <:Any, INDEX}) = getastuple(obj, lens)
Setfield.get(obj, lens::Union{BatchLens{names, <:Any, KEY},
                              BatchLens{names, <:Any, PROPERTY}}) where names =
    NamedTuple{names}(getastuple(obj, lens))

BatchLens{<:Any, objacc, valueacc}(names::Vararg{Symbol}) where {objacc, valueacc} =
    BatchLens{names, objacc, valueacc}()

Base.show(io::IO, lens::BatchLens{names, PROPERTY}) where names =
    print_apply(io, Prefixed(prefixof(BatchLens), lenstypenameof(lens)), names)


text:
function compare(a::MyMessage, b::MyMessage; tol::Float64=1e-14)
    a.field1 == b.field1 && isapprox(a.field2, b.field2, atol=tol)
end

function handle_file(lcmlog::LCMLog; N=1)
    for i in 1 : N
        handle(lcmlog) || break
    end
    nothing
end

function make_test_lcmlog()
    # record a temporary log file
    lcmlogdir = joinpath(@__DIR__, "testdata", "testlog.lcm")

    # recreate the messages locally for comparison with those in the test log file
    msg1 = MyMessage(23, 1.234)
    msg2 = MyMessage(24, 2.345)

    lcmlog = LCMLog(lcmlogdir)
    subscribe(lcmlog, "CHANNEL_1", (channel, data) -> @test(compare(decode(data, MyMessage), msg1)))
    subscribe(lcmlog, "CHANNEL_2", (channel, msg) -> @test(compare(msg, msg2)), MyMessage)
    lcmlog
end

@testset "LCM log overrun" begin
    lcmlog = make_test_lcmlog()
    handle_file(lcmlog, N=100)
    close(lcmlog)
end

@testset "Read until end" begin
    lcmlog = make_test_lcmlog()
    while handle(lcmlog); end
    close(lcmlog)
end

# Test for issue #56
include(joinpath(dirname(@__FILE__), "lcmtypes", "image_metadata_t.jl"))
@testset "Encode/decode issue #56" begin
    function handleData(channel, msg::image_metadata_t, msgs)
        @show msg
        push!(msgs, msg)
    end

    msgs = []
    lcmlogdir = joinpath(dirname(@__FILE__),"testdata","image_metadata_log.lcm")
    lc = LCMLog(lcmlogdir)
    LCMCore.subscribe(lc, "CHANNEL_1", (c, d) -> handleData(c, d, msgs), image_metadata_t)

    # Run while there is data
    while handle(lc)
    end
    close(lc)
    # Assert - should be 2 messages
    @test length(msgs) == 2
end

## Code used to create the testlog.lcm LCM log file used in this test
# lcm = LCM()
# msg1 = MyMessage(23, 1.234)
# msg2 = MyMessage(24, 2.345)
# # run `lcm-logger testlog.lcm`
# publish(lcm, "CHANNEL_1", msg1)
# publish(lcm, "CHANNEL_2", msg2)
# # terminate lcm-logger
# close(lcm)


text:
module PickAxisArrays
export PickAxisArray

struct PickAxisArray{T,N,AT<:AbstractArray,P,PERM} <: AbstractArray{T,N}
    parent::AT
end

struct Window
    i::Int
    pre::Int
    after::Int
end

function PickAxisArray(parent, indmask, perm = nothing)
    f = findall(isequal(true), indmask)
    f2 = findall(isequal(Colon()), indmask)
    f3 = findall(i -> isa(i, Tuple{Int,Int}), indmask)
    o = sort([f; f2; f3])
    o = isempty(f2) ? o : replace(o, map(i -> i => Colon(), f2)...)
    o = isempty(f3) ? o : replace(o, map(i -> i => Window(i, indmask[i]...), f3)...)
    nsub = 0
    for i = 1:length(o)
        if o[i] isa Colon
            nsub += 1
        elseif o[i] isa Window
            o[i] = Window(o[i].i - nsub, o[i].pre, o[i].after)
        else
            o[i] = o[i] - nsub
        end
    end
    if perm !== nothing
        length(perm) != length(f2) + length(f3) && error("Not a valid permutation")
        perm = (perm...,)
    end
    im = (o...,)
    T = anycol(im) ? Array{eltype(parent),ncol(im)} : eltype(parent)
    nd = length(f)
    PickAxisArray{T,nd,typeof(parent),im,perm}(parent)
end
indmask(::PickAxisArray{<:Any,<:Any,<:Any,i}) where {i} = i
getind(i, j) = i[j]
getind(i, j::Colon) = j
function getind(i, j::Window)
    c = i[j.i]
    c-j.pre:c+j.after
end
function Base.size(p::PickAxisArray)
    im = indmask(p)
    ii2 = filter(i->!isa(i,Union{Colon,Window}),im)
    permsize(p,map(i->size(p.parent,i),ii2))
end

permsize(::PickAxisArray{<:Any,<:Any,<:Any,<:Any,P}, x) where {P} = x[P]
permsize(::PickAxisArray{<:Any,<:Any,<:Any,<:Any,nothing}, x) = x
permout(::PickAxisArray{<:Any,<:Any,<:Any,<:Any,P}, x) where {P} = permutedims(x, P)
permout(::PickAxisArray{<:Any,<:Any,<:Any,<:Any,nothing}, x) = x
function Base.view(p::PickAxisArray, i::Integer...)
    inew = map(j -> getind(i, j), indmask(p))
    r = permout(p, view(p.parent, inew...))
    r
end
function Base.getindex(p::PickAxisArray, i::Integer...)
    inew = map(j -> getind(i, j), indmask(p))
    permout(p, getindex(p.parent, inew...))
end
anycol(::Tuple{}) = false
anycol(t::Tuple) = anycol(first(t), Base.tail(t))
anycol(::Union{Colon,Window}, t::Tuple) = true
anycol(i, ::Tuple{}) = false
anycol(::Union{Colon,Window}, ::Tuple{}) = true
anycol(i, t::Tuple) = anycol(first(t), Base.tail(t))

ncol(t::Tuple) = ncol(first(t), Base.tail(t), 0)
ncol(::Union{Colon,Window}, t::Tuple, n) = ncol(first(t), Base.tail(t), n + 1)
ncol(i, ::Tuple{}, n) = n
ncol(::Union{Colon,Window}, ::Tuple{}, n) = n + 1
ncol(i, t::Tuple, n) = ncol(first(t), Base.tail(t), n)

Base.getindex(p::PickAxisArray, i::CartesianIndex) = p[i.I...]
Base.view(p::PickAxisArray, i::CartesianIndex) = view(p,i.I...)
end


text:
import PlutoSliderServer
import PlutoSliderServer.MoreAnalysis
import PlutoSliderServer.Pluto

using Test

file = joinpath(@__DIR__, "parallelpaths4.jl")

newpath = tempname()
Pluto.readwrite(file, newpath)

notebook = Pluto.load_notebook(newpath)
# Run pluto's analysis. This is like opening the notebook, without actually running it
# s = Pluto.ServerSession()
# Pluto.update_save_run!(s, notebook, notebook.cells; run)
notebook.topology = Pluto.updated_topology(notebook.topology, notebook, notebook.cells)

# bound_variables = (map(notebook.cells) do cell
#     MoreAnalysis.find_bound_variables(cell.parsedcode)
# end)

# @show bound_variables

connections = MoreAnalysis.bound_variable_connections_graph(notebook)
# @show connections

@test !isempty(connections)
wanted_connections = Dict(
    :x         => [:y, :x],
    :y         => [:y, :x],
    :show_dogs => [:show_dogs],
    :b         => [:b],
    :c         => [:c],
    :five1     => [:five1],
    :five2     => [:five2],
    :six1      => [:six2, :six1],
    :six2      => [:six3, :six2, :six1],
    :six3      => [:six3, :six2],
    :cool1     => [:cool1, :cool2],
    :cool2     => [:cool1, :cool2],
    :world     => [:world],
    :boring    => [:boring],
)

transform(d) = Dict(k => sort(v) for (k, v) in d)

@test transform(connections) == transform(wanted_connections)


text:
# Writer the byte b in w.
function write(w::WritableFile, b::UInt8)
    write(w, Ref(b))
end

# Read and return a byte from f. Throws EOFError if there is no more byte to read.
function read(f::ReadableFile, ::Type{UInt8})
    # This function needs to be fast because readbytes, readstring, etc.
    # uses it. Avoid function calls when possible.
    b = Vector{UInt8}(undef, 1)
    c = read(f, b)
    c[1]
end


text:
# test SE3 kernel and MMD

using ApproxManifoldProducts
using Manifolds
using Test

##

@testset "Test isapprox function on basic SpecialEuclidean(3)" begin

##

# FIXME hacking the manifolds here, needs consolidation
# a = kde!(randn(6,50))
# b = kde!(randn(6,50))
# c_ = randn(6,50)
# c_[1,:] .+= 50
# c = kde!(c_)

M = SpecialEuclidean(3)
u0 = ProductRepr(zeros(3),[1 0 0; 0 1 0; 0 0 1.0])
ϵ = identity_element(M, u0)
N = 50

pts1 = [exp(M, ϵ, hat(M, ϵ, [0.5*randn(3);0.1*randn(3)])) for i in 1:N]
pts2 = [exp(M, ϵ, hat(M, ϵ, [0.5*randn(3);0.1*randn(3)])) for i in 1:N]
pts3 = [exp(M, ϵ, hat(M, ϵ, [0.5*randn()+50;0.5*randn(2);0.1*randn(3)])) for i in 1:N]

##

ret = mmd(M, pts1, pts2)
@test ret < 1

ret = mmd(M, pts1, pts3)
@test 0.001 < ret

##

A = ManifoldKernelDensity(M, pts1)
B = ManifoldKernelDensity(M, pts2)
C = ManifoldKernelDensity(M, pts3)

@test isapprox(A, B)
@test !isapprox(A, C)

##

show(A)

##

end

#


text:
using NetworkFlows
using Base.Test

## Test link.jl ##
#import NetworkFlows.zero_to_one!
@test zero_to_one!([(0,1),(1,2)]) == [(1,2),(2,3)]
@test zero_to_one!([(3,1),(1,2)]) != [(3,1),(1,2)]
@test zero_to_one!([(0,1,1.),(1,2,2.)]) == [(1,2,1.),(2,3,2.)]
@test zero_to_one!([(3,1,0.),(1,2,-2.1)]) != [(3,1,0.),(1,2,-2.1)]

import NetworkFlows.Arc, NetworkFlows.ghostArc, NetworkFlows.simpleArc
@test Arc(0,1,0.).sym == ghostArc(1).sym
@test Arc(0,1,0.).head == ghostArc(1).head
@test Arc(0,1,0.).cap == ghostArc(1).cap
@test Arc(0,2,1.).sym == simpleArc(2).sym
@test Arc(0,2,1.).head == simpleArc(2).head
@test Arc(0,2,1.).cap == simpleArc(2).cap

## Test network.jl
import NetworkFlows.Network, NetworkFlows.findArc
g1 = Network(1,2,[1,2,2],[Arc(0,2,1.)])
g2 = Network(2,1,[1,2,2],[simpleArc(2)])
@test g1.source == g2.target
@test g2.source == g1.target
@test g1.tails == g2.tails
@test length(g1.links) == length(g2.links)
@test findArc(g1,1,2) == findArc(g2,1,2)

## Test io.jl
#import NetworkFlows.print
edges1 = [(0,1,5.),(0,2,7.),(0,5,1.75),(0,6,1.3),(1,3,2.),(2,1,1.),
  (2,3,11.),(2,4,8.),(3,5,6.),(4,5,3.),(6,5,1.25)]
zero_to_one!(edges1)
g3 = Network(edges1, true, 1, 6)
print(g3)

import NetworkFlows.orderedDIMACS
orderedDIMACS(pwd()*"/../ioexample/unordereddimacs",
  pwd()*"/../ioexample/ordereddimacs")
gdimacs = Network(pwd()*"/../ioexample/ordereddimacs",:DIMACS)
gcsv = Network(pwd()*"/../ioexample/csv",:CSV,(1,6))

## Test search.jl
import NetworkFlows.bfs
println("BFS: ",bfs(g3,:Path))

## Test flow.jl
# import NetworkFlows.edmondsKarp, NetworkFlows.connectivity
# import NetworkFlows.kishimoto
@test edmondsKarp(g3)[1] == 12
@test connectivity(g3) == 4
@test edmondsKarp(gdimacs)[1] == 15
@test edmondsKarp(gcsv)[1] == 15
edges2 =
[(1,2,10.),(1,3,20.),(1,4,40.),(1,5,80.),(2,6,1.),(3,6,2.),(4,6,4.),(5,6,8.)]
g4 = Network(edges2,true,1,6)
@test kishimoto(g4,1)[1] == edmondsKarp(g4)[1]
@test kishimoto(g4,2)[1] == 14.
@test kishimoto(g4,3)[1] == 9.
@test kishimoto(g4,4)[1] == 4.

## Test cut.jl
#import NetworkFlows.mincut
@test mincut(g3)[1] == edmondsKarp(g3)[1]
@test mincut(g3)[2] != edmondsKarp(g3)[2]

# Test extmrflow
# import NetworkFlows.breakingPoints
@test breakingPoints(g4) ==
  [(0.,0.,4),(1.0,4.0,3),(2.0,7.0,2),(4.0,11.0,1),(8.0,15.0,0)]
edges3 =
[(1,2,10.),(1,3,20.),(1,4,40.),(1,5,80.),(2,6,2.),(3,6,2.),(4,6,2.),(5,6,8.)]
g5 = Network(edges3,true,1,6)
@test breakingPoints(g5) ==
  [(0.0,0.0,4),(-1.0,-1.0,3),(-1.0,-1.0,2),(2.0,8.0,1),(4.0,14.0,0)]

  # Test mla.jl
  # import NetworkFlows.mixedMLA, NetworkFlows.successMLA
  edges4 = [(1,2,2.),(1,3,3.),(1,4,5.),(2,5,2.),(3,5,3.),(4,5,5.),
    (5,6,3.),(5,7,3.),(5,8,3.),(6,9,3.),(7,9,3.),(8,9,3.),
    (9,10,7.),(9,11,7.),(10,12,7.),(11,12,7.)]
  g6 = Network(edges4,true,1,12)
  @test mixedMLA(g6) == [(4.5,6.0),(9.0,9.0)]
  @test successMLA(g6) == 0.5


text:
#GracePlot: Lightweight units
#-------------------------------------------------------------------------------

#==TODO:
use proper units/quantity module, as described in:
https://github.com/ma-laforge/testcases/blob/master/units_test/units_test.jl
==#

#==Constants
===============================================================================#
const INCH_IN_METERS = .0254
const INCH_IN_POINTS = 72

#==Type definitions
===============================================================================#

abstract AbstractQuantity
val(x::AbstractQuantity) = x.v

abstract AbstractLength <: AbstractQuantity

#Typographic (DTP) "point":
immutable TPoint <: AbstractLength
	v::Number
end

immutable Inch <: AbstractLength
	v::Number
end

immutable Meter <: AbstractLength
	v::Number
end

#==Conversion functions:
===============================================================================#
Base.convert(::Type{Inch}, x::Meter) = Inch(val(x)*(1/INCH_IN_METERS))
Base.convert(::Type{Inch}, x::TPoint) = Inch(val(x)*(1/INCH_IN_POINTS))
Base.convert(::Type{Meter}, x::Inch) = Meter(val(x)*INCH_IN_METERS)
Base.convert(::Type{TPoint}, x::Inch) = TPoint(val(x)*INCH_IN_POINTS)

#Indirect conversions:
Base.convert(::Type{Meter}, x::TPoint) = convert(Meter, convert(Inch, x))
Base.convert(::Type{TPoint}, x::Meter) = convert(TPoint, convert(Inch, x))

TPoint(x::TPoint) = x
Meter(x::Meter) = x
Inch(x::Inch) = x
TPoint(x) = convert(TPoint, x)
Meter(x) = convert(Meter, x)
Inch(x) = convert(Inch, x)

#end


text:
using Test

@testset "Functions" begin
    include("functions.jl")
end
@testset "Sets" begin
    include("sets.jl")
end
@testset "Constraints" begin
    include("constraints.jl")
end
@testset "Model" begin
    include("model.jl")
end
@testset "Universal Fallback" begin
    include("universalfallback.jl")
end
@testset "Parser" begin
    include("parser.jl")
end
@testset "Mock Optimizer" begin
    include("mockoptimizer.jl")
end
@testset "Caching Optimizer" begin
    include("cachingoptimizer.jl")
end
@testset "Copy" begin
    include("copy.jl")
end


text:
include("posteriorquantileinterval.jl")

data = rand(1,100)

calcquantileint(data, 5, 95)

data = rand(100,4)

calcquantileint(data)

data = rand(100)

calcquantileint(data)

max_stars_repo_name: SamuelWiqvist/abcmcmcrickermodel

max_issues_repo_name: SamuelWiqvist/abcmcmcrickermodel

max_forks_repo_name: SamuelWiqvist/abcmcmcrickermodel


text:
export networkModel

type networkModel
  dt::Float64
  groups::Dict
  parameters::AbstractString
  equations::AbstractString
  initialization::AbstractString
  connection::AbstractString
  spike::AbstractString
  reset::AbstractString
  synapse::AbstractString
  plasticity::AbstractString
  record::Array
  nCores::Int64
  solver::AbstractString
  synctime::Float64
  saveVars::Array
  postProcs::AbstractString
  networkModel()=new(0.1,Dict(),"","","","","","","","",[],1,"Euler1",0.0,[],"")
end


text:
abstract type CountingBasedProbabilityEstimator <: ProbabilitiesEstimator end
export CountOccurrences
import DelayEmbeddings: AbstractDataset

"""
    CountOccurrences  <: CountingBasedProbabilityEstimator

A probabilities/entropy estimator based on straight-forward counting of distinct elements in
a univariate time series or multivariate dataset. From these counts, construct histograms.
Sum-normalize histograms to obtain probability distributions.
"""
struct CountOccurrences <: ProbabilitiesEstimator end

function probabilities(x::AbstractDataset, est::CountOccurrences)
    probabilities(x)
end


text:
# Kmer Iterator
# =============
#
# Iterator over all k-mers in a biological sequence.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Iterate through every k-mer in a nucleotide sequence
struct EachKmerIterator{T<:Kmer,S<:Sequence}
    seq::S
    step::Int
    start::Int
end

"""
    each(::Type{Kmer{T,k}}, seq::Sequence[, step=1])

Initialize an iterator over all k-mers in a sequence `seq` skipping ambiguous
nucleotides without changing the reading frame.

# Arguments
* `Kmer{T,k}`: k-mer type to enumerate.
* `seq`: a nucleotide sequence.
* `step=1`: the number of positions between iterated k-mers

# Examples
```
# iterate over DNA codons
for (pos, codon) in each(DNAKmer{3}, dna"ATCCTANAGNTACT", 3)
    @show pos, codon
end
```
"""
function each(::Type{Kmer{T,K}}, seq::Sequence, step::Integer=1) where {T,K}
    if eltype(seq) ∉ (DNA, RNA)
        throw(ArgumentError("element type must be either DNA or RNA nucleotide"))
    elseif !(0 ≤ K ≤ 32)
        throw(ArgumentError("k-mer length must be between 0 and 32"))
    elseif step < 1
        throw(ArgumentError("step size must be positive"))
    end
    return EachKmerIterator{Kmer{T,K},typeof(seq)}(seq, step, 1)
end

eachkmer(seq::BioSequence{A}, K::Integer, step::Integer=1) where {A<:DNAAlphabet} = each(DNAKmer{Int(K)}, seq, step)
eachkmer(seq::BioSequence{A}, K::Integer, step::Integer=1) where {A<:RNAAlphabet} = each(RNAKmer{Int(K)}, seq, step)
eachkmer(seq::ReferenceSequence, K::Integer, step::Integer=1) = each(DNAKmer{Int(K)}, seq, step)

Base.eltype(::Type{EachKmerIterator{T,S}}) where {T,S} = Tuple{Int,T}

function Base.IteratorSize(::Type{EachKmerIterator{T,S}}) where {T,S}
    return Base.SizeUnknown()
end

Base.iterate(it::EachKmerIterator) = iterate(it, (it.start, UInt64(0)))

function Base.iterate(
        it::EachKmerIterator{T},
        state::Tuple{Int, UInt64}) where {T}

    k = kmersize(T)
    pos, kmer = state
    isok = true

    # faster path: return the next overlapping kmer if possible
    if it.step < k && pos + k - 1 ≤ lastindex(it.seq) && pos != it.start
        offset = k - it.step
        if it.step == 1
            nt = inbounds_getindex(it.seq, pos+offset)
            kmer = kmer << 2 | trailing_zeros(nt)
            isok &= iscertain(nt)
        else
            for i in 1:it.step
                nt = inbounds_getindex(it.seq, pos+i-1+offset)
                kmer = kmer << 2 | trailing_zeros(nt)
                isok &= iscertain(nt)
            end
        end
        if isok
            return (pos, T(kmer)), (pos+it.step, kmer)
        end
    end

    while pos + k - 1 ≤ lastindex(it.seq)
        kmer, ok = extract_kmer_impl(it.seq, pos, k)
        if ok
            return (pos, T(kmer)), (pos+it.step, kmer)
        end
        pos += it.step
    end

    return nothing
end

function extract_kmer_impl(seq, from, k)
    kmer::UInt64 = 0
    isok = true
    for i in 1:k
        nt = inbounds_getindex(seq, from+i-1)
        kmer = kmer << 2 | trailing_zeros(nt)
        isok &= iscertain(nt)
    end
    return kmer, isok
end


text:
## Point in Polyhedron Algorithm
#  By Lane, Magedson and Rarick

using SparseArrays;
using LinearAlgebraicRepresentation, ViewerGL;
Lar = LinearAlgebraicRepresentation;
GL =  ViewerGL;
using AlphaStructures;

"""
	pointInPolyhedron3D(
		p::Array{Float64,1},
		V::Lar.Points,
		EV::Lar.Cells,
		FV::Lar.Cells
	)::Int64
Evaluate if `p` belongs to the complex defined by `[V, EV, FV]`.

This method computes if a given point `p` is inner, outer or edger with respect
to the complex defined by the sets `V, EV, FV`. Respectivelly it return
the integer value `+1`, `-1` or `0`.

_Note._ If something weired happens then `-3` is returned and the value for the
solid angle related to `p` is shown.
"""
function pointInPolyhedron3D(
		p::Array{Float64,1},
		V::Lar.Points,
		EV::Lar.Cells,
		FV::Lar.Cells
	)::Int64

	# Faces triangulation is performed (if needed)
	if max([length(σ) for σ in FV]...) > 3
		cop_EV = Lar.coboundary_0(EV::Lar.Cells);
		cop_EW = convert(Lar.ChainOp, cop_EV);
		cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
    	triFV = Lar.triangulate(convert(Lar.Points, V'), [cop_EV, cop_FE]);
	else
		triFV = [[σ] for σ in FV]
	end

	# Solid Angle Initialization
	Ω = 0.0;

	# For the triangulation of each face ...
    for face_complex in triFV
		# ... compute the axis face
        axis = Lar.cross(
			V[:, face_complex[1][2]] .- V[:, face_complex[1][1]],
			V[:, face_complex[1][3]] .- V[:, face_complex[1][1]]
		);
		# If the point lies on a face plan ...
		if Lar.dot(axis, V[:, face_complex[1][1]]) == Lar.dot(axis, p)
			# ... If for some face ...
			face_ref = ∪(face_complex...);
			for Σ in face_complex
				points = projectTo2D([Σ p], Σ);
				position = pointInPolyhedron2D(
					points[:, 1],   points[:, [2 3 4]],   [[1 2], [2 3], [3 1]]
				);
				# ... the point is inner with respect to that face
				#     than it is inner to the polyhedron;
				if position == 1
					return (0, face_ref);
				# ... if it is on an edge of that polyhedron
				#     than it is either inner to that face or to that edge;
				elseif position ∋ 0;
					if position[2] ∈ EV
						return position;
					else
						return (0, face_ref);
					end
				end
			end
			# ... however, if no faces contains it, then the contribute of
			#     that face to the solid angle is equal to zero.
		# ... else if the point is not coplanar to a face ...
		else
			# ... the normal to the face is used to determine if `p` is inner
			#     (i.e. it is opposite to the side of the axis) or it is outer
			isInner = Lar.dot(axis, V[:, face_complex[1][1]]) < Lar.dot(axis, p);
        	for Σ in face_complex
				@show isInner (-1)^!isInner
				@show simplexSolidAngle(p, V[:, Σ])/π
	            Ω = Ω + (-1)^isInner * simplexSolidAngle(p, V[:, Σ])/π;
	        end
		end
    end
    if isapprox(abs(Ω), 4.0; atol=1e-1)
        return 1;
    elseif isapprox(abs(Ω), 2.0; atol=1e-1)
        return 0;
    elseif isapprox(abs(Ω), 0.0; atol=1e-1)
        return -1;
    else
        @show Ω;
		return -3;
    end

end

"""
    simplexSolidAngle(v::Int64, σ::Lar.Points)::Float64
Evaluate the Solid Angle for the `v`-th vertex of the simplex `σ`.
`σ` must be a simplex, therefore, since the dimension is three, then `σ` must be
made of four vertices (non collinear).
"""
function simplexSolidAngle(v::Int64, σ::Lar.Points)::Float64
    @assert size(σ) == (3, 4)
		"ERROR: the simplex σ must be made of four 3D points.";
    p = σ[:, v];
    Σ = σ[:, setdiff(1:n, v)]
    simplexSolidAngle(p, Σ);
end

"""
    simplexSolidAngle(p::Array{Float64,1}, Σ::Lar.Points)::Float64
Evaluate the Solid Angle of `p` with respect to `Σ`.
`Σ ∪ p` must be a simplex, therefore, since 3 is the dimension, then `Σ` must be
made of `3` vertices (non collinear).

#ToDo check collinearity
"""
function simplexSolidAngle(p::Array{Float64,1}, Σ::Lar.Points)::Float64
    @assert size(Σ) == (3, 3) "ERROR: Σ must be made of three 3D points."
    @assert length(p) == 3 "ERROR: p must be a 3D point."

	# If two base vertices are equilocated then the solid angle is null.
	if |([Σ[:, i] == Σ[:, j] for i = 1:3 for j=i+1:3]...)
		return 0.;
	end

	# If `p` is equilocated with a Σ point then the solid angle is not defined.
	if |([p == Σ[:, i] for i = 1 : 3]...)
		println("The solid angle for $p with respect to $Σ is not defined.");
		return NaN;
	end

	# If `p` is coplanar with the base point then the angle is
	#  - 0 if `p` is not in the triangle defined by `Σ`
	#  - 2π otherwise
	if areCoplanar(p, Σ[:, 1], Σ[:, 2], Σ[:, 3])
		print("$p is coplanar with $Σ ");
		points = projectTo2D([Σ p], Σ);
		if pointInPolyhedron2D(points[:, 1], points[:, [2 3 4]], [[1 2], [2 3], [3 1]]);
			println("and it is inside of it.");
			return 2π;
		else
			println("but it is outside of it.");
			return 0.;
		end
	end

	# If `p` is collinear with two base point then the solid angle is null
#	if |([areCollinear(p, Σ[:, i], Σ[:, j]) for i = 1 : 3 for j = i+1 : 3]...)
#		println("$p is collinar with two points of $Σ.");
#		return 0;
#	end

	# Evaluate vectors from p to Σ vertices
    A = Σ[:, 1] .- p;
    B = Σ[:, 2] .- p;
    C = Σ[:, 3] .- p;
	# Evaluate magnitudes
	a = Lar.norm(A);
	b = Lar.norm(B);
	c = Lar.norm(C);
	tan05Ω = Lar.dot(A, Lar.cross(B, C))/(
		a*b*c + Lar.dot(A, B)*c + Lar.dot(B, C)*a + Lar.dot(C, A)*b
	);
    return abs(2 * atan(tan05Ω));
end

"""
	areCoplanar(p, a, b, c)::Bool
Determine if `p` belogns to the plane determined by `a`, `b` and `c`.

Namely it evluates the normal of the plane defined by `pa` and `pb` and checks
if it is normal to the `pc` vector.
"""
function areCoplanar(
		p::Array{Float64,1},
		a::Array{Float64,1},
		b::Array{Float64,1},
		c::Array{Float64,1}
	)::Bool
	return Lar.dot(p.-a, Lar.cross(p.-b, p.-c)) == 0;
end

"""
	pointInPolyhedron2D(
		p::Array{Float64,1},
		V::Lar.Points,
		EV::Lar.Cells,
		CHECKS = true
	)::Int64
Evaluate if `p` belongs to the 2D complex defined by `[V, EV]`.

This method computes if a given point `p` is inner, outer or edger with respect
to the complex defined by the sets `V, EV`. Respectivelly it return
the integer value `+1`, `-1` or `0`.
The optional argument `CHECK` could be set to false in order to avoid checks
over the data structures.

_Note._ If something weired happens then `-3` is returned and the value for the
inner angle related to `p` is shown.
"""
function pointInPolyhedron2D(
		p::Array{Float64,1},
		V::Lar.Points,
		EV::Lar.Cells;
		CHECKS=true
	)::Int64

	ω = 0.0;

	if CHECKS
		@assert size(V, 1) == 2
			"ERROR: V expect to have 2 coordinates but $size(V, 2) where found";
		for i = 1 : size(V, 2)
			@assert length([e for e in EV if i ∈ e]) % 2 == 0
				"ERROR: each point must have an even number of incident edges.";
		end
	end

	if |([isequal(V[:, i], p) for i = 1 : size(V, 2)]...)
		return 0;
	end

	Vvec = V .- p;

	for e in EV
		α = rot2Dangle(Vvec[:, e[1]], Vvec[:, e[2]]) / π;
		if isapprox(α, 1.0)
			return 0;
		elseif α < 1
			ω = ω + α;
		else
			ω = ω + α - 2.;
		end
	end

	if     isapprox(ω, 2.; atol=1e-1) return +1;
	elseif isapprox(ω, 0.; atol=1e-1) return -1; end
	@show "Evaluated angle of $(ω * 180) degrees.";
	return -3;
end

"""
	rot2Dangle(
		x::Array{Float64,1}, y::Array{Float64,1}
		degree = false
	)::Float64
Evaluates radiant anticlockwise rotation angle between two points `x` and `y`.

This method evaluates the angle between `x` and `y` with respect to the origin.
The angle is computed in anticlockwise order and it is outputted either
in radiants (by default) or in degree (if specified via optional argument).

# Examples
```juliadoctest
julia> rot2Dangle([1.; -1.], [1., 1]) / π
1.0

julia> rot2Dangle([1.; -1.], [1., 1], degree = true)
90.0

julia> rot2Dangle([1.; 1.], [1., -1])/π
1.5

julia> rot2Dangle([1.; 1.], [1., -1], degree = true)
270.0
```

"""
function rot2Dangle(
		x::Array{Float64,1}, y::Array{Float64,1};
		degree = false
	)::Float64

	Δθ = rot2Dangle(y, degree = degree) - rot2Dangle(x, degree = degree);
	if degree
		return (Δθ + 360) % 360;
	end
	return (Δθ + 2 * π) % (2 * π);
end

"""
	rot2Dangle(x::Array{Float64,1}; degree = false)::Float64
Evaluates radiant anticlockwise rotation angle of `x` by the positive `x`-axis.

This method evaluates the angle between the positive `x`-axis and the point `x`
with respect to the origin.
The angle is computed in anticlockwise order and it is outputted either
in radiants (by default) or in degree (if specified via optional argument).

# Examples
```juliadoctest
julia> rot2Dangle([1.; 0.]) / π
0.0

julia> rot2Dangle([1.; -1]) / π
1.75

julia> rot2Dangle([1.; -1.], degree = true)
315.0
```
"""
function rot2Dangle(x::Array{Float64,1}; degree = false)::Float64

	if isequal(x, [0.;0.])
		println("Cannot evaluate the angle for (0,0)");
		return NaN;
	end

	abs_α = acos(Lar.dot(x, [1.;0.])/(Lar.norm(x) * Lar.norm([1.;0.])));
	is_negative = x[2]<0;
	α = (-1)^is_negative * abs_α + is_negative * 2 * π;
	if degree
		return α / π * 180;
	end
	return α;
end

"""
	projectTo2D(V::Lar.Points, Σ::Lar.Points)::Lar.Points

Projects the points `V` over the plane defined by `Σ`
"""
function projectTo2D(V::Lar.Points, Σ::Lar.Points; CHECKS = true)::Lar.Points

	dim, m = size(V);

	if CHECKS
		@assert dim == 3 "ERROR: points are not 3D.";
		@assert isequal(size(Σ), (3, 3)) "ERROR: Σ not well defined.";
	end

	O = Σ[:, 1];
	xaxis = Σ[:, 2] .- O;
	yaxis = Σ[:, 3] .- O;
	n = Lar.cross(xaxis, yaxis);
	xaxis = xaxis ./ sqrt(sum(xaxis.^2));
	yaxis = Lar.cross(xaxis, n);
	yaxis = yaxis ./ sqrt(sum(yaxis.^2));

	#pr = [[Lar.dot(V[:, i], xaxis); Lar.dot(V[:, i], yaxis)] for i in 1 : 4];
	xcord = [Lar.dot(V[:, i], xaxis) for i in 1 : m];
	ycord = [Lar.dot(V[:, i], yaxis) for i in 1 : m];

	return [xcord ycord]';
end