text:
@testset "BroadcastStyle" begin
    A = typeof(AxisArray(ones(2,2)))
    B = typeof(ones(2,2))
    @test Base.BroadcastStyle(A) isa AxisIndices.AxisArrayStyle
    SA = Base.BroadcastStyle(A)
    SB = Base.BroadcastStyle(B)
    @test Base.BroadcastStyle(SA, SA) isa AxisIndices.AxisArrayStyle
    @test Base.BroadcastStyle(SA, SB) isa AxisIndices.AxisArrayStyle
    @test Base.BroadcastStyle(SB, SA) isa AxisIndices.AxisArrayStyle
end

@testset "combine" begin
    #=
    @testset "broadcast_axis" begin
        @test @inferred(broadcast_axis(Axis(1:2), Axis(1:2)))       isa Axis{Int64,Int64,UnitRange{Int64},Base.OneTo{Int64}}
        @test @inferred(broadcast_axis(Axis(1:2), SimpleAxis(1:2))) isa Axis{Int64,Int64,UnitRange{Int64},UnitRange{Int64}}
        @test @inferred(broadcast_axis(Axis(1:2), 1:2))             isa Axis{Int64,Int64,UnitRange{Int64},UnitRange{Int64}}

        @test @inferred(broadcast_axis(SimpleAxis(1:2), SimpleAxis(1:2))) isa SimpleAxis{Int,UnitRange{Int}}
        @test @inferred(broadcast_axis(SimpleAxis(1:2), 1:2))             isa SimpleAxis{Int,UnitRange{Int}}
        @test @inferred(broadcast_axis(SimpleAxis(1:2), Axis(1:2)))       isa Axis{Int64,Int64,UnitRange{Int64},UnitRange{Int64}}

        @test @inferred(broadcast_axis(1:2, SimpleAxis(1:2))) isa SimpleAxis{Int,UnitRange{Int}}
        @test @inferred(broadcast_axis(1:2, 1:2))             isa UnitRange{Int}
        @test @inferred(broadcast_axis(1:2, Axis(1:2)))       isa Axis{Int64,Int64,UnitRange{Int64},UnitRange{Int64}}

        @test @inferred(broadcast_axis(1:2, Symbol.(1:2))) isa Vector{Symbol}
        @test @inferred(broadcast_axis(Symbol.(1:2), 1:2)) isa Vector{Symbol}

        @test @inferred(broadcast_axis(1:2, string.(1:2))) isa Vector{String}
        @test @inferred(broadcast_axis(string.(1:2), 1:2)) isa Vector{String}

        #@test_throws ErrorException("No method available for combining keys of type Int64 and String.") combine_keys(12, "x")

        @test @inferred(broadcast_axis(1:2, SimpleAxis(1:2))) isa SimpleAxis{Int,UnitRange{Int}}
    end

    @test @inferred(broadcast_axis(SimpleAxis(1:2), SimpleAxis(1:2))) == SimpleAxis(1:2)
    =#

    @testset "Multiple broadcasts" begin
        @test (AxisArray(ones(2)) .- AxisArray(ones(2))) .^ 2 == zeros(2)
    end

    @test @inferred(Base.Broadcast.broadcast_shape((1:10,), (1:10, 1:10), (1:10,))) == (1:10, 1:10)
    b1 = Broadcast.broadcast_shape(
        axes(CartesianIndices((1,))),
        axes(CartesianIndices((3, 2, 2))),
        axes(CartesianIndices((3, 2, 2)))
    )
    b2 = Broadcast.broadcast_shape(
        axes(CartesianAxes((1,))),
        axes(CartesianIndices((3, 2, 2))),
        axes(CartesianAxes((3, 2, 2)))
    )
    @test b1 == b2

    b1 = Broadcast.broadcast_shape(
        axes(LinearIndices((1,))),
        axes(LinearIndices((3, 2, 2))),
        axes(LinearIndices((3, 2, 2)))
    )
    b2 = Broadcast.broadcast_shape(
        axes(LinearAxes((1,))),
        axes(LinearIndices((3, 2, 2))),
        axes(LinearAxes((3, 2, 2)))
    )
    @test b1 == b2

    @test_throws DimensionMismatch Broadcast.broadcast_shape(axes(LinearAxes((3, 2, 2))), axes(LinearAxes((3, 2, 5))))

    @testset "combine_axes" begin
        @test Broadcast.combine_axes(CartesianIndices((1,)), CartesianIndices((3, 2, 2)), CartesianIndices((3, 2, 2))) ==
                @inferred(Broadcast.combine_axes(CartesianAxes((1,)), CartesianIndices((3, 2, 2)), CartesianAxes((3, 2, 2))))
        @test Broadcast.combine_axes(LinearIndices((1,)), LinearIndices((3, 2, 2)), LinearIndices((3, 2, 2))) ==
                @inferred(Broadcast.combine_axes(LinearAxes((1,)), CartesianAxes((3, 2, 2)), CartesianAxes((3, 2, 2))))

        cartinds = CartesianIndices((2, 2));
        cartaxes = CartesianAxes((2:3, 3:4));
        @test keys.(Broadcast.combine_axes(cartaxes, cartaxes, cartaxes)) == (2:3, 3:4)
        @test keys.(Broadcast.combine_axes(cartaxes, cartaxes, cartinds)) == (2:3, 3:4)
        @test keys.(Broadcast.combine_axes(cartaxes, cartinds, cartaxes)) == (2:3, 3:4)
        @test keys.(Broadcast.combine_axes(cartinds, cartaxes, cartaxes)) == (2:3, 3:4)
    end
end

