text:
using Unitful
using UnitfulAssets
using Test
using Decimals
using FixedPointDecimals

using UnitfulAssets: @asset

module TestCash
    using Unitful
    using UnitfulAssets
    using UnitfulAssets: @asset

    @asset Cash TSP TestPataca
end

const localpromotion = Unitful.promotion # needed for the new currency dimension
Unitful.register(TestCash) # needed for new Units
merge!(Unitful.promotion, localpromotion) # only needed with new dimensions

module TestCashSymbol
    using Unitful
    using UnitfulAssets

    @unit ¥ "¥" YuanSign 1.0u"CNY" true
end

Unitful.register(TestCashSymbol) # needed for new Units

# create exchange market from a fixer.io json file
fixer_exchmkt = Dict(
    "2020-11-01" => UnitfulAssets.generate_exchmkt(
        joinpath(@__DIR__, "exchmkt_data", "2020-11-01_fixer.json"),
        Val(:fixer)
    )
)
# create exchange market from a currencylayer json file
currencylayer_exchmkt = Dict(
    "2020-11-25" => UnitfulAssets.generate_exchmkt(
        joinpath(@__DIR__, "exchmkt_data", "2020-11-25_currencylayer.json"),
        Val(:currencylayer)
    )
)
# rates on Nov 27, 2020
exch_mkt_27nov2020 = generate_exchmkt([
    ("EUR","USD") => 1.19536, ("USD","EUR") => 0.836570,
    ("EUR","GBP") => 1.11268, ("GBP","EUR") => 0.898734,
    ("USD","CAD") => 1.29849, ("CAD","USD") => 0.770125,
    ("USD","BRL") => 5.33897, ("BRL","USD") => 0.187302
])

# rates with rational numbers
exch_mkt_from_dict_and_rationals = generate_exchmkt(Dict([
    ("EUR","USD") => 119536//100000, ("USD","EUR") => 836570//1000000
]))

# rates with decimals
exch_mkt_from_dict_and_decimals = generate_exchmkt(Dict([
    ("EUR","USD") => Decimal(1.19536), ("USD","EUR") => Decimal(0.836570)
]))

test_exch_mkt = generate_exchmkt(("EUR","TSP") => 1234567.89)

# rates to test broadcasting
BRLGBP_timeseries = Dict(
    "2011-01-01" => generate_exchmkt(("BRL","GBP") => 0.38585),
    "2012-01-01" => generate_exchmkt(("BRL","GBP") => 0.34587),
    "2013-01-01" => generate_exchmkt(("BRL","GBP") => 0.29998),
    "2014-01-01" => generate_exchmkt(("BRL","GBP") => 0.25562),
    "2015-01-02" => generate_exchmkt(("BRL","GBP") => 0.24153),
    "2016-01-03" => generate_exchmkt(("BRL","GBP") => 0.17093),
    "2017-01-02" => generate_exchmkt(("BRL","GBP") => 0.24888),
    "2018-01-02" => generate_exchmkt(("BRL","GBP") => 0.22569),
    "2019-01-04" => generate_exchmkt(("BRL","GBP") => 0.21082),
    "2020-01-04" => generate_exchmkt(("BRL","GBP") => 0.18784)
)

@testset "Assets" begin
    # conversions
    @test Unitful.Quantity(1,u"EUR") == 1u"EUR"
    @test Unitful.unit(1u"BRL") == u"BRL"
    @test Unitful.unit(1u"TSP") == u"TSP"
    @test UnitfulAssets.is_asset("USD")
    @test UnitfulAssets.is_asset("TSP")
    @test !UnitfulAssets.is_asset("ABC")
    @test typeof(UnitfulAssets.@asset Cash AAA TripleAs) <: Unitful.FreeUnits
end

@testset "Arithmetics" begin
    @test 1u"EUR" + 2u"EUR" === 3u"EUR"
    @test 1.0u"USD" + 2.0u"USD" === 3.0u"USD"
    @test (1//2)u"CAD" + (3//2)u"CAD" === (2//1)u"CAD"
    @test Decimal(1.5)u"BRL" + Decimal(2.5)u"BRL" == Decimal(4.0)u"BRL"
    @test 13u"EUR"/2u"BRL" === 6.5u"EUR/BRL"
    @test 10u"AUD"/2u"AUD" === 5.0
    @test FixedDecimal{Int,2}(7.21)u"BRL" / FixedDecimal{Int,2}(2)u"kg" === FixedDecimal{Int,2}(3.60)u"BRL/kg"
    @test_throws Unitful.DimensionError (1u"BRL" + 1u"EUR")
end

@testset "Exchanges" begin
    @test uconvert(u"BRL", 1u"EUR", fixer_exchmkt["2020-11-01"]) == 6.685598u"BRL"
    @test uconvert(u"kBRL", 2u"MEUR", fixer_exchmkt["2020-11-01"]) == 13371.196u"kBRL"
    @test uconvert(u"EUR", 1u"BRL", fixer_exchmkt["2020-11-01"], mode=-1) == 0.149575251159283u"EUR"
    @test uconvert(u"CAD", 1u"USD", currencylayer_exchmkt["2020-11-25"]) == 1.30045u"CAD"
    @test uconvert(u"BRL", 2u"hUSD", currencylayer_exchmkt["2020-11-25"]) == 1064.8188u"BRL"
    @test uconvert(u"USD", 1u"EUR", exch_mkt_27nov2020) == 1.19536u"USD"
    @test uconvert(u"EUR", 1u"USD", exch_mkt_27nov2020, mode=-1) == 0.8365680631776201u"EUR"
    @test uconvert(u"CAD", 1u"EUR", exch_mkt_27nov2020, mode=2) == 1.5521630063999998u"CAD"
    @test uconvert(u"EUR", 1u"CAD", exch_mkt_27nov2020, mode=-2) == 0.6442622301116068u"EUR"
    @test uconvert(u"BRL", 1000u"CAD", exch_mkt_27nov2020, mode=2) == 4111.674271249999u"BRL"
    @test uconvert(u"BRL", 1000u"CAD", exch_mkt_27nov2020, mode=-2) == 4111.6768608248185u"BRL"
    @test uconvert.(u"BRL", 1000u"GBP", values(BRLGBP_timeseries), mode=-1) ≈ [2591.68, 2891.26, 4018.00, 4743.38, 5850.35, 3333.56, 4140.27, 3912.06, 4430.86, 5323.68]u"BRL" (atol=0.01u"BRL")
    @test uconvert(u"EUR", Decimal(1.5)u"USD", exch_mkt_from_dict_and_decimals) == Decimal(0, 1254855, -6) * u"EUR"
    @test uconvert(u"EUR", Decimal(1.5)u"USD", exch_mkt_from_dict_and_decimals, mode=-1) == Decimal(0, 125485209476643019676, -20)u"EUR"
    @test uconvert(u"EUR", (3//2)u"USD", exch_mkt_from_dict_and_rationals) == (250971//200000)u"EUR"
    @test uconvert(u"EUR", (3//2)u"USD", exch_mkt_from_dict_and_rationals, mode=-1) == (9375//7471)u"EUR"
    @test uconvert(u"TSP", 1u"EUR", test_exch_mkt) == 1234567.89u"TSP"
    @test_throws ArgumentError uconvert(u"EUR", 1u"BRL", fixer_exchmkt["2020-11-01"])
    @test_throws ArgumentError uconvert(u"CAD", 1u"BRL", fixer_exchmkt["2020-11-01"], mode=2)
    @test_throws ArgumentError uconvert(u"EUR", 1u"CAD", exch_mkt_27nov2020)
    @test_throws ArgumentError uconvert(u"m",1u"km", exch_mkt_27nov2020)
    @test_throws ArgumentError generate_exchmkt(("EUR","USD") => Decimal(-1.0))
end

@testset "Symbols" begin
    @test uconvert(u"€", 1u"EUR") == 1u"€"
    @test uconvert(u"USdollar", 1u"USD") == 1u"USdollar"
    @test uconvert(u"CAdollar", 1u"CAD") == 1u"CAdollar"
    @test uconvert(u"BRL", 1u"Real") == 1u"BRL"
    @test uconvert(u"£", 1u"GBP") == 1u"£"
    @test uconvert(u"¥", 1u"CNY") == 1u"¥"
end

nothing


text:
# with z-proxy. refactor/remove this later.
function evalitpproxymixture(u_rad, Es::Vector{zCompoundFIDType{T,SST}};
    w::Vector{T} = ones(T, length(Es)))::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out = zero(Complex{T})

    for n = 1:length(Es)
        out += w[n]*evalitpproxycompound(u_rad, Es[n])
    end

    return out
end

# with z-proxy.
function evalitpproxycompound(u_rad, A::zCompoundFIDType{T,SST})::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out_sys = evalzitpproxysys(A.zs, A.core.gs, u_rad, A.core.ss_params)

    out_singlets = evalzsinglets(u_rad, A.core.d_singlets,
    A.core.αs_singlets, A.core.Ωs_singlets, A.core.λ0,
    A.core.κs_λ_singlets, A.zs_singlets)

    return out_sys + out_singlets
end

function evalzsinglets(u_rad::T, d::Vector{T},
    αs_singlets::Vector{T}, Ωs_singlets, λ0::T, λ_multipliers::Vector{T},
    zs_singlets::Vector{Complex{T}}) where T <: Real

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for i = 1:length(αs_singlets)
        τ = u_rad - d[i]

        λ = λ0*λ_multipliers[i]
        Ω = Ωs_singlets[i]
        out += zs_singlets[i]*αs_singlets[i]/(λ+im*(τ-Ω))
    end
    return out
end

function evalzitpproxysys(zs::Vector{Vector{Complex{T}}}, gs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysFIDType1{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ

    @assert length(d) == length(gs)

    out = zero(Complex{T})

    #u_rad = 2*π*u
    for i = 1:length(gs)
        r = u_rad - d[i]

        for k = 1:length(gs[i])
            out += zs[i][k]*gs[i][k](r, κs_λ[i])
        end
    end

    return out
end

function evalκitpproxysys(zs::Vector{Vector{Complex{T}}}, gs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysFIDType2{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ

    @assert length(d) == length(gs)

    out = zero(Complex{T})

    #u_rad = 2*π*u
    for i = 1:length(gs)

        for k = 1:length(gs[i])
            r = u_rad - d[i][k]

            out += zs[i][k]*gs[i][k](r, κs_λ[i][k])
        end
    end

    return out
end

ext: jl


text:
using Makie
using GeometryBasics
using FractalTools

# Construct 2d mesh
msh2 = GeometryBasics.mesh(Tesselation(Rect2D(0., 0., 10., 10.), (50, 50)))

# Define functions, nodes and observables
f1(x, y) = sin(x) * sin(y)
f2(x, y) = sin(x) * cos(y)
f3(x, y) = cos(x) * sin(y)
f4(x, y) = cos(x) * cos(y)
func = Node{Any}(f1)
msh3 = lift(func) do f
    pnts = [Point(pnt[1], pnt[2], f(pnt[1], pnt[2])) for pnt in msh2.position]
    fcs = faces(msh2)
    GeometryBasics.Mesh(pnts, fcs)
end

# Consruct a figure
fig = Figure()

# Insert a menu
menu = Menu(fig, options=zip(["1", "2", "3", "4"], [f1, f2, f3, f4]))
on(menu.selection) do s
    func[] = s
    update_cam!(ls.scene, Vec3(10.), Vec3(0.))
    center!(ls.scene)
end
gd = fig[1, 1] = vgrid!(Label(fig, "Choose function"), menu,tellheight=false)

# Insert a mehs plot
ls = fig[1, 2] = LScene(fig, scenekw=(camera=cam3d!, raw=false), tellheight=false, tellwidth=false)
foo = mesh!(ls.scene, msh3, color=first.(msh3[].position))
wireframe!(ls.scene, msh3, linewidth=3)

# Display figure
display(fig)


text:
module PhysVectors

import Base: convert
import LinearAlgebra

abstract type PhysVector end
abstract type PositionVector <: PhysVector end
abstract type DisplacementVector <: PhysVector end
abstract type MomentumVector <: PhysVector end

"""
  Generic three vector, used where you don't want the type system to restrict operations.
"""
struct XYZVector <: PhysVector
    x::Real
    y::Real
    z::Real
end
Base.show(io::IO, v::XYZVector) = print(io, "[x: ",v.x,", y:",v.y,", z:",v.z,"]")
Base.:+(v1::XYZVector, v2::XYZVector) = XYZVector(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z)
Base.:-(v1::XYZVector, v2::XYZVector) = XYZVector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z)

"""
  Generic four vector, used where you don't want the type system to restrict operations.
"""
struct XYZTVector <: PhysVector
    x::Real
    y::Real
    z::Real
    t::Real
end
Base.show(io::IO, v::XYZTVector) = print(io, "[x: ",v.x,", y: ",v.y,", z:",v.z,", t:",v.t,"]")
Base.:+(v1::XYZTVector, v2::XYZTVector) = XYZTVector(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.t + v2.t)
Base.:-(v1::XYZTVector, v2::XYZTVector) = XYZTVector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.t - v2.t)

"""
    Three vector describing a spatial displacement.

    Displacements can be summed, subtracted, scaled, and rotated.
"""
struct Displacement3Vector <: DisplacementVector
    dx::Real
    dy::Real
    dz::Real
end
Base.show(io::IO, v::Displacement3Vector) = print(io, "[disp: ",v.dx,",",v.dy,",",v.dz,"]")
Base.:+(v1::Displacement3Vector, v2::Displacement3Vector) = Displacement3Vector(v1.dx + v2.dx, v1.dy + v2.dy, v1.dz + v2.dz)
Base.:-(v1::Displacement3Vector, v2::Displacement3Vector) = Displacement3Vector(v1.dx - v2.dx, v1.dy - v2.dy, v1.dz - v2.dz)

"""
    Four vector describing a spacetime displacement.

    Four displacements can be boosted, summed, subtracted, and their spatial components can be scaled and rotated.
"""
struct Displacement4Vector <: DisplacementVector
    dx::Real
    dy::Real
    dz::Real
    dt::Real
end
Base.show(io::IO, v::Displacement4Vector) = print(io, "[disp: ",v.dx,",",v.dy,",",v.dz,",",v.dt,"]")
Base.:+(v1::Displacement4Vector, v2::Displacement4Vector) = Displacement4Vector(v1.dx + v2.dx, v1.dy + v2.dy, v1.dz + v2.dz, v1.dt + v2.dt)
Base.:-(v1::Displacement4Vector, v2::Displacement4Vector) = Displacement4Vector(v1.dx - v2.dx, v1.dy - v2.dy, v1.dz - v2.dz, v1.dt - v2.dt)

"""
    Three vector describing a spatial position.

    Can be handled by a restricted set of operators, e.g. you cannot add positions, but you can subtract them, producing a displacement vector.

    Positions can be rotated and translated (with a displacement vector).
"""
struct Position3Vector <: PositionVector
    x::Real
    y::Real
    z::Real
end
Base.show(io::IO, v::Position3Vector) = print(io, "[pos: ",v.x,",",v.y,",",v.z,"]")

Base.:-(v1::Position3Vector, v2::Position3Vector) = Displacement3Vector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z)
Base.:+(v1::Position3Vector, v2::Displacement3Vector) = Position3Vector(v1.x + v2.dx, v1.y + v2.dy, v1.z + v2.dz)
Base.:-(v1::Position3Vector, v2::Displacement3Vector) = Position3Vector(v1.x - v2.dx, v1.y - v2.dy, v1.z - v2.dz)

"""
    Four vector describing a spacetime position.

    Can be handled by a restricted set of operators, e.g. you cannot add positions, but you can subtract them, producing a displacement vector.

    Four positions can be boosted, and translated, their spatial components can be rotated.
"""
struct Position4Vector <: PositionVector
    x::Real
    y::Real
    z::Real
    t::Real
end
Base.show(io::IO, v::Position4Vector) = print(io, "[pos: ",v.x,",",v.y,",",v.z,",",v.t,"]")

Base.:-(v1::Position4Vector, v2::Position4Vector) = Displacement4Vector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.t - v2.t)
Base.:+(v1::Position4Vector, v2::Displacement4Vector) = Position4Vector(v1.x + v2.dx, v1.y + v2.dy, v1.z + v2.dz, v1.t + v2.dt)
Base.:-(v1::Position4Vector, v2::Displacement4Vector) = Position4Vector(v1.x - v2.dx, v1.y - v2.dy, v1.z - v2.dz, v1.t - v2.dt)

"""
    Three vector describing a momentum.

    Three momenta can be summed, subtracted, scaled, and rotated.
"""
struct Momentum3Vector <: MomentumVector
    px::Real
    py::Real
    pz::Real
end
Base.show(io::IO, v::Momentum3Vector) = print(io, "[mom: ",v.px,",",v.py,",",v.pz,"]")
Base.:+(v1::Momentum3Vector, v2::Momentum3Vector) = Momentum3Vector(v1.dx + v2.dx, v1.dy + v2.dy, v1.dz + v2.dz)
Base.:-(v1::Momentum3Vector, v2::Momentum3Vector) = Momentum3Vector(v1.dx - v2.dx, v1.dy - v2.dy, v1.dz - v2.dz)

"""
    Four vector describing a four momentum.

    Four momenta can be boosted, summed, subtracted.
"""
struct Momentum4Vector <: MomentumVector
    px::Real
    py::Real
    pz::Real
    e::Real
end
Base.show(io::IO, v::Momentum4Vector) = print(io, "[mom: ",v.px,",",v.py,",",v.pz,",",v.e,"]")
Base.:+(v1::Momentum4Vector, v2::Momentum4Vector) = Momentum4Vector(v1.dx + v2.dx, v1.dy + v2.dy, v1.dz + v2.dz, v1.dt + v2.dt)
Base.:-(v1::Momentum4Vector, v2::Momentum4Vector) = Momentum4Vector(v1.dx - v2.dx, v1.dy - v2.dy, v1.dz - v2.dz, v1.dt - v2.dt)

"""
    Internal function for getting the spatial component of a fourvector while respecting the vector type hierarchy
"""
function _threecomp(v::XYZTVector)
    XYZVector(v.x,v.y,v.z)
end
_threecomp(v::Position4Vector) = Position3Vector(v.x,v.y,v.z)
_threecomp(v::Displacement4Vector) = Displacement3Vector(v.dx,v.dy,v.dz)
_threecomp(v::Momentum4Vector) = Momentum3Vector(v.px,v.py,v.pz)

"""
    Internal function for array-izing typed vectors to use julia BLAS
"""
function _toarr(v::XYZVector)
    [v.x,v.y,v.z]
end
_toarr(v::XYZTVector) = [v.x,v.y,v.z,v.t]
_toarr(v::Position3Vector) = [v.x,v.y,v.z]
_toarr(v::Position4Vector) = [v.x,v.y,v.z,v.t]
_toarr(v::Displacement3Vector) = [v.dx,v.dy,v.dz]
_toarr(v::Displacement4Vector) = [v.dx,v.dy,v.dz,v.dt]
_toarr(v::Momentum3Vector) = [v.px,v.py,v.pz]
_toarr(v::Momentum4Vector) = [v.px,v.py,v.pz,v.e]

ThreeVectors = Union{
    XYZVector,
    Position3Vector,
    Displacement3Vector,
    Momentum3Vector
}

LinearAlgebra.:⋅(v1::ThreeVectors, v2::ThreeVectors) = _toarr(v1) ⋅ _toarr(v2)
LinearAlgebra.:×(v1::ThreeVectors, v2::ThreeVectors) = _toarr(v1) × _toarr(v2)

mag2(v::ThreeVectors) = LinearAlgebra.norm(_toarr(v))
mag(v::ThreeVectors) = sqrt(mag2(v))

DirectionVectors = Union{
    XYZVector,
    Displacement3Vector,
    Momentum3Vector
}

Base.convert(::Type{XYZVector}, x::Displacement3Vector) = XYZVector(x.dx, x.dy, x.dz)
Base.convert(::Type{XYZVector}, x::Momentum3Vector) = XYZVector(x.px, x.py, x.pz)
Base.convert(::Type{XYZTVector}, x::Displacement4Vector) = XYZTVector(x.dx, x.dy, x.dz, x.dt)
Base.convert(::Type{XYZTVector}, x::Momentum4Vector) = XYZTVector(x.px, x.py, x.pz, x.e)

"""
    Build a three rotation matrix around the x axis
"""
function rotxmat(ang)
    c = cos(ang)
    s = sin(ang)
    [  1  0  0  ;
       0  c -s  ;
       0  s  c  ]
end

"""
    Build a three rotation matrix around the y axis
"""
function rotymat(ang)
    c = cos(ang)
    s = sin(ang)
    [  c  0  s  ;
       0  1  0  ;
      -s  0  c  ]
end

"""
    Build a three rotation matrix around the z axis
"""
function rotzmat(ang)
    c = cos(ang)
    s = sin(ang)
    [  c -s  0  ;
       s  c  0  ;
       0  0  1  ]
end

"""
    Rotate a three vector by rotation matrix, respects the physics vector type
"""
function rotate(v::T, matrix::Array{R,2}) where {T<:ThreeVectors, R<:Real}
    rotated = matrix * _toarr(v)
    T(rotated[1], rotated[2], rotated[3])
end

"""
    Rotate a three vector around x by an angle, respects the physics vector type
"""
rotateaboutx(v::T, ang::Real) where {T<:ThreeVectors} = rotate(v, rotxmat(ang))
"""
    Rotate a three vector around y by an angle, respects the physics vector type
"""
rotateabouty(v::T, ang::Real) where {T<:ThreeVectors} = rotate(v, rotymat(ang))
"""
    Rotate a three vector around z by an angle, respects the physics vector type
"""
rotateaboutz(v::T, ang::Real) where {T<:ThreeVectors} = rotate(v, rotzmat(ang))

"""
    Rotate a three vector by roll around x, then by pitch around y, then by yaw around z, respects the physics vector type
"""
function rotatexyz(v::T, angyaw::Real, angpitch::Real, angroll::Real)  where {T<:ThreeVectors}
    rotate(v, rotzmat(angyaw) * rotymat(angpitch) * rotxmat(angroll))
end

"""
    Build a three rotation matrix around an arbitrary three axis
"""
function rotationmat(axis::DirectionVectors, ang::Real)
    # From https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
    v = _toarr(unit(axis))

    c = cos(ang)
    s = sin(ang)
    diagc = LinearAlgebra.Diagonal([c,c,c])
    crossprodmat = [     0 -v[3]  v[2] ;
                      v[3]     0 -v[1] ;
                     -v[2]  v[1]     0 ]
    outerprod = v * v'

    diagc + s * crossprodmat + (1-c) * outerprod
end

"""
    Rotate a three vector by roll around x, then by pitch around y, then by yaw around z, respects the physics vector type
"""
rotateaboutaxis(v::T, axis::DirectionVectors, ang::Real) where {T<:ThreeVectors} = rotate(v, rotationmat(axis, ang))

ScaleableVectors = Union{
    XYZVector,
    Displacement3Vector,
    Momentum3Vector
}

"""
    Scales a three vector by a different scale along each principal axis.
"""
function scale(v::T, xs::Real, ys::Real, zs::Real) where {T<:ScaleableVectors}
    scaled = [xs, ys, zs] .* _toarr(v)
    T(scaled[1], scaled[2], scaled[3])
end

"""
    Scales a three vector by a single scalefactor
"""
scale(v::ScaleableVectors, s::Real) = scale(v,s,s,s)

"""
    Normalizes a vector to length 1
"""
unit(v::ScaleableVectors) = scale(v, 1.0/mag(v))

function Base.:*(s::Real, v::T) where {T<:ScaleableVectors}
    scale(v, s)
end
Base.:*(v::T, s::Real) where {T<:ScaleableVectors} = s * v
Base.:/(v::T, s::Real) where {T<:ScaleableVectors} = (1.0/s) * v
Base.:-(v::T) where {T<:ScaleableVectors} = -1.0 * v

FourVectors = Union{
    XYZTVector,
    Position4Vector,
    Displacement4Vector,
    Momentum4Vector
}

function LinearAlgebra.:⋅(v1::FourVectors, v2::FourVectors)
    a1 = _toarr(v1)
    a2 = _toarr(v2)
    (a1[4] * a2[4]) - LinearAlgebra.:⋅(a1[1:3], a2[1:3])
end

threemag2(v::FourVectors) = mag2(_threecomp(v))
threemag(v::FourVectors) = sqrt(threemag2(v))

function fourmag2(v::FourVectors)
    a = _toarr(v)
    a[4]^2 - (a[1]^2 + a[2]^2 + a[3]^2)
end

function fourmag(v::FourVectors)
    sqrt(fourmag2(v))
end

BoostableVectors = Union{
    XYZTVector,
    Position4Vector,
    Momentum4Vector
}

function boost(r::T, dir::DirectionVectors, β::Real) where {T<:BoostableVectors}
    # From https://en.wikipedia.org/wiki/Lorentz_transformation#Proper_transformations
    n = Base.convert(XYZVector, unit(dir))
    γ = 1.0/(sqrt(1.0 - β^2))

    boosted =
    [ 1+(γ-1)*(n.x^2)    (γ-1)*n.x*n.y    (γ-1)*n.x*n.z  -γ*β*n.x;
        (γ-1)*n.y*n.z  1+(γ-1)*(n.y^2)    (γ-1)*n.y*n.z  -γ*β*n.y;
        (γ-1)*n.x*n.z    (γ-1)*n.y*n.z  1+(γ-1)*(n.z^2)  -γ*β*n.z;
             -γ*β*n.x         -γ*β*n.y         -γ*β*n.z         γ ] * _toarr(r)
    T(boosted[1],boosted[2],boosted[3],boosted[4])
end

function boost(r::T, v::XYZVector) where {T<:BoostableVectors}
    boost(r, v, mag(v))
end

β(v::Momentum4Vector) = threemag(v)/v.e
γ(v::Momentum4Vector) = 1.0/sqrt(1.0 - β(v)^2)

function getboostvector(v::Momentum4Vector)
  β(v) * unit(_threecomp(v))
end

mass2(v::Momentum4Vector) = fourmag2(v)
mass(v::Momentum4Vector) = fourmag(v)

end

max_stars_repo_name: luketpickering/nujuphys

max_issues_repo_name: luketpickering/nujuphys

max_forks_repo_name: luketpickering/nujuphys


text:
module Types

import Base: ==

mutable struct PkgMirror
    name::String
    url::String
end

==(a::PkgMirror, b::PkgMirror) = a.name == b.name && a.url == b.url

mutable struct Finalizer end
FINALIZER = Finalizer()

end


text:
module NearestNeighborModels

# ===================================================================
# IMPORTS
import MLJModelInterface
import MLJModelInterface: @mlj_model, metadata_model, metadata_pkg,
    Table, Continuous, Count, Finite, OrderedFactor, Multiclass
import NearestNeighbors
import StatsBase
import Tables

using Distances
using FillArrays
using LinearAlgebra
using Statistics

# ===================================================================
## EXPORTS
export list_kernels, ColumnTable, DictTable

# Export KNN models
# KNN models are exported automatically by `@mjl_model`

# Re-Export Distance Metrics from `Distances.jl`
export Euclidean, Cityblock, Minkowski, Chebyshev, Hamming, WeightedEuclidean,
    WeightedCityblock, WeightedMinkowski

# Export KNN Kernels
export DualU, DualD, Dudani, Fibonacci, Inverse, ISquared, KNNKernel, Macleod, Rank,
    ReciprocalRank, UDK, Uniform, UserDefinedKernel, Zavreal

# ===================================================================
## CONSTANTS
const Vec{T} = AbstractVector{T}
const Mat{T} = AbstractMatrix{T}
const Arr{T, N} = AbstractArray{T, N}
const ColumnTable =  Tables.ColumnTable
const DictTable = Tables.DictColumns
const MultiUnivariateFinite = Union{DictTable, ColumnTable}

# Define constants for easy referencing of packages
const MMI = MLJModelInterface
const NN = NearestNeighbors
const PKG = "NearestNeighborModels"

# Definitions of model descriptions for use in model doc-strings.
const KNNRegressorDescription = """
    K-Nearest Neighbors regressor: predicts the response associated with a new point
    by taking an weighted average of the response of the K-nearest points.
    """

const KNNClassifierDescription = """
    K-Nearest Neighbors classifier: predicts the class associated with a new point
    by taking a vote over the classes of the K-nearest points.
    """

const KNNCoreFields = """
    * `K::Int=5` : number of neighbors
    * `algorithm::Symbol = :kdtree` : one of `(:kdtree, :brutetree, :balltree)`
    * `metric::Metric = Euclidean()` : a `Metric` object for the distance between points
    * `leafsize::Int = 10` : at what number of points to stop splitting the tree
    * `reorder::Bool = true` : if true puts points close in distance close in memory
    * `weights::KNNKernel=Uniform()` : kernel used in assigning weights to the
        k-nearest neighbors for each observation. An instance of one of the types in
        `list_kernels()`. User-defined weighting functions can be passed by wrapping the
        function in a `UDF` kernel. If sample weights `w` are passed during machine
        construction e.g `machine(model, X, y, w)` then the weight assigned to each
        neighbor is the product of the `KNNKernel` generated weight and the corresponding
        neighbor sample weight.

    """

const SeeAlso = """
    See also the
    [package documentation](https://github.com/KristofferC/NearestNeighbors.jl).
    For more information about the kernels see the paper by Geler et.al
    [Comparison of different weighting schemes for the kNN classifier
    on time-series data]
    (https://perun.pmf.uns.ac.rs/radovanovic/publications/2016-kais-knn-weighting.pdf).
    """

const MultitargetKNNClassifierFields = """
    ## Keywords Parameters

    $KNNCoreFields
    * `output_type::Type{<:MultiUnivariateFinite}=DictTable` : One of
       (`ColumnTable`, `DictTable`). The type of table type to use for predictions.
       Setting to `ColumnTable` might improve performance for narrow tables while setting to
       `DictTable` improves performance for wide tables.

    $SeeAlso

    """

const KNNFields = """
    ## Keywords Parameters

    $KNNCoreFields

    $SeeAlso

    """

# ===================================================================
# Includes
include("utils.jl")
include("kernels.jl")
include("models.jl")

# ===================================================================
# List of all models interfaced
const MODELS = (
    KNNClassifier, KNNRegressor, MultitargetKNNRegressor, MultitargetKNNClassifier
)

# ====================================================================
# PKG_METADATA
metadata_pkg.(
    MODELS,
    name = "NearestNeighborModels",
    uuid = "6f286f6a-111f-5878-ab1e-185364afe411",
    url = "https://github.com/alan-turing-institute/NearestNeighborModels.jl",
    license = "MIT",
    julia = true,
    is_wrapper = false
)

end # module


text:
export ClassEncoding, BinaryClassEncoding, MultinomialClassEncoding
export OneHotClassEncoding
export labelencode, labeldecode, groupindices
export nclasses, labels, classDistribution

import MLBase.labelencode
import MLBase.labeldecode
import MLBase.groupindices
import Base.show

# ==========================================================================

abstract ClassEncoding

abstract BinaryClassEncoding <: ClassEncoding
abstract MultinomialClassEncoding <: ClassEncoding

# ==========================================================================

module ClassEncodings

export ZeroOneClassEncoding, SignedClassEncoding, MultivalueClassEncoding, OneOfKClassEncoding, OneHotClassEncoding

using MLBase
using ..SupervisedLearning.ClassEncoding
using ..SupervisedLearning.BinaryClassEncoding
using ..SupervisedLearning.MultinomialClassEncoding

immutable ZeroOneClassEncoding{T} <: BinaryClassEncoding
  labelmap::LabelMap{T}

  function ZeroOneClassEncoding(labelmap::LabelMap{T})
    numLabels = length(labelmap.vs)
    if numLabels != 2
      throw(ArgumentError("The given target vector must have exactly two classes"))
    end
    new(labelmap)
  end
end

ZeroOneClassEncoding{T}(lm::LabelMap{T}) =
  ZeroOneClassEncoding{T}(lm)

ZeroOneClassEncoding{T}(targets::Vector{T}) =
  ZeroOneClassEncoding{T}(labelmap(targets))

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

immutable SignedClassEncoding{T} <: BinaryClassEncoding
  labelmap::LabelMap{T}

  function SignedClassEncoding(labelmap::LabelMap{T})
    numLabels = length(labelmap.vs)
    if numLabels != 2
      throw(ArgumentError("The given target vector must have exactly two classes"))
    end
    new(labelmap)
  end
end

SignedClassEncoding{T}(lm::LabelMap{T}) =
  SignedClassEncoding{T}(lm)

SignedClassEncoding{T}(targets::Vector{T}) =
  SignedClassEncoding{T}(labelmap(targets))

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

immutable MultivalueClassEncoding{T} <: MultinomialClassEncoding
  labelmap::LabelMap{T}
  nlabels::Int
  zeroBased::Bool

  function MultivalueClassEncoding(labelmap::LabelMap{T}, zeroBased = false)
    numLabels = length(labelmap.vs)
    if numLabels < 2
      throw(ArgumentError("The given target vector has less than two classes"))
    end
    new(labelmap, numLabels, zeroBased)
  end
end

MultivalueClassEncoding{T}(lm::LabelMap{T}; zero_based = false) =
  MultivalueClassEncoding{T}(lm, zero_based)

MultivalueClassEncoding{T}(targets::Vector{T}; zero_based = false) =
  MultivalueClassEncoding{T}(labelmap(targets), zero_based)

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

immutable OneOfKClassEncoding{T} <: MultinomialClassEncoding
  labelmap::LabelMap{T}
  nlabels::Int

  function OneOfKClassEncoding(labelmap::LabelMap{T})
    numLabels = length(labelmap.vs)
    if numLabels < 2
      throw(ArgumentError("The given target vector has less than two classes"))
    end
    new(labelmap, numLabels)
  end
end

OneOfKClassEncoding{T}(lm::LabelMap{T}) =
  OneOfKClassEncoding{T}(lm)

OneOfKClassEncoding{T}(targets::Vector{T}) =
  OneOfKClassEncoding{T}(labelmap(targets))

end

# ==========================================================================

@reexport using SupervisedLearning.ClassEncodings

typealias OneHotClassEncoding OneOfKClassEncoding

# ==========================================================================

nclasses(ce::BinaryClassEncoding) = 2
nclasses(ce::MultinomialClassEncoding) = ce.nlabels

labels{C<:ClassEncoding}(ce::C) = ce.labelmap.vs

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

function groupindices{T}(classEncoding::ClassEncoding, targets::Vector{T})
  groupindices(classEncoding.labelmap, targets)
end

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

function classDistribution{T}(labelmap::LabelMap{T}, targets::Vector{T})
  labelmap.vs, map(length, groupindices(labelmap, targets))
end

function classDistribution{T}(classEncoding::ClassEncoding, targets::Vector{T})
  classEncoding.labelmap.vs, map(length, groupindices(classEncoding, targets))
end

# ==========================================================================

function labelencode{T}(classEncoding::ZeroOneClassEncoding{T}, targets::Vector{T})
  indicies = labelencode(classEncoding.labelmap, targets)
  float(indicies - 1)
end

function labeldecode{T}(classEncoding::ZeroOneClassEncoding{T}, values::AbstractVector{Float64})
  indicies = round(Integer, values + 1)
  labeldecode(classEncoding.labelmap, indicies)
end

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

function labelencode{T}(classEncoding::SignedClassEncoding{T}, targets::Vector{T})
  indicies = labelencode(classEncoding.labelmap, targets)
  2(indicies - 1.5)
end

function labeldecode{T}(classEncoding::SignedClassEncoding{T}, values::AbstractVector{Float64})
  indicies = round(Integer, (values / 2.) + 1.5)
  labeldecode(classEncoding.labelmap, indicies)
end

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

function labelencode{T}(classEncoding::MultivalueClassEncoding{T}, targets::Vector{T})
  labelencode(classEncoding.labelmap, targets) - classEncoding.zeroBased*1.
end

function labeldecode{T}(classEncoding::MultivalueClassEncoding{T}, values::AbstractVector{Float64})
  indicies = round(Integer, values + classEncoding.zeroBased*1.)
  labeldecode(classEncoding.labelmap, indicies)
end

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

function labelencode{T}(classEncoding::OneOfKClassEncoding{T}, targets::Vector{T})
  indicies = labelencode(classEncoding.labelmap, targets)
  convert(Matrix{Float64}, indicatormat(indicies)) # this doesn't work if the indexseq is broken (e.g. [1,2,5])
end

function labeldecode{T}(classEncoding::OneOfKClassEncoding{T}, values::AbstractMatrix{Float64})
  numLabels = classEncoding.nlabels
  indicies = map(safeRound, values' * collect(1:numLabels))
  labeldecode(classEncoding.labelmap, indicies)
end

# ==========================================================================

function getLabelString{T}(labelmap::LabelMap{T})
  labels = labelmap.vs
  c = length(labels)
  if c > 10
    labels = labels[1:10]
    labelString = string(join(labels, ", "), ", ... [TRUNC]")
  else
    labelString = join(labels, ", ")
  end
end

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

function show{T}(io::IO, classEncoding::ZeroOneClassEncoding{T})
  labelString = getLabelString(classEncoding.labelmap)
  print(io,
        """
        ZeroOneClassEncoding (Binary) to {0, 1}
          .labelmap  ...  encoding for: {$labelString}""")
end

function show{T}(io::IO, classEncoding::SignedClassEncoding{T})
  labelString = getLabelString(classEncoding.labelmap)
  print(io,
        """
        SignedClassEncoding (Binary) to {-1, 1}
          .labelmap  ...  encoding for: {$labelString}""")
end

function show{T}(io::IO, classEncoding::MultivalueClassEncoding{T})
  c = classEncoding.nlabels
  zB = classEncoding.zeroBased
  labelString = getLabelString(classEncoding.labelmap)
  print(io,
        """
        MultivalueClassEncoding (Multinomial) to range $((1:c) - zB*1)
          .nlabels   ...  $c classes
          .labelmap  ...  encoding for: {$labelString}""")
end

function show{T}(io::IO, classEncoding::OneOfKClassEncoding{T})
  c = classEncoding.nlabels
  labelString = getLabelString(classEncoding.labelmap)
  print(io,
        """
        OneOfKClassEncoding (Multinomial) to one-out-of-$c hot-vector
          .nlabels   ...  $c classes
          .labelmap  ...  encoding for: {$labelString}""")
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
for a in ("", "t", "c")
  for b in ("", "t", "c")
    local sym = a == b == "" ? :* : Symbol("A", a, "_mul_B", b)
    if !isdefined(Base, sym) continue end
    q = QuoteNode(sym)

    @eval begin
      function Base.$(sym)(x::ANode, y::ANode)
        @assert size(x, $(a == "" ? 2 : 1)) == size(y, $(b == "" ? 1 : 2))
        ANode($q, (x, y), (size(x, $(a == "" ? 1 : 2)), size(y, $(b == "" ? 2 : 1))))
      end

      function do_forward!(n::ANode{$q}, out, in1, in2)
        BLAS.gemm!($(a == "" ? 'N' : 'T'), $(b == "" ? 'N' : 'T'), one(eltype(out)), in1, in2, zero(eltype(out)), out)
      end

      gradient_node(n::ANode{$q}, ::Type{Val{1}}, bp::ANode) = $(a == "" ? (b=="" ? :(bp*n.input[2].') : :(bp*n.input[2])) : (b=="" ? :(n.input[2]*bp.') : :(n.input[2].'*bp.')))

      gradient_node(n::ANode{$q}, ::Type{Val{2}}, bp::ANode) = $(b == "" ? (a=="" ? :(n.input[1].'*bp) : :(n.input[1]*bp)) : (a=="" ? :(bp.'*n.input[1]) : :(bp.'*n.input[1].')))
    end
  end
end


text:
# Business Calendars (adapted from Ito.jl and BusinessDays.jl)
using Base.Dates

abstract BusinessCalendar

abstract WesternCalendar <: BusinessCalendar
abstract OrthodoxCalendar <: BusinessCalendar

# target calendar
type TargetCalendar <: BusinessCalendar end

# for simply moving foward and backward in time
type NullCalendar <: BusinessCalendar end

type JointCalendar{B <: BusinessCalendar, C <: BusinessCalendar} <: BusinessCalendar
  cal1::B
  cal2::C
end

# US Calendars
abstract UnitedStatesCalendar <: WesternCalendar

type USSettlementCalendar <: UnitedStatesCalendar; end
type USNYSECalendar <: UnitedStatesCalendar; end
type USNERCCalendar <: UnitedStatesCalendar; end
type USGovernmentBondCalendar <: UnitedStatesCalendar; end

# UK Calendars
abstract UnitedKingdomCalendar <: WesternCalendar

type UKSettlementCalendar <: UnitedKingdomCalendar end
type UKLSECalendar <: UnitedKingdomCalendar end
type UKLMECalendar <: UnitedKingdomCalendar end

abstract BusinessDayConvention
type Unadjusted <: BusinessDayConvention end
type ModifiedFollowing <: BusinessDayConvention end
type Following <: BusinessDayConvention end

# easter functions
function easter_rata(y::Int)

  c::Int64
	e::Int64
	p::Int64

   # Algo R only works after 1582
   if y < 1582
        # Are you using this? Send me a postcard!
        error("Year cannot be less than 1582. Provided: $(y).")
   end

	# Century
   c = div( y , 100) + 1

   # Shifted Epact
   e = mod(14 + 11*(mod(y, 19)) - div(3*c, 4) + div(5+8*c, 25), 30)

   # Adjust Epact
   if (e == 0) || ((e == 1) && ( 10 < mod(y, 19) ))
   	e += 1
   end

   # Paschal Moon
   p = Date(y, 4, 19).instant.periods.value - e

   # Easter: locate the Sunday after the Paschal Moon
   return p + 7 - mod(p, 7)
end

# Returns Date
function easter_date(y::Int)
	# Compute the gregorian date for Rata Die number
     return Date(Dates.rata2datetime( easter_rata(y) ))
end

# calendar functions
# advance{B <: BusinessDayConvention}(time_period::Day, cal::NullCalendar, dt::Date, ::B) = dt += time_period
# advance{B <: BusinessDayConvention}(time_period::Union{Week, Month, Year}, cal::NullCalendar, dt::Date, ::B) = dt += time_period

function advance(days::Day, cal::BusinessCalendar, dt::Date, biz_conv::BusinessDayConvention = Following())
  n = Int(days)
  if n > 0
    while n > 0
      dt += Day(1)
      while !is_business_day(cal, dt)
        dt += Day(1)
      end
      n -= 1
    end
  else
    while (n < 0)
      dt -= Day(1)
      while !is_business_day(cal, dt)
        dt -= Day(1)
      end
      n += 1
    end
  end

  return dt
end

function advance(time_period::Union{Week, Month, Year}, cal::BusinessCalendar, dt::Date, biz_conv::BusinessDayConvention = Following())
  dt += time_period
  return adjust(cal, biz_conv, dt)
end

is_business_day(cal::NullCalendar, ::Date) = true

function is_business_day(cal::BusinessCalendar, dt::Date)
  if dayofweek(dt) in [6, 7] || is_holiday(cal, dt)
    return false
  else
    return true
  end
end

# In the United States, if a holiday falls on Saturday, it's observed on the preceding Friday.
# If it falls on Sunday, it's observed on the next Monday.
function adjustweekendholidayUS(dt::Date)
	if dayofweek(dt) == 6
		return dt - Dates.Day(1)
	end

	if dayofweek(dt) == 7
		return dt + Dates.Day(1)
	end

	return dt
end

is_holiday(joint::JointCalendar, dt::Date) = is_holiday(joint.cal1, dt) || is_holiday(joint.cal2, dt)

function is_holiday(::USSettlementCalendar , dt::Date)

	const yy = year(dt)
	const mm = month(dt)
	const dd = day(dt)

	if (
			# New Year's Day
			adjustweekendholidayUS(Date(yy, 1, 1)) == dt
			||
			# New Year's Day on the previous year when 1st Jan is Saturday
			(mm == 12 &&  dd == 31 && dayofweek(dt) == Friday)
			||
			# Birthday of Martin Luther King, Jr.
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 1)
			||
			# Washington's Birthday
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 2)
			||
			# Memorial Day is the last Monday in May
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == daysofweekinmonth(dt) && mm == 5)
			||
			# Independence Day
			adjustweekendholidayUS(Date(yy, 7, 4)) == dt
			||
			# Labor Day is the first Monday in September
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == 1 && mm == 9)
			||
			# Columbus Day is the second Monday in October
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == 2 && mm == 10)
			||
			# Veterans Day
			adjustweekendholidayUS(Date(yy, 11, 11)) == dt
			||
			# Thanksgiving Day is the fourth Thursday in November
			(dayofweek(dt) == 4 && dayofweekofmonth(dt) == 4 && mm == 11)
			||
			# Christmas
			adjustweekendholidayUS(Date(yy, 12, 25)) == dt
		)
		return true
	end

	return false
end

function is_holiday(::USGovernmentBondCalendar, dt::Date)
  const yy = year(dt)
	const mm = month(dt)
	const dd = day(dt)
	if (
			# New Year's Day
			adjustweekendholidayUS(Date(yy, 1, 1)) == dt
			||
			# New Year's Day on the previous year when 1st Jan is Saturday
			(mm == 12 &&  dd == 31 && dayofweek(dt) == Friday)
			||
			# Birthday of Martin Luther King, Jr.
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 1)
			||
			# Washington's Birthday
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 2)
			||
      # Good Friday
      easter_date(yy) - Day(2) == dt
      ||
			# Memorial Day is the last Monday in May
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == daysofweekinmonth(dt) && mm == 5)
			||
			# Independence Day
			adjustweekendholidayUS(Date(yy, 7, 4)) == dt
			||
			# Labor Day is the first Monday in September
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == 1 && mm == 9)
			||
			# Columbus Day is the second Monday in October
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == 2 && mm == 10)
			||
			# Veterans Day
			adjustweekendholidayUS(Date(yy, 11, 11)) == dt
			||
			# Thanksgiving Day is the fourth Thursday in November
			(dayofweek(dt) == 4 && dayofweekofmonth(dt) == 4 && mm == 11)
			||
			# Christmas
			adjustweekendholidayUS(Date(yy, 12, 25)) == dt
		)
		return true
	end

	return false