@testset "Broadcasting" begin
    x = Axis(1:10, 1:10)
    y = SimpleAxis(1:10)
    z = 1:10

    @test @inferred(broadcasted(bstyle, -, 2, x)) ==
          @inferred(broadcasted(bstyle, -, 2, y)) ==
          broadcasted(bstyle, -, 2, z)

    @test @inferred(broadcasted(bstyle, -, x, 2)) ==
          @inferred(broadcasted(bstyle, -, y, 2)) ==
          broadcasted(bstyle, -, z, 2)

    @test @inferred(broadcasted(bstyle, +, 2, x)) ==
          @inferred(broadcasted(bstyle, +, 2, y)) ==
          broadcasted(bstyle, +, 2, z)
    @test isa(broadcasted(bstyle, +, 2, x), Axis)
    @test isa(broadcasted(bstyle, +, 2, y), SimpleAxis)

    @test @inferred(broadcasted(bstyle, +, x, 2)) ==
          @inferred(broadcasted(bstyle, +, y, 2)) ==
          broadcasted(bstyle, +, z, 2)
    @test isa(broadcasted(bstyle, +, x, 2), Axis)
    @test isa(broadcasted(bstyle, +, y, 2), SimpleAxis)

    @test @inferred(broadcasted(bstyle, *, 2, x)) ==
          @inferred(broadcasted(bstyle, *, 2, y)) ==
          broadcasted(bstyle, *, 2, z)

    @test @inferred(broadcasted(bstyle, *, x, 2)) ==
          @inferred(broadcasted(bstyle, *, y, 2)) ==
          broadcasted(bstyle, *, z, 2)

    @test broadcasted(bstyle, -, x, 1//2) == 1//2:19//2
    @test broadcasted(bstyle, -, y, 1//2) == 1//2:19//2
end

@testset "Binary broadcasting operations (.+)" begin

    #= TODO moved to docs
    @testset "standard case" begin
        a = AxisArray(ones(3), (2:4,))
        @test a .+ a == 2ones(3)
        @test keys(axes(a .+ a, 1)) == 2:4

        @test a .+ a.+ a == 3ones(3)
        @test keys(axes(a .+ a .+ a, 1)) == 2:4

        @test (a .= 0 .* a .+ 7) == [7, 7, 7]
    end
    =#

    #= TODO Need to explain how order matters
    @testset "Order matters" begin
        x = AxisArray(ones(3, 5), (:, nothing))
        y = AxisArray(ones(3, 5), (nothing, :y))

        lhs = x .+ y
        rhs = y .+ x
        @test dimnames(lhs) == (:x, :y) == dimnames(rhs)
        @test lhs == 2ones(3, 5) == rhs
    end
    =#

   @testset "broadcasting" begin
        v = AxisArray(zeros(3,), (2:4,))
        m = AxisArray(ones(3, 3), (2:4, 3:5))
        s = 0

        @test @inferred(v .+ m) == ones(3, 3) == @inferred(m .+ v)
        @test @inferred(s .+ m) == ones(3, 3) == @inferred(m .+ s)
        @test @inferred(s .+ v .+ m) == ones(3, 3) == @inferred(m .+ s .+ v)

        @test keys.(axes(v .+ m)) == (2:4, 3:5) == keys.(axes(m .+ v))
        @test keys.(axes(s .+ m)) == (2:4, 3:5) == keys.(axes(m .+ s))
        @test keys.(axes(s .+ v .+ m)) == (2:4, 3:5) == keys.(axes(m .+ s .+ v))
    end
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using LinearAlgebra: qr, cond

function get_V0(train)
   # get list of atomic numbers
   Zs = AtomicNumber[]
   for (at, E) in train
      Zs = unique( [Zs; at.Z] )
   end

   # setup lsq system for E0s
   A = zeros(length(train), length(Zs))
   y = zeros(length(train))
   for (it, (at, E)) in enumerate(train)
      y[it] = E
      for (iz, z) in enumerate(Zs)
         A[it, iz] = length(findall(at.Z .== z))
      end
   end
   E0s = A \ y
   @info("  E0s = $(E0s)")
   syms = chemical_symbol.(Zs)
   return JuLIP.OneBody([syms[i] => E0s[i] for i = 1:length(Zs)]...)
end

function trainset(Vref, Ntrain; kwargs...)
   train = []
   for n = 1:Ntrain
      at = rand(Vref; kwargs...)
      push!(train, (at = at, E = energy(Vref, at), F = forces(Vref, at)))
   end
   # remove E0s
   V0 = get_V0(train)
   train = [ (at = at, E = E - energy(V0, at), F = F)
               for (at, E,  F) in train ]
   return train
end

function lsq(train, basis; verbose=true, wE = 1.0, wF = 1.0)
   @info("lsq info")
   nobs = sum( 1+3*length(t.at) for t in train )
   @info("  nobs = $(nobs); nbasis = $(length(basis))")
   A = zeros(nobs, length(basis))
   y = zeros(nobs)
   irow = 0
   for (at, E, F) in train
      # add energy to the lsq system
      irow += 1
      y[irow] = wE * E / length(at)
      A[irow, :] = wE * energy(basis, at) / length(at)

      # add forces to the lsq system
      nf = 3*length(at)
      y[(irow+1):(irow+nf)] = wF * mat(F)[:]
      Fb = forces(basis, at)
      for ib = 1:length(basis)
         A[(irow+1):(irow+nf), ib] = wF * mat(Fb[ib])[:]
      end
      irow += nf
   end
   qrF = qr(A)
   c = qrF \ y
   relrmse = norm(A * c - y) / norm(y)
   @info("   cond(R) = $(cond(qrF.R)); relrmse = $(relrmse)")
   return c, relrmse
end

ext: jl


text:
include("utilities.jl")
include("mcmc.jl")
include("hmc.jl")
# include("sampling.jl")


text:
module SparseOps

using LinearAlgebra
using SparseArrays
using StaticArrays

using ..Defs

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

export ID
"""
Tagged ID (or index) for various IDVector and SparseOp

This type is a pure ordinal; it cannot be used for counting.
"""
struct ID{Tag}
    id::Int
end

Base.Int(i::ID) = i.id
Base.convert(::Type{Int}, i::ID) = Int(i)

const subscripts = ["₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉"]
Base.show(io::IO, i::ID{Tag}) where {Tag} = print(io, i.id, subscripts[Tag + 1])

Base.:(==)(i::ID{Tag}, j::ID{Tag}) where {Tag} = i.id == j.id
Base.:(==)(::ID, ::ID) = false
Base.:(<)(i::ID{Tag}, j::ID{Tag}) where {Tag} = i.id < j.id
Base.isless(i::ID{Tag}, j::ID{Tag}) where {Tag} = isless(i.id, j.id)
function Base.hash(i::ID{Tag}, h::UInt) where {Tag}
    return hash(0x4e3ca100, hash(Tag, hash(i.id, h)))
end

Base.zero(::Type{ID{Tag}}) where {Tag} = ID{Tag}(0)
Base.:+(i::ID{Tag}, j::Int) where {Tag} = ID{Tag}(i.id + j)
Base.:-(i::ID{Tag}, j::Int) where {Tag} = ID{Tag}(i.id - j)
Base.:-(i::ID{Tag}, j::ID{Tag}) where {Tag} = i.id - j.id

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

export IDVector
struct IDVector{Tag,T}
    vec::Vector{T}
end

IDVector{Tag}(xs::Vector{T}) where {Tag,T} = IDVector{Tag,T}(xs)
IDVector{Tag,T}() where {Tag,T} = IDVector{Tag}(T[])

function Base.show(io::IO, xs::IDVector{Tag}) where {Tag}
    return print(io, xs.vec, subscripts[Tag + 1])
end

Base.:(==)(xs::IDVector{Tag}, ys::IDVector{Tag}) where {Tag} = xs.vec == ys.vec
Base.:(==)(::IDVector, ::IDVector) = false
function Base.isless(xs::IDVector{Tag}, y::IDVector{Tag}) where {Tag}
    return isless(xs.vec, ys.vec)
end
function Base.hash(xs::IDVector{Tag}, h::UInt) where {Tag}
    return hash(0x934f6884, hash(Tag, hash(xs.vec, h)))
end

Base.axes(xs::IDVector{Tag}, d) where {Tag} = ID{Tag}(1):ID{Tag}(length(xs))
Base.axes(xs::IDVector) = (axes(xs, 1),)
Base.eltype(::Type{<:IDVector{<:Any,T}}) where {T} = T
Base.first(xs::IDVector) = first(xs.vec)
Base.@propagate_inbounds function Base.getindex(xs::IDVector{Tag},
                                                i::ID{Tag}) where {Tag}
    return getindex(xs.vec, i.id)
end
Base.@propagate_inbounds function Base.getindex(xs::IDVector{Tag2},
                                                is::IDVector{Tag1,ID{Tag2}}) where {Tag1,
                                                                                Tag2}
    return IDVector{Tag1}([xs[i] for i in is])
end
Base.@propagate_inbounds function Base.getindex(xs::IDVector{Tag,T},
                                                is::SVector{N,ID{Tag}}) where {T,
                                                                               N,
                                                                               Tag}
    return SVector{N,T}(Tuple(xs[i] for i in is))
end
Base.keys(xs::IDVector) = axes(xs, 1)
Base.last(xs::IDVector) = last(xs.vec)
Base.length(xs::IDVector) = length(xs.vec)
Base.ndims(::IDVector) = 1
Base.ndims(::Type{<:IDVector}) = 1
Base.@propagate_inbounds function Base.setindex!(xs::IDVector{Tag}, val,
                                                 i::ID{Tag}) where {Tag}
    return setindex!(xs.vec, val, i.id)
end
Base.size(xs::IDVector) = (length(xs),)
Base.size(xs::IDVector, d) = d == 1 ? length(xs) : 1

Base.pop!(xs::IDVector) = pop!(xs.vec)
Base.popfirst!(xs::IDVector) = popfirst!(xs.vec)
Base.push!(xs::IDVector, vals...) = (push!(xs.vec, vals...); xs)
Base.pushfirst!(xs::IDVector, vals...) = (pushfirst!(xs.vec, vals...); xs)

Base.copy(xs::IDVector{Tag,T}) where {Tag,T} = IDVector{Tag,T}(copy(xs.vec))
Base.iterate(xs::IDVector, state...) = iterate(xs.vec, state...)
Base.map(f, xs::IDVector{Tag}) where {Tag} = IDVector{Tag}(map(f, xs.vec))
function Base.map(f, xs::IDVector{Tag}, yss::IDVector{Tag}...) where {Tag}
    return IDVector{Tag}(map(f, xs.vec, (ys.vec for ys in yss)...))
end
Base.mapreduce(f, op, xs::IDVector; kws...) = mapreduce(f, op, xs.vec; kws...)
Base.reduce(f, xs::IDVector; kws...) = reduce(f, xs.vec; kws...)

Base.:+(xs::IDVector{Tag}) where {Tag} = IDVector{Tag}(+xs.vec)
Base.:-(xs::IDVector{Tag}) where {Tag} = IDVector{Tag}(-xs.vec)
function Base.:+(xs::IDVector{Tag}, ys::IDVector{Tag}) where {Tag}
    return IDVector{Tag}(xs.vec + ys.vec)
end
function Base.:-(xs::IDVector{Tag}, ys::IDVector{Tag}) where {Tag}
    return IDVector{Tag}(xs.vec - ys.vec)
end
Base.:*(a::Number, xs::IDVector{Tag}) where {Tag} = IDVector{Tag}(a * xs.vec)
Base.:\(a::Number, xs::IDVector{Tag}) where {Tag} = IDVector{Tag}(a \ xs.vec)
Base.:*(xs::IDVector{Tag}, a::Number) where {Tag} = IDVector{Tag}(xs.vec * a)
Base.:/(xs::IDVector{Tag}, a::Number) where {Tag} = IDVector{Tag}(xs.vec / a)

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

# Sparse matrix constructors

export MakeSparse
struct MakeSparse{Tv,Iv}
    m::Iv
    n::Iv
    combine::Any
    I::Vector{Iv}
    J::Vector{Iv}
    V::Vector{Tv}
    function MakeSparse{Tv,Iv}(m::Integer, n::Integer,
                               combine=nothing) where {Tv,Iv}
        return new{Tv,Iv}(m, n, combine, Vector{Iv}(), Vector{Iv}(),
                          Vector{Tv}())
    end
    function MakeSparse{Tv}(m::Integer, n::Integer, combine=nothing) where {Tv}
        return MakeSparse{Tv,Int}(m, n, combine)
    end
end

function SparseArrays.sparse(ms::MakeSparse)
    combines = ms.combine ≡ nothing ? () : (combine,)
    return sparse(ms.I, ms.J, ms.V, ms.m, ms.n, combines...)
end

Base.size(ms::MakeSparse) = (ms.m, ms.n)

Base.@propagate_inbounds function Base.setindex!(ms::MakeSparse, v, i::Integer,
                                                 j::Integer)
    @assert 1 ≤ i ≤ ms.m
    @assert 1 ≤ j ≤ ms.n
    push!(ms.I, i)
    push!(ms.J, j)
    push!(ms.V, v)
    return v
end

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

# Sparse matrix iterators

export sparse_column
"""
All nonzero row indices and values for a column ∈ a sparse matrix
"""
function sparse_column(A::SparseMatrixCSC{T,I}, col::Integer) where {T,I}
    return SparseMatrixCSCColumn(Val(:RowVal), I, A, col)
end

# TODO: write this as `keys(sparse_column(...))`
export sparse_column_rows
"""
All row indices of nonzero values for a column ∈ a sparse matrix
"""
function sparse_column_rows(A::SparseMatrixCSC{T,I}, col::Integer) where {T,I}
    return SparseMatrixCSCColumn(Val(:Row), I, A, col)
end

export sparse_column_values
"""
All nonzero values for a column ∈ a sparse matrix
"""
function sparse_column_values(A::SparseMatrixCSC{T,I}, col::Integer) where {T,I}
    return SparseMatrixCSCColumn(Val(:Val), I, A, col)
end

struct SparseMatrixCSCColumn{F,Idx,T,I} <: AbstractVector{T}
    rowvals::Vector{I}
    nonzeros::Vector{T}
    nzrange::UnitRange{Int}
    function SparseMatrixCSCColumn(::F, ::Type{Idx}, A::SparseMatrixCSC{T,I},
                                   col::Integer) where {F,Idx,T,I}
        return new{F,Idx,T,I}(rowvals(A), nonzeros(A), nzrange(A, col))
    end
end

function process(iter::SparseMatrixCSCColumn{Val{:RowVal},Idx}, ind) where {Idx}
    return Idx(iter.rowvals[ind]), iter.nonzeros[ind]
end
function process(iter::SparseMatrixCSCColumn{Val{:Row},Idx}, ind) where {Idx}
    return Idx(iter.rowvals[ind])
end
process(iter::SparseMatrixCSCColumn{Val{:Val}}, ind) = iter.nonzeros[ind]

function Base.eltype(::SparseMatrixCSCColumn{Val{:RowVal},Idx,T,I}) where {Idx,
                                                                           T,I}
    return Tuple{Idx,T}
end
Base.eltype(::SparseMatrixCSCColumn{Val{:Row},Idx,T,I}) where {Idx,T,I} = Idx
Base.eltype(::SparseMatrixCSCColumn{Val{:Val},Idx,T}) where {Idx,T} = T
Base.first(iter::SparseMatrixCSCColumn) = process(iter, first(iter.nzrange))
Base.@propagate_inbounds function Base.getindex(iter::SparseMatrixCSCColumn, i)
    return process(iter, iter.nzrange[i])
end
Base.keys(iter::SparseMatrixCSCColumn) = iter.nzrange
Base.last(iter::SparseMatrixCSCColumn) = process(iter, last(iter.nzrange))
Base.length(iter::SparseMatrixCSCColumn) = length(iter.nzrange)
Base.size(iter::SparseMatrixCSCColumn, i) = length(iter)
Base.size(iter::SparseMatrixCSCColumn) = (length(iter),)

function Base.iterate(iter::SparseMatrixCSCColumn, state...)
    ind_next = iterate(iter.nzrange, state...)
    ind_next ≡ nothing && return nothing
    ind, next = ind_next
    return process(iter, ind), next
end

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

export SparseOp
"""
A type-tagged sparse operator
"""
struct SparseOp{Tag1,Tag2,T}
    op::SparseMatrixCSC{T,Int}

    function SparseOp{Tag1,Tag2,T}(op::SparseMatrixCSC{T,Int}) where {Tag1,Tag2,
                                                                      T}
        return new{Tag1,Tag2,T}(op)
    end
end

function SparseOp{Tag1,Tag2}(op::SparseMatrixCSC{T,Int}) where {Tag1,Tag2,T}
    return SparseOp{Tag1,Tag2,T}(op)
end

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

# I/O

function Base.show(io::IO, A::SparseOp{Tag1,Tag2,T}) where {Tag1,Tag2,T}
    # Convert to CSR -- this is expensive!
    println(io, "$(size(A.op,1))×$(size(A.op,2)) SparseOp{$Tag1,$Tag2,$T} ",
            "with $(nnz(A.op)) stored entries:")
    println(io, "  col: rows")
    # Aop′ = permutedims(A.op)
    Aop′ = A.op
    for j in 1:size(Aop′, 2)
        print(io, "  [$j]: ")
        if !T.mutable && sizeof(T) == 0
            # Immutable types with size 0 have always the same value,
            # hence we don't need to output it
            print(io, "[")
            didoutput = false
            for i in sparse_column_rows(Aop′, j)
                didoutput && print(io, ", ")
                didoutput = true
                print(io, "$i")
            end
            print(io, "]")
        else
            didoutput = false
            for (i, v) in sparse_column(Aop′, j)
                didoutput && print(io, ", ")
                didoutput = true
                print(io, "[$i:$v]")
            end
        end
        println(io)
    end
end

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

# Comparisons

function Base.:(==)(A::SparseOp{Tag1,Tag2},
                    B::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    A ≡ B && return true
    return A.op == B.op
end
function Base.:(<)(A::SparseOp{Tag1,Tag2},
                   B::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return A.op < B.op
end
function Base.isequal(A::SparseOp{Tag1,Tag2},
                      B::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return isequal(A.op, B.op)
end
function Base.hash(A::SparseOp{Tag1,Tag2}, h::UInt) where {Tag1,Tag2}
    return hash(0xf9ecec86, hash(Tag2, hash(Tag1, hash(A.op, h))))
end

function Base.rand(::Type{SparseOp{Tag1,Tag2,T}}, nrows::Integer,
                   ncols::Integer) where {Tag1,Tag2,T}
    p = clamp(4 / min(nrows, ncols), 0, 1)
    return SparseOp{Tag1,Tag2}(sprand(T, nrows, ncols, p))
end

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

# Collection

Base.eltype(::Type{<:SparseOp{<:Any,<:Any,T}}) where {T} = T
Base.isempty(A::SparseOp) = isempty(A.op)
Base.iterate(A::SparseOp, state...) = iterate(A.op, state...)
Base.length(A::SparseOp) = length(A.op)
# function Base.map(f::F, A::SparseOp{Tag1,Tag2},
#                   Bs::SparseOp{Tag1,Tag2}...) where {F,Tag1,Tag2}
#     if nnz(A.op) == 0
#         U = typeof(f(one(eltype(A))))
#         return zero(SparseOp{Tag1,Tag2,U}, size(A, 1), size(A, 2))
#     end
#     SparseOp{Tag1,Tag2}(map(f, A.op, map(B -> B.op, Bs)...))
# end
function Base.map(f, A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    U = typeof(f(one(eltype(A))))
    Bop = SparseMatrixCSC{U,Int}(A.op.m, A.op.n, A.op.colptr, A.op.rowval,
                                 map(f, A.op.nzval))
    return SparseOp{Tag1,Tag2}(Bop)
end
function Base.reduce(f, A::SparseOp{Tag1,Tag2}, Bs::SparseOp{Tag1,Tag2}...;
                     kw...) where {Tag1,Tag2}
    return reduce(f, A.op, map(B -> B.op, Bs)...; kw...)
end

function sparse_column(A::SparseOp{Tag1,Tag2}, col::ID{Tag2}) where {Tag1,Tag2}
    return SparseMatrixCSCColumn(Val(:RowVal), ID{Tag1}, A.op, col.id)
end
function sparse_column_rows(A::SparseOp{Tag1,Tag2},
                            col::ID{Tag2}) where {Tag1,Tag2}
    return SparseMatrixCSCColumn(Val(:Row), ID{Tag1}, A.op, col.id)
end
function sparse_column_values(A::SparseOp{Tag1,Tag2},
                              col::ID{Tag2}) where {Tag1,Tag2}
    return SparseMatrixCSCColumn(Val(:Val), ID{Tag1}, A.op, col.id)
end

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

# Operators are an abstract matrix

Base.IndexStyle(::Type{<:SparseOp}) = IndexCartsian()
Base.axes(A::SparseOp) = ntuple(dir -> axes(A.op, dir), ndims(A))
@inline function Base.axes(A::SparseOp{Tag1,Tag2}, dir) where {Tag1,Tag2}
    dir == 1 && return ID{Tag1}(1):ID{Tag1}(size(A, 1))
    dir == 2 && return ID{Tag2}(1):ID{Tag2}(size(A, 2))
    return throw(DimensionMismatch())
end
Base.eachindex(A::SparseOp) = CartesianIndices(axes(A))
Base.@propagate_inbounds function Base.getindex(A::SparseOp, inds...)
    return throw(ArgumentError())
end
Base.@propagate_inbounds function Base.getindex(A::SparseOp{Tag1,Tag2},
                                                i::ID{Tag1},
                                                j::ID{Tag2}) where {Tag1,Tag2}
    return getindex(A.op, i.id, j.id)
end
Base.ndims(::SparseOp) = 2
Base.size(A::SparseOp) = size(A.op)
Base.size(A::SparseOp, dir) = size(A.op, dir)

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

# Vector space

function Base.zero(::Type{SparseOp{Tag1,Tag2,T}}, nrows::Integer,
                   ncols::Integer) where {Tag1,Tag2,T}
    return SparseOp{Tag1,Tag2}(spzeros(T, nrows, ncols))
end
Base.iszero(A::SparseOp) = iszero(A.op)

function Base.:+(A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(+A.op)
end
function Base.:-(A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(-A.op)
end

function Base.:+(A::SparseOp{Tag1,Tag2},
                 B::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(A.op + B.op)
end
function Base.:-(A::SparseOp{Tag1,Tag2},
                 B::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(A.op - B.op)
end

function Base.:*(a::Number, A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(a * A.op)
end
function Base.:\(a::Number, A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(a \ A.op)
end
function Base.:*(A::SparseOp{Tag1,Tag2}, a::Number) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(A.op * a)
end
function Base.:/(A::SparseOp{Tag1,Tag2}, a::Number) where {Tag1,Tag2}
    return SparseOp{Tag1,Tag2}(A.op / a)
end

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

# Algebra

function Base.one(::Type{SparseOp{Tag,Tag,T}}, nrows::Integer,
                  ncols::Integer) where {Tag,T}
    @assert nrows == ncols
    return SparseOp{Tag,Tag}(sparse(one(T) * I, nrows, ncols))
end
Base.isone(A::SparseOp{Tag,Tag}) where {Tag} = A.op == I

function Base.:*(A::SparseOp{Tag1,Tag2},
                 B::SparseOp{Tag2,Tag3}) where {Tag1,Tag2,Tag3}
    Cop = A.op * B.op
    if eltype(Cop) <: Integer
        dropzeros!(Cop)
    end
    return SparseOp{Tag1,Tag3}(Cop)
end

function Base.adjoint(A::SparseOp{Tag1,Tag2}) where {Tag1,Tag2}
    return SparseOp{Tag2,Tag1}(permutedims(A.op))
end

function Base.:*(A::SparseOp{Tag1,Tag2}, x::IDVector{Tag2}) where {Tag1,Tag2}
    return IDVector{Tag1}(A.op * x.vec)
end
function Base.:\(A::SparseOp{Tag1,Tag2}, x::IDVector{Tag1}) where {Tag1,Tag2}
    return IDVector{Tag2}(A.op \ x.vec)
end

end


text:
abstract type ELFRel{H <: ELFHandle} end
abstract type ELFRela{H <: ELFHandle} end

@io struct ELFRel32{H <: ELFHandle} <: ELFRel{H}
    r_offset::UInt32
    r_info::UInt32
end

@io struct ELFRela32{H <: ELFHandle} <: ELFRela{H}
    r_offset::UInt32
    r_info::UInt32
    r_addend::Int32
end

@io struct ELFRel64{H <: ELFHandle} <: ELFRel{H}
    r_offset::UInt64
    r_info::UInt64
end

@io struct ELFRela64{H <: ELFHandle} <: ELFRela{H}
    r_offset::UInt64
    r_info::UInt64
    r_addend::Int64
end

# Access to relocations
struct Relocations{T <: ELFRel, S <: SectionRef}
    sec::S
end
function Relocations(sec::SectionRef)
    is64 = isa(sec.handle.file,ELF64.File)
    isRela = sec.header.sh_type == SHT_RELA
    Relocations{is64 ? (isRela ? ELF64.Rela : ELF64.Rel) : (isRela ? ELF32.Rela : ELF32.Rel),typeof(sec)}(sec)
end

struct RelocationRef{T <: ELFRel} <: ObjFileBase.RelocationRef{ELFHandle}
    h::ELFHandle
    reloc::T
end

deref(x::RelocationRef) = x.reloc

entrysize{T}(s::Relocations{T}) = sizeof(T)
endof{T}(s::Relocations{T}) = div(s.sec.header.sh_size,entrysize(s))
length(r::Relocations) = endof(r)
function getindex{T}(s::Relocations{T},n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    offset = sectionoffset(s.sec) + (n-1)*entrysize(s)
    seek(s.sec.handle,offset)
    RelocationRef{T}(s.sec.handle,unpack(s.sec.handle, T))
end

start(s::Relocations) = 1
done(s::Relocations,n) = n > length(s)
next(s::Relocations,n) = (x=s[n];(x,n+1))

# Utilities for extracting information from relocation entries
r_sym(rel::Union{ELF32.Rel,ELF32.Rela}) = rel.r_info >> 8
r_type(rel::Union{ELF32.Rel,ELF32.Rela}) = rel.r_info % UInt8
r_sym(rel::Union{ELF64.Rel,ELF64.Rela}) = rel.r_info >> 32
r_type(rel::Union{ELF64.Rel,ELF64.Rela}) = rel.r_info % UInt32
addend(rel::Union{ELF32.Rela,ELF64.Rela}) = rel.r_addend
addend(rel::Union{ELF32.Rel,ELF64.Rel}) = 0

r_sym(x) = r_sym(deref(x))
r_type(x) = r_type(deref(x))
addend(x) = addend(deref(x))

# The relocation to apply
struct RelocToApply
    value::UInt64
    size::UInt8
end

function compute_value(h, symbols, reloc, secs)
    # Here one could to external symbol lookups, etc, but I'm not interested
    UInt(symbolvalue(symbols[r_sym(reloc) + 1], secs)), deref(symbols[r_sym(reloc) + 1]).st_value
end

# Apply relocations in `buffer`. `h` should be the buffer being relocated
function relocate!(buffer, h; LOI = nothing)
    sects = Sections(h)
    rel_sects = filter(sects) do sec
        sht = deref(sec).sh_type
        sht == SHT_RELA || sht == SHT_REL
    end
    for sec in rel_sects
        symbols = Symbols(link_sec(sec))
        relocated_sec = info_sec(sec)
        for reloc in Relocations(sec)
            if r_type(reloc) == 0
                continue
            end
            Value, STValue = compute_value(h, symbols, reloc, sects)
            if LOI === nothing
                SectionLoadAddress = deref(relocated_sec).sh_addr
            else
                SectionLoadAddress = getSectionLoadAddress(LOI,relocated_sec)
            end
            rta = compute_relocation(h, reloc, SectionLoadAddress, Value, STValue)
            seek(buffer, deref(relocated_sec).sh_offset + deref(reloc).r_offset)
            write(buffer, rta.size == 8 ? rta.value :
                rta.size == 4 ? convert(UInt32,rta.value) :
                rta.size == 2 ? convert(UInt16,rta.value) :
                rta.size == 1 ? convert(UInt8,rta.value) :
                error("Unsupported Relocation Size"))
        end
    end
end

function compute_relocation(h, reloc, SectionLoadAddress, Value, STValue)
    if h.file.header.e_machine == EM_X86_64
        return compute_X86_64(reloc, SectionLoadAddress, Value, STValue)
    else
        error("Unknown Machine!")
    end
end

function compute_X86_64(reloc, SectionLoadAddress, Value, STValue)
    # The distinction between Value and STValue deserves some explanation.
    # The added difficulty stems from the fact that we're also relocating
    # relocatable object whose sections need not be contiguous in memory,
    # which is an assumption that is implicitly made in the ELF relocation
    # rules. Instead, we need to confront the possibility that the location
    # at which we're applying the relocation and `Value` are contained in
    # sections with different slides. To counteract this, `Value` is the full
    # load address of the target. It is to be used in computations that demand
    # comparison to other absolute addresses. `STValue` on the other hand is the
    # symbol's `st_value`

    # Any error generated by the conversions below indicates a bad
    # value (the conversions are checked)
    kind = r_type(reloc)
    Offset = deref(reloc).r_offset
    Addend = addend(reloc)
    Value += Addend
    STValue += Addend
    if kind == R_X86_64_64
        return RelocToApply(Value, 8)
    elseif kind == R_X86_64_32
        return RelocToApply(convert(UInt32,STValue), 4)
    elseif kind == R_X86_64_32S
        return RelocToApply(convert(Int32,reinterpret(Int,STValue)), 4)
    else
        PCOffset = Value - UInt(SectionLoadAddress + Offset)
        if kind == R_X86_64_PC64
            return RelocToApply(PCOffset, 8)
        elseif kind == R_X86_64_PC32
            return RelocToApply(convert(Int32,PCOffset), 4)
        elseif kind == R_X86_64_GLOB_DAT || kind == R_X86_64_JUMP_SLOT
            return RelocToApply(Value, 8)
        elseif kind == R_X86_64_RELATIVE
            return RelocToApply(Addend #= + LoadAddress =#, 8)
        else
            error("Unknown relocation ($(R_X86_64[kind]))")
        end
    end
end


text:
function tdvp(
  order::TDVPOrder, solver, PH, time_step::Number, psi::MPS; current_time=0.0, kwargs...
)
  orderings = ITensorTDVP.orderings(order)
  sub_time_steps = ITensorTDVP.sub_time_steps(order)
  sub_time_steps *= time_step
  global info
  for substep in 1:length(sub_time_steps)
    psi, PH, info = tdvp(
      orderings[substep], solver, PH, sub_time_steps[substep], psi; current_time, kwargs...
    )
    current_time += sub_time_steps[substep]
  end
  return psi, PH, info
end

function tdvp(direction::Base.Ordering, solver, PH, time_step::Number, psi::MPS; kwargs...)
  PH = copy(PH)
  psi = copy(psi)

  if length(psi) == 1
    error(
      "`tdvp` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.",
    )
  end

  nsite::Int = get(kwargs, :nsite, 2)
  reverse_step::Bool = get(kwargs, :reverse_step, true)
  normalize::Bool = get(kwargs, :normalize, false)
  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
  svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
  observer = get(kwargs, :observer!, NoObserver())
  outputlevel = get(kwargs, :outputlevel, 0)
  sw = get(kwargs, :sweep, 1)
  current_time = get(kwargs, :current_time, 0.0)

  maxdim::Integer = get(kwargs, :maxdim, typemax(Int))
  mindim::Integer = get(kwargs, :mindim, 1)
  cutoff::Real = get(kwargs, :cutoff, 1E-16)
  noise::Real = get(kwargs, :noise, 0.0)

  N = length(psi)
  set_nsite!(PH, nsite)
  if direction == Base.Forward
    if !isortho(psi) || orthocenter(psi) != 1
      orthogonalize!(psi, 1)
    end
    @assert isortho(psi) && orthocenter(psi) == 1
    position!(PH, psi, 1)
  elseif direction == Base.Reverse
    if !isortho(psi) || orthocenter(psi) != N - nsite + 1
      orthogonalize!(psi, N - nsite + 1)
    end
    @assert(isortho(psi) && (orthocenter(psi) == N - nsite + 1))
    position!(PH, psi, N - nsite + 1)
  end

  maxtruncerr = 0.0
  for (b, ha) in sweepnext(N; ncenter=nsite)
    # unidirectional (half-)sweeps only, skip over the other direction
    if direction == Base.Forward && ha == 2
      continue
    elseif direction == Base.Reverse && ha == 1
      continue
    end

    # Do 'forwards' evolution step
    set_nsite!(PH, nsite)
    position!(PH, psi, b)
    if nsite == 1
      phi1 = psi[b]
    elseif nsite == 2
      phi1 = psi[b] * psi[b + 1]
    end
    phi1, info = solver(PH, time_step, phi1; current_time, outputlevel)

    current_time += time_step

    normalize && (phi1 /= norm(phi1))

    spec = nothing
    if nsite == 1
      psi[b] = phi1
    elseif nsite == 2
      ortho = ha == 1 ? "left" : "right"

      drho = nothing
      if noise > 0.0 && ha == 1
        drho = noise * noiseterm(PH, phi, ortho)
      end

      spec = replacebond!(
        psi,
        b,
        phi1;
        maxdim,
        mindim,
        cutoff,
        eigen_perturbation=drho,
        ortho=ortho,
        normalize,
        which_decomp,
        svd_alg,
      )
      maxtruncerr = max(maxtruncerr, spec.truncerr)
    end

    #
    # Do backwards evolution step
    #
    if reverse_step && (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1)
      b1 = (ha == 1 ? b + 1 : b)
      Δ = (ha == 1 ? +1 : -1)
      if nsite == 2
        phi0 = psi[b1]
      elseif nsite == 1
        uinds = uniqueinds(phi1, psi[b + Δ])
        U, S, V = svd(phi1, uinds)
        psi[b] = U
        phi0 = S * V
        if ha == 1
          ITensors.setleftlim!(psi, b)
        elseif ha == 2
          ITensors.setrightlim!(psi, b)
        end
      end

      set_nsite!(PH, nsite - 1)
      position!(PH, psi, b1)

      phi0, info = solver(PH, -time_step, phi0; current_time, outputlevel)

      current_time -= time_step

      normalize && (phi0 ./= norm(phi0))

      if nsite == 2
        psi[b1] = phi0
      elseif nsite == 1
        psi[b + Δ] = phi0 * psi[b + Δ]
        if ha == 1
          ITensors.setrightlim!(psi, b + Δ + 1)
        elseif ha == 2
          ITensors.setleftlim!(psi, b + Δ - 1)
        end
      end
      set_nsite!(PH, nsite)
    end

    if outputlevel >= 2
      @printf("Sweep %d, half %d, bond (%d,%d) \n", sw, ha, b, b + 1)
      print("  Truncated using")
      @printf(" cutoff=%.1E", cutoff)
      @printf(" maxdim=%.1E", maxdim)
      print(" mindim=", mindim)
      print(" current_time=", round(current_time; digits=3))
      println()
      if spec != nothing
        @printf(
          "  Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b))
        )
      end
      flush(stdout)
    end

    half_sweep_is_done = ((b == 1 && ha == 2) || (b == N && ha == 1))
    update!(
      observer;
      psi,
      bond=b,
      sweep=sw,
      half_sweep=ha,
      spec,
      outputlevel,
      half_sweep_is_done,
      current_time,
    )
  end

  # Just to be sure:
  normalize && normalize!(psi)

  return psi, PH, TDVPInfo(maxtruncerr)
end


text:
"""

- `⁺pᵢ` : Particle i at next step
- `pᵢ`  : Particle i
- `pⱼ`  : Particle j
- `ϵ`   : Softening parameter
"""
function evaluate_Euler!(⁺aᵢ, pᵢ, pⱼ, ϵ)
    r = pⱼ.r .- pᵢ.r
    r⁻³ = (norm(r)^2 + ϵ^2)^(-3/2)

    @. ⁺aᵢ += G * pⱼ.m * r⁻³ * r
end

function evaluate_Euler!(ps, ϵ)
    for pᵢ in ps
        pᵢ.a .= 0.
        for pⱼ in ps
            pᵢ != pⱼ && evaluate_Euler!(pᵢ.a, pᵢ, pⱼ, ϵ)
        end
    end
end

function update!(p::AbstractParticle, Δt)
    @. p.v += p.a * Δt
    @. p.r += p.v * Δt
end

function update!(ps, Δt)
    for p in ps
        update!(p, Δt)
    end
end

function run_Euler(ps, params_sim)
    @unpack ϵ, Δt, t_end = params_sim

    times = (0:Δt:t_end)
    ts, rs, vs, as = prep_snapshot(ps, params_sim)

    for (i, t) in enumerate(times)
        evaluate_Euler!(ps, ϵ)
        save_snapshot!(ts, rs, vs, as, i, t, ps, params_sim)
        update!(ps, Δt)
    end
    ts, rs, vs, as
end

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
#=
This file is part of Pseudospectra.jl.

Julia implementation
Copyright (c) 2017 Ralph A. Smith

Portions derived from EigTool:
 Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
 of the University of Oxford, and the EigTool Developers. All rights reserved.
 EigTool is maintained on GitHub:  https://github.com/eigtool

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    transient_bestlb(ps_data, ftype, it_range) => pts,bnds,sel_pt

compute lower bound for transient growth of a linear system

`ftype` may be `:exp` or `:powers` to indicate the sort of transient.

This is based on Trefethen & Embree, Spectra and Pseudospectra, chs. 15, 16.
"""
function transient_bestlb(ps_data::PSAStruct,ftype::Symbol,it_range;
                          dk=0.25, method=:new)
    A = ps_data.matrix
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    x,y,Z = zoom.x, zoom.y, zoom.Z
    badresult = [NaN],[NaN],[NaN]
    if ftype == :exp
        R = 1.0 ./ Z # resolvent norm
        ωA = 0.5 * maximum(eigvals(A+A')) # numerical abscissa, Eq. 14.2seq.
        if method == :old
            xcrit = 0
        else
            xcrit = ωA
        end
        xtmp = [ifelse(xi<=xcrit,NaN,xi) for xi in x]
        if all(isnan.(xtmp))
            @warn("No points in portrait reach threshold [Re(z)>" *
                 "$xcrit]: expand grid and recompute")
            return badresult
        end
        X = xtmp .* ones(1,length(y))
        t_pts = collect(it_range)
        nt = length(t_pts)
        bnd = zeros(nt)
        sel_pt = Matrix{eltype(x)}(undef,nt,2)
        for it in eachindex(t_pts)
            t = t_pts[it]
            eat = exp.(t.*X)
            if method == :old
                bnds = eat ./ (1 .+ (eat .- 1) ./ (X .* R)) # Eq. 15.11
            elseif method == :new
                # Eq. 15.20, using K as defined at top of Thm. 15.5
                bnds = eat .- (eat .- exp(ωA*t)) ./ ((X .- ωA) .* R)
            else
                # combine the two
                bnds1 = eat ./ (1 .+ (eat .- 1) ./ (X .* R))
                bnds2 = eat .- (eat .- exp(ωA*t)) ./ ((X .- ωA) .* R)
                bnds = [max(x1,x2) for (x1,x2) in zip(bnds1,bnds2)]
            end
            if all(isnan.(bnds))
                @warn("no meaningful results obtained")
                return badresult
            else
                bnds = map(x->isnan(x) ? 0 : x,bnds)
                m,idx = findmax(vec(bnds))
            end
            j,i = map(x->x+1,divrem(idx-1,size(X,1)))
            bnd[it] = m
            sel_pt[it,:] = [x[i],y[j]]
        end
        data_pts = t_pts
    elseif ftype == :powers
        R = 1.0 ./ Z
        X = x .* ones(1,length(y))
        Y = ones(length(x),1) .* y'
        pts = abs.(X .+ 1im * Y)
        anrm = norm(A)
        if method == :old
            rcrit = 1
        else
            rcrit = anrm
        end
        pts = map(x->ifelse(x<=rcrit,NaN,x), pts)
        if all(isnan.(pts))
            @warn("No points achieve threshold [|z| > $rcrit]; "
                 * "expand grid and recompute")
            return badresult
        end
        k_pts = collect(it_range[1]:dk:it_range[end]) # collect(it_range)
        nt = length(k_pts)
        bnd = zeros(nt)
        sel_pt = Matrix{eltype(x)}(nt,2)
        nA = norm(A)
        for it in eachindex(k_pts)
            k = k_pts[it]
            pk = pts .^k
            if method == :old
                bnds = pk ./ (1+(pk-1)./((pts-1) .* (pts .* R -1))) # Eq. 16.15
            elseif method == :new
                bnds = pk .- (pk .- nA^k) ./ ((pts -nA).*R) # Eq. 16.24
            else
                bnds1 = pk ./ (1+(pk-1)./((pts-1) .* (pts .* R -1))) # Eq. 16.15
                bnds2 = pk .- (pk .- nA^k) ./ ((pts -nA).*R) # Eq. 16.24
                bnds = [max(x1,x2) for (x1,x2) in zip(bnds1,bnds2)]
            end
            m,idx = findmax(vec(bnds))
            j,i = map(x->x+1,divrem(idx-1,size(X,1)))
            bnd[it] = m
            sel_pt[it,:] =  [x[i],y[j]]
        end
        data_pts = k_pts
    else
        throw(ArgumentError("incomprehensible ftype=$ftype"))
    end
    data_pts,bnd,sel_pt
end


text:
#
# Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
# Copyright 2015 Oliver Heimlich (oheim@posteo.de)
#
# Original author: Marco Nehmeier (unit tests in libieeep1788)
# Converted into portable ITL format by Oliver Heimlich with minor corrections.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#Language imports

#Test library imports
using FactCheck

#Arithmetic library imports
using ValidatedNumerics

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_isCommonInterval_test") do
    @fact iscommon(Interval(-27.0, -27.0)) --> true
    @fact iscommon(Interval(-27.0, 0.0)) --> true
    @fact iscommon(Interval(0.0, 0.0)) --> true
    @fact iscommon(Interval(-0.0, -0.0)) --> true
    @fact iscommon(Interval(-0.0, 0.0)) --> true
    @fact iscommon(Interval(0.0, -0.0)) --> true
    @fact iscommon(Interval(5.0, 12.4)) --> true
    @fact iscommon(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023)) --> true
    @fact iscommon(entireinterval(Float64)) --> false
    @fact iscommon(∅) --> false
    @fact iscommon(Interval(-Inf, 0.0)) --> false
    @fact iscommon(Interval(0.0, Inf)) --> false
end

facts("minimal_isCommonInterval_dec_test") do
    @fact iscommon(DecoratedInterval(Interval(-27.0, -27.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(-27.0, 0.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(0.0, 0.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0.0, -0.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0.0, 0.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(0.0, -0.0), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(5.0, 12.4), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023), com)) --> true
    @fact iscommon(DecoratedInterval(Interval(-27.0, -27.0), trv)) --> true
    @fact iscommon(DecoratedInterval(Interval(-27.0, 0.0), def)) --> true
    @fact iscommon(DecoratedInterval(Interval(0.0, 0.0), dac)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0.0, 0.0), def)) --> true
    @fact iscommon(DecoratedInterval(Interval(0.0, -0.0), dac)) --> true
    @fact iscommon(DecoratedInterval(Interval(5.0, 12.4), def)) --> true
    @fact iscommon(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023), trv)) --> true
    @fact iscommon(DecoratedInterval(entireinterval(Float64), dac)) --> false
    @fact iscommon(DecoratedInterval(∅, trv)) --> false
    @fact iscommon(DecoratedInterval(∅, trv)) --> false
    @fact iscommon(DecoratedInterval(Interval(-Inf, 0.0), trv)) --> false
    @fact iscommon(DecoratedInterval(Interval(0.0, Inf), def)) --> false
end

facts("minimal_isSingleton_test") do
    @fact isthin(Interval(-27.0, -27.0)) --> true
    @fact isthin(Interval(-2.0, -2.0)) --> true
    @fact isthin(Interval(12.0, 12.0)) --> true
    @fact isthin(Interval(17.1, 17.1)) --> true
    @fact isthin(Interval(-0.0, -0.0)) --> true
    @fact isthin(Interval(0.0, 0.0)) --> true
    @fact isthin(Interval(-0.0, 0.0)) --> true
    @fact isthin(Interval(0.0, -0.0)) --> true
    @fact isthin(∅) --> false
    @fact isthin(entireinterval(Float64)) --> false
    @fact isthin(Interval(-1.0, 0.0)) --> false
    @fact isthin(Interval(-1.0, -0.5)) --> false
    @fact isthin(Interval(1.0, 2.0)) --> false
    @fact isthin(Interval(-Inf, -0x1.fffffffffffffp1023)) --> false
    @fact isthin(Interval(-1.0, Inf)) --> false
end

facts("minimal_isSingleton_dec_test") do
    @fact isthin(DecoratedInterval(Interval(-27.0, -27.0), def)) --> true
    @fact isthin(DecoratedInterval(Interval(-2.0, -2.0), trv)) --> true
    @fact isthin(DecoratedInterval(Interval(12.0, 12.0), dac)) --> true
    @fact isthin(DecoratedInterval(Interval(17.1, 17.1), com)) --> true
    @fact isthin(DecoratedInterval(Interval(-0.0, -0.0), def)) --> true
    @fact isthin(DecoratedInterval(Interval(0.0, 0.0), com)) --> true
    @fact isthin(DecoratedInterval(Interval(-0.0, 0.0), def)) --> true
    @fact isthin(DecoratedInterval(Interval(0.0, -0.0), dac)) --> true
    @fact isthin(DecoratedInterval(∅, trv)) --> false
    @fact isthin(DecoratedInterval(∅, trv)) --> false
    @fact isthin(DecoratedInterval(entireinterval(Float64), def)) --> false
    @fact isthin(DecoratedInterval(Interval(-1.0, 0.0), dac)) --> false
    @fact isthin(DecoratedInterval(Interval(-1.0, -0.5), com)) --> false
    @fact isthin(DecoratedInterval(Interval(1.0, 2.0), def)) --> false
    @fact isthin(DecoratedInterval(Interval(-Inf, -0x1.fffffffffffffp1023), dac)) --> false
    @fact isthin(DecoratedInterval(Interval(-1.0, Inf), trv)) --> false
end

facts("minimal_isMember_test") do
    @fact -27.0 ∈ Interval(-27.0, -27.0) --> true
    @fact in(-27.0, Interval(-27.0, -27.0)) --> true
    @fact -27.0 ∈ Interval(-27.0, 0.0) --> true
    @fact in(-27.0, Interval(-27.0, 0.0)) --> true
    @fact -7.0 ∈ Interval(-27.0, 0.0) --> true
    @fact in(-7.0, Interval(-27.0, 0.0)) --> true
    @fact 0.0 ∈ Interval(-27.0, 0.0) --> true
    @fact in(0.0, Interval(-27.0, 0.0)) --> true
    @fact -0.0 ∈ Interval(0.0, 0.0) --> true
    @fact in(-0.0, Interval(0.0, 0.0)) --> true
    @fact 0.0 ∈ Interval(0.0, 0.0) --> true
    @fact in(0.0, Interval(0.0, 0.0)) --> true
    @fact 0.0 ∈ Interval(-0.0, -0.0) --> true
    @fact in(0.0, Interval(-0.0, -0.0)) --> true
    @fact 0.0 ∈ Interval(-0.0, 0.0) --> true
    @fact in(0.0, Interval(-0.0, 0.0)) --> true
    @fact 0.0 ∈ Interval(0.0, -0.0) --> true
    @fact in(0.0, Interval(0.0, -0.0)) --> true
    @fact 5.0 ∈ Interval(5.0, 12.4) --> true
    @fact in(5.0, Interval(5.0, 12.4)) --> true
    @fact 6.3 ∈ Interval(5.0, 12.4) --> true
    @fact in(6.3, Interval(5.0, 12.4)) --> true
    @fact 12.4 ∈ Interval(5.0, 12.4) --> true
    @fact in(12.4, Interval(5.0, 12.4)) --> true
    @fact 0.0 ∈ entireinterval(Float64) --> true
    @fact in(0.0, entireinterval(Float64)) --> true
    @fact 5.0 ∈ entireinterval(Float64) --> true
    @fact in(5.0, entireinterval(Float64)) --> true
    @fact 6.3 ∈ entireinterval(Float64) --> true
    @fact in(6.3, entireinterval(Float64)) --> true
    @fact 12.4 ∈ entireinterval(Float64) --> true
    @fact in(12.4, entireinterval(Float64)) --> true
    @fact -71.0 ∈ Interval(-27.0, 0.0) --> false
    @fact in(-71.0, Interval(-27.0, 0.0)) --> false
    @fact 0.1 ∈ Interval(-27.0, 0.0) --> false
    @fact in(0.1, Interval(-27.0, 0.0)) --> false
    @fact -0.01 ∈ Interval(0.0, 0.0) --> false
    @fact in(-0.01, Interval(0.0, 0.0)) --> false
    @fact 0.000001 ∈ Interval(0.0, 0.0) --> false
    @fact in(0.000001, Interval(0.0, 0.0)) --> false
    @fact 111110.0 ∈ Interval(-0.0, -0.0) --> false
    @fact in(111110.0, Interval(-0.0, -0.0)) --> false
    @fact 4.9 ∈ Interval(5.0, 12.4) --> false
    @fact in(4.9, Interval(5.0, 12.4)) --> false
    @fact -6.3 ∈ Interval(5.0, 12.4) --> false
    @fact in(-6.3, Interval(5.0, 12.4)) --> false
    @fact 0.0 ∈ ∅ --> false
    @fact in(0.0, ∅) --> false
    @fact -4535.3 ∈ ∅ --> false
    @fact in(-4535.3, ∅) --> false
    @fact -Inf ∈ ∅ --> false
    @fact in(-Inf, ∅) --> false
    @fact Inf ∈ ∅ --> false
    @fact in(Inf, ∅) --> false
    @fact -Inf ∈ entireinterval(Float64) --> false
    @fact in(-Inf, entireinterval(Float64)) --> false
    @fact Inf ∈ entireinterval(Float64) --> false
    @fact in(Inf, entireinterval(Float64)) --> false
end

facts("minimal_isMember_dec_test") do
    @fact -27.0 ∈ DecoratedInterval(Interval(-27.0, -27.0), trv) --> true
    @fact in(-27.0, DecoratedInterval(Interval(-27.0, -27.0), trv)) --> true
    @fact -27.0 ∈ DecoratedInterval(Interval(-27.0, 0.0), def) --> true
    @fact in(-27.0, DecoratedInterval(Interval(-27.0, 0.0), def)) --> true
    @fact -7.0 ∈ DecoratedInterval(Interval(-27.0, 0.0), dac) --> true
    @fact in(-7.0, DecoratedInterval(Interval(-27.0, 0.0), dac)) --> true
    @fact 0.0 ∈ DecoratedInterval(Interval(-27.0, 0.0), com) --> true
    @fact in(0.0, DecoratedInterval(Interval(-27.0, 0.0), com)) --> true
    @fact -0.0 ∈ DecoratedInterval(Interval(0.0, 0.0), trv) --> true
    @fact in(-0.0, DecoratedInterval(Interval(0.0, 0.0), trv)) --> true
    @fact 0.0 ∈ DecoratedInterval(Interval(0.0, 0.0), def) --> true
    @fact in(0.0, DecoratedInterval(Interval(0.0, 0.0), def)) --> true
    @fact 0.0 ∈ DecoratedInterval(Interval(-0.0, -0.0), dac) --> true
    @fact in(0.0, DecoratedInterval(Interval(-0.0, -0.0), dac)) --> true
    @fact 0.0 ∈ DecoratedInterval(Interval(-0.0, 0.0), com) --> true
    @fact in(0.0, DecoratedInterval(Interval(-0.0, 0.0), com)) --> true
    @fact 0.0 ∈ DecoratedInterval(Interval(0.0, -0.0), trv) --> true
    @fact in(0.0, DecoratedInterval(Interval(0.0, -0.0), trv)) --> true
    @fact 5.0 ∈ DecoratedInterval(Interval(5.0, 12.4), def) --> true
    @fact in(5.0, DecoratedInterval(Interval(5.0, 12.4), def)) --> true
    @fact 6.3 ∈ DecoratedInterval(Interval(5.0, 12.4), dac) --> true
    @fact in(6.3, DecoratedInterval(Interval(5.0, 12.4), dac)) --> true
    @fact 12.4 ∈ DecoratedInterval(Interval(5.0, 12.4), com) --> true
    @fact in(12.4, DecoratedInterval(Interval(5.0, 12.4), com)) --> true
    @fact 0.0 ∈ DecoratedInterval(entireinterval(Float64), trv) --> true
    @fact in(0.0, DecoratedInterval(entireinterval(Float64), trv)) --> true
    @fact 5.0 ∈ DecoratedInterval(entireinterval(Float64), def) --> true
    @fact in(5.0, DecoratedInterval(entireinterval(Float64), def)) --> true
    @fact 6.3 ∈ DecoratedInterval(entireinterval(Float64), dac) --> true
    @fact in(6.3, DecoratedInterval(entireinterval(Float64), dac)) --> true
    @fact 12.4 ∈ DecoratedInterval(entireinterval(Float64), trv) --> true
    @fact in(12.4, DecoratedInterval(entireinterval(Float64), trv)) --> true
    @fact -71.0 ∈ DecoratedInterval(Interval(-27.0, 0.0), trv) --> false
    @fact in(-71.0, DecoratedInterval(Interval(-27.0, 0.0), trv)) --> false
    @fact 0.1 ∈ DecoratedInterval(Interval(-27.0, 0.0), def) --> false
    @fact in(0.1, DecoratedInterval(Interval(-27.0, 0.0), def)) --> false
    @fact -0.01 ∈ DecoratedInterval(Interval(0.0, 0.0), dac) --> false
    @fact in(-0.01, DecoratedInterval(Interval(0.0, 0.0), dac)) --> false
    @fact 0.000001 ∈ DecoratedInterval(Interval(0.0, 0.0), com) --> false
    @fact in(0.000001, DecoratedInterval(Interval(0.0, 0.0), com)) --> false
    @fact 111110.0 ∈ DecoratedInterval(Interval(-0.0, -0.0), trv) --> false
    @fact in(111110.0, DecoratedInterval(Interval(-0.0, -0.0), trv)) --> false
    @fact 4.9 ∈ DecoratedInterval(Interval(5.0, 12.4), def) --> false
    @fact in(4.9, DecoratedInterval(Interval(5.0, 12.4), def)) --> false
    @fact -6.3 ∈ DecoratedInterval(Interval(5.0, 12.4), dac) --> false
    @fact in(-6.3, DecoratedInterval(Interval(5.0, 12.4), dac)) --> false
    @fact 0.0 ∈ DecoratedInterval(∅, trv) --> false
    @fact in(0.0, DecoratedInterval(∅, trv)) --> false
    @fact 0.0 ∈ DecoratedInterval(∅, trv) --> false
    @fact in(0.0, DecoratedInterval(∅, trv)) --> false
    @fact -4535.3 ∈ DecoratedInterval(∅, trv) --> false
    @fact in(-4535.3, DecoratedInterval(∅, trv)) --> false
    @fact -4535.3 ∈ DecoratedInterval(∅, trv) --> false
    @fact in(-4535.3, DecoratedInterval(∅, trv)) --> false
    @fact -Inf ∈ DecoratedInterval(∅, trv) --> false
    @fact in(-Inf, DecoratedInterval(∅, trv)) --> false
    @fact -Inf ∈ DecoratedInterval(∅, trv) --> false
    @fact in(-Inf, DecoratedInterval(∅, trv)) --> false
    @fact Inf ∈ DecoratedInterval(∅, trv) --> false
    @fact in(Inf, DecoratedInterval(∅, trv)) --> false
    @fact Inf ∈ DecoratedInterval(∅, trv) --> false
    @fact in(Inf, DecoratedInterval(∅, trv)) --> false
    @fact -Inf ∈ DecoratedInterval(entireinterval(Float64), trv) --> false
    @fact in(-Inf, DecoratedInterval(entireinterval(Float64), trv)) --> false
    @fact Inf ∈ DecoratedInterval(entireinterval(Float64), def) --> false
    @fact in(Inf, DecoratedInterval(entireinterval(Float64), def)) --> false
end
# FactCheck.exitstatus()


text:
# This file was generated, do not modify it. # hide
clf = machine(KNNClassifier(K=3), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

misclassification_rate(ŷ, y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
# v0.6.0

function countto(n::Integer)
    i = zero(n)
    println("Counting...")
    while i < n
        i += 1
    end
    println("Done!")
end

@time countto(10 ^ 5)
@time countto(10 ^ 10)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module OscarCI

import GitHub
import JSON
import TOML

using ..Helpers

# these are not reexported from the main module but should be easily
# accessible for the ci runner by using OscarDevTools.OscarCI

export github_auth, github_repo, github_repo_exists, find_branch,
       parse_meta, ci_matrix, github_json,
       parse_job, job_meta_env, job_pkgs, github_env_runtests,
       github_env_run_doctests

######
### defaults for julia-version, os and branches

const default_os = [ "ubuntu-latest" ]
const default_julia = [ "~1.6.0-0" ]
const default_branches = [ "<matching>", "release" ]

### end defaults
######

global gh_auth = nothing

function github_auth(;token::AbstractString="")
   if !isempty(token)
      global gh_auth = GitHub.authenticate(token)
   elseif haskey(ENV,"GITHUB_TOKEN")
      global gh_auth = GitHub.authenticate(ENV["GITHUB_TOKEN"])
   else
      @info "Using anonymous github auth"
      global gh_auth = GitHub.AnonymousAuth()
   end
end

function github_repo(pkg::AbstractString; fork=nothing)
   isnothing(gh_auth) && github_auth()
   return GitHub.repo(pkg_url(pkg; full=false, fork=fork); auth=gh_auth)
end

function github_repo_exists(name::AbstractString, org::AbstractString)
   try
      github_repo(name; fork=org)
      return true
   catch
      return false
   end
end

parse_meta(file::AbstractString) = TOML.parsefile(file)

# we try to find a matching branch in the main repo for that pkg
# or in a given fork, or via a full 'url#branch'
function find_branch(pkg::AbstractString, branch::AbstractString; fork=nothing)
   isnothing(gh_auth) && github_auth()
   if startswith(branch, "https://")
      (pkg_fork, pkg_branch, _) = pkg_parsebranch(pkg, branch)
      return find_branch(pkg, pkg_branch; fork=pkg_fork)
   end
   @info "locating branch '$branch' for '$pkg'" * (isnothing(fork) ? "" : " (in fork '$fork')")
   if !isnothing(fork)
      try
         GitHub.reference(github_repo(pkg; fork=fork),
                                 "heads/$branch";
                                 auth=gh_auth)
         @info "  -> found in fork: '$fork'"
         return (pkg_url(pkg; full=true, fork=fork), branch, fork)
      catch
         @info "  -- not found in fork: '$fork'"
      end
   end
   if branch != "master"
      try
         GitHub.reference(github_repo(pkg), "heads/$branch"; auth=gh_auth)
         @info "  -> found in main repo"
         return (pkg_url(pkg; full=true), branch, nothing)
      catch
         @info "  -- not found in main repo"
      end
      @warn "  ** branch '$branch' not found, using default 'master' branch"
   end
   # fallback to master branch in default repo
   return (pkg_url(pkg; full=true), "master", nothing)
end

# generate a dict describing branches, includes, os and julia versions
# for use in github actions
function ci_matrix(meta::Dict{String,Any}; pr=0, fork=nothing, active_repo=nothing)
   isnothing(gh_auth) && github_auth()

   matrix = Dict{String,Any}(meta["env"])
   active_pkg = pkg_from_repo(active_repo)

   # assign defaults if unset
   get!(matrix,"os",default_os)
   get!(matrix,"julia-version",default_julia)
   # matrix["branches"] will be replaced later via global_branches
   global_branches = get(matrix,"branches",copy(default_branches))

   pr_branch = ""
   if pr > 0 && !isnothing(active_pkg) && isnothing(fork)
      @info "fetching $active_pkg PR #$pr."
      ghpr = GitHub.pull_request(github_repo(active_pkg), pr; auth=gh_auth)
      # check if this comes from a fork
      if ghpr.head.ref != "master"
         if ghpr.head.repo.full_name != pkg_url(active_pkg; full=false)
            fork = ghpr.head.user.login
         end
         pr_branch = ghpr.head.ref
         # replace '<matching>' with pr_branch
         global_branches[global_branches.=="<matching>"] .= pr_branch
      else
         @warn "PR branch name 'master' cannot be used for branch autodetection"
      end
      # TODO: we might even look into pr.body and parse the branch from there?
      # e.g.: look for such a line
      # SomePkg.jl: otherUser/SomePkg.jl#branchname
   end

   global_axis_pkgs = []

   # for each package lookup branch with the same name in fork and main repo
   for (pkg,pkgmeta) in meta["pkgs"]
      # adjustments for first version of toml files
      # keep for compat reasons for now
      if isa(pkgmeta, Array)
         @info "legacy: mapping meta to branches"
         meta["pkgs"][pkg] = Dict("branches" => pkgmeta,
                                  "test" => pkg == "Oscar",
                                  "testoptions" => [])
         pkgmeta = meta["pkgs"][pkg]
      end

      totest = get!(pkgmeta, "test", false)
      testopts = get!(pkgmeta, "testoptions", [])

      # ignore currently active repo
      pkg == active_pkg && continue

      # add pkgs without 'branches' entry to the global pkg list
      # and dont create a separate axis
      if !haskey(pkgmeta,"branches")
         push!(global_axis_pkgs,pkg)
         continue
      end

      branches = pkgmeta["branches"]

      if !isempty(pr_branch) && "<matching>" in branches
         (url, branch, pkg_fork) = find_branch(pkg, pr_branch; fork=fork)
         # replace '<matching>' with pr_branch
         if !isnothing(pkg_fork) || !in(branch,branches)
            branches[branches.=="<matching>"] .= isnothing(pkg_fork) ?
                                                 branch : "$url#$branch"
         end
      end
      # remove matching branch specifier if it cannot be found
      filter!(p->p!="<matching>",branches)
      if !isempty(branches)
         matrix[pkg] = [Dict("name" => pkg_parsebranch(pkg,branch)[3],
                             "branch" => branch,
                             "test" => totest,
                             "options" => testopts)
                           for branch in branches]
      end
   end
   if !isempty(global_axis_pkgs)
      pkgs = Dict(pkg => Dict("test" => meta["pkgs"][pkg]["test"],
                              "options" => meta["pkgs"][pkg]["testoptions"])
                  for pkg in global_axis_pkgs )
      namestr = "["*join(global_axis_pkgs,",")*"]"
      branchdicts = []
      for branch in global_branches
         if branch == "<matching>"
            @warn "'<matching>' specified but no PR branch could be found"
            continue
         end
         push!(branchdicts,Dict("name" => "$namestr#$branch",
                                "branch" => branch,
                                "pkgs" => deepcopy(pkgs)))
         # we need to record the fork-url for the matching branch if necessary
         if branch == pr_branch
            namestr_branches = []
            branchfound = false
            for (pkgname, pkgmeta) in pkgs
               (url, pkg_branch, pkg_fork) = find_branch(pkgname, branch; fork=fork)
               push!(namestr_branches,"$pkgname#$pkg_branch")
               if !isnothing(pkg_fork)
                  branchdicts[end]["pkgs"][pkgname]["branch"] = "$url#$pkg_branch"
                  namestr_branches[end] = "$pkgname@$pkg_fork#$pkg_branch"
                  branchfound = true
               elseif pkg_branch != "master"
                  branchdicts[end]["pkgs"][pkgname]["branch"] = "$pkg_branch"
                  branchfound = true
               else
                  # put master in dict, which will be kept only if it is not
                  # already in the global branch list
                  # (to make sure <matching> is removed)
                  branchdicts[end]["pkgs"][pkgname]["branch"] = "$pkg_branch"
               end
            end
            branchdicts[end]["name"] = "matching: ["*join(namestr_branches,",")*"]"
            if !branchfound && "master" in global_branches
               # no matching branch found and master already exists
               pop!(branchdicts)
            end
         end
      end
      matrix["branches"] = branchdicts
   end

   # add includes for custom configurations
   matrix["include"] = []
   if haskey(meta, "include")
      for (name,inc) in meta["include"]
         named_include = Dict()
         for (key,val) in inc
            if key in ("os","julia-version")
               named_include[key] = val
            else
               if val == "<matching>"
                  (url, branch, _) = find_branch(key, pr_branch; fork=fork)
                  val = "$url#$branch"
               end
               named_include[key] = Dict{String,Any}(
                                                     "name" => "$(pkg_parsebranch(key,val)[3])",
                                                     "branch" => val
                                                    )
               named_include[key]["test"] = haskey(meta["pkgs"],key) ?
               meta["pkgs"][key]["test"] : false
               named_include[key]["options"] = haskey(meta["pkgs"],key) ?
               meta["pkgs"][key]["testoptions"] : []
            end
         end
         push!(matrix["include"],named_include)
      end
   end
   return matrix
end

# this allows setting a github output variable 'matrix'
# which we can then use as input for the matrix-strategy
github_json(github_matrix::Dict{String,Any}) =
   "::set-output name=matrix::" * JSON.json(github_matrix)

# extract some data from the matrix context of one job

function job_meta(job_json::AbstractString)
   job_dict = JSON.parse(job_json)
   meta = Dict{String,Any}(filter(p->(!in(p.first, ("os","julia-version","branches"))), job_dict))
   if haskey(job_dict,"branches")
      global_branch = job_dict["branches"]["branch"]
      for (pkgname, pkgmeta) in job_dict["branches"]["pkgs"]
         if !haskey(meta,pkgname)
            meta[pkgname] = Dict("branch" => get(pkgmeta,"branch",
                                                 global_branch),
                                 "test" => pkgmeta["test"],
                                 "options" => pkgmeta["options"],
                                 "name" => "$pkgname#$global_branch" )
         end
      end
   end
   return meta
end

job_meta_env(var::AbstractString) = job_meta(ENV[var])

job_pkgs(job::Dict) = Dict{String,Any}(pkg => val["branch"] for (pkg,val) in job)

# keep this for now for older OscarCI.yml files
parse_job(job_json::AbstractString) = job_pkgs(job_meta(job_json))

# generate one long julia command that tests all packages with test==true
# and passes any test_args to Pkg.test
# the testcommand is then stored as an env variable an run in the next step
# (using the gitub actions `>> $GITHUB_ENV style`
# doing the tests separately via github actions syntax seems complicated
function github_env_runtests(job::Dict; varname::String, filename::String)
   testcmd = ["using Pkg;"]
   for (pkg, param) in job
      if get(param, "test", false)
         if length(get(param, "options", [])) > 0
            push!(testcmd, """Pkg.test("$pkg"; test_args=$(string.(param["options"])));""")
         else
            push!(testcmd, """Pkg.test("$pkg");""")
         end
      end
   end
   open(filename, "a") do io
      println(io, "$varname=", join(testcmd))
   end
end

function github_env_run_doctests(job::Dict; varname::String, filename::String)
   testcmd = ["using Pkg;", "Pkg.add(\"Documenter\");", "Pkg.add(\"DocumenterCitations\");", "Pkg.add(\"DocumenterMarkdown\"); using Documenter;"]
   for (pkg, param) in job
      if pkg == "Oscar"
         push!(testcmd, """using Oscar; DocMeta.setdocmeta!(Oscar, :DocTestSetup, :(using Oscar, Oscar.Graphs); recursive = true); doctest(Oscar);""")
      elseif pkg != "Polymake"
         push!(testcmd, """using $pkg; DocMeta.setdocmeta!($pkg, :DocTestSetup, :(using $pkg); recursive = true); doctest($pkg);""")
      end
   end
   open(filename, "a") do io
      println(io, "$varname=", join(testcmd))
   end
end

end


text:
abstract LTISystem

+(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...)


text:
using NLopt
using Plots
using SparseArrays

struct ModelParameters
    γ; # RRA
    μ; # risk premium
    σ; # vol of risky asset
    ν; # vol of int rate
    ρ; # discount rate
    η; # correlation between risky return and risk free rate
    rgrid; # grid for interest rate
    rmin;  # min and max interest rates in the grid
    rmax;  # -//-
    Δr; # step in interest rates grid
    rbar; # [For the geometric OU]
    θ; # [For the geometric OU]
end

function solve_optimal_risky_share(A, params)
    # This function takes a vector of A (investment opportuinities)
    # for a given value of interest rate and a struct of model
    # parameters. It then numerically solves for optimal risky
    # asset allocation.
    rgrid = params.rgrid
    γ = params.γ
    σ = params.σ
    μ = params.μ
    η = params.η
    ν = params.ν;

    # Modifying for non-zero correlation
    αopt=zeros(size(A))
    for i = 1:size(rgrid)[1]
        ri = rgrid[i]
        Ai = A[i]

        if i == 1
            Ailower = A[1]
            Aihigher = A[2]
        elseif i == I
            Ailower = A[I-1]
            Aihigher = A[I]
        else
            Ailower = A[i-1]
            Aihigher = A[i+1]
        end

        # Forming function to solve numerically
        function to_minimize(x::Vector, grad::Vector)
            ((ri + x[1]*μ - 0.5*x[1]^2*σ^2)^(-γ) * (μ - x[1]*σ^2) + Ai*(1-γ)*x[1]*σ^2 + (Aihigher - Ailower)*σ*ri*ν*η/(2*Δr))^2
        end

        function positive_consumption_constraint(x::Vector, grad::Vector)
            -(ri + x[1]*μ - 0.5*x[1]^2*σ^2)
        end

        opt = Opt(:LN_COBYLA, 1)
        lower_bounds!(opt, [-5])
        upper_bounds!(opt, [Inf])
        ftol_abs!(opt, 1e-20)

        min_objective!(opt, to_minimize)
        inequality_constraint!(opt, (x,g) -> positive_consumption_constraint(x,g), 1e-8)

        (minf, minx, ret) = optimize(opt, [0.7])

        αopt[i] = minx[1]
    end

    return αopt
end

function form_matrix_B(α, params)
    # This function forms matrix B that is used in the iteration
    # step to solve for vector of investment opportunities A
    # given a vector of optimal risky shares

    # Getting parameters:
    Δr = params.Δr;
    γ = params.γ;
    ν = params.ν;
    σ = params.σ;
    γ = params.γ;
    η = params.η;
    I = size(params.rgrid)[1];
    rbar = params.rbar;
    θ = params.θ;

    # Modifying for nonzero correlation between risky and risk free asset.
    # Log risk free rate follows a random walk (implies a positive drift for
    # the level of interest rate). Using central differences
    xvec = ν^2 * rgrid.^2/(2 *Δr^2) .- (1-γ)*α.*rgrid*σ*ν*η/(2*Δr) .- rgrid .* 0.5 * ν^2 *0.5/Δr;
    yvec = 0.5*(1-γ)^2 .* α.^2 .* σ^2 .- rgrid.^2 .*ν^2/Δr^2;
    zvec = ν^2 * rgrid.^2/(2 *Δr^2) .+ (1-γ)*α.*rgrid*σ*ν*η/(2*Δr) .+ rgrid .* 0.5 * ν^2 *0.5/Δr;
    v1 = xvec[1];
    vI = zvec[I];
    yvec[1] = yvec[1] + v1;
    yvec[I] = yvec[I] + vI;

    # Forming matrix B using
    B = spdiagm(-1 => xvec[2:I], 0=>yvec, 1=>zvec[1:(I-1)])
    return B
end

function calculate_utility_flow(α, params)
    μ = params.μ;
    γ = params.γ;
    σ = params.σ;
    rgrid = params.rgrid;
    return (rgrid .+ α.*μ .- 0.5.*α.^2 .* σ^2).^(1-γ)
end

function update_value_function(Acurrent, params)
    # 1. Given current A vector calculate optimal risky share
    αopt = solve_optimal_risky_share(Acurrent, params)

    # 2. Given optimal risky share calculate utility flow
    u = calculate_utility_flow(αopt, params)

    # 3. Given optimal alpha calculate "transition" matrix B
    B = form_matrix_B(αopt, params)

    # 4. Given old A vector, flow utility and "transition"
    # matrix update A vector
    I = size(params.rgrid)[1];
    ρ = params.ρ;
    Δ = 1;
    γ = params.γ;

    Aupdated = ((1/Δ + ρ) .* spdiagm(0 => ones(I)) .- B) \ (u .+ Acurrent./Δ);

    return Aupdated, αopt, u
end

function iterate_value_function(params, A_init, tol)
    A_prev = A_init;
    i = 0;
    dist = Inf
    while dist > tol
        i += 1;
        A_next,_ = update_value_function(A_prev, params);

        # Calculating the difference between value functions:
        dist = sum(abs.(A_next .- A_prev));
        A_prev = A_next;
        if i % 1000 == 0
            @show i
            @show dist
            print("\n")
        end
    end

    # Iterate one more step to get all parameters:
    return update_value_function(A_prev, params)
end

# Start with a guess of negative exponential declining A
rmin = 0;
rmax = 1;
I = 1001;
rgrid = range(rmin, stop = rmax, length = I);
Δr = (rmax - rmin)/(I-1);

I = size(rgrid)[1];
A_init = 10000 * exp.(-rgrid);

# params = ModelParameters(4, 0.04, 0.18, 0.01, 0.05, -0.25, rgrid,rmin,rmax,Δr)
params = ModelParameters(3, 0.04, 0.18, 0.01, 0.075, -0.0, rgrid,rmin,rmax,Δr,0.1,1)

Anew, αopt, u = update_value_function(A_init, params)
Plots.plot(Anew)
Plots.plot(αopt)
Plots.plot(u)

Anew, αopt, u = update_value_function(Anew, params)
Plots.plot(rgrid,Anew)
Plots.plot(rgrid,αopt)
Plots.plot(u*(1-params.γ))

# 1. Given current A vector calculate optimal risky share
αopt = solve_optimal_risky_share(A_init, params)

# 2. Given optimal risky share calculate utility flow
u = calculate_utility_flow(αopt, params)

# 3. Given optimal alpha calculate "transition" matrix B
B = form_matrix_B(αopt, params)
sum(B)

mathematicaColor1 = RGBA(0.37,0.51,0.71,1)
mathematicaColor2 = RGBA(0.88,0.61,0.14,1)
mathematicaColor3 = RGBA(0.56,0.69,0.19,1)

# Solving for different value of correlation:
params1 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0.99, rgrid,rmin,rmax,Δr,0.02,1)
params2 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0.0, rgrid,rmin,rmax,Δr,0.02,1)
params3 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, -0.99, rgrid,rmin,rmax,Δr,0.02,1)

A1, α1, flow1 = iterate_value_function(params1, A_init, 1e-6)
A2, α2, flow2 = iterate_value_function(params2, A_init, 1e-6)
A3, α3, flow3 = iterate_value_function(params3, A_init, 1e-6)

Plots.plot(rgrid, α1, label = string("Corr = ", params1.η),
    xaxis = ("Risk Free Rate", (0,0.1)), yaxis = "Risky Share",legend=:bottomleft, lw = 2.75,grid=false,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.plot!(rgrid, α2, label = string("Corr = ", params2.η), lw = 2.75, color=mathematicaColor2)
Plots.plot!(rgrid, α3, label = string("Corr = ", params3.η), lw = 2.75, color=mathematicaColor3)
Plots.savefig("../figures/share_dynamic_vary_corr.pdf")

# Plotting hedging demand: difference between α under non-zero correlation and
# α under zero correlation.

# Solving for different value of Risk Premium:
params4 = ModelParameters(3, 0.02, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)
params5 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)
params6 = ModelParameters(3, 0.10, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)

A4, α4, flow4 = iterate_value_function(params4, A_init, 1e-4)
A5, α5, flow5 = iterate_value_function(params5, A_init, 1e-4)
A6, α6, flow6 = iterate_value_function(params6, A_init, 1e-7)

Plots.plot(rgrid, α4, label = string("RP = ", params4.μ*100,"%"),
    xaxis = ("Risk Free Rate", (0,0.1)), yaxis = "Risky Share",legend=:bottomleft, lw = 2.75,grid=false,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.plot!(rgrid, α5, label = string("RP = ", params5.μ*100,"%"), lw = 2.75, color=mathematicaColor2)
Plots.plot!(rgrid, α6, label = string("RP = ", params6.μ*100,"%"), lw = 2.75, color=mathematicaColor3)
Plots.savefig("../figures/share_dynamic_vary_rp.pdf")

# Solving for different value of Discount Rate:
params7 = ModelParameters(3, 0.06, 0.18, 0.01, 0.050, 0, rgrid,rmin,rmax,Δr,0.02,1)
params8 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)
params9 = ModelParameters(3, 0.06, 0.18, 0.01, 0.100, 0, rgrid,rmin,rmax,Δr,0.02,1)

A7, α7, flow7 = iterate_value_function(params7, A_init, 1e-7)
A8, α8, flow8 = iterate_value_function(params8, A_init, 1e-7)
A9, α9, flow9 = iterate_value_function(params9, A_init, 1e-7)

Plots.plot(rgrid, α7, label = string("Discount Rate  = ", params7.ρ*100,"%"),
    xaxis = ("Risk Free Rate", (0,0.1)), yaxis = "Risky Share",legend=:bottomleft, lw = 2.75,grid=false,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.plot!(rgrid, α8, label = string("Discount Rate  = ", params8.ρ*100,"%"), lw = 2.75, color=mathematicaColor2)
Plots.plot!(rgrid, α9, label = string("Discount Rate  = ", params9.ρ*100,"%"), lw = 2.75, color=mathematicaColor3)
Plots.savefig("../figures/share_dynamic_vary_discount_rate.pdf")

# Solving for different value of Risk Aversion:
params10 = ModelParameters(2, 0.06, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)
params11 = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)
params12 = ModelParameters(4, 0.06, 0.18, 0.01, 0.075, 0, rgrid,rmin,rmax,Δr,0.02,1)

A10, α10, flow10 = iterate_value_function(params10, A_init, 1e-7)
A11, α11, flow11 = iterate_value_function(params11, A_init, 1e-7)
A12, α12, flow12 = iterate_value_function(params12, A_init, 1e-7)

Plots.plot(rgrid, α10, label = string("RRA = ", params10.γ),
    xaxis = ("Risk Free Rate", (0,0.1)), yaxis = "Risky Share",legend=:bottomleft, lw = 2.75,grid=false,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.plot!(rgrid, α11, label = string("RRA = ", params11.γ), lw = 2.75, color=mathematicaColor2)
Plots.plot!(rgrid, α12, label = string("RRA = ", params12.γ), lw = 2.75, color=mathematicaColor3)
Plots.savefig("../figures/share_dynamic_vary_rra.pdf")

################################################################
# Numerically calculating drift of log consumption

params = ModelParameters(3, 0.06, 0.18, 0.01, 0.075, 0.0, rgrid,rmin,rmax,Δr,0.1,1)
A, α, flow=iterate_value_function(params, A_init, 1e-6)
γ = params.γ
σ = params.σ
μ = params.μ
ν = params.ν
ρ = params.ρ

dα = (α[2:(I-1)] .- α[1:(I-2)])./Δr;
ddα = (α[3:I] - 2*α[2:(I-1)] + α[1:(I-2)])./Δr;
α = α[2:(I-1)]
r = rgrid[2:(I-1)]
f = r .+ α * μ .- 0.5*α.^2 * σ^2;
df = 1 .+ dα * μ .- α .* dα * σ^2;
ddf = ddα * μ .- (dα.^2 + α .* ddα) * σ^2;
dlogc = 0.5*ν^2 * r .* (df .* f .+ ddf .* f - df.^2)./(f.^2);

Plots.plot(r,dlogc,xaxis = ("Risk Free Rate", (0,0.25)),label="",
    yaxis = "Log Consumption Drift",lw = 2.75,grid=false,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.savefig("../figures/dynamic_model_log_consumption_drift.pdf")

# Calculating instantenous volatility of log(c):
voldlogc = ν * r .* df./f

Plots.plot(r,dlogc,xaxis = ("Risk Free Rate", (0,0.25)),label="Drift",
    yaxis = "Log Consumption Drift",lw = 2.75,grid=false,legend=:topleft,
    yguidefontsize=12,guidefontsize=12,legendfontsize=12,tickfontsize=12, color=mathematicaColor1)
Plots.plot!(r,voldlogc,lw = 2.75, color=mathematicaColor2, label = "Volatility")


text:
using KeplerTools
using ProfileView

include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
stime1, stime2 = 0.,                    852. * 6. * 3600.
ftime1, ftime2 = 151. * 6. * 3600.,     453. * 6. * 3600.
startorb = Orbit(kerbin.eqradius + 100000, kerbin)
endorb   = Orbit(duna.eqradius   + 100000, duna)

ProfileView.@profview pc = Porkchop(startorb, endorb, stime1, stime2, ftime1, ftime2; npts=100)
ProfileView.@profview qpc = fastPorkchop(startorb, endorb, stime1, stime2, ftime1, ftime2; npts=100)


text:
function test_most_popular()
    println("-- Testing MostPopular recommender")

    data = DataAccessor(sparse([1 2 3; 4 5 0]))
    recommender = MostPopular(data)
    build!(recommender)
    @test ranking(recommender, 1, 1) == 2.0
    @test ranking(recommender, 1, 3) == 1.0

    n_user, n_item = 5, 10
    events = [Event(1, 2, 1), Event(3, 2, 1), Event(2, 6, 4)]
    data = DataAccessor(events, n_user, n_item)
    recommender = MostPopular(data)
    build!(recommender)
    @test ranking(recommender, 1, 1) == 0.0
    @test ranking(recommender, 1, 2) == 2.0
    @test ranking(recommender, 1, 6) == 1.0
end

test_most_popular()