end

function is_holiday(::USNYSECalendar, dt::Date)
  const yy = year(dt)
	const mm = month(dt)
	const dd = day(dt)
	if (
			# New Year's Day
			adjustweekendholidayUS(Date(yy, 1, 1)) == dt
			||
			# New Year's Day on the previous year when 1st Jan is Saturday
			(mm == 12 &&  dd == 31 && dayofweek(dt) == Friday)
			||
      # Birthday of Martin Luther King, Jr.
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 1)
			||
			# Washington's Birthday
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) ==3 && mm == 2)
			||
      # Good Friday
      easter_date(yy) - Day(2) == dt
      ||
			# Memorial Day is the last Monday in May
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == daysofweekinmonth(dt) && mm == 5)
			||
			# Independence Day
			adjustweekendholidayUS(Date(yy, 7, 4)) == dt
			||
			# Labor Day is the first Monday in September
			(dayofweek(dt) == 1 && dayofweekofmonth(dt) == 1 && mm == 9)
			||
			# Thanksgiving Day is the fourth Thursday in November
			(dayofweek(dt) == 4 && dayofweekofmonth(dt) == 4 && mm == 11)
			||
			# Christmas
			adjustweekendholidayUS(Date(yy, 12, 25)) == dt
		)
		return true
	end

  # Special Closings
  if (
    # Hurricane Sandy
    (yy == 2012 && mm == 10 && dd in (29, 30))
    ||
    # President Ford's funeral
    (yy == 2007 && mm == 1 && dd == 2)
    ||
    # President Reagan's funeral
    (yy == 2004 && mm == 6 && dd == 11)
    ||
    # 9/11
    (yy == 2001 && mm == 9 && dd in (11, 12, 13, 14))
    )
    return true
  end

	return false
end

function adjustweekendholidayUK(dt::Date)

  if dayofweek(dt) == 6
		return dt + Day(2)
	end

  if dayofweek(dt) == 7
		return dt + Day(1)
	end

	return dt
end

## UK Calendar functions
function is_holiday(c::Union{UKSettlementCalendar, UKLSECalendar}, dt::Date)
  const yy = year(dt)
	const mm = month(dt)
	const dd = day(dt)

	if (
    # New Year's Day
    adjustweekendholidayUK(Date(yy, 1, 1)) == dt
    ||
    # Good Friday
    easter_date(yy) - Day(2) == dt
    ||
    # Easter MONDAY
    easter_date(yy) + Day(1) == dt
    ||
    # first MONDAY of May (Early May Bank Holiday)
    (dayofweek(dt) == 1 && dayofweekofmonth(dt) == 1 && mm == 5)
    ||
    # last MONDAY of MAY (Spring Bank Holiday)
    (dayofweek(dt) == 1 && dayofweekofmonth(dt) == daysofweekinmonth(dt) && mm == 5)
    ||
    # last MONDAY of August (Summer Bank Holiday)
    (dayofweek(dt) == 1 && dayofweekofmonth(dt) == daysofweekinmonth(dt) && mm == 8)
    ||
    # Christmas (possibly moved to MONDAY or Tuesday)
    adjustweekendholidayUK(Date(yy, 12, 25)) == dt
    ||
    # Boxing Day (possibly moved to MONDAY or TUESDAY)
    adjustweekendholidayUK(adjustweekendholidayUK(Date(yy, 12, 25)) + Day(1)) == dt
    )
    return true
  end

  # Fixed holidays
  if (
    # Substitute date for Spring Bank Holiday
		(dt == Date(2012, 06, 04))
		||
		# Diamond Jubilee of Queen Elizabeth II.
		(dt == Date(2012, 06, 05))
		||
		# Golden Jubilee of Queen Elizabeth II.
		(dt == Date(2002, 06, 03))
		||
		# Substitute date for Spring Bank Holiday
		(dt == Date(2002, 06, 04))
		||
		# Wedding of Prince William and Catherine Middleton
		(dt == Date(2011, 04, 29))
    )
    return true
  end

  return false
end

function is_holiday(::TargetCalendar, dt::Date)
  const yy = year(dt)
	const mm = month(dt)
	const dd = day(dt)
  const easter_sun = easter_date(yy)

  if (
    # New Years Day
    (mm == 1 && dd == 1)
    ||
    # Good Friday
    easter_sun - Day(2) == dt
    ||
    # Easter Monday
    easter_sun + Day(1) == dt
    ||
    # Int'l Labor Day
    (mm == 5 && dd == 1)
    ||
    # Christmas
    (mm == 12 && dd == 25)
    ||
    # Day of Goodwill
    (mm == 12 && dd == 26)
    )
    return true
  else
    return false
  end
end

# adjustments
adjust(::BusinessCalendar, ::Unadjusted, d::Date) = d

function adjust(cal::BusinessCalendar, ::Union{ModifiedFollowing, Following}, d::Date)
  while !is_business_day(cal, d)
    d += Day(1)
  end

  return d
end

adjust(cal::BusinessCalendar, d::Date) = adjust(cal, Following(), d)


text:
abstract type AbstractExogenous end

abstract type AbstractProcess <: AbstractExogenous end
abstract type DiscreteProcess <: AbstractProcess end
abstract type ContinuousProcess <: AbstractProcess end
abstract type IIDExogenous <: ContinuousProcess end

abstract type AbstractDiscretizedProcess <: DiscreteProcess end

# this is a bit crude but not performance critical, for now
function node(::Type{Point}, dp::DiscreteProcess, i::Int)
    v = node(dp, i)
    SVector{length(v)}(v)
end
function inode(::Type{Point}, dp::DiscreteProcess, i::Int, j::Int)
    v = inode(dp, i, j)
    SVector{length(v)}(v)
end

###
### Discretized process
###

# date-t grid has a known structure
mutable struct DiscretizedProcess{TG<:Grid} <: AbstractDiscretizedProcess
    grid::TG
    integration_nodes::Array{Matrix{Float64},1}
    integration_weights::Array{Vector{Float64},1}
end

n_nodes(dp::DiscretizedProcess) = n_nodes(dp.grid)
node(dp::DiscretizedProcess, i) = node(dp.grid, i)
n_inodes(dp::DiscretizedProcess, i::Int) = size(dp.integration_nodes[i], 1)
inode(dp::DiscretizedProcess, i::Int, j::Int) = dp.integration_nodes[i][j, :]
iweight(dp::DiscretizedProcess, i::Int, j::Int) = dp.integration_weights[i][j]

function Product(gdp1::DiscretizedProcess, gdp2::DiscretizedProcess)
  In = [ gridmake(gdp1.integration_nodes[i], gdp2.integration_nodes[j] ) for i = 1:n_nodes(gdp1)  for j = 1:n_nodes(gdp2) ]
  Iw =[kron(gdp1.integration_weights[i], gdp2.integration_weights[j] ) for i = 1:n_nodes(gdp1)  for j = 1:n_nodes(gdp2) ]
  N=length(gdp1.grid.min)+length(gdp2.grid.min)
  i_grid  = Product(gdp1.grid, gdp2.grid)

  return DiscretizedProcess(i_grid, In, Iw)
end

# date-t grid is unstructured
mutable struct DiscreteMarkovProcess <: AbstractDiscretizedProcess
    grid::UnstructuredGrid
    transitions::Matrix{Float64}
    values::Matrix{Float64}
    i0::Int
end

discretize(::Type{DiscreteMarkovProcess}, mp::DiscreteMarkovProcess) = mp
discretize(mp::DiscreteMarkovProcess) = mp

DiscreteMarkovProcess(transitions::Matrix{Float64}, values::Matrix{Float64}) =
    DiscreteMarkovProcess(UnstructuredGrid{size(values,2)}(values), transitions, values)

DiscreteMarkovProcess(grid::UnstructuredGrid, transitions::Matrix{Float64}, values::Matrix{Float64}) =
    DiscreteMarkovProcess(grid, transitions, values, 1)

n_nodes(dp::DiscreteMarkovProcess) = size(dp.values, 1)
n_inodes(dp::DiscreteMarkovProcess, i::Int) = size(dp.values, 1)
inode(dp::DiscreteMarkovProcess, i::Int, j::Int) = dp.values[j, :]
iweight(dp::DiscreteMarkovProcess, i::Int, j::Int) = dp.transitions[i, j]
node(dp::DiscreteMarkovProcess, i) = dp.values[i, :]
default_index(dp::DiscreteMarkovProcess) = dp.i0

function MarkovProduct(mc1::DiscreteMarkovProcess, mc2::DiscreteMarkovProcess)
    Q = gridmake(mc1.values, mc2.values)
    P = fkron(mc1.transitions, mc2.transitions)
    return DiscreteMarkovProcess(P, Q)
end

function MarkovProduct(mcs::DiscreteMarkovProcess...)
    reduce(MarkovProduct, mcs)
end

# date-t grid is empty

mutable struct DiscretizedIIDProcess <: AbstractDiscretizedProcess
    # nodes::Matrix{Float64}
    grid::EmptyGrid
    integration_nodes::Matrix{Float64}
    integration_weights::Vector{Float64}
end

DiscretizedIIDProcess(x, w) = DiscretizedIIDProcess(EmptyGrid(), x, w)

n_nodes(dp::DiscretizedIIDProcess) = 0
n_inodes(dp::DiscretizedIIDProcess, i::Int) = size(dp.integration_nodes, 1)
inode(dp::DiscretizedIIDProcess, i::Int, j::Int) = dp.integration_nodes[j, :]
iweight(dp::DiscretizedIIDProcess, i::Int, j::Int) = dp.integration_weights[j]
node(dip::DiscretizedIIDProcess, i) = zeros(n_inodes(dip, 1))

# Normal law

struct MvNormal <: IIDExogenous
    mu::Vector{Float64}
    Sigma::Matrix{Float64}
end

MvNormal(Sigma::Matrix{Float64}) = MvNormal(zeros(size(Sigma, 1)), Sigma)
MvNormal(sigma::Float64) = MvNormal(reshape([sigma], 1, 1))

function discretize(mvn::MvNormal)
    n = fill(5, size(mvn.mu))
    x, w = QE.qnwnorm(n, mvn.mu, mvn.Sigma)
    DiscretizedIIDProcess(x, w)
end

function Base.rand(mvn::MvNormal, args...)
    dist = Distributions.MvNormal(mvn.mu, mvn.Sigma)
    return rand(dist, args...)
end

function simulate(mvn::MvNormal, N::Integer, T::Int, e0::Vector{Float64}; stochastic::Bool=true)
    dist = QE.MVNSampler(mvn.mu, mvn.Sigma)
    d = length(mvn.mu)
    out = zeros(d, N, T)
    for i in 1:N
        out[:, i, 1] = e0
    end
    if stochastic
        for t in 2:T
            out[:, :, t] = rand(dist, N)
        end
    end
    AxisArray(out, Axis{:V}(1:d), Axis{:N}(1:N), Axis{:T}(1:T))
end

function simulate(mvn::MvNormal, N::Integer, T::Int; stochastic::Bool=true)
    simulate(mvn, N, T, mvn.mu; stochastic=stochastic)
end

function response(mvn::MvNormal, e1::AbstractVector; T::Int=40)
    d = length(mvn.mu)
    out = zeros(d, T)
    out[:, 2] = e1
    return out
end

function simulate(process::DiscreteMarkovProcess, N::Int, T::Int, i0::Int)
    mc_qe = QE.MarkovChain(process.transitions)
    inds = Array{Int}(undef, T, N)
    QE.simulate_indices!(inds, mc_qe, init=i0)
    AxisArray(inds, Axis{:T}(1:T),  Axis{:N}(1:N))
end

function simulate(process::DiscreteMarkovProcess, N::Int, T::Int, m0::AbstractVector{Float64})
    # try to find index in process.values. IF we do, then simulate using
    # the corresponding row index. If we don't, then throw an error
    for i in 1:size(process.values, 1)
        if process.values[i, :] == m0
            return simulate(process, N, T, i)
        end
    end
    error("Couldn't find the vector `m0` in process.values. "*
          "Try passing an integer `i0` instead")
end

function simulate_values(process::DiscreteMarkovProcess, N::Int, T::Int, i0::Int)
    inds = simulate(process, N, T, i0)
    n_values = size(process.values, 2)
    out_values = Array{Float64}(n_values, N, T)
    for n in 1:N, t in 1:T
        out_values[:, n, t] = process.values[inds[t, n], :]
    end
    AxisArray(out_values, Axis{:n}(1:n_values), Axis{:T}(1:T), Axis{:N}(1:N))
end

# VAR 1

mutable struct VAR1 <: ContinuousProcess
    mu::Array{Float64,1}
    R::Array{Float64,2}
    Sigma::Array{Float64,2}
end

VAR1(;rho::Float64=0.0, Sigma=ones(1,1)) = VAR1(rho, Sigma)

function VAR1(R::Array{Float64,2}, Sigma::Array{Float64,2})
    M = zeros(size(R, 1))
    @assert size(R)==size(Sigma)
    return VAR1(M, R, Sigma)
end

function VAR1(rho::Array{Float64,1}, Sigma::Array{Float64,2})
    R = diagm(rho)
    return VAR1(R, Sigma)
end

function VAR1(rho::Float64, Sigma::Array{Float64,2})
    p = size(Sigma,1)
    R = Matrix(rho*I, p, p)
    return VAR1(R, Sigma)
end

function discretize(var::VAR1)
    d = size(var.R, 1)
    n_states = ones(Int, d)*5
    n_integration = ones(Int, d)*5
    dis = discretize(var, n_states, n_integration)
    return dis
end

function discretize(var::VAR1, n_states::Array{Int,1}, n_integration::Array{Int,1}; n_std::Int=2)
    R = var.R
    M = var.mu
    Sigma = var.Sigma
    S = QE.solve_discrete_lyapunov(R,Sigma)  # asymptotic variance
    sig = diag(S)
    min = var.mu - n_std*sqrt.(sig)
    max = var.mu + n_std*sqrt.(sig)
    grid = CartesianGrid{length(min)}(min,max,n_states)
    # discretize innovations
    x,w = QE.qnwnorm(n_integration, zeros(size(var.Sigma,1)), var.Sigma)
    integration_nodes = [ cat([(M + R*(node(grid, i)-M) + x[j,:])' for j=1:size(x,1)]...; dims=1) for i in 1:n_nodes(grid)]
    integration_weights = [w for i in 1:n_nodes(grid)]
    return DiscretizedProcess(grid, integration_nodes, integration_weights)
end

discretize(::Type{DiscretizedProcess}, var::VAR1; args...) = discretize(var; args...)

function discretize(::Type{DiscreteMarkovProcess}, var::VAR1; N::Int=3)

    # it would be good to have a special type of VAR1 process
    # which has a scalar autoregressive term
    R =  var.R
    d = size(R, 1)
    ρ = R[1, 1]
    @assert maximum(abs, R-Matrix(ρ*I,d,d))<1e-16

    sigma = var.Sigma

    if size(var.Sigma, 1) == 1
        mc_qe = QE.rouwenhorst(N, ρ, sqrt(sigma[1]))
        return DiscreteMarkovProcess(mc_qe.p, appenddim(collect(mc_qe.state_values)))
    end

    L = chol(sigma) # sigma = L'*L
    NN = fill(N, d) # default: same number of nodes in each dimension

    components = [QE.rouwenhorst(N, ρ, 1.0) for N in NN]
    mc_components = [DiscreteMarkovProcess(mc.p, appenddim(collect(mc.state_values))) for mc in components]
    mc_prod = MarkovProduct(mc_components...)
    mc_prod.values = mc_prod.values*L'
    return mc_prod
end

"""
```julia
simulate(var::VAR1, N::Int, T::Int, x0::Vector{Float64}
         stochastic::Bool=true, irf::Bool=false, e0::Vector{Float64}=zeros(0))
```

Construct `N` simulated paths of length `T` the vector autoregressive process
`var`, starting each simulated path from `x0`.

If `irf` is `true`, then `e0` must be a vector of the same length as `var.mu`
and `e0` will be used as the initial shocks.

If `stochastic` is `true`, construct a sequence of shocks by drawing `N` paths
of length `T` from the distribution `Normal(0, var.Sigma)` (where `0` is a
vector of zeros of the appropriate length). Otherwise, if `stochastic` is
false, set all shocks to 0 (unless `irf` is true -- see above).

The output is an `AxisArray` contining variables, paths, and time on the 3 axes
"""
function simulate(var::VAR1, N::Int, T::Int, x0::Vector{Float64};
                  stochastic::Bool=true, irf::Bool=false, e0::Vector{Float64}=zeros(0))

    n = size(var.mu, 1)
    XN = zeros(n, N, T)
    E = zeros(n, N, T)

    if stochastic
        dist = QE.MVNSampler(zeros(n), var.Sigma)
        for jj in 1:N
            E[:, jj, :] = rand(dist, T)
        end
    end

    if irf
        E[:, 1, :] = repeat(e0,1,T)
    end

    # Initial conditions
    for i in 1:N
        XN[:, i, 1] = x0
        for  ii in 1:T-1
            XN[:, i, ii+1] = var.mu+var.R*(XN[:, i, ii]-var.mu)+E[:, i, ii]
        end
    end
    AxisArray(XN, Axis{:V}(1:n), Axis{:N}(1:N), Axis{:T}(1:T))
end

function simulate(var::VAR1, N::Int, T::Int; kwargs...)
    return simulate(var, N, T, var.mu; kwargs...)
end

function response(var::VAR1, x0::AbstractVector,
                  e1::AbstractVector; T::Int=40)
    simulate(var, 1, T, x0; stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1, e1::AbstractVector; T::Int=40)
    simulate(var, 1, T; stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1, x0::AbstractVector, index_s::Int; T::Int=40)
    e1 = zeros(size(var.mu, 1))
    Impulse = sqrt.(diag(var.Sigma)[index_s])
    e1[index_s] = Impulse
    simulate(var, 1, T, x0; stochastic=false, irf=true, e0=e1)[1, :, :]
end

function response(var::VAR1; T::Int=40)
    e1 = sqrt.(diag(var.Sigma))
    simulate(var, 1, T; stochastic=false, irf=true, e0=e1)[1, :, :]
end

function ErgodDist(var::VAR1, N::Int, T::Int)

    # Simulate for T period, N times
    #  #    - simulate a VAR1 process for T periods N times#
    n = size(var.mu, 1)
    # srand(123) # Setting the seed
    d = MvNormal(var.Sigma)
    VAR_process = simulate(var, N, T)
    E = VAR_process[2]
    XN = VAR_process[1]
    # Computing moments
    Mean_sim = mean(squeeze(mean(VAR_process, 3), 2), 2)
    #Computing the std of simulted Processes (Covariance matrix)
    E1_d = (E[:, :, 1] - repeat(mean(E[:, :, 1], 2), inner=[1, T]))
    Sigma_sim = mean(diag(cov(E1_d[:, :])))  # mean of covariances across simulations

    # Autocorrelation matrix
    X_d = zeros(N, T, 2)
    X_d0 = zeros(N, T-1, 2)
    X_d1 = zeros(N, T-1, 2)
    R_sim = zeros(2, 2)

    for ii in 1:2
        X_d[:, :, ii] = (XN[:, :, 2] - repeat(mean(XN[:, :, ii], 2), inner=[1, T]))
        X_d0[:, :, ii] = X_d[:, 1:end-1, ii]
        X_d1[:, :, ii] = (XN[:, 2:end, ii] - repeat(mean(XN[:, 2:end, ii], 2), inner=[1, T-1]))
        R_sim[ii, ii] = mean(diag( cov(X_d0[:, :, ii], X_d1[:, :, ii] )/sqrt.(var(X_d0[:, :, ii]))/sqrt.(var(X_d1[:, :, ii]))))
        if ii == 2
            R_sim[ii-1, ii] = mean(diag(cov(X_d0[:, :, ii-1], X_d1[:, :, ii])/sqrt.(var(X_d0[:, :, ii-1]))/sqrt.(var(X_d1[:, :, ii]))))
            R_sim[ii, ii-1] = mean(diag(cov(X_d0[:, :, ii], X_d1[:, :, ii-1])/sqrt.(var(X_d0[:, :, ii]))/sqrt.(var(X_d1[:, :, ii-1]))))
        end
    end

    return Mean_sim, Sigma_sim, R_sim
end

#### ProductProcess

mutable struct ProductProcess <: AbstractProcess
    process_1::AbstractProcess
    process_2::AbstractProcess
end

ProductProcess(p) = p

function discretize(::Type{DiscreteMarkovProcess}, pp::ProductProcess; opt1=Dict(), opt2=Dict())
    p1 = discretize(DiscreteMarkovProcess, pp.process_1; opt1...)
    p2 = discretize(DiscreteMarkovProcess, pp.process_2; opt2...)
    return MarkovProduct(p1,p2)
end

function discretize(::Type{DiscretizedProcess}, pp::ProductProcess; opt1=Dict(), opt2=Dict())
  p1 = discretize(DiscretizedProcess, pp.process_1; opt1...)
  p2 = discretize(DiscretizedProcess, pp.process_2; opt2...)
  return Product(p1,p2)
end

function discretize(pp::ProductProcess; kwargs...)
    return discretize(DiscreteMarkovProcess, pp; kwargs...)
end

## special processes (could be implemented as types later on)

function DeathProcess(mu::Float64)
    values = [0.0 1.0;]'
    transitions = [(1-mu) mu; 0 1]
    DiscreteMarkovProcess(transitions, values)
end
function PoissonProcess(mu::Float64, K::Int)
    values = (0:K)[:,:]*1.0
    transitions = zeros(K+1, K+1)
    for i=1:K
        transitions[i,i] = 1-mu
        transitions[i,i+1] = mu
    end
    transitions[K+1,K+1] = 1
    DiscreteMarkovProcess(transitions, values)
end
function AgingProcess(mu::Float64, K::Int)
    values = zeros(K+1,2)
    values[:,1] = 0:K
    values[1,2] = 1
    transitions = zeros(K+1, K+1)
    transitions[1,1] = 1
    for i=2:K
        transitions[i,i+1] = (1-mu)
        transitions[i,1] = mu
    end
    transitions[end,1] = 1
    dp = DiscreteMarkovProcess(transitions, values)
    dp.i0 = 2
    dp
end

get_integration_nodes(dprocess::Dolo.AbstractDiscretizedProcess, i::Int)= [(iweight(dprocess,i,j), inode(dprocess,i,j), j) for j in 1:n_inodes(dprocess,i) if iweight(dprocess,i,j)!=0]

get_integration_nodes(::typeof(Point), dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=[(iweight(dprocess,i,j), inode(Point,dprocess,i,j), j) for j in 1:n_inodes(dprocess,i) if iweight(dprocess,i,j)!=0]

# type unstable
#
# get_integration_nodes(dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=Iterators.filter( x -> (x[1]!=0), ((iweight(dprocess,i,j), inode(dprocess,i,j), j) for j in 1:n_inodes(dprocess,i)) )
#
# get_integration_nodes(::typeof(Point), dprocess::Dolo.AbstractDiscretizedProcess, i::Int)=Iterators.filter( x -> (x[1]!=0), ((iweight(dprocess,i,j), inode(Point,dprocess,i,j), j) for j in 1:n_inodes(dprocess,i)) )

# compatibility names
const AR1 = VAR1
const MarkovChain = DiscreteMarkovProcess
const Normal = MvNormal
const GDP = DiscretizedProcess

MarkovChain(;transitions=ones(1,1), values=[range(1,size(transitions,1))...]) = MarkovChain(transitions, values)


text: const N = 10

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
using MPI
using ITensors
using ITensorParallel
using Random

MPI.Init()

ITensors.BLAS.set_num_threads(1)
ITensors.Strided.disable_threads()
ITensors.disable_threaded_blocksparse()

seed = 1234
Random.seed!(seed)

mpi_projmpo = true

Nx, Ny = 8, 1
N = Nx * Ny

sites = siteinds("S=1/2", N; conserve_qns=true)

lattice = square_lattice(Nx, Ny; yperiodic=false)

function heisenberg_hamiltonian(lattice)
  opsum = OpSum()
  for b in lattice
    opsum += 0.5, "S+", b.s1, "S-", b.s2
    opsum += 0.5, "S-", b.s1, "S+", b.s2
    opsum += "Sz", b.s1, "Sz", b.s2
  end
  return opsum
end
opsum = heisenberg_hamiltonian(lattice)

# Define how the Hamiltonian will be partitioned
# into a sum of sub-Hamiltonians.
function in_partition(sites::Tuple{Int,Int}, p, nparts)
  i, j = sites
  return p == mod1(i, nparts)
end

# Number of partitions
nparts = 4

nprocs = MPI.Comm_size(MPI.COMM_WORLD)
nparts_per_proc = nparts ÷ nprocs
opsums = collect(
  Iterators.partition(partition(opsum, nparts; in_partition=in_partition), nparts_per_proc)
)

PH = if mpi_projmpo
  n = MPI.Comm_rank(MPI.COMM_WORLD) + 1
  MPISum(ThreadedProjMPOSum([MPO(opsum, sites) for opsum in opsums[n]]))
else
  ProjMPOSum(H)
end

state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
psi0 = randomMPS(sites, state, 10)

sweeps = Sweeps(10)
setmaxdim!(sweeps, 20, 60, 100, 100, 200, 400, 800)
setcutoff!(sweeps, 1E-10)
@show sweeps

energy, psi = dmrg(PH, psi0, sweeps)
@show energy

MPI.Finalize()


text:
"""
Tim Holy
High performance computing 1

"""

value = add2(rand(5, 5))

matrix = add2([rand(5, 5), rand(5, 5)])
value2=add2(matrix)
add2([1.0, 2.0])
methods(add2)
"""
doing different things when passed different types of input

"""

m = @which add2(arr = [1.0, 2.0])

m

"""package

method analysis
there's a function in it  called method instance (using MethodAnalysis)
"""
using MethodAnalysis
#pipe , tim ,Vector{Core.MethodInstance}
methodinstances(@which add2([1, 2]))
"""
 MethodInstance for add2(::Matrix{Float64})
 MethodInstance for add2(::Vector{Matrix{Float64}})
 MethodInstance for add2(::Vector{Int64})
hypothesize: these method instances (of the method)
there's only one method
but julia itself would generate multiple instances
(of that method) - for handling different argument types
we get two MethodInstances from a single Method
form of specialization
julia does Automatically

you can
can specialize methods
 can create methods
 append(me:update) methods
 specific types
 have multiple methods  (for a single function)

 automatic specialization
 generatie by the compiler  - compiler specialization

"""
"""
passing a novel Type

let's make another , indexable type

    we should generate another method instance

julia when adds an int & Float64
it will return a Float64
"""

dataType1 = typeof(add2((1, 2, 0)))
dataType2 = typeof(add2((1, 2.0)))

"""
let's call this MethodInstances again
"""
#Base.method_instances()
methodinstances()

@code

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
#module ModiaBase

using .Instantiation
using .ModelElaboration

using Unitful
using .StructuralTransform
using .Synchronous: sample, Clock, previous, hold, positive, positiveChange, positiveEdge

Var(; args...) = Variable(; args...)

Float(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Float64, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

Boolean(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Bool, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

Integ(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Int, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

Str(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, String, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

#=
Float(; args...) = Var(T=Float64; args...)
Float0(; args...) = Float(size=(); args...)

Boolean(; args...) = Var(T=Bool; args...)
Boolean0(; args...) = Boolean(size=(); args...)

Integ(; args...) = Var(T=Int; args...)
Integ0(; args...) = Integ(size=(); args...)

Str(; args...) = Var(T=String; args...)
Str0(; args...) = Str(size=(); args...)
=#

Parameter(; args...) = Variable(variability=parameter; args...)
Parameter(value; args...) = Variable(variability=parameter, value=value; args...)

Par(; args...) = Variable(variability=parameter; args...)
Par(value; args...) = Variable(variability=parameter, value=value; args...)

const undefined=nothing

#end


text:
using FourVectors
using FourVectors.StaticArrays
using Test

@test FourVector(1,2,3; t=4) == FourVector(MVector(1,2,3,4.0))
@test FourVector(1,2,3; t=4) == Particle(p=[1,2,3], E=4.0)

@testset "creation - $Particle" for Particle in [Particle, SParticle]
    @test Particle(0,0,3; msq=4^2) == Particle(p=[0,0,3], E=5.0)

    @test_throws ArgumentError Particle(1,2,3)
    @test_throws ArgumentError Particle(1,2,3; msq=4^2, E=3.3)

    @test mass(Particle(1,2,3; msq=4)) ≈ 2

    @test collect(FourVector(1,2,3; t=4)) == [1,2,3,4]

    x = Particle(1.1,2.2,3.0; msq=4)

    newp = [1.1,1.1,1.1]
    if x isa FourVector{Float64, <:MVector}
        x.p .= newp
    else
        x = Particle(1.1, 1.1, 1.1; msq=4)
    end
    @test invmasssq(x) == x.E^2 - sum(abs2,newp)

    @test x.E == x.T == x[0] == x[4]
    @test x.Px == x.X == x.x == x.px
    @test x.Py == x.Y == x.y == x.py
    @test x.Pz == x.Z == x.z == x.pz
    @test psq(x) == sum(abs2, x[1:3])
    @test_throws ArgumentError x.SomethingElse

    @test length(x) == 4
    @test mass(x * 1000) ≈ 1000*mass(x)

    @test sphericalangles(Particle(0,4,3; msq=4^2)) == (3/5,π/2)
    @test boostfactor(Particle(4,0,0; msq=3^2)) == 5/3

    if x isa FourVector{Float64, <:MVector}
        x1,x2 = x[[1,2]]
        x[[1,2]] .= x[[2,1]]
        @test x[[1,2]] == [x2,x1]
    end
end


text:
using LinearAlgebra: normalize, cross

"""
    turtle

A graphics turtle meant to move in all 3-Dimensions. It contains the following fields:

* `pos` : The coordinates of the turtle.
* `direction` : The Roll Axis - direction in which the turtle is facing.
* `normal` : The Yaw axis - direction normal to the turtle's body.
* `penBlock` : The block which will be used to draw.
* `stepSize` : The size of a turtle's step.
* `penDown` : Activation state of `penBlock`.

The default constructor initializes the `turtle` at the player's position,roll axis is the positive x-direction,
with the yaw axis pointing downwards. `penBlock` is gold by default and is activated with the `stepSize` = 0.5
"""
mutable struct turtle
    pos::Tuple{Float64, Float64, Float64}
    direction::Array{Float64, 1} # Direction turtle is facing, Roll axis
    normal::Array{Float64, 1} # Direction normal to the turtle's body, Yaw axis
    penBlock::PiCraft.Block
    penWidth::Float64
    penScale::Float64
    penDown::Bool

    turtle(;pos = getPos(), direction = [1.0; 0.0; 0.0], normal = [0.0; -1.0; 0.0], penBlock = Block(41), penWidth = 1.0, penScale = 1.0, penDown = true) = new(pos, direction, normal, penBlock, penWidth, penScale, penDown)
end

"""
    move(t::turtle, s::Int)

Move the turtle `t` forward `s` units.
"""
function move(t::turtle, s::Real)
    t.penDown && drawLine(t.pos, Tuple(t.pos .+ s*t.direction), t.penBlock, t.penWidth)
    t.pos = Tuple(t.pos .+ s*t.direction)
end

"""
    pitch(t::turtle, θ::Real)

Pitch the turtle `t` by `θ` degrees.
"""
function pitch(t::turtle, θ::Real)
    yawAxis = cross(t.normal, t.direction)
    t.normal, t.direction = normalize(t.normal*cosd(θ) + cross(yawAxis, t.normal)*sind(θ)), normalize(t.direction*cosd(θ) + cross(yawAxis, t.direction)*sind(θ))
end

"""
    yaw(t::turtle, θ::Real)

Yaw the turtle `t` by `θ` degrees.
"""
function yaw(t::turtle, θ::Real)
    t.direction = normalize(t.direction*cosd(θ) + cross(t.normal, t.direction)*sind(θ))
end

"""
    roll(t::turtle, θ::Real)

Roll the turtle `t` by `θ` degrees.
"""
function roll(t::turtle, θ::Real)
    t.normal = normalize(t.normal*cosd(θ) + cross(t.direction, t.normal)*sind(θ))
end


text:
module ELF

using StructIO

# Bring in ObjectFile definitions
using ObjectFile
importall ObjectFile

# Load in imported C #define constants
include("constants.jl")

# Start to bring in concrete types, in the order they're needed
include("ELFInternal.jl")
include("ELFHeader.jl")
include("ELFHandle.jl")
include("ELFSection.jl")
include("ELFStrTab.jl")
include("ELFSegment.jl")
include("ELFSymbol.jl")
include("ELFDynEntry.jl")
include("ELFDynamicLink.jl")

# Holdovers from ObjFileBase that I haven't cleaned up yet
#include("ELFRelocation.jl")
#include("ELFDebug.jl")

end # module ELF


text:
include("abstract.jl")

type MITgcmDatas <: DataFiles
  Dir :: ASCIIString
  nx :: Int64
  ny :: Int64
  nz :: Int64
  dxspacing :: Float64
  dyspacing :: Float64
  dt :: Float64
  files :: Array{ASCIIString,1}
end

function MITgcmDatas(Dir; nx=200,ny=200, nz=60,
  dxspacing = 100.0, dyspacing = 100.0, dt = 30.0)
  files = map(ASCIIString,readdir(Dir))
  return MITgcmDatas(Dir, nx,ny,nz, dxspacing, dyspacing, dt,files)
end

export MITgcmDatas


text:
###### update ####################

function update{T,N,D}(ta::TimeArray{T,N,D}, tstamp::D, val::Vector{T})
    if tstamp < maximum(ta.timestamp)
        error("only appending operations supported")
    else
        t    = vcat(ta.timestamp, tstamp)
        vals = vcat(ta.values, val')
    end
    TimeArray(t, vals, ta.colnames, ta.meta)
end

function update{T,N,D}(ta::TimeArray{T,N,D}, tstamp::D, val::T)
    if tstamp < maximum(ta.timestamp)
        error("only appending operations supported")
    else
        t    = vcat(ta.timestamp, tstamp)
        vals = vcat(ta.values, val)
    end
    TimeArray(t, vals, ta.colnames, ta.meta)
end

###### rename ####################

function rename{T,N,D}(ta::TimeArray{T,N,D}, colnames::Vector)
    TimeArray(ta.timestamp, ta.values, colnames, ta.meta)
end

function rename{T,N,D}(ta::TimeArray{T,N,D}, colnames::String)
    TimeArray(ta.timestamp, ta.values, [colnames], ta.meta)
end


text:
# # Automatic differentiation through spherical harmonic transforms
# This example finds a positive value of $\lambda$ in:
# ```math
# f(r) = \sin[\lambda (k\cdot r)],
# ```
# for some $k,r\in\mathbb{S}^2$ such that $\int_{\mathbb{S}^2} f^2 {\rm\,d}\Omega = 1$.
# We do this by using derivative information through:
# ```math
# \dfrac{\partial f}{\partial \lambda} = (k\cdot r) \cos[\lambda (k\cdot r)].
# ```

using FastTransforms, LinearAlgebra

# The colatitudinal grid (mod $\pi$):
N = 15
θ = (0.5:N-0.5)/N

# The longitudinal grid (mod $\pi$):
M = 2*N-1
φ = (0:M-1)*2/M

# We precompute a spherical harmonic--Fourier plan:
P = plan_sph2fourier(Float64, N)

# And an FFTW Fourier analysis plan on $\mathbb{S}^2$:
PA = plan_sph_analysis(Float64, N, M)

# Our choice of $k$ and angular parametrization of $r$:
k = [2/7, 3/7, 6/7]
r = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]

# Our initial guess for $\lambda$:
λ = 1.0

# Then we run Newton iteration and grab an espresso:
for _ in 1:7
    F = [sin(λ*(k⋅r(θ,φ))) for θ in θ, φ in φ]
    Fλ = [(k⋅r(θ,φ))*cos(λ*(k⋅r(θ,φ))) for θ in θ, φ in φ]
    U = P\(PA*F)
    Uλ = P\(PA*Fλ)
    global λ = λ - (norm(U)^2-1)/(2*sum(U.*Uλ))
    println("λ: $(rpad(λ, 18)) and the 2-norm: $(rpad(norm(U), 18))")
end


text:
#
# Copyright (c) 2020, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above
#      copyright notice, this list of conditions and the following
#      disclaimer in the documentation and/or other materials provided
#      with the distribution.
# * Neither the name of Intel Corporation nor the names of its
#      contributors may be used to endorse or promote products
#      derived from this software without specific prior written
#      permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

#*******************************************************************
#
# NAME:    nstream
#
# PURPOSE: To compute memory bandwidth when adding a vector of a given
#          number of double precision values to the scalar multiple of
#          another vector of the same length, and storing the result in
#          a third vector.
#
# USAGE:   The program takes as input the number
#          of iterations to loop over the triad vectors, the length of the
#          vectors, and the offset between vectors
#
#          <progname> <# iterations> <vector length> <offset>
#
#          The output consists of diagnostics to make sure the
#          algorithm worked, and of timing statistics.
#
# NOTES:   Bandwidth is determined as the number of words read, plus the
#          number of words written, times the size of the words, divided
#          by the execution time. For a vector length of N, the total
#          number of words read and written is 4*N*sizeof(double).
#
#
# HISTORY: This code is loosely based on the Stream benchmark by John
#          McCalpin, but does not follow all the Stream rules. Hence,
#          reported results should not be associated with Stream in
#          external publications
#
#          Converted to Python by Jeff Hammond, October 2017.
#
# *******************************************************************

# ********************************************************************
# read and test input parameters
# ********************************************************************

function do_initialize(A, B, C, N)
    for i in 1:N
        @inbounds A[i] = 0.0
        @inbounds B[i] = 2.0
        @inbounds C[i] = 2.0
    end
end

function do_nstream(A, B, C, s, N)
    for i in 1:N
        @inbounds A[i] += B[i] + s * C[i]
    end
end

function do_norm(A, N)
    asum = 0.0
    for i in 1:N
        @inbounds asum += abs(A[i])
    end
    return asum
end

function main()
    println("Parallel Research Kernels version")
    println("Julia STREAM triad: A = B + scalar * C")

    if length(ARGS) != 2
        println("argument count = ", length(ARGS))
        println("Usage: ./nstream <# iterations> <vector length>")
        exit(1)
    end

    argv = map(x->parse(Int64,x),ARGS)

    # iterations
    iterations = argv[1]
    if iterations < 1
        println("ERROR: iterations must be >= 1")
        exit(2)
    end

    # vector length
    vlength = argv[2]
    if vlength < 1
        println("ERROR: length must be >= 1")
        exit(3)
    end

    println("Number of iterations     = ", iterations)
    println("Vector length            = ", vlength)

    # ********************************************************************
    # ** Allocate space for the input and transpose matrix
    # ********************************************************************

    A = zeros(Float64,vlength)
    B = zeros(Float64,vlength)
    C = zeros(Float64,vlength)
    precompile(do_initialize, (Array{Float64,1}, Int64))
    do_initialize(A, B, C, vlength)

    # precompile hot functions to smooth performance measurement
    precompile(do_nstream, (Array{Float64,1}, Array{Float64,1}, Array{Float64,1}, Float64, Int64))

    scalar = 3.0

    t0 = time_ns()

    for k in 0:iterations
        if k==0
            t0 = time_ns()
        end
        do_nstream(A, B, C, scalar, vlength)
    end

    t1 = time_ns()
    nstream_time = (t1 - t0) * 1.e-9

    # ********************************************************************
    # ** Analyze and output results.
    # ********************************************************************

    ar = 0.0
    br = 2.0
    cr = 2.0
    for k in 0:iterations
        ar += br + scalar * cr
    end

    ar *= vlength

    precompile(do_norm, (Array{Float64,1}, Int64))
    asum = do_norm(A, vlength)

    epsilon = 1.e-8
    if abs(ar-asum)/asum < epsilon
        println("Solution validates")
        avgtime = nstream_time/iterations
        nbytes = 4.0 * vlength * sizeof(Float64)
        println("Rate (MB/s): ",1.e-6*nbytes/avgtime, " Avg time (s): ", avgtime)
    else
        println("Failed Validation on output array");
        println("        Expected checksum: ",ar);
        println("        Observed checksum: ",asum);
        println("ERROR: solution did not validate")
        exit(1)
    end
end

main()

max_stars_repo_name: hattom/Kernels

max_issues_repo_name: hattom/Kernels

max_forks_repo_name: hattom/Kernels


text:
using DashBootstrapComponents

nav = html_div([
    dbc_nav([
        dbc_navlink("Internal link", href = "/l/components/nav"),
        dbc_navlink("External link", href = "https://github.com"),
        dbc_navlink("External relative", href = "/l/components/nav", external_link = true),
        dbc_navlink("Button", id = "button-link", n_clicks = 0),
    ]),
    html_br(),
    html_p(id = "button-clicks"),
]);

callback!(app, Output("button-clicks", "children"), Input("button-link", "n_clicks")) do n
    return "Button clicked $n times"
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using LinearAlgebra: norm
using Random: MersenneTwister
using LoopVectorization: @avxt
using JuMP
using OSQP

export generatePerturbation, MCMCTrajectory, MinimumHoleSize, detectTransitions, minimumdistance, FormNetwork, GenerateConstraints

defaultRNG = MersenneTwister()

"""
    generatePerturbation([ rng,] samples, σ, windowlength)

Generates a random perturbation trajectory.

The trajectory is made of a random gaussian vector, filtered
through convolution with normalized Hann window of a given length.

### Arguments
 - `rng`          -- Random number generator (Optional)
 - `samples`      -- The number of timesteps for the perturbation
 - `σ`            -- The standard deviation of any given Gaussian sample
 - `windowlength` -- The length of the Hann window

 ### Returns
 A filtered Gaussian vector of length `samples`

 ### Notes
 Performs the convolution directly, assumes ``windowlength << samples``
"""
function generatePerturbation(rng, samples, σ, windowlength)

    perturbation = σ * randn(rng, samples + windowlength)
    window = sin.((1:windowlength)/(windowlength+1) * π).^2

    # Rescale the window to preserve variance
    window ./= norm(window)

    return [sum(window .* perturbation[i:i+windowlength-1])
                                                    for i in 1:samples]
end

export generatePerturbation_v

"""
    generatePerturbation_v([ rng,] samples, σ, windowlength)

Generates a random perturbation trajectory.

The trajectory is made of a random gaussian vector, filtered
through convolution with normalized Hann window of a given length.
Equivalent to generatePerturbation but uses LoopVectorization.

### Arguments
 - `rng`          -- Random number generator (Optional)
 - `samples`      -- The number of timesteps for the perturbation
 - `σ`            -- The standard deviation of any given Gaussian sample
 - `windowlength` -- The length of the Hann window

 ### Returns
 A filtered Gaussian vector of length `samples`

 ### Notes
 Performs the convolution directly, assumes ``windowlength << samples``
"""
function generatePerturbation_v(rng,samples, σ, windowlength)

    perturbation = σ * randn(rng, samples + windowlength)
    window = sin.((1:windowlength)/(windowlength+1) * π).^2

    # Rescale the window to preserve variance
    window ./= norm(window)

    # Allocate Memory
    out = zeros(samples)
    @avxt for i in 1:samples
        for j in 1:windowlength
            out[i] += window[j] * perturbation[i+j-1]
        end
    end

    return out
end

generatePerturbation(samples, σ, windowlength) =
        generatePerturbation(defaultRNG, samples, σ, windowlength)

generatePerturbation_v(samples, σ, windowlength) =
        generatePerturbation_v(defaultRNG, samples, σ, windowlength)

"""
    MCMCTrajectory([rng,] fₓ, x₀, α, stepsize, samples, σ, windowlength)

Generates a trajectory through the state × perturbation space based
on a noisy gradient descent process using the provided gradient.

The trajectory evolves according to ``x_{t+1} = x_t + α(uₜ - fₓ(x_t))``,
where ``uₜ`` represents the random perturbation.

### Arguments
 - `rng`          -- Random number generator (Optional)
 - `fₓ`           -- Gradient of the function of interest
 - `x₀`           -- Initial Point
 - `α`            -- Gradient Descent step size
 - `samples`      -- Number of timesteps for the trajectory
 - `σ`            -- Standard deviation of the Gaussian perturbation
 - `windowlength` -- Length of the Hann window for the perturbation process

### Returns
A trajectory in the state × perturbation space with states in the first
half of the rows, and perturbations in the second half.
Additionally, returns the gradients along the path.

### Notes
Performs the convolution directly, assumes ``windowlength << samples``
"""
function MCMCTrajectory(rng, fₓ, x₀, α, samples, σ, windowlength)
    N = length(x₀)

    # Initialize the state vector
    state = zeros(2*length(x₀),samples)
    state[1:N,1]     .= x₀
    @inbounds for i in 1:N
        state[N+i,:] .= generatePerturbation_v(rng, samples,σ,windowlength)
    end
    gradients = zeros(N, samples)

    diff = zeros(N)
    @inbounds for i in 2:samples
        xₜ = @view state[1:N,i-1]
        uₜ = @view state[N+1:end,i]
        g  = @view gradients[:,i-1]
        g .= fₓ(xₜ)
        diff .= uₜ .- g
        state[1:N,i] .= xₜ .+  α .* diff
    end
    g = @view gradients[:,end]
    g .= fₓ(state[1:N,end])

    return (state, gradients)
end

MCMCTrajectory(fₓ, x₀, α, samples, σ, windowlength) =
        MCMCTrajectory(defaultRNG, fₓ, x₀, α, samples, σ, windowlength)

"""
    resampleSpatialTrajectory(state, δ)

Resamples the trajectory by removing steps closer than δ in parameters.

### Arguments
 - `state`  --  A matrix where columns represent N samples in the
                M-dimensional parameter-control space
 - `δ`      --  Minimum step distance

### Returns
A vector of N booleans with false if the sample should be removed.
"""
function resampleSpatialTrajectory(state, δ)
    M, N = size(state)
    param_range = 1:Int(M/2)
    keep = ones(Bool,N)
    dist = 0
    for i in 2:N
        dist += norm(state[param_range,i] - state[param_range,i-1])
        if  dist < δ
            keep[i] = false
        else
            dist = 0
        end
    end
    return keep
end

"""
    GenerateConstraints(gradients, directions)

Generates the constraints used to compute the minimum spacing between
the forward path and reverse path in the augmented space.

### Arguments
 - `gradients`  -- Gradients along the trajectory
 - `directions` -- normalized step direction along the trajectory

### Returns
    (A1, A2, b1, b2)

Two matrices `A1, A2` and two vector `b1, b2` such that min_{u,v} ||u - v||
s.t. A1 u ≥ b1, A2 v ≥ b2 results in the minimum path discrepancy using the
threshold approximation.
"""
function GenerateConstraints(gradients,directions)

    # Split into two halves
    M, N = size(gradients)
    mid = Int(ceil(N/2))

    A1 = directions[:,1:mid-1]
    A2 = -directions[:,mid+1:end]

    b1 = zeros(mid-1)
    b2 = zeros(mid-1)
    @inbounds @simd for i in 1:mid-1
        b1[i] = directions[:,i]' * gradients[:,i]
        b2[i] = -directions[:,mid+i]' * gradients[:,mid+i]
    end

    return (A1', A2', b1, b2)
end

"""
    SolveMinProb(A1, A2, b1, b2, λ[, prune = true])

Helper function solving min_{u,v} ||u - v||^2 + λ(||u||^2 + ||v||^2)
s.t. A1 u ≥ b1, A2 v ≥ b2 with optional pruning based on requiring the
path to not curve too much. Uses JuMP and OSQP for now, could be changed later.

### Arguments
 - `A1, A2, b1, b2, λ`  -- The parameters of optimization problem
 - `prune`              -- Whether or not to prune the constraints

### Returns
||u-v|| for the minimizing `u,v`
"""
function SolveMinProb(A1, A2, b1, b2, λ, prune)
    M1, N1 = size(A1)
    M2, N2 = size(A2)

    # Prune
    if prune
        for i in reverse(1:M1-1)
            if A1[i,:]' * A1[end,:] < .5
                A1 = A1[i+1:end,:]
                b1 = b1[i+1:end]
                break
            end
        end
        for i in 2:M2
            if A2[i,:]' * A2[1,:] < .5
                A2 = A2[1:i-1,:]
                b2 = b2[1:i-1]
                break
            end
        end
    end
    M1, N1 = size(A1)
    M2, N2 = size(A2)

    # Flag for monotonic --> no constraints = zero distance
    if M1 == 0 || M2 == 0
        return 0
    end

    model = Model(OSQP.Optimizer)
    set_silent(model)
    @variable(model, x[1:N1])
    @variable(model, y[1:N2])
    @constraint(model, A1 * x .>= b1)
    @constraint(model, A2 * y .>= b2)
    @objective(model, Min, x' * x - 2 * y' * x + y' * y + λ*(x' * x + y' * y))
    optimize!(model)

    x_v = value.(x)
    y_v = value.(y)

    v_comp = x_v' * x_v - 2 * y_v' * x_v + y_v' * y_v

    if isinf(v_comp) || isnan(v_comp) || v_comp < 0
        return 0
    else
        return sqrt(v_comp)
    end
end

SolveMinProb(A1, A2, b1, b2, λ) = SolveMinProb(A1, A2, b1, b2, λ, true)

"""
    MinimumHoleSize(gradients, directions, λ)

Given a matrix representing the gradients along a path,
returns the minimum size of the geometric hole for trajectories
along that path.

### Arguments
 - `gradients`  -- M×N matrix of N timesteps of gradients
 - `directions` -- M×N matrix of N timesteps of step directions
 - `λ`          -- ℓ2 regularization perameter

### Returns
The list of the gradients
"""
function MinimumHoleSize(gradients, directions, λ)
    A1, A2, b1, b2 = GenerateConstraints(gradients, directions)
    return SolveMinProb(A1, A2, b1, b2, λ)
end

"""
    detectTransitions(discrepancies, threshold, pathlength)

Detects the transitions by applying a simple threshold
"""
function detectTransitions(discrepancies, threshold, pathlength)
    unstable = Int.(findall(discrepancies .> threshold) .+ floor(pathlength/2))
    stable = Int.(findall(discrepancies .<= threshold) .+ floor(pathlength/2))
    return unstable, stable
end

"""
    minimumdistance(A,B)

Returns the minimum distance between two sets of points.
Each column corresponds to a point.
"""
function minimumdistance(A,B)
    M, N1 = size(A)
    ~, N2 = size(B)

    min_dist = Inf * ones(N1)

    Threads.@threads for i in 1:N1
        tmp = zeros(M)
        @inbounds for j in 1:N2
            d = euclideandist_fast!(tmp,A,B,i,j)
            min_dist[i] = min_dist[i]<d ? min_dist[i] : d
        end
    end

    return sqrt(minimum(min_dist[1:N1]))
end

"""
    euclideandist_fast!(tmp,A,B,i,j)

Optimized routine for Euclidean distance between two matrix columns
"""
function euclideandist_fast!(tmp,A,B,i,j)
    tmp .= view(A,:,i) .- view(B,:,j)
    return sum(abs2,tmp)
end

"""
    FormNetwork(dataset, transitions, threshold)

Takes a dataset of positions in the parameter space and the transition labels,
generates the network structure.

### Arguments
 - `dataset`        -- The trajectory points
 - `transitions`    -- Vector of indicies categorized as transitions
 - `threshold`      -- Threshold for merging two clusters

### Returns
    (edges, means)

 - `edges`  -- Array of tuples representing connections between nodes
 - `means`  -- Average of points labeled for each node
"""
function FormNetwork(dataset, transitions, threshold)

    # Begin by splitting the dataset along the transitions
    # These form the initial clusters.
    clusterindices = [Array{Int}(1:transitions[1]-1)]
    last_ind = transitions[1]
    for i in 2:length(transitions)
        if transitions[i] > last_ind+2
            push!(clusterindices, last_ind+1:transitions[i]-1)
        end
        last_ind = transitions[i]
    end
    if transitions[end] < size(dataset)[2]
        push!(clusterindices, transitions[end]+1:size(dataset)[2])
    end

    # Next, cluster the clusters by a threshold, recording indices
    clusters_included = [[i] for i in 1:length(clusterindices)]
    ind1 = 1
    while ind1 < length(clusterindices)
        ind2 = ind1+1
        no_merge = true
        while ind2 <= length(clusterindices)
            cluster1 = @view dataset[:,clusterindices[ind1]]
            cluster2 = @view dataset[:,clusterindices[ind2]]
            dist = minimumdistance(cluster2,cluster1)

            # If clusters are close enough, merge
            if dist < threshold
                append!(clusterindices[ind1], clusterindices[ind2])
                append!(clusters_included[ind1],clusters_included[ind2])
                deleteat!(clusterindices,ind2)
                deleteat!(clusters_included,ind2)
                no_merge = false
            else
                ind2 += 1
            end
        end
        if no_merge
            ind1 += 1
        end
    end

    # Using sequential connections, form the graph (assuming sparsity)
    edges = []
    for i in 1:length(clusters_included)
        for j in i+1:length(clusters_included)
            # Check for inclusion of sequential initial clusters
            offsets = [i1 - i2  for i1 in clusters_included[i]
                                for i2 in clusters_included[j]]
            connection = findfirst((x)-> (abs(x)==1), offsets)
            if !isnothing(connection)
                push!(edges, [i,j])
            end
        end
    end

    # Finally, as a coarse approximation, compute the means of the clusters
    means = zeros(size(dataset)[1],length(clusterindices))
    for i in 1:length(clusterindices)
        means[:,i] = sum(dataset[:,clusterindices[i]], dims=2)/length(clusterindices[i])
    end

    return (edges, means)
end


text:
#!/usr/bin/env

using Luxor

using Test

using Random
Random.seed!(42)

function test_circle_intersections(fname)
    Drawing(2000, 2000, fname)
    origin()
    setopacity(0.995)
    setline(0.1)
    tiles = Tiler(2000, 2000, 5, 5)
    setmode("darken")
    for (pos, n) in tiles
            @layer begin
                translate(pos)
                circle1 = (O, rand(30:tiles.tilewidth/3))
                circle2 = (O + (rand(40:70), rand(-20:20)), rand(40:tiles.tilewidth/3))
                flag, C, D = intersectioncirclecircle(circle1[1], circle1[2], circle2[1], circle2[2])
                if flag
                    sethue("slategray3")
                    circle(circle1..., :fill)
                    sethue("orange3")
                    circle(circle2..., :fill)
                    sethue("black")
                    circle(circle1..., :stroke)
                    circle(circle2..., :stroke)
                    sethue("red")
                    circle.([C, D], 4, :fill)
                else
                    sethue("black")
                    circle.([circle1[1], circle2[1]], [circle1[2], circle2[2]], :stroke)
                end
            end
    end
    @test finish() == true
end

fname = "intersection_circle_circle_test.png"
test_circle_intersections(fname)
println("...finished test: output in $(fname)")


text:
@with_kw struct DetectAndAvoidMDP <: MDP{Array{Float32},Float32}
    ddh_max::Float64 = 1.0 # vertical acceleration limit [m/s²]
    collision_threshold::Float64 = 50.0 # collision threshold [m]
    reward_collision::Float64 = -100.0 # reward obtained if collision occurs
    reward_change::Float64 = -1 # reward obtained if action changes
    actions::Vector{Float64} = [-5.0, 0.0, 5.0] # Available actions to command for vertical rate [m/s]
    px = DiscreteNonParametric([2.0, 0.0, -2.0], [0.25, 0.5, 0.25])

    h0_dist = Distributions.Uniform(-200, 200)
    dh0_dist = Distributions.Uniform(-10, 10)
    a_prev0 = 0.0
    τ0 = 40
end

function POMDPs.gen(mdp::DetectAndAvoidMDP, s, a, x, rng::AbstractRNG=Random.GLOBAL_RNG)
    t = transition(mdp, s, a, x)
    (sp=rand(t), r=reward(mdp, s, a))
end

function POMDPs.transition(mdp::DetectAndAvoidMDP, s, a, x, rng::AbstractRNG=Random.GLOBAL_RNG)
    a = x == 0 ? 0.0 : a # COC if don't detect

    h, dh, a_prev, τ = s

    # Update the dynamics
    h = h + dh
    #if a != 0.0
    if abs(a - dh) < mdp.ddh_max
        dh += a - dh
    else
        dh += sign(a - dh) * mdp.ddh_max
    end
    #end
    a_prev = a
    τ = max(τ - 1.0, -1.0)
    SparseCat([Float32[h, dh+x, a_prev, τ] for x in mdp.px.support], mdp.px.p)
end

function POMDPs.reward(mdp::DetectAndAvoidMDP, s, a)
    h, dh, a_prev, τ = s

    r = 0.0
    if isfailure(mdp, s)
        # We collided
        r += mdp.reward_collision
    end
    if a != a_prev
        # We changed our action
        r += mdp.reward_change
    end
    r
end

POMDPs.convert_s(::Type{Array{Float32}}, v::V where {V<:AbstractVector{Float64}}, ::DetectAndAvoidMDP) = Float32.(v)
POMDPs.convert_s(::Type{V} where {V<:AbstractVector{Float64}}, s::Array{Float32}, ::DetectAndAvoidMDP) = Float64.(s)

function POMDPs.initialstate(mdp::DetectAndAvoidMDP)
    ImplicitDistribution((rng) -> Float32[rand(mdp.h0_dist), rand(mdp.dh0_dist), mdp.a_prev0, mdp.τ0])
end

POMDPs.actions(mdp::DetectAndAvoidMDP) = mdp.actions
POMDPs.actionindex(mdp::DetectAndAvoidMDP, a) = findfirst(mdp.actions .== a)

disturbanceindex(mdp::DetectAndAvoidMDP, x) = findfirst(mdp.px.support .== x)
disturbances(mdp::DetectAndAvoidMDP) = mdp.px.support

function isfailure(mdp::DetectAndAvoidMDP, s)
    h, dh, a_prev, τ = s
    abs(h) < mdp.collision_threshold && abs(τ) < eps()
end

function POMDPs.isterminal(mdp::DetectAndAvoidMDP, s)
    h, dh, a_prev, τ = s
    τ < 0.0
end

POMDPs.discount(mdp::DetectAndAvoidMDP) = 0.99

## Here is a solver that gives the optimal policy
struct OptimalDetectAndAvoidPolicy <: Policy
    𝒜
    grid
    Q
end

function OptimalDetectAndAvoidPolicy(mdp::DetectAndAvoidMDP, hs=range(-200, 200, length=21), dhs=range(-10, 10, length=21), τs=range(0, 40, length=41))
    grid = RectangleGrid(hs, dhs, actions(mdp), τs)

    𝒮 = [[h, dh, a_prev, τ] for h in hs, dh in dhs, a_prev in actions(mdp), τ in τs]

    # State value function
    U = zeros(length(𝒮))

    # State-action value function
    Q = [zeros(length(𝒮)) for a in actions(mdp)]

    # Solve with backwards induction value iteration
    for (si, s) in enumerate(𝒮)
        for (ai, a) in enumerate(actions(mdp))
            Tsa = transition(mdp, s, a, 1.0)
            Q[ai][si] = reward(mdp, s, a)
            Q[ai][si] += sum(isterminal(mdp, s′) ? 0.0 : Tsa.probs[j] * GridInterpolations.interpolate(grid, U, vec(s′)) for (j, s′) in enumerate(Tsa.vals))
        end
        U[si] = maximum(q[si] for q in Q)
    end
    return OptimalDetectAndAvoidPolicy(actions(mdp), grid, Q)
end

function POMDPs.action(policy::OptimalDetectAndAvoidPolicy, s)
    a_best = first(policy.𝒜)
    q_best = -Inf
    for (a, q) in zip(policy.𝒜, policy.Q)
        q = GridInterpolations.interpolate(policy.grid, q, s)
        if q > q_best
            a_best, q_best = a, q
        end
    end
    return a_best
end


text:
module AWSMetadata

using ..AWSMetadataUtilities: _get_aws_sdk_js_files, _get_service_and_version,
    _generate_low_level_definitions, _generate_high_level_definitions
using Base64
using GitHub
using HTTP
using JSON
using OrderedCollections: OrderedDict

const services_path = joinpath(@__DIR__, "AWSServices.jl")

"""
    parse_aws_metadata()

Generate low and high level wrappers for each AWS Service based on their definitions in the
[aws-sdk-js GitHub Repository](https://github.com/aws/aws-sdk-js/tree/master/apis).

Low level wrappers are written into `src/AWSServices.jl`, while high level wrappers API
wrappers are written into their respective files in `src/services/{service}.jl`.
"""
function parse_aws_metadata()
    auth = GitHub.authenticate(ENV["GITHUB_AUTH"])
    repo_name = "aws/aws-sdk-js"

    files = _get_aws_sdk_js_files(repo_name, auth)

    _generate_low_level_wrappers(files, repo_name, auth)
    _generate_high_level_wrapper(files, repo_name, auth)
end

function _generate_low_level_wrappers(services::AbstractArray, repo_name::String, auth::GitHub.OAuth2)
    service_definitions = sort(_generate_low_level_definitions(services, repo_name, auth))

    template = """
    # This file is auto-generated by AWSMetadata.jl

    module AWSServices

    using AWS
    using OrderedCollections: LittleDict

    $(join(service_definitions, "\n"))

    end
    """

    open(services_path, "w") do f
        print(f, template)
    end

    return template
end

function _generate_high_level_wrapper(services::AbstractArray{<:AbstractDict}, repo_name::String, auth::GitHub.OAuth2)
    for service in services
        service_name = service["path"]
        @info "Generating high level wrapper for $service_name"

        service_blob = blob(repo_name, service["sha"]; auth=auth)
        service = JSON.parse(String(base64decode(service_blob.content)))

        service_name = lowercase(service["metadata"]["serviceId"])
        service_name = replace(service_name, ' ' => '_')
        operations = service["operations"]
        shapes = service["shapes"]

        protocol = service["metadata"]["protocol"]

        operations = sort!(_generate_high_level_definitions(service_name, protocol, operations, shapes))

        service_path = joinpath(@__DIR__, "services", "$service_name.jl")
        open(service_path, "w") do f
            println(f, """
                    # This file is auto-generated by AWSMetadata.jl
                    using AWS
                    using AWS.AWSServices: $service_name
                    using AWS.Compat
                    using AWS.UUIDs
                    """)
            join(f, operations, "\n")
        end
    end
end

end