text:
## Shootouts

type Shootout
  setups::Vector{Dict{Symbol,Any}}
  times#::Vector{Float64}
  errors#::Vector{uType}
  effs#::Vector{Float64} # Efficiencies
  effratios#::Matrix{uEltype}
  solutions
  names::Vector{String}
  N::Int
  bestidx::Int
  winner::String
end

type ShootoutSet
  shootouts::Vector{Shootout}
  probs#::Vector{DEProblem}
  probaux#::Vector{Dict{Symbol,Any}}
  N::Int
  winners::Vector{String}
end

function ode_shootout(prob::AbstractODEProblem,setups;appxsol=nothing,numruns=20,names=nothing,error_estimate=:final,kwargs...)
  N = length(setups)
  errors = Vector{Float64}(N)
  solutions = Vector{AbstractODESolution}(N)
  effs = Vector{Float64}(N)
  times = Vector{Float64}(N)
  effratios = Matrix{Float64}(N,N)
  timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
  dense_errors = error_estimate ∈ DENSE_ERRORS
  if names == nothing
    names = [string(setups[i][:alg]) for i=1:N]
  end
  for i in eachindex(setups)
    sol = solve(prob,setups[i][:alg];timeseries_errors=timeseries_errors,
    dense_errors = dense_errors,kwargs...,setups[i]...) # Compile and get result
    sol = solve(prob,setups[i][:alg],sol[:],sol.t,sol.k;timeseries_errors=timeseries_errors,
    dense_errors = dense_errors,kwargs...,setups[i]...) # Compile and get result
    t = @elapsed for j in 1:numruns
      sol = solve(prob,setups[i][:alg],sol[:],sol.t,sol.k;kwargs...,setups[i]...)
    end
    if appxsol != nothing
      errsol = appxtrue(sol,appxsol)
      errors[i] = errsol.errors[error_estimate]
      solutions[i] = errsol
    else
      errors[i] = sol.errors[error_estimate]
      solutions[i] = sol
    end
    effs[i] = 1/(errors[i]*t)
    t = t/numruns
    times[i] = t
  end
  for j in 1:N, i in 1:N
    effratios[i,j] = effs[i]/effs[j]
  end
  bestidx = find((y)->y==maximum(effs),effs)[1]; winner = names[bestidx]
  return Shootout(setups,times,errors,effs,effratios,solutions,names,N,bestidx,winner)
end

function ode_shootoutset{T<:AbstractODEProblem}(probs::Vector{T},setups;probaux=nothing,numruns=20,names=nothing,kwargs...)
  N = length(probs)
  shootouts = Vector{Shootout}(N)
  winners = Vector{String}(N)
  if names == nothing
    names = [string(setups[i][:alg]) for i=1:length(setups)]
  end
  if probaux == nothing
    probaux = Vector{Dict{Symbol,Any}}(N)
    for i in 1:N
      probaux[i] = Dict{Symbol,Any}()
    end
  end
  for i in eachindex(probs)
    shootouts[i] = ode_shootout(probs[i],setups;numruns=numruns,names=names,kwargs...,probaux[i]...)
    winners[i] = shootouts[i].winner
  end
  return ShootoutSet(shootouts,probs,probaux,N,winners)
end

length(shoot::Shootout) = shoot.N
Base.size(shoot::Shootout) = length(shoot)
Base.endof(shoot::Shootout) = length(shoot)
Base.getindex(shoot::Shootout,i::Int) = shoot.effs[i]
Base.getindex(shoot::Shootout,::Colon) = shoot.effs

function print(io::IO, shoot::Shootout)
  println(io,"Names: $(shoot.names), Winner: $(shoot.winner)")
  println(io,"Efficiencies: $(shoot.effs)")
  println(io,"EffRatios: $(shoot.effratios[shoot.bestidx,:])")
  println(io,"Times: $(shoot.times)")
  println(io,"Errors: $(shoot.errors)")
end

function show(io::IO, shoot::Shootout)
  println(io,"Winner: $(shoot.winner)")
  println(io,"EffRatios: $(shoot.effratios[shoot.bestidx,:])")
end

length(set::ShootoutSet) = set.N
Base.size(set::ShootoutSet) = length(set)
Base.endof(set::ShootoutSet) = length(set)
Base.getindex(set::ShootoutSet,i::Int) = set.shootouts[i]
Base.getindex(set::ShootoutSet,::Colon) = set.shootouts

function print(io::IO, set::ShootoutSet)
  println(io,"ShootoutSet of $(set.N) shootouts")
  println(io,"Winners: $(set.winners)")
end

function show(io::IO, set::ShootoutSet)
  println(io,"ShootoutSet of $(set.N) shootouts ")
end

## WorkPrecisions

type WorkPrecision
  prob
  abstols
  reltols
  errors
  times
  name
  N::Int
end

type WorkPrecisionSet
  wps::Vector{WorkPrecision}
  N::Int
  abstols
  reltols
  prob
  setups
  names
end

function ode_workprecision(prob::AbstractODEProblem,alg,abstols,reltols;name=nothing,numruns=20,appxsol=nothing,error_estimate=:final,kwargs...)
  N = length(abstols)
  errors = Vector{Float64}(N)
  times = Vector{Float64}(N)
  if name == nothing
    name = "WP-Alg"
  end
  timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
  dense_errors = error_estimate ∈ DENSE_ERRORS
  for i in 1:N
    sol = solve(prob,alg;kwargs...,abstol=abstols[i],
    reltol=reltols[i],timeseries_errors=timeseries_errors,
    dense_errors = dense_errors) # Compile and get result
    sol = solve(prob,alg,sol[:],sol.t,sol.k;kwargs...,abstol=abstols[i],
    reltol=reltols[i],timeseries_errors=timeseries_errors,
    dense_errors = dense_errors) # Compile and get result
    t = @elapsed for j in 1:numruns
      sol = solve(prob,alg,sol[:],sol.t,sol.k;kwargs...,abstol=abstols[i],
      reltol=reltols[i],timeseries_errors=timeseries_errors,
      dense_errors = dense_errors)
    end
    t = t/numruns

    if appxsol != nothing
      errsol = appxtrue(sol,appxsol)
      errors[i] = errsol.errors[error_estimate]
    else
      errors[i] = sol.errors[error_estimate]
    end
    times[i] = t
  end
  return WorkPrecision(prob,abstols,reltols,errors,times,name,N)
end

function ode_workprecision_set(prob::AbstractODEProblem,abstols,reltols,setups;numruns=20,names=nothing,appxsol=nothing,kwargs...)
  N = length(setups)
  wps = Vector{WorkPrecision}(N)
  if names == nothing
    names = [string(setups[i][:alg]) for i=1:length(setups)]
  end
  for i in 1:N
    wps[i] = ode_workprecision(prob,setups[i][:alg],abstols,reltols;numruns=numruns,appxsol=appxsol,name=names[i],kwargs...,setups[i]...)
  end
  return WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names)
end

length(wp::WorkPrecision) = wp.N
Base.size(wp::WorkPrecision) = length(wp)
Base.endof(wp::WorkPrecision) = length(wp)
Base.getindex(wp::WorkPrecision,i::Int) = wp.times[i]
Base.getindex(wp::WorkPrecision,::Colon) = wp.times

function print(io::IO, wp::WorkPrecision)
  println(io,"Name: $(wp.name)")
  println(io,"Times: $(wp.times)")
  println(io,"Errors: $(wp.errors)")
end

function show(io::IO, wp::WorkPrecision)
  println(io,"Name: $(wp.name)")
  println(io,"Times: $(wp.times)")
  println(io,"Errors: $(wp.errors)")
end

length(wp_set::WorkPrecisionSet) = wp_set.N
Base.size(wp_set::WorkPrecisionSet) = length(wp_set)
Base.endof(wp_set::WorkPrecisionSet) = length(wp_set)
Base.getindex(wp_set::WorkPrecisionSet,i::Int) = wp_set.wps[i]
Base.getindex(wp_set::WorkPrecisionSet,::Colon) = wp_set.wps

function print(io::IO, wp_set::WorkPrecisionSet)
  println(io,"WorkPrecisionSet of $(wp_set.N) wps")
  println(io,"Names: $(wp_set.names)")
end

function show(io::IO, wp_set::WorkPrecisionSet)
  println(io,"WorkPrecisionSet of $(wp_set.N) wps")
end

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
module ModelingTest
using OnlineStats, GLM, Base.Test

@testset "Modeling" begin
@testset "LinReg" begin
    n, p = 10000, 10
    x = randn(n, p)
    β = collect(1.:p)
    y = x * β + randn(n)
    o = LinReg(x, y, intercept = false)
    @test coef(o) ≈ x \ y
    o2 = LinReg(10, intercept = false)
    fit!(o2, x[1:500, :], y[1:500], 500)
    fit!(o2, x[501:1000, :], y[501:1000], 1)
    @test_approx_eq_eps coef(o) coef(o2) .5

    # vs. GLM
    l = lm(x, y)
    @test predict(l, x) ≈ predict(o, x)
    xi = randn(p)
    @test predict(l, xi') ≈ collect(predict(o, xi))
    @test loss(o, x, y) ≈ 0.5 * mean(abs2(y - predict(o, x)))
    ltab = coeftable(l)
    otab = coeftable(o)

    @test_approx_eq_eps vcov(o)[1:5]    vcov(l)[1:5]    .1
    @test_approx_eq_eps stderr(o)[1:5]  stderr(l)[1:5]  .1
    @test_approx_eq_eps confint(o)[1:5] confint(l)[1:5] .1

    @test coef(LinReg(10, intercept = false)) == zeros(10)
    @test coef(LinReg(10)) == zeros(11)

    o = LinReg(10)
    @test coef(o) == value(o)
    predict(o, randn(100, 10)) == zeros(100)
end
@testset "QuantRegMM" begin
    n, p = 10000, 10
    x = randn(n, p)
    β = collect(1.:p)
    y = x * β + randn(n)

    o = QuantRegMM(x, y)
    fit!(o, x, y, 10)

    @test coef(o) == value(o)
    @test value(o) == o.β
end
@testset "BiasVector / BiasMatrix" begin
    x = randn(100, 10)
    y = randn(100)
    yb = BiasVector(y)
    xb = BiasMatrix(x)

    @test length(yb)    == 101
    @test size(yb)      == (101,)
    @test length(xb)    == 100 * 10 + 100
    @test size(xb)      == (100, 11)
    @test yb[101]       == 1
    @test xb[1, 11]     == 1
    xb[1, 1] = 2.0
    yb[1] = 2.0
    @test xb[1, 1]      == 2.0
    @test yb[1]         == 2.0
end
@testset "TwoWayInteractionVector/Matrix" begin
    x = rand(5)
    v = TwoWayInteractionVector(x)
    @test v[5]  == x[5]
    @test v[6]  == x[1] * x[2]
    @test v[7]  == x[1] * x[3]
    @test v[8]  == x[1] * x[4]
    @test v[9]  == x[1] * x[5]
    @test v[10] == x[2] * x[3]

    x = rand(10, 5)
    m = TwoWayInteractionMatrix(x)
    @test m[1, 1]   == x[1, 1]
    @test m[1, 5]   == x[1, 5]
    @test m[1, 6]   == x[1, 1] * x[1, 2]
    @test m[1, 10]  == x[1, 2] * x[1, 3]
end
end

end#module

ext: jl


text:
module Linear1DTests
println("Testing Linear interpolation in 1D...")
using Interpolations, Base.Test

f(x) = sin((x-3)*2pi/9 - 1)
xmax = 10

A = Float64[f(x) for x in 1:xmax]

## ExtrapError

itp1 = Interpolation(A, Linear(OnGrid()), ExtrapError())

for x in 1:.2:xmax
    @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
end

@test_throws BoundsError itp1[-3]
@test_throws BoundsError itp1[xmax+.1]

## ExtrapNaN

itp2 = Interpolation(A, Linear(OnGrid()), ExtrapNaN())

for x in 1:.2:xmax
    @test_approx_eq_eps f(x) itp2[x] abs(.1*f(x))
end

xlo, xhi = itp2[.9], itp2[xmax+.2]
@test isnan(xlo)
@test isnan(xhi)

## ExtrapConstant

itp3 = Interpolation(A, Linear(OnGrid()), ExtrapConstant())
for x in 3.1:.2:4.3
    @test_approx_eq_eps f(x) itp3[x] abs(.1*f(x))
end

xlo, xhi = itp3[.9], itp3[xmax+.2]
@test xlo == A[1]
@test xhi == A[end]

# Rational element types
A = Rational{Int}[x//10+2 for x in 1:10]
itp = Interpolation(A, Linear(OnGrid()), ExtrapNaN())
@test itp[11//10] == (11//10)//10+2

end

module Linear2DTests
println("Testing Linear interpolation in 2D...")
using Interpolations, Base.Test

f(x,y) = sin(x/10)*cos(y/6)

xmax, ymax = 30,10

A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]

itp1 = Interpolation(A, Linear(OnGrid()), ExtrapError())

for x in 2.1:.2:xmax-1, y in 1.9:.2:ymax-.9
    @test_approx_eq_eps f(x,y) itp1[x,y] .01
end

end


text:
using DelimitedFiles

# This is a function to simulate a "chaos game"
function chaos_game(n::Int, shape_points)

    # Initializing the output array and the initial point
    output_points = zeros(n,2)
    point = [rand(), rand()]

    for i = 1:n
        output_points[i,:] .= point
        point = 0.5*(rand(shape_points) .+ point)
    end

    return output_points

end

# This will generate a Sierpinski triangle with a chaos game of n points for an
# initial triangle with three points on the vertices of an equilateral triangle:
#     A = (0.0, 0.0)
#     B = (0.5, sqrt(0.75))
#     C = (1.0, 0.0)
# It will output the file sierpinski.dat, which can be plotted after
shape_points = [[0.0, 0.0],
                [0.5, sqrt(0.75)],
                [1.0, 0.0]]
output_points = chaos_game(10000, shape_points)
writedlm("sierpinski.dat", output_points)


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

bpfnative = Pkg.PackageSpec(path = root_directory)
Pkg.develop(bpfnative)
Pkg.precompile()
Pkg.test("Dagger")


text:
using LinearAlgebra: checksquare
using LinearAlgebra.BLAS: gemv, gemv!, gemm!, trsm!, axpy!, ger!

#####
##### `svd`
#####

function rrule(::typeof(svd), X::AbstractMatrix{<:Real})
    F = svd(X)
    ∂X = Rule() do Ȳ::NamedTuple{(:U,:S,:V)}
        svd_rev(F, Ȳ.U, Ȳ.S, Ȳ.V)
    end
    return F, ∂X
end

function rrule(::typeof(getproperty), F::SVD, x::Symbol)
    if x === :U
        rule = Ȳ->(U=Ȳ, S=zero(F.S), V=zero(F.V))
    elseif x === :S
        rule = Ȳ->(U=zero(F.U), S=Ȳ, V=zero(F.V))
    elseif x === :V
        rule = Ȳ->(U=zero(F.U), S=zero(F.S), V=Ȳ)
    elseif x === :Vt
        # TODO: This could be made to work, but it'd be a pain
        throw(ArgumentError("Vt is unsupported; use V and transpose the result"))
    end
    update = (X̄::NamedTuple{(:U,:S,:V)}, Ȳ)->_update!(X̄, rule(Ȳ), x)
    return getproperty(F, x), (Rule(rule, update), DNERule())
end

function svd_rev(USV::SVD, Ū::AbstractMatrix, s̄::AbstractVector, V̄::AbstractMatrix)
    # Note: assuming a thin factorization, i.e. svd(A, full=false), which is the default
    U = USV.U
    s = USV.S
    V = USV.V
    Vt = USV.Vt

    k = length(s)
    T = eltype(s)
    F = T[i == j ? 1 : inv(@inbounds s[j]^2 - s[i]^2) for i = 1:k, j = 1:k]

    # We do a lot of matrix operations here, so we'll try to be memory-friendly and do
    # as many of the computations in-place as possible. Benchmarking shows that the in-
    # place functions here are significantly faster than their out-of-place, naively
    # implemented counterparts, and allocate no additional memory.
    Ut = U'
    FUᵀŪ = _mulsubtrans!(Ut*Ū, F)  # F .* (UᵀŪ - ŪᵀU)
    FVᵀV̄ = _mulsubtrans!(Vt*V̄, F)  # F .* (VᵀV̄ - V̄ᵀV)
    ImUUᵀ = _eyesubx!(U*Ut)        # I - UUᵀ
    ImVVᵀ = _eyesubx!(V*Vt)        # I - VVᵀ

    S = Diagonal(s)
    S̄ = Diagonal(s̄)

    Ā = _add!(U * FUᵀŪ * S, ImUUᵀ * (Ū / S)) * Vt
    _add!(Ā, U * S̄ * Vt)
    _add!(Ā, U * _add!(S * FVᵀV̄ * Vt, (S \ V̄') * ImVVᵀ))

    return Ā
end

#####
##### `cholesky`
#####

function rrule(::typeof(cholesky), X::AbstractMatrix{<:Real})
    F = cholesky(X)
    ∂X = Rule(Ȳ->chol_blocked_rev(Matrix(Ȳ), Matrix(F.U), 25, true))
    return F, ∂X
end

function rrule(::typeof(getproperty), F::Cholesky, x::Symbol)
    if x === :U
        if F.uplo === 'U'
            ∂F = Ȳ->UpperTriangular(Ȳ)
        else
            ∂F = Ȳ->LowerTriangular(Ȳ')
        end
    elseif x === :L
        if F.uplo === 'L'
            ∂F = Ȳ->LowerTriangular(Ȳ)
        else
            ∂F = Ȳ->UpperTriangular(Ȳ')
        end
    end
    return getproperty(F, x), (Rule(∂F), DNERule())
end

# See "Differentiation of the Cholesky decomposition" (Murray 2016), pages 5-9 in particular,
# for derivations. Here we're implementing the algorithms and their transposes.

"""
    level2partition(A::AbstractMatrix, j::Integer, upper::Bool)

Returns views to various bits of the lower triangle of `A` according to the
`level2partition` procedure defined in [1] if `upper` is `false`. If `upper` is `true` then
the transposed views are returned from the upper triangle of `A`.

[1]: "Differentiation of the Cholesky decomposition", Murray 2016
"""
function level2partition(A::AbstractMatrix, j::Integer, upper::Bool)
    n = checksquare(A)
    @boundscheck checkbounds(1:n, j)
    if upper
        r = view(A, 1:j-1, j)
        d = view(A, j, j)
        B = view(A, 1:j-1, j+1:n)
        c = view(A, j, j+1:n)
    else
        r = view(A, j, 1:j-1)
        d = view(A, j, j)
        B = view(A, j+1:n, 1:j-1)
        c = view(A, j+1:n, j)
    end
    return r, d, B, c
end

"""
    level3partition(A::AbstractMatrix, j::Integer, k::Integer, upper::Bool)

Returns views to various bits of the lower triangle of `A` according to the
`level3partition` procedure defined in [1] if `upper` is `false`. If `upper` is `true` then
the transposed views are returned from the upper triangle of `A`.

[1]: "Differentiation of the Cholesky decomposition", Murray 2016
"""
function level3partition(A::AbstractMatrix, j::Integer, k::Integer, upper::Bool)
    n = checksquare(A)
    @boundscheck checkbounds(1:n, j)
    if upper
        R = view(A, 1:j-1, j:k)
        D = view(A, j:k, j:k)
        B = view(A, 1:j-1, k+1:n)
        C = view(A, j:k, k+1:n)
    else
        R = view(A, j:k, 1:j-1)
        D = view(A, j:k, j:k)
        B = view(A, k+1:n, 1:j-1)
        C = view(A, k+1:n, j:k)
    end
    return R, D, B, C
end

"""
    chol_unblocked_rev!(Ā::AbstractMatrix, L::AbstractMatrix, upper::Bool)

Compute the reverse-mode sensitivities of the Cholesky factorization in an unblocked manner.
If `upper` is `false`, then the sensitivites are computed from and stored in the lower triangle
of `Ā` and `L` respectively. If `upper` is `true` then they are computed and stored in the
upper triangles. If at input `upper` is `false` and `tril(Ā) = L̄`, at output
`tril(Ā) = tril(Σ̄)`, where `Σ = LLᵀ`. Analogously, if at input `upper` is `true` and
`triu(Ā) = triu(Ū)`, at output `triu(Ā) = triu(Σ̄)` where `Σ = UᵀU`.
"""
function chol_unblocked_rev!(Σ̄::AbstractMatrix{T}, L::AbstractMatrix{T}, upper::Bool) where T<:Real
    n = checksquare(Σ̄)
    j = n
    @inbounds for _ in 1:n
        r, d, B, c = level2partition(L, j, upper)
        r̄, d̄, B̄, c̄ = level2partition(Σ̄, j, upper)

        # d̄ <- d̄ - c'c̄ / d.
        d̄[1] -= dot(c, c̄) / d[1]

        # [d̄ c̄'] <- [d̄ c̄'] / d.
        d̄ ./= d
        c̄ ./= d

        # r̄ <- r̄ - [d̄ c̄'] [r' B']'.
        r̄ = axpy!(-Σ̄[j,j], r, r̄)
        r̄ = gemv!(upper ? 'n' : 'T', -one(T), B, c̄, one(T), r̄)

        # B̄ <- B̄ - c̄ r.
        B̄ = upper ? ger!(-one(T), r, c̄, B̄) : ger!(-one(T), c̄, r, B̄)
        d̄ ./= 2
        j -= 1
    end
    return (upper ? triu! : tril!)(Σ̄)
end

function chol_unblocked_rev(Σ̄::AbstractMatrix, L::AbstractMatrix, upper::Bool)
    return chol_unblocked_rev!(copy(Σ̄), L, upper)
end

"""
    chol_blocked_rev!(Σ̄::AbstractMatrix, L::AbstractMatrix, nb::Integer, upper::Bool)

Compute the sensitivities of the Cholesky factorization using a blocked, cache-friendly
procedure. `Σ̄` are the sensitivities of `L`, and will be transformed into the sensitivities
of `Σ`, where `Σ = LLᵀ`. `nb` is the block size to use. If the upper triangle has been used
to represent the factorization, that is `Σ = UᵀU` where `U := Lᵀ`, then this should be
indicated by passing `upper = true`.
"""
function chol_blocked_rev!(Σ̄::AbstractMatrix{T}, L::AbstractMatrix{T}, nb::Integer, upper::Bool) where T<:Real
    n = checksquare(Σ̄)
    tmp = Matrix{T}(undef, nb, nb)
    k = n
    if upper
        @inbounds for _ in 1:nb:n
            j = max(1, k - nb + 1)
            R, D, B, C = level3partition(L, j, k, true)
            R̄, D̄, B̄, C̄ = level3partition(Σ̄, j, k, true)

            C̄ = trsm!('L', 'U', 'N', 'N', one(T), D, C̄)
            gemm!('N', 'N', -one(T), R, C̄, one(T), B̄)
            gemm!('N', 'T', -one(T), C, C̄, one(T), D̄)
            chol_unblocked_rev!(D̄, D, true)
            gemm!('N', 'T', -one(T), B, C̄, one(T), R̄)
            if size(D̄, 1) == nb
                tmp = axpy!(one(T), D̄, transpose!(tmp, D̄))
                gemm!('N', 'N', -one(T), R, tmp, one(T), R̄)
            else
                gemm!('N', 'N', -one(T), R, D̄ + D̄', one(T), R̄)
            end

            k -= nb
        end
        return triu!(Σ̄)
    else
        @inbounds for _ in 1:nb:n
            j = max(1, k - nb + 1)
            R, D, B, C = level3partition(L, j, k, false)
            R̄, D̄, B̄, C̄ = level3partition(Σ̄, j, k, false)

            C̄ = trsm!('R', 'L', 'N', 'N', one(T), D, C̄)
            gemm!('N', 'N', -one(T), C̄, R, one(T), B̄)
            gemm!('T', 'N', -one(T), C̄, C, one(T), D̄)
            chol_unblocked_rev!(D̄, D, false)
            gemm!('T', 'N', -one(T), C̄, B, one(T), R̄)
            if size(D̄, 1) == nb
                tmp = axpy!(one(T), D̄, transpose!(tmp, D̄))
                gemm!('N', 'N', -one(T), tmp, R, one(T), R̄)
            else
                gemm!('N', 'N', -one(T), D̄ + D̄', R, one(T), R̄)
            end

            k -= nb
        end
        return tril!(Σ̄)
    end
end

function chol_blocked_rev(Σ̄::AbstractMatrix, L::AbstractMatrix, nb::Integer, upper::Bool)
    return chol_blocked_rev!(copy(Σ̄), L, nb, upper)
end


text:
using ReachabilityAnalysis
using ReachabilityAnalysis: solve
using JuMP, Ipopt
using Plots
using LaTeXStrings

# IVP parameters
A = [0.0 1; -4π 0]
X0 = BallInf([0.0, 10.0], 0.1)
X = Universe(2)

# discretization parameters
δ = 0.01
b = [1.0, -1.0]
B = hcat(b)
u1 = -3.0

# varying input domains
Us_raw = [
          Singleton([100 * u1]),
          Singleton([u1]),
          BallInf([0.0], abs(u1)),
          BallInf([u1], 1.0)
         ]

Us = [
      (ConstantInput(linear_map(B, Us_raw[1])), "{$(Us_raw[1].element[1])}"),
      (ConstantInput(linear_map(B, Us_raw[2])), "{$u1}"),
      (ConstantInput(linear_map(B, Us_raw[3])), "□(0, $(abs(u1)))"),
      (ConstantInput(linear_map(B, Us_raw[4])), "□($u1, 0.1)")
     ]

# methods
algs = [
        (FirstOrderZonotope(), "Zonotope"),
        (CorrectionHull(order=4), "Correction hull"),
        (FirstOrder(), "First order"),
        (ForwardBackward(solver=Ipopt.Optimizer()), "Forward/backward"),
        (Forward(), "Forward")
       ]

# results
U2Ω0s = []
for (U, name_U) in Us
    local Ω0s = []
    local P = @ivp(x' = A*x + u, x ∈ X, u ∈ U, x(0) ∈ X0)
    for (alg, name_alg) in algs
        local Pd = discretize(P, δ, alg)
        local Ω0 = Pd.x0
        push!(Ω0s, Ω0)
    end
    push!(U2Ω0s, Ω0s)
end

# plot sets
colors = [:orange, :green, :blue, :purple, :red, :black]
plots = []
for (i, (U_raw, (U, name_U))) in enumerate(zip(Us_raw, Us))
    local fig = plot(leg=:bottomleft, xlab="U = B·$name_U")
    for (j, (alg, name_alg)) in enumerate(algs)
        local Ω0 = U2Ω0s[i][j]
        if name_alg == "Forward/backward"
            Ω0 = overapproximate(Ω0, PolarDirections(30))  # workaround because σ is not available
        end
        # j+1 because d/dt is missing
        plot!(Ω0, alpha=0, color=colors[j+1], linecolor=colors[j+1], linealpha=1, lw=2, lab="Ω₀ $(name_alg)")
    end

    local P = @ivp(x' = A*x + u, x ∈ X, u ∈ U, x(0) ∈ X0)
    sol = solve(P, δ, alg=LGG09(δ=δ/10, dirs=collect(PolarDirections(40))))
    alpha = i == 1 ? 0.2 : 0.05  # first plot should use a darker gray
    plot!(sol, vars=(1, 2), color=:gray, alpha=alpha, linealpha=0, lab="")

    savefig("experiment_U_$i.pdf")
    push!(plots, fig)
end

return plots


text:
reservedWords = Dict{String, Any}()
reservedWords["SentenceType"] = Set{AbstractString}()
push!(reservedWords["SentenceType"], "is")
push!(reservedWords["SentenceType"], "react")
push!(reservedWords["SentenceType"], "uptake")
push!(reservedWords["SentenceType"], "secrete")
push!(reservedWords["SentenceType"], "bind")
push!(reservedWords["SentenceType"], "unbind")
push!(reservedWords["SentenceType"], "phosphorylate")
push!(reservedWords["SentenceType"], "dephosphorylate")
push!(reservedWords["SentenceType"], "induce")
push!(reservedWords["SentenceType"], "activate")
push!(reservedWords["SentenceType"], "repress")
push!(reservedWords["SentenceType"], "inhibit")
push!(reservedWords["SentenceType"], "catalyze")
reservedWords["the"] = "negligible"
reservedWords["The"] = "negligible"
reservedWords["cell"] = "negligible"
reservedWords["together"] = "negligible"
reservedWords["to"] = "negligible"
reservedWords["form"] = "negligible"
reservedWords["expression"] = "negligible"
reservedWords["of"] = "negligible"
reservedWords["transcription"] = "negligible"
reservedWords["translation"] = "negligible"
reservedWords["conversion"] = "negligible"
reservedWords["is"] = "is"
reservedWords["was"] = "is"
reservedWords["are"] = "is"
reservedWords["were"] = "is"
reservedWords["react"] = "react"
reservedWords["reacted"] = "react"
reservedWords["reacting"] = "react"
reservedWords["reacts"] = "react"
reservedWords["into"] = "into"
reservedWords["irreversible"] = "irreversible"
reservedWords["irreversibly"] = "irreversible"
reservedWords["uptake"] = "uptake"
reservedWords["uptakes"] = "uptake"
reservedWords["uptaking"] = "uptake"
reservedWords["uptook"] = "uptake"
reservedWords["via"] = "via"
reservedWords["through"] = "via"
reservedWords["the"] = "negligible"
reservedWords["cell"] = "negligible"
reservedWords["secrete"] = "secrete"
reservedWords["secretes"] = "secrete"
reservedWords["secreting"] = "secrete"
reservedWords["secreted"] = "secrete"
reservedWords["via"] = "via"
reservedWords["through"] = "via"
reservedWords["the"] = "negligible"
reservedWords["cell"] = "negligible"
reservedWords["bind"] = "bind"
reservedWords["bound"] = "bind"
reservedWords["binds"] = "bind"
reservedWords["binding"] = "bind"
reservedWords["associate"] = "bind"
reservedWords["associates"] = "bind"
reservedWords["associated"] = "bind"
reservedWords["associating"] = "bind"
reservedWords["complex"] = "bind"
reservedWords["complexes"] = "bind"
reservedWords["complexed"] = "bind"
reservedWords["complexing"] = "bind"
reservedWords["into"] = "into"
reservedWords["together"] = "negligible"
reservedWords["unbind"] = "unbind"
reservedWords["unbinds"] = "unbind"
reservedWords["unbound"] = "unbind"
reservedWords["unbinding"] = "unbind"
reservedWords["dissociate"] = "unbind"
reservedWords["dissociates"] = "unbind"
reservedWords["dissociated"] = "unbind"
reservedWords["dissociating"] = "unbind"
reservedWords["into"] = "into"
reservedWords["phosphorylate"] = "phosphorylate"
reservedWords["phosphorylates"] = "phosphorylate"
reservedWords["phosphorylated"] = "phosphorylate"
reservedWords["phosphorylating"] = "phosphorylate"
reservedWords["reversibly"] = "reversible"
reservedWords["reversible"] = "reversible"
reservedWords["at"] = "at"
reservedWords["site"] = "site"
reservedWords["dephosphorylate"] = "dephosphorylate"
reservedWords["dephosphorylates"] = "dephosphorylate"
reservedWords["dephosphorylated"] = "dephosphorylate"
reservedWords["dephosphorylating"] = "dephosphorylate"
reservedWords["at"] = "at"
reservedWords["site"] = "site"
reservedWords["activate"] = "activate"
reservedWords["activates"] = "activate"
reservedWords["activating"] = "activate"
reservedWords["activated"] = "activate"
reservedWords["promote"] = "induce"
reservedWords["promotes"] = "induce"
reservedWords["promoting"] = "induce"
reservedWords["promoted"] = "induce"
reservedWords["induce"] = "induce"
reservedWords["induces"] = "induce"
reservedWords["inducing"] = "induce"
reservedWords["induced"] = "induce"
reservedWords["upregulate"] = "activate"
reservedWords["upregulates"] = "activate"
reservedWords["upregulated"] = "activate"
reservedWords["upregulating"] = "activate"
reservedWords["the"] = "negligible"
reservedWords["expression"] = "negligible"
reservedWords["of"] = "negligible"
reservedWords["transcription"] = "negligible"
reservedWords["translation"] = "negligible"
reservedWords["inhibit"] = "inhibit"
reservedWords["inhibits"] = "inhibit"
reservedWords["inhibiting"] = "inhibit"
reservedWords["inhibited"] = "inhibit"
reservedWords["repress"] = "repress"
reservedWords["represses"] = "repress"
reservedWords["repressing"] = "repress"
reservedWords["repressed"] = "repress"
reservedWords["downregulate"] = "repress"
reservedWords["downregulates"] = "repress"
reservedWords["downregulated"] = "repress"
reservedWords["downregulating"] = "repress"
reservedWords["the"] = "negligible"
reservedWords["expression"] = "negligible"
reservedWords["of"] = "negligible"
reservedWords["transcription"] = "negligible"
reservedWords["translation"] = "negligible"
reservedWords["catalyze"] = "catalyze"
reservedWords["catalyzes"] = "catalyze"
reservedWords["catalyzed"] = "catalyze"
reservedWords["catalyzing"] = "catalyze"
reservedWords["catalyse"] = "catalyze"
reservedWords["catalyses"] = "catalyze"
reservedWords["catalysed"] = "catalyze"
reservedWords["catalysing"] = "catalyze"
reservedWords["reversible"] = "reversible"
reservedWords["into"] = "into"
reservedWords["the"] = "negligible"
reservedWords["conversion"] = "negligible"
reservedWords["of"] = "negligible"


text:
import Flux:gpu
using CuArrays

gpu(x::SubArray) = CuArray{Float32}(x)

const is_using_gpu = true


text:
const ϵ = 1.0e-6

function decibel(value::Real)
    return 10 * log10(value)
end

# Test odd-length signal
@testset "summarize PSD" begin
    N = 90000
    fs = 100.
    num_segments = 13
    smoothing_width_factor = 2.
    T = Float64
    data = Array{T, 2}(undef, num_segments, N)
    for i in 1:num_segments
        for j in 1:N
            @inbounds data[i, j] = i
        end
    end

    expect_mean = decibel(7.)
    expect_min = decibel(1.)
    expect_max = decibel(13.)
    expect_median = decibel(7.)
    result, center_periods = summarize_psd(data, fs, smoothing_width_factor)
    @test eltype(result) == T
    @test eltype(center_periods) == T
    @test size(result[1, :]) == size(center_periods)
    @test result[1, 1] - expect_mean <= ϵ
    @test result[2, 1] - expect_min <= ϵ
    @test result[3, 1] - expect_max <= ϵ
    @test result[4, 1] - expect_median <= ϵ
end

# Test the lenght and type of frequency and period arrays
@testset "get frequency and periods array" begin
    N = 900 # 900 seconds, i.e. 15 minutes
    fs = 100.
    smoothing_width_factor = 2.
    T = Float64

    left_freqs, right_freqs, center_periods = get_freqs_and_periods(fs, N, smoothing_width_factor)

    expect_length = 98
    @test eltype(left_freqs) == T
    @test eltype(right_freqs) == T
    @test eltype(center_periods) == T
    @test length(left_freqs) == expect_length
    @test length(right_freqs) == expect_length
    @test length(center_periods) == expect_length
end


text:
using Printf
using StringEncodings
using Base.Filesystem

mutable struct  CSVinfo
     fieldname ::String
     signalStringName::String
     keyFieldNames::String
     fieldNameDict::Dict
     Area2FiberLen
     # PosNames
     function CSVinfo(;fieldname::String ="定位界标1	定位界标2	xxxx主机IP	通道号	定位界标1光程	定位界标2光程	子防区编号	防区起点GPS经度坐标	防区起点GPS纬度坐标	防区终点GPS经度坐标	防区终点GPS纬度坐标	左右岸	地名信息	界标之间围栏长度	管理处	分局	分局编号	管理处编号	类型	主机编号	防区编号	通道号	归属IP",
          signalStringName::String="定位界标1光程",
          keyFieldNames::String ="定位界标1光程	定位界标2光程	地名信息	管理处编号	防区编号	通道号	归属IP	定位界标1	定位界标2"
           )
          Area2FiberLen=[ ]
          new(fieldname,signalStringName,keyFieldNames ,Dict(),Area2FiberLen)
     end

end

#
function initDict!(CSVinfo1::CSVinfo)
     CSVinfo1.fieldNameDict = Dict()
     ts = split(CSVinfo1.fieldname,"\t" )
     for t = 1: length( ts )
          CSVinfo1.fieldNameDict[ ts[t] ] = 0
     end
     # return fieldNameDict
end

function readinCSVFile!(CSVinfo1::CSVinfo;csvfile="../Jl-Aux/glbCfg/org/全局-xxxx主机-防区信息.csv")
     initDict!(CSVinfo1 )
     fields= split( CSVinfo1.keyFieldNames,"\t")
     fieldlen =length( fields)
     # csvfileAbsname = csvfile

     csvfileAbsname =joinpath(pwd(), csvfile)
     fin =nothing
     while true
          println( csvfileAbsname)
          fin = open(csvfileAbsname,"r" )
         # @info "here ",fin
          i1 = read(fin,UInt8); #@info "here ", i1
          i2= read(fin,UInt8);#@info "here ", i2
          i3 = read(fin,UInt8);#@info "here ", i3
          if i1 ==  UInt8(239) &&  i2 == UInt8(187) && i3 == UInt8(191)
               # println( "UTF8 -BOM file ")
               break
          else
               close(fin )
               fin = open(csvfileAbsname, enc"GBK", "r");
               # println( "GBK file ")
               break
          end
     end
     CSVinfo1.Area2FiberLen =[]
     for line in readlines(fin)
          ts = split(line,"," )
          if occursin(CSVinfo1.signalStringName , line) ## 首行
               for t = 1: length( ts )
                    CSVinfo1.fieldNameDict[ ts[t] ] = t
               end
          else
               a=[]

               for i =1:fieldlen
                    field = fields[i]
                    idx = CSVinfo1.fieldNameDict[ field]
                    val  = ts[ idx ]
                    append!(a,[val])
               end
               append!(CSVinfo1.Area2FiberLen,[ a ])
          end
     end
end
function show(CSVinfo1::CSVinfo )
     fields= split( CSVinfo1.keyFieldNames,"\t")
     fieldlen =length( fields)
     for  row  in CSVinfo1.Area2FiberLen
          for i =1:fieldlen
               fieldname = fields[ i ]
               print("$fieldname:$(row[i]) ,")
          end
          println()
     end
end

function getrow(CSVinfo1,归属IP,strchnid,strAreaID )
     for rowobj in CSVinfo1.Area2FiberLen
     #"定位界标1光程	定位界标2光程	地名信息	管理处编号	防区编号	通道号	归属IP"
          if rowobj[ 7 ] == 归属IP && rowobj[ 6 ] == strchnid
               # println( rowobj[ 9 ][ end-2:end]  )
               if rowobj[ 5 ] == strAreaID
                    return rowobj
               elseif rowobj[ 9 ][ end-2:end] == strAreaID[2:end]  # 末尾的标牌是该管理处
                    return rowobj
               end
          end
     end
  return

end

function  outfile(CSVinfo1::CSVinfo;归属IP)
 tgrFolder = Base.Filesystem.dirname( csvfile1  )
 for chnid = 1:2
   filename ="defence-" * 归属IP *"Sheet "*string( chnid)*".csv"
   tgrfullname = Base.Filesystem.joinpath(tgrFolder, filename )
   fobj = open( tgrfullname ,enc"GBK","w")
   write(fobj ,"编号,通道,起点光程,防区号,防区描述,使能\n")
fiberpos1 =""
idx =0
地名信息 =""
  for AreaID = 1:10000
                idx = AreaID
	strAreaID = @sprintf("_%03d", AreaID )
	rowDictObj = getrow(CSVinfo1,归属IP,string( chnid),strAreaID )
	if  rowDictObj  == nothing
		break
	end
	# if rowDictObj[ 3 ]   == nothing
	# 	地名信息 =""
	# else
		地名信息 = rowDictObj[ 3 ]
     # end
     if rowDictObj[ 5 ] == strAreaID
          fiberpos1 = rowDictObj[ 1 ]
          rowStr = @sprintf("%04d,%d,%s,_%03d,%s,1\n",idx ,chnid, fiberpos1 ,  AreaID  , 地名信息  )
     else
          fiberpos1 = rowDictObj[ 2 ]
          rowStr = @sprintf("%04d,%d,%s,_%03d,%s,1\n",idx ,chnid, fiberpos1 ,  AreaID  , ""  )
     end
	write( fobj , rowStr )
  end

   close(fobj  )
end
end

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

csvfile1 = "../Jl-Aux/glbCfg/org/全局-xxxx主机-防区信息.csv"
CSVinfoObj = CSVinfo()
readinCSVFile!(CSVinfoObj ,csvfile=csvfile1);
#show( CSVinfoObj )
IPs= "10.17.171.41 10.17.174.41 10.17.178.41 10.17.36.41 10.17.39.41 10.17.46.41 10.17.5.41 10.18.132.41 10.18.135.41 10.18.136.41 10.19.138.41 10.19.139.41 10.19.141.41 10.19.164.41"
for  归属IP1 in split(IPs," ")
	#归属IP1= "10.19." * string(A4) * ".41"
                println(  归属IP1 )
	outfile( CSVinfoObj ,归属IP= 归属IP1)
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
module Minibatch

export VectorBatch, MaskedBatch, ϵ

using NNlib
using Flux

################
# NN functions #
################

const ϵ = 1e-7

function NNlib.softmax(x::AbstractArray, axis=1)
    out = exp.(x .- maximum(x, axis))
    out ./= sum(out, axis) .+ ϵ
    return out
end

#################
# AxisInfo type #
#################

# TODO non-varlen axes should have (optional?) static dimension
const AxisInfo{N} = NTuple{N, Bool}

###############
# Batch types #
###############

# currently it's intentionally not <:AbstractArray
# because that provides a lot of magic I want to avoid
abstract type AbstractBatch{T, A} end
axisinfo(::AbstractBatch{T, A}) where {T, A} = A

struct VectorBatch{T, A} <: AbstractBatch{T, A}
    data::Vector{T}
end
VectorBatch(data::Vector{T}, axes::AxisInfo) where {T} = VectorBatch{T, axes}(data)

Base.length(b::VectorBatch) = length(b.data)
Base.:(==)(a::VectorBatch{T1, A}, b::VectorBatch{T2, A}) where {T1, T2, A} = a.data == b.data
Base.isapprox(a::VectorBatch{T1, A}, b::VectorBatch{T2, A}) where {T1, T2, A} = a.data ≈ b.data

# TODO unsure if it should be mutable
mutable struct MaskedBatch{T, A} <: AbstractBatch{T, A}
    data
    mask
end
function MaskedBatch(data::AbstractVector, axes::AxisInfo)
    dims = (a ? maximum(size(v, d) for v in data) : size(first(data), d) for (d, a) in enumerate(axes))
    maskdims = (a ? maximum(size(v, d) for v in data) : 1 for (d, a) in enumerate(axes))
    batch = fill!(similar(first(data), dims..., length(data)), 0)
    mask = fill!(similar(first(data), maskdims..., length(data)), 0)
    for (i, x) in enumerate(data)
        setindex!(batch, x, indices(x)..., i)
        mask[(1:(a ? size(x, d) : 1) for (d, a) in enumerate(axes))..., i] .= 1
    end
    return MaskedBatch{typeof(first(data)), axes}(batch, mask)
end
function sizes(b::MaskedBatch{T, A}) where {T, A}
    bs = length(b)
    mask = Array{Int}(b.mask)
    return vcat((a ? sum(mask, d)[(1 for _ in A)..., :]' :
                 fill(size(b.data, d), bs)' for (d, a) in enumerate(A))...)
end
function VectorBatch(b::MaskedBatch{T, A}) where {T, A}
    s = sizes(b)
    xs = [b.data[(1:n for n in s[:, i])..., i] for i in 1:length(b)]
    return VectorBatch{T, A}(xs)
end

MaskedBatch(b::VectorBatch{T, A}) where {T, A} = MaskedBatch(b.data, A)
Base.length(b::MaskedBatch) = last(size(b.mask))
Base.:(==)(a::MaskedBatch{T1, A}, b::MaskedBatch{T2, A}) where {T1, T2, A} = a.data == b.data && a.mask == b.mask
Base.isapprox(a::MaskedBatch{T1, A}, b::MaskedBatch{T2, A}) where {T1, T2, A} = a.data ≈ b.data && a.mask == b.mask

# TODO also CatBatch using CatArrays.jl?

####################################
# Methods/overdubs for batch types #
####################################

function _samebatchsizes(xs::Vararg{AbstractBatch})
    bs = length(first(xs))
    all(length(x) == bs for x in xs) || error("batch size mismatch")
    return bs
end
function _vbcall(f, T, xs...)
    bs = _samebatchsizes((x for x in xs if x isa T)...)
    return T([f((x isa T ? x.data[i] : x for x in xs)...) for i in 1:bs])
end
Base.:+(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:VectorBatch = _vbcall(+, T, xs...)
Base.:-(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:VectorBatch = _vbcall(-, T, xs...)
Base.:*(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:VectorBatch = _vbcall(*, T, xs...)
Base.:/(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:VectorBatch = _vbcall(/, T, xs...)
Base.getindex(x::T, inds...) where T<:VectorBatch = _vbcall(getindex, T, x, inds...)
Base.broadcast(f, xs::Vararg{Union{T, AbstractArray, Number}}) where T<:VectorBatch = _vbcall(broadcast, T, f, xs...)
function Base.transpose(b::VectorBatch{T, A}) where {T, A}
    axes = tuple(A[2], A[1], A[3:end]...)
    return VectorBatch{T, axes}([transpose(x) for x in b.data])
end
softmax(x::T, axis=1) where T<:VectorBatch = _vbcall(softmax, T, x, axis)

# _vbcall is unsafe for contractions because it doesn't check that axes
# which must have static dimension actually do at the type level. We
# could fix that by including explicit checks, as in the methods below,
# but even when _vbcall wouldn't break.
function Base.:*(a::VectorBatch{T, A}, b::AbstractVector) where {T<:AbstractMatrix, A}
    A[2] == false || error("cannot contract axes with static and dynamic dimension")
    return VectorBatch{T, (A[1],)}([a.data[i] * b for i in 1:length(a)])
end

function Base.:*(a::VectorBatch{T, A}, b::VectorBatch{T, B}) where {T<:AbstractMatrix, A, B}
    A[2] == B[1] || error("cannot contract axes with static and dynamic dimension")
    bs = _samebatchsizes(a, b)
    return VectorBatch{T, (A[1], B[2])}([a.data[i] * b.data[i] for i in 1:bs])
end

function Base.:*(a::VectorBatch{T1, A}, b::VectorBatch{T2, B}) where {T1<:AbstractMatrix, T2<:AbstractVector, A, B}
    A[2] == B[1] || error("cannot contract axes with static and dynamic dimension")
    bs = _samebatchsizes(a, b)
    return VectorBatch{T, (A[1],)}([a.data[i] * b.data[i] for i in 1:bs])
end

_rezero!(b::MaskedBatch) = b.data .*= b.mask
# TODO this should really be promote_containersize
function _samemasks(xs::Vararg{<:MaskedBatch{T2, B}}) where {T2, B}
    fs = first(xs).mask
    if any(B)
        all(x.mask == fs for x in xs) || error("MaskedBatch size mismatch")
    end
    return fs
end
function Base.broadcast(f, xs::Vararg{Union{T, AbstractArray, Number}}) where T<:MaskedBatch
    mask = _samemasks((x for x in xs if x isa T)...)
    res = T(broadcast(f, (x isa T ? x.data : x for x in xs)...), mask)
    _rezero!(res)
    return res
end
Base.:+(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:MaskedBatch = broadcast(+, xs...)
Base.:-(xs::Vararg{Union{T, AbstractArray, Number}}) where T<:MaskedBatch = broadcast(-, xs...)

# contractions between axes with static dimension

function Base.dot(a::AbstractVector, b::MaskedBatch{T, (false,)}) where {T}
    return MaskedBatch{T, ()}(b.data'*a, b.mask[1])
end
function Base.dot(a::MaskedBatch{T, (false,)}, b::AbstractVector) where {T}
    return MaskedBatch{T, ()}(a.data'*b, a.mask[1])
end

function Base.:*(a::AbstractMatrix, b::MaskedBatch{T, (false,)}) where {T}
    return MaskedBatch{T, (false,)}(a * b.data, b.mask)
end
function Base.:*(a::MaskedBatch{T, A}, b::AbstractVector) where {T<:AbstractMatrix, A}
    A[2] == false || error("cannot contract axes with static and dynamic dimension")
    s1, s2, bs = size(a.data)
    data = reshape(reshape(a.data, s1, s2*bs) * b, s1, bs)
    return MaskedBatch{T, (A[1],)}(data, a.mask[:, 1, :])
end

function Base.:*(a::AbstractMatrix, b::MaskedBatch{T, B}) where {T<:AbstractMatrix, B}
    B[1] == false || error("cannot contract axes with static and dynamic dimension")
    s1, s2, bs = size(b.data)
    data = reshape(a * reshape(b.data, s1, s2*bs), size(a, 1), s2, bs)
    return MaskedBatch{T, B}(data, b.mask)
end
function Base.:*(a::MaskedBatch{T, A}, b::AbstractMatrix) where {T<:AbstractMatrix, A}
    A[2] == false || error("cannot contract axes with static and dynamic dimension")
    s1, s2, bs = size(a.data)
    amat = reshape(permutedims(a.data, [2, 3]), s1*bs, s2)
    data = permutedims(reshape(amat * b, s1, bs, size(b, 2)), [2, 3])
    return MaskedBatch{T, A}(data, a.mask)
end

# contractions between axes with dynamic dimension

# TODO should things work when dynamic dimension = 0?
function Base.dot(a::MaskedBatch{T, A}, b::MaskedBatch{T, B}) where {T<:AbstractVector, A, B}
    data = sum(a.data .* b.data, 1)[1, :] # TODO replace with batchedmul when fast
    return MaskedBatch{T, ()}(data, a.mask[:, 1])
end

# the methods below use batched matrix-matrix and matrix-vector products;
# these are provided as BLAS extensions by cuBLAS and MKL but aren't
# wrapped in Base; here we implement sequential fallbacks but TODO these
# functions should actually be in another package and have those fast methods
function batchedmul(a::AbstractArray{T, 3}, b::AbstractArray{T, 2}) where {T}
    (bs = size(a, 3)) == size(b, 2) || error("batch size mismatch")
    res = similar(a, size(a, 1), bs)
    for i in 1:bs
        res[:, i] = a[:, :, i] * b[:, i]
    end
    return res
end
function batchedmul(a::AbstractArray{T, 3}, b::AbstractArray{T, 3}) where {T}
    (bs = size(a, 3)) == size(b, 3) || error("batch size mismatch")
    res = similar(a, size(a, 1), size(b, 2), bs)
    for i in 1:bs
        res[:, :, i] = a[:, :, i] * b[:, :, i]
    end
    return res
end

function Base.:*(a::MaskedBatch{T1, A}, b::MaskedBatch{T2, B}) where {T1<:AbstractMatrix, T2<:AbstractVector, A, B}
    A[2] == B[1] || error("cannot contract axes with static and dynamic dimension")
    data = batchedmul(a.data, b.data)
    mask = batchedmul(a.mask[:, 1:1, :], b.mask[1:1, :])
    return MaskedBatch{T2, (A[1],)}(data, mask)
end

function Base.:*(a::MaskedBatch{T, A}, b::MaskedBatch{T, B}) where {T<:AbstractMatrix, A, B}
    A[2] == B[1] || error("cannot contract axes with static and dynamic dimension")
    data = batchedmul(a.data, b.data)
    mask = batchedmul(a.mask[:, 1:1, :], b.mask[1:1, :, :])
    return MaskedBatch{T, (A[1], B[2])}(data, mask)
end

function NNlib.softmax(x::T, axis=1) where T<:MaskedBatch
    data = exp.(x.data .- maximum(x.data, axis)) .* x.mask
    data ./= sum(data, axis) .+ ϵ
    return T(data, x.mask)
end

function Base.transpose(x::MaskedBatch{T, A}) where {T, A}
    dims = collect(1:ndims(x.data))
    dims[1:2] = [2, 1]
    data = permutedims(x.data, dims)
    mask = permutedims(x.mask, dims)
    axes = tuple(A[2], A[1], A[3:end]...)
    return MaskedBatch{T, axes}(data, mask)
end

# generalize this
function Base.getindex(a::AbstractMatrix, ::Colon, b::MaskedBatch{T, B}) where {T<:AbstractVector, B}
    data = b.data .+ (1 .- b.mask)
    data = a[:, data]
    mask = reshape(b.mask, 1, size(b.mask)...)
    axes = tuple(false, B...)
    # TODO it's a nice coincidence that this typeof works here; it won't in other methods
    return MaskedBatch{typeof(a), axes}(data .* mask, mask)
end

function Base.mapslices(f, x::MaskedBatch{T, A}, dim::Int) where {T, A}
    A[dim] && error("reduction along dynamic dimension not implemented")
    return MaskedBatch{T, A}(mapslices(f, x.data, dim) .* x.mask, x.mask)
end

Base.mean(x::MaskedBatch, dim::Int) = mapslices(mean, x, dim)
Base.std(x::MaskedBatch, dim::Int) = mapslices(std, x, dim)
Base.sum(x::MaskedBatch, dim::Int) = mapslices(sum, x, dim)
Base.prod(x::MaskedBatch, dim::Int) = mapslices(prod, x, dim)

function Flux.normalise(x::MaskedBatch{<:AbstractVecOrMat})
  μ′ = mean(x, 1)
  σ′ = std(x, 1)
  return (x .- μ′) ./ (σ′ .+ ϵ)
end

end # module


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    mcismol, mcismolsize,
    mcesmol, mcesmolsize

"""
    mcismol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...
        ) -> Tuple{Dict{Int,Int},Symbol}

Compute maximum common induced substructure (MCIS) of mol1 and mol2.

## Keyword arguments

- connected(Bool): if true, apply connected MCS constraint.
- topological(Bool): if true, apply topological constraint.
- diameter(Int): distance cutoff for topological constraint.
- tolerance(Int): distance mismatch tolerance for topological constraint.
- timeout(Int): abort calculation and return suboptimal results if the execution
time has reached the given value (default=60, in seconds).
- targetsize(Int): abort calculation and return suboptimal result so far if the
given mcs size achieved.
"""
function mcismol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...)
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return findmcis(
        mol1, mol2, nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

mcismolsize(mol1, mol2; kwargs...) = length(mcismol(mol1, mol2; kwargs...)[1])

"""
    mcesmol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...
        ) -> Tuple{Dict{Int,Int},Symbol}

Compute maximum common edge induced substructure (MCES) of mol1 and mol2.
"""
function mcesmol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...)
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return findmces(
        mol1, mol2, nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

mcesmolsize(mol1, mol2; kwargs...) = length(mcesmol(mol1, mol2; kwargs...)[1])

# TODO: subgraphview


text:
depths = parse.(Int, readlines("day01.in"))
println("Q1: ", count(>(0), diff(depths)))
avgs = depths[1:end-2] .+ depths[2:end-1] .+ depths[3:end]
println("Q2: ", count(>(0), diff(avgs)))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
export splitstrip

function splitstrip{T<:AbstractString}(line::T, splitter::String)
    convert(Array{T,1}, split(strip(line), splitter))
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()

##
using CSV
using DataFrames
using Pipe
using DataStructures
using StatsBase
using ProgressMeter
using Random
Random.seed!(2002261988307380348)

##
using Conda
Conda.pip_interop(true)
Conda.pip("install", "ete3")

ENV["PYTHON"] = ""
Pkg.build("PyCall")
using PyCall

ete3 = pyimport("ete3")

##

try
    mkdir("../data")
catch e
end

try
    mkdir("../data/wals")
catch e
end

languagesF = "../data/wals/languages.csv"
valsF = "../data/wals/values.csv"
paramsF = "../data/wals/parameters.csv"
codesF = "../data/wals/codes.csv"

!(isfile(languagesF) && isfile(valsF) && isfile(paramsF)) && begin
    try
        mkdir("tmp")
    catch e
    end

    wals2020 = download(
        "https://github.com/cldf-datasets/wals/archive/v2020.zip",
        "tmp/wals2020.zip",
    )
    run(`unzip $wals2020 -d tmp/`)
    cp("tmp/wals-2020/cldf/languages.csv", languagesF, force = true)
    cp("tmp/wals-2020/cldf/values.csv", valsF, force = true)
    cp("tmp/wals-2020/cldf/parameters.csv", paramsF, force = true)
    cp("tmp/wals-2020/cldf/codes.csv", codesF, force = true)
    rm("tmp", recursive = true)
end

##

languages = CSV.read(languagesF, DataFrame)

vals = CSV.read(valsF, DataFrame)

params = CSV.read(paramsF, DataFrame)

codes = CSV.read(codesF, DataFrame)
##

woFeatures = ["87A", "86A", "85A", "88A", "89A", "83A", "90A", "82A"]

data = unstack(
    (@pipe vals |>
           filter(x -> x.Parameter_ID ∈ woFeatures, _) |>
           select(_, [:Language_ID, :Parameter_ID, :Value])),
    :Language_ID,
    :Parameter_ID,
    :Value,
)

##

filter!(x -> x.Parameter_ID ∈ woFeatures, codes)

select!(codes, [:ID, :Parameter_ID, :Name, :Number])

filter!(x -> x.Number ∈ [1,2], codes)
##

for i in 1:size(data,1), j in 2:size(data,2)
    v = data[i,j]
    if !ismissing(v) && v > 2
        data[i,j] = missing
    end
end

##
nValues = vec(8 .- mapslices(x -> sum(ismissing.(x)), Array(data), dims=2))
insertcols!(data, 10, :nValues => nValues)

sort!(data, :nValues, rev=true)

##

dropmissing!(languages, :Glottocode)

data = innerjoin(
    data,
    select(languages, [:ID, :Glottocode]),
    on = :Language_ID => :ID,
)

unique!(data, :Glottocode)

filter!(x -> x.nValues >= 6, data)

##

try
    mkdir("../data")
catch e
end

try
    mkdir("../data/asjp")
catch e
end

languagesF = "../data/asjp/languages.csv"

!isfile(languagesF) && begin
    try
        mkdir("tmp")
    catch e
    end

    asjp18 = download(
        "https://zenodo.org/record/3835952/files/lexibank/asjp-v18.zip",
        "tmp/asjp-v18.zip",
    )
    run(`unzip $asjp18 -d tmp/`)
    cp("tmp/lexibank-asjp-fb8987f/cldf/languages.csv", languagesF, force = true)
    rm("tmp", recursive = true)
end
asjp = CSV.read("../data/asjp/languages.csv", DataFrame)

##

longnames =
    [replace(r.classification_wals * "." * r.Name, "-" => "_") for r in eachrow(asjp)]

insertcols!(asjp, 2, :longname => longnames)

data = innerjoin(
    dropmissing(select(asjp, [:Glottocode, :longname])),
    data,
    on = :Glottocode,
)

##

asjp18ClusteredF = download(
    "https://osf.io/tdma5/download",
    "../data/asjp18Clustered.csv",
)

asjp18Clustered = CSV.read(asjp18ClusteredF, DataFrame)

asjp18Clustered[:,:glot_fam] = replace.(asjp18Clustered.glot_fam, "'" => "_")

longnames = [
    replace(join([r.wls_fam, r.wls_gen, r.doculect], "."), "-" => "_") for
    r in eachrow(asjp18Clustered)
]

insertcols!(asjp18Clustered, 2, :longname => longnames)
##
ln2g = unique(select(asjp18Clustered, [:longname, :glot_fam]))
data = innerjoin(data, ln2g, on=:longname)

##

taxa = unique(data.longname)
taxa = [x for x in taxa if x in asjp18Clustered.longname]

asjpCC = filter(x -> x.longname ∈ taxa, asjp18Clustered)

data = filter(x -> x.longname ∈ taxa, data)

fDict = Dict(
    "82A" => "VS",
    "83A" => "VO",
    "85A" => "PN",
    "86A" => "NG",
    "87A" => "NA",
    "88A" => "ND",
    "89A" => "NNum",
    "90A" => "NRc"
)

rename!(data, fDict)

select!(data, [:longname, :glot_fam, :VS, :VO, :PN, :NG, :NA, :ND, :NNum, :NRc])

CSV.write("../data/charMtx.csv", data)

##

features = names(data)[3:end]

fPairs = [
    join([f1, f2], "-") for (i, f1) in enumerate(features) for
    (j, f2) in enumerate(features) if i < j
]

open("../data/fpairs.txt", "w") do file
    for fp in fPairs
        write(file, fp)
        write(file, "\n")
    end
end
##

famFreqs = sort(combine(groupby(data, :glot_fam), nrow), :nrow, rev=true)

open("../data/families.csv", "w") do file
    write(file, join(famFreqs.glot_fam, "\n"))
    write(file, "\n")
end

glot2Plus = filter(x -> x.nrow > 1, famFreqs).glot_fam

isolates = filter(x -> x.nrow == 1, famFreqs).glot_fam

# open("../data/isolates.csv", "w") do file
#     write(file, join(isolates, "\n"))
#     write(file, "\n")
# end
#

CSV.write("../data/famFrequencies.csv", famFreqs)

##

codingDict = Dict(
    (missing, missing) => "-",
    (missing, 1) => "-",
    (missing, 2) => "-",
    (1, missing) => "-",
    (2, missing) => "-",
    (1,1) => "a",
    (1,2) => "b",
    (2,1) => "c",
    (2,2) => "d",
)

##

pairMtx = DataFrame(taxon = taxa)
for fp in fPairs
    f1, f2 = Symbol.(split(fp, "-"))
    insertcols!(
        pairMtx,
        fp => [codingDict[x] for x in zip(data[:, f1], data[:, f2])],
    )
end

CSV.write("../data/fpairMtx.csv", pairMtx)

##

try
    mkdir("../data/posteriorTrees")
catch e
end

for fm in isolates
    l = first(filter(x -> x.glot_fam == fm, data).longname)
    nex = "($l:.01, dummy:.01);"
    open("../data/posteriorTrees/" * fm * ".posterior.tree", "w") do file
        write(file, nex)
    end
end

##

worldGlotF = download("https://osf.io/jyvgt/download", "../data/world_fullGlot.tre")

glot = ete3.Tree(worldGlotF)

glot.prune(taxa)

##

concepts = unique(asjpCC.concept)

sounds = first.(sort(unique(split(join(asjpCC.simplified), ""))))

df = @pipe asjpCC |>
      select(_, [:longname, :concept, :simplified]) |>
      groupby(_, [:longname, :concept]) |>
      combine(_, :simplified => join => :words) |>
      unstack(_, :longname, :concept, :words)

scDict = Dict()
@showprogress for (i,l) in enumerate(taxa), (j,c) in enumerate(concepts), s in sounds
    if ismissing(df[i, j+1])
        scDict[l,c,s] = "-"
    else
        scDict[l,c,s] = Int(s ∈ df[i, j+1])
    end
end

##

asjpCC = filter(x -> x.longname ∈ taxa, asjp18Clustered)

try
    mkdir("../data/asjpNex/")
catch e
end

try
    mkdir("../data/asjpNex/output")
catch e
end

glot3 = filter(x->x.nrow>2, famFreqs).glot_fam

open("../data/glot3.txt", "w") do file
    write(file, join(glot3, "\n")*"\n")
end

##

ccc = unique(asjpCC.cClass)

lc = Dict()
for l in taxa, c in concepts
    lc[l,c] = false
end

for (l,c) in zip(asjpCC.longname, asjpCC.concept)
    lc[l,c] = true
end

cc2c = Dict(zip(asjpCC.cClass, asjpCC.concept))

ccDict = Dict()
@showprogress for l in taxa, cc in ccc
    if lc[l, cc2c[cc]]
        ccDict[l,cc] = "0"
    else
        ccDict[l,cc] = "-"
    end
end

for (l, cc) in zip(asjpCC.longname, asjpCC.cClass) |> unique
    ccDict[l, cc] = "1"
end

ccMtx = Array{String}(undef, size(taxa,1), size(ccc,1))
@showprogress for (i,l) in enumerate(taxa), (j,cc) in enumerate(ccc)
    ccMtx[i,j] = ccDict[l,cc]
end
##

l2glot = Dict(zip(asjpCC.longname, asjpCC.glot_fam))

famV = [l2glot[l] for l in taxa]
##

scChar_ = []
@showprogress for c in concepts
    sChar = Array{String}(undef, length(taxa), length(sounds))
    for (i,l) in enumerate(taxa), (j,s) in enumerate(sounds)
        sChar[i,j] = string(scDict[l,c,s])
    end
    push!(scChar_, sChar)
end

scMtx = hcat(scChar_...)

##

@showprogress for fm in glot3
    fmTaxa = taxa[famV .== fm]
    ccChar = ccMtx[famV .== fm,:]
    scChar = scMtx[famV .== fm,:]
    charMtx = hcat(ccChar, scChar)
    pad = maximum(length.(fmTaxa))+5
    nex = """
#NEXUS

BEGIN DATA;
DIMENSIONS ntax=$(length(fmTaxa)) NCHAR=$(size(charMtx,2));
FORMAT DATATYPE=restriction GAP=? MISSING=- interleave=yes;
MATRIX

"""
    for (i,l) in enumerate(fmTaxa)
        nex
        nex *= rpad(l, pad)
        nex *= join(charMtx[i,:])
        nex *= "\n"
    end
    nex *= """

;

END;
"""
    open("../data/asjpNex/"*fm*".nex", "w") do file
        write(file, nex)
    end
end
##

try
    mkdir("revbayes")
catch e
end

glot2 = filter(x -> x.nrow==2, famFreqs).glot_fam

@showprogress for fm in glot2
    fmTaxa = taxa[famV .== fm]
    ccChar = ccMtx[famV .== fm,:]
    scChar = scMtx[famV .== fm,:]
    charMtx = hcat(ccChar, scChar)

    pad = maximum(length.(fmTaxa))+5
    nex = """
#NEXUS

BEGIN DATA;
DIMENSIONS ntax=$(length(fmTaxa)) NCHAR=$(size(charMtx,2));
FORMAT DATATYPE=restriction GAP=? MISSING=- interleave=yes;
MATRIX

"""
    for (i,l) in enumerate(fmTaxa)
        nex
        nex *= rpad(l, pad)
        nex *= join(charMtx[i,:])
        nex *= "\n"
    end
    nex *= """

;

END;
"""

    open("../data/asjpNex/"*fm*".nex", "w") do file
        write(file, nex)
    end
rb = """
family = "$fm"
source("../phylogeny.Rev")
"""
    open("revbayes/$(fm).Rev", "w") do file
        write(file, rb)
    end
end

##

geoData = innerjoin(
    data[:, [:longname, :glot_fam]],
    asjp[:, [:longname, :Longitude, :Latitude]],
    on=:longname
)

CSV.write("../data/geoData.csv", geoData)

max_stars_repo_name: gerhardJaeger/phylogeneticTypology

max_issues_repo_name: gerhardJaeger/phylogeneticTypology

max_forks_repo_name: gerhardJaeger/phylogeneticTypology


text:
export ArrayObservable, from

import Base: ==
import Base: show

abstract type Scalarness end

struct Scalar              <: Scalarness end
struct NonScalar           <: Scalarness end
struct UndefinedScalarness <: Scalarness end

scalarness(::Type)                   = UndefinedScalarness()
scalarness(::Type{<:Number})         = Scalar()
scalarness(::Type{<:Char})           = Scalar()
scalarness(::Type{<:AbstractArray})  = NonScalar()
scalarness(::Type{<:Tuple})          = NonScalar()
scalarness(::Type{<:AbstractString}) = NonScalar()

as_array(x::T) where T = as_array(scalarness(T), x)

as_array(::Scalar, x)              = [ x ]
as_array(::NonScalar, x)           = collect(x)
as_array(::UndefinedScalarness, x) = error("Value of type $(typeof(x)) has undefined scalarness type. \nConsider implement scalarness(::Type{<:$(typeof(x))}).")

"""
    ArrayObservable{D, H}(values::Vector{D}, scheduler::H) where { D, H }

ArrayObservable wraps a regular Julia array into an observable. Uses scheduler object to schedule messages delivery.

# Constructor arguments
- `values`: array of values to be wrapped
- `scheduler`: Scheduler-like object

See also: [`Subscribable`](@ref), [`from`](@ref)
"""
@subscribable struct ArrayObservable{D, H} <: ScheduledSubscribable{D}
    values    :: Vector{D}
    scheduler :: H
end

getscheduler(observable::ArrayObservable) = observable.scheduler

function on_subscribe!(observable::ArrayObservable, actor, scheduler)
    for value in observable.values
        next!(actor, value, scheduler)
    end
    complete!(actor, scheduler)
    return voidTeardown
end

"""
    from(x; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }
    from(a::Vector{D}; scheduler::H = AsapScheduler()) where { D, H <: AbstractScheduler }

Creation operator for the `ArrayObservable` that emits either a single value if x has a `Scalar` trait specification or a collection of values if x has a `NonScalar` trait specification.
Throws an ErrorException if x has `UndefinedScalarness` trait type. To specify scalarness for arbitrary type T some can implement an additional method
for `scalarness(::Type{<:MyType})` function and to specify scalarness behavior. Optionally accepts custom scheduler-like object to schedule messages delivery.

# Arguments
- `x`: an object to be wrapped into array of values
- `scheduler`: optional, scheduler-like object

For an object `x` to be a valid input for a from operator it must implement
`Rocket.scalarness(::Type{ <: T })` method which should return either `Rocket.Scalar` or
`Rocket.NonScalar` objects. In first case `from` operator will treat `x` as a single scalar value and will wrap
it into a vector while in the second case `from` operator will convert `x` object into an array using `collect` function.

For arbitrary iterable objects consider using `iterable` creation operator.

# Note
`from` operators creates a copy of `x` using `collect` on a given object.

# Examples

```jldoctest
using Rocket

source = from([ 0, 1, 2 ])
subscribe!(source, logger())
;

# output

[LogActor] Data: 0
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed

```

```jldoctest
using Rocket

source = from(( 0, 1, 2 ))
subscribe!(source, logger())
;

# output

[LogActor] Data: 0
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed

```

```jldoctest
using Rocket

source = from(0)
subscribe!(source, logger())
;

# output

[LogActor] Data: 0
[LogActor] Completed

```

```jldoctest
using Rocket

source = from("Hello, world!")
subscribe!(source, logger())
;

# output

[LogActor] Data: H
[LogActor] Data: e
[LogActor] Data: l
[LogActor] Data: l
[LogActor] Data: o
[LogActor] Data: ,
[LogActor] Data:
[LogActor] Data: w
[LogActor] Data: o
[LogActor] Data: r
[LogActor] Data: l
[LogActor] Data: d
[LogActor] Data: !
[LogActor] Completed

```

See also: [`ArrayObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref), [`iterable`](@ref)
"""
from(x; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }               = from(as_array(x); scheduler = scheduler)
from(a::Vector{D}; scheduler::H = AsapScheduler()) where { D, H <: AbstractScheduler } = ArrayObservable{D, H}(copy(a), scheduler)

Base.:(==)(left::ArrayObservable{D, H},  right::ArrayObservable{D, H}) where { D, H } = left.values == right.values
Base.:(==)(left::ArrayObservable,        right::ArrayObservable)                      = false

Base.show(io::IO, ::ArrayObservable{D, H}) where { D, H } = print(io, "ArrayObservable($D, $H)")


text:
# This file is based on LinearAlgebra/src/bidiag.jl, a part of Julia. License is MIT: https://julialang.org/license

# Bidiagonal matrices
struct Bidiagonal{T,DV<:AbstractVector{T},EV<:AbstractVector{T}} <: AbstractBandedMatrix{T}
    dv::DV      # diagonal
    ev::EV      # sub/super diagonal
    uplo::Char # upper bidiagonal ('U') or lower ('L')
    function Bidiagonal{T,DV,EV}(dv, ev, uplo::AbstractChar) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}}
        require_one_based_indexing(dv, ev)
        if length(ev) != max(length(dv)-1, 0)
            throw(DimensionMismatch("length of diagonal vector is $(length(dv)), length of off-diagonal vector is $(length(ev))"))
        end
        new{T,DV,EV}(dv, ev, uplo)
    end
end
function Bidiagonal{T,DV,EV}(dv, ev, uplo::Symbol) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}}
    Bidiagonal{T,DV,EV}(dv, ev, char_uplo(uplo))
end
function Bidiagonal{T}(dv::AbstractVector, ev::AbstractVector, uplo::Union{Symbol,AbstractChar}) where {T}
    Bidiagonal(convert(AbstractVector{T}, dv)::AbstractVector{T},
               convert(AbstractVector{T}, ev)::AbstractVector{T},
               uplo)
end

"""
    Bidiagonal(dv::V, ev::V, uplo::Symbol) where V <: AbstractVector

Constructs an upper (`uplo=:U`) or lower (`uplo=:L`) bidiagonal matrix using the
given diagonal (`dv`) and off-diagonal (`ev`) vectors. The result is of type `Bidiagonal`
and provides efficient specialized linear solvers, but may be converted into a regular
matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short). The length of `ev`
must be one less than the length of `dv`.

# Examples
```jldoctest
julia> dv = [1, 2, 3, 4]
4-element Vector{Int64}:
 1
 2
 3
 4

julia> ev = [7, 8, 9]
3-element Vector{Int64}:
 7
 8
 9

julia> Bu = Bidiagonal(dv, ev, :U) # ev is on the first superdiagonal
4×4 Bidiagonal{Int64, Vector{Int64}}:
 1  7  ⋅  ⋅
 ⋅  2  8  ⋅
 ⋅  ⋅  3  9
 ⋅  ⋅  ⋅  4

julia> Bl = Bidiagonal(dv, ev, :L) # ev is on the first subdiagonaldataonal
4×4 Bidiagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅  ⋅
 7  2  ⋅  ⋅
 ⋅  8  3  ⋅
 ⋅  ⋅  9  4
```
"""
Bidiagonal(dv::DV, ev::EV, uplo::Symbol) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} =
    Bidiagonal{T,DV,EV}(dv, ev, char_uplo(uplo))
Bidiagonal(dv::DV, ev::EV, uplo::AbstractChar) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} =
    Bidiagonal{T,DV,EV}(dv, ev, uplo)
function Bidiagonal(dv::AbstractVector{T}, ev::AbstractVector{V}, ulo) where {T,V}
    TV = promote_type(T,V)
    Bidiagonal(convert(AbstractVector{TV}, dv), convert(AbstractVector{TV}, ev), ulo)
end

"""
    Bidiagonal(A, uplo::Symbol)

Construct a `Bidiagonal` matrix from the main diagonal of `A` and
its first super- (if `uplo=:U`) or sub-diagonal (if `uplo=:L`).

# Examples
```jldoctest
julia> A = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]
4×4 Matrix{Int64}:
 1  1  1  1
 2  2  2  2
 3  3  3  3
 4  4  4  4

julia> Bidiagonal(A, :U) # contains the main diagonal and first superdiagonal of A
4×4 Bidiagonal{Int64, Vector{Int64}}:
 1  1  ⋅  ⋅
 ⋅  2  2  ⋅
 ⋅  ⋅  3  3
 ⋅  ⋅  ⋅  4

julia> Bidiagonal(A, :L) # contains the main diagonal and first subdiagonaldataonal of A
4×4 Bidiagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅  ⋅
 2  2  ⋅  ⋅
 ⋅  3  3  ⋅
 ⋅  ⋅  4  4
```
"""
function Bidiagonal(A::AbstractMatrix, uplo::Symbol)
    Bidiagonal(diag(A, 0), diag(A, uplo === :U ? 1 : -1), uplo)
end

Bidiagonal(A::Bidiagonal) = A
Bidiagonal{T}(A::Bidiagonal{T}) where {T} = A
Bidiagonal{T}(A::Bidiagonal) where {T} = Bidiagonal{T}(A.dv, A.ev, A.uplo)

function getindex(A::Bidiagonal{T}, i::Integer, j::Integer) where T
    if !((1 <= i <= size(A,2)) && (1 <= j <= size(A,2)))
        throw(BoundsError(A,(i,j)))
    end
    if i == j
        return A.dv[i]
    elseif A.uplo == 'U' && (i == j - 1)
        return A.ev[i]
    elseif A.uplo == 'L' && (i == j + 1)
        return A.ev[j]
    else
        return zero(T)
    end
end

function setindex!(A::Bidiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.dv[i] = x
    elseif A.uplo == 'U' && (i == j - 1)
        @inbounds A.ev[i] = x
    elseif A.uplo == 'L' && (i == j + 1)
        @inbounds A.ev[j] = x
    elseif !iszero(x)
        throw(ArgumentError(string("cannot set entry ($i, $j) off the ",
            "$(istriu(A) ? "upper" : "lower") bidiagonal band to a nonzero value ($x)")))
    end
    return x
end

#Converting from Bidiagonal to dense Matrix
function Matrix{T}(A::Bidiagonal) where T
    n = size(A, 1)
    B = zeros(T, n, n)
    if n == 0
        return B
    end
    for i = 1:n - 1
        B[i,i] = A.dv[i]
        if A.uplo == 'U'
            B[i, i + 1] = A.ev[i]
        else
            B[i + 1, i] = A.ev[i]
        end
    end
    B[n,n] = A.dv[n]
    return B
end
Matrix(A::Bidiagonal{T}) where {T} = Matrix{T}(A)
Array(A::Bidiagonal) = Matrix(A)
promote_rule(::Type{Matrix{T}}, ::Type{<:Bidiagonal{S}}) where {T,S} =
    @isdefined(T) && @isdefined(S) ? Matrix{promote_type(T,S)} : Matrix
promote_rule(::Type{Matrix}, ::Type{<:Bidiagonal}) = Matrix

# Converting to LinearAlgebra
convert(::Type{LinearAlgebra.Bidiagonal}, B::Bidiagonal) = LinearAlgebra.Bidiagonal(B.dv, B.ev, B.uplo)

#Converting from Bidiagonal to Tridiagonal
function Tridiagonal{T}(A::Bidiagonal) where T
    dv = convert(AbstractVector{T}, A.dv)
    ev = convert(AbstractVector{T}, A.ev)
    z = fill!(similar(ev), zero(T))
    A.uplo == 'U' ? Tridiagonal(z, dv, ev) : Tridiagonal(ev, dv, z)
end
promote_rule(::Type{<:Tridiagonal{T}}, ::Type{<:Bidiagonal{S}}) where {T,S} =
    @isdefined(T) && @isdefined(S) ? Tridiagonal{promote_type(T,S)} : Tridiagonal
promote_rule(::Type{<:Tridiagonal}, ::Type{<:Bidiagonal}) = Tridiagonal

# When asked to convert Bidiagonal to AbstractMatrix{T}, preserve structure by converting to Bidiagonal{T} <: AbstractMatrix{T}
AbstractMatrix{T}(A::Bidiagonal) where {T} = convert(Bidiagonal{T}, A)

convert(T::Type{<:Bidiagonal}, m::AbstractMatrix) = m isa T ? m : T(m)

# For B<:Bidiagonal, similar(B[, neweltype]) should yield a Bidiagonal matrix.
# On the other hand, similar(B, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(B::Bidiagonal, ::Type{T}) where {T} = Bidiagonal(similar(B.dv, T), similar(B.ev, T), B.uplo)
# The method below is moved to SparseArrays for now
# similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

####################
# Generic routines #
####################

function show(io::IO, M::Bidiagonal)
    # TODO: make this readable and one-line
    summary(io, M); println(io, ":")
    print(io, " diag:")
    print_matrix(io, (M.dv)')
    print(io, M.uplo == 'U' ? "\n super:" : "\n sub:")
    print_matrix(io, (M.ev)')
end

size(M::Bidiagonal) = (length(M.dv), length(M.dv))
function size(M::Bidiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d <= 2
        return length(M.dv)
    else
        return 1
    end
end

#Elementary operations
for func in (:conj, :copy, :real, :imag)
    @eval ($func)(M::Bidiagonal) = Bidiagonal(($func)(M.dv), ($func)(M.ev), M.uplo)
end

adjoint(B::Bidiagonal) = Adjoint(B)
transpose(B::Bidiagonal) = Transpose(B)
adjoint(B::Bidiagonal{<:Real}) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? :L : :U)
transpose(B::Bidiagonal{<:Number}) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? :L : :U)
function Base.copy(aB::Adjoint{<:Any,<:Bidiagonal})
    B = aB.parent
    return Bidiagonal(map(x -> copy.(adjoint.(x)), (B.dv, B.ev))..., B.uplo == 'U' ? :L : :U)
end
function Base.copy(tB::Transpose{<:Any,<:Bidiagonal})
    B = tB.parent
    return Bidiagonal(map(x -> copy.(transpose.(x)), (B.dv, B.ev))..., B.uplo == 'U' ? :L : :U)
end

iszero(M::Bidiagonal) = iszero(M.dv) && iszero(M.ev)
isone(M::Bidiagonal) = all(isone, M.dv) && iszero(M.ev)
function istriu(M::Bidiagonal, k::Integer=0)
    if M.uplo == 'U'
        if k <= 0
            return true
        elseif k == 1
            return iszero(M.dv)
        else # k >= 2
            return iszero(M.dv) && iszero(M.ev)
        end
    else # M.uplo == 'L'
        if k <= -1
            return true
        elseif k == 0
            return iszero(M.ev)
        else # k >= 1
            return iszero(M.ev) && iszero(M.dv)
        end
    end
end
function istril(M::Bidiagonal, k::Integer=0)
    if M.uplo == 'U'
        if k >= 1
            return true
        elseif k == 0
            return iszero(M.ev)
        else # k <= -1
            return iszero(M.ev) && iszero(M.dv)
        end
    else # M.uplo == 'L'
        if k >= 0
            return true
        elseif k == -1
            return iszero(M.dv)
        else # k <= -2
            return iszero(M.dv) && iszero(M.ev)
        end
    end
end
isdiag(M::Bidiagonal) = iszero(M.ev)

function tril!(M::Bidiagonal, k::Integer=0)
    n = length(M.dv)
    if !(-n - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
    elseif M.uplo == 'U' && k < 0
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif k < -1
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif M.uplo == 'U' && k == 0
        fill!(M.ev,0)
    elseif M.uplo == 'L' && k == -1
        fill!(M.dv,0)
    end
    return M
end

function triu!(M::Bidiagonal, k::Integer=0)
    n = length(M.dv)
    if !(-n + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least",
            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
    elseif M.uplo == 'L' && k > 0
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif k > 1
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif M.uplo == 'L' && k == 0
        fill!(M.ev,0)
    elseif M.uplo == 'U' && k == 1
        fill!(M.dv,0)
    end
    return M
end

function diag(M::Bidiagonal, n::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of n
    if n == 0
        return copyto!(similar(M.dv, length(M.dv)), M.dv)
    elseif (n == 1 && M.uplo == 'U') ||  (n == -1 && M.uplo == 'L')
        return copyto!(similar(M.ev, length(M.ev)), M.ev)
    elseif -size(M,1) <= n <= size(M,1)
        return fill!(similar(M.dv, size(M,1)-abs(n)), 0)
    else
        throw(ArgumentError(string("requested diagonal, $n, must be at least $(-size(M, 1)) ",
            "and at most $(size(M, 2)) for an $(size(M, 1))-by-$(size(M, 2)) matrix")))
    end
end

function +(A::Bidiagonal, B::Bidiagonal)
    if A.uplo == B.uplo || length(A.dv) == 0
        Bidiagonal(A.dv+B.dv, A.ev+B.ev, A.uplo)
    else
        newdv = A.dv+B.dv
        Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(B.ev), newdv, typeof(newdv)(A.ev)) : (typeof(newdv)(A.ev), newdv, typeof(newdv)(B.ev)))...)
    end
end

function -(A::Bidiagonal, B::Bidiagonal)
    if A.uplo == B.uplo || length(A.dv) == 0
        Bidiagonal(A.dv-B.dv, A.ev-B.ev, A.uplo)
    else
        newdv = A.dv-B.dv
        Tridiagonal((A.uplo == 'U' ? (typeof(newdv)(-B.ev), newdv, typeof(newdv)(A.ev)) : (typeof(newdv)(A.ev), newdv, typeof(newdv)(-B.ev)))...)
    end
end

-(A::Bidiagonal)=Bidiagonal(-A.dv,-A.ev,A.uplo)
*(A::Bidiagonal, B::Number) = Bidiagonal(A.dv*B, A.ev*B, A.uplo)
*(B::Number, A::Bidiagonal) = A*B
/(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.uplo)
\(B::Number, A::Bidiagonal) = A/B

function ==(A::Bidiagonal, B::Bidiagonal)
    if A.uplo == B.uplo
        return A.dv == B.dv && A.ev == B.ev
    else
        return iszero(A.ev) && iszero(B.ev) && A.dv == B.dv
    end
end

function check_A_mul_B!_sizes(C, A, B)
    require_one_based_indexing(C)
    require_one_based_indexing(A)
    require_one_based_indexing(B)
    nA, mA = size(A)
    nB, mB = size(B)
    nC, mC = size(C)
    if nA != nC
        throw(DimensionMismatch("sizes size(A)=$(size(A)) and size(C) = $(size(C)) must match at first entry."))
    elseif mA != nB
        throw(DimensionMismatch("second entry of size(A)=$(size(A)) and first entry of size(B) = $(size(B)) must match."))
    elseif mB != mC
        throw(DimensionMismatch("sizes size(B)=$(size(B)) and size(C) = $(size(C)) must match at first second entry."))
    end
end

# function to get the internally stored vectors for Bidiagonal and [Sym]Tridiagonal
# to avoid allocations in A_mul_B_td! below (#24324, #24578)
_diag(A::Tridiagonal, k) = k == -1 ? A.dl : k == 0 ? A.d : A.du
_diag(A::SymTridiagonal, k) = k == 0 ? A.dv : A.ev
function _diag(A::Bidiagonal, k)
    if k == 0
        return A.dv
    elseif (A.uplo == 'L' && k == -1) || (A.uplo == 'U' && k == 1)
        return A.ev
    else
        return diag(A, k)
    end
end

factorize(A::Bidiagonal) = A

Base._sum(A::Bidiagonal, ::Colon) = sum(A.dv) + sum(A.ev)
function Base._sum(A::Bidiagonal, dims::Integer)
    res = Base.reducedim_initarray(A, dims, zero(eltype(A)))
    n = length(A.dv)
    if n == 0
        # Just to be sure. This shouldn't happen since there is a check whether
        # length(A.dv) == length(A.ev) + 1 in the constructor.
        return res
    elseif n == 1
        res[1] = A.dv[1]
        return res
    end
    @inbounds begin
        if (dims == 1 && A.uplo == 'U') || (dims == 2 && A.uplo == 'L')
            res[1] = A.dv[1]
            for i = 2:length(A.dv)
                res[i] = A.ev[i-1] + A.dv[i]
            end
        elseif (dims == 1 && A.uplo == 'L') || (dims == 2 && A.uplo == 'U')
            for i = 1:length(A.dv)-1
                res[i] = A.ev[i] + A.dv[i]
            end
            res[end] = A.dv[end]
        elseif dims >= 3
            if A.uplo == 'U'
                for i = 1:length(A.dv)-1
                    res[i,i]   = A.dv[i]
                    res[i,i+1] = A.ev[i]
                end
            else
                for i = 1:length(A.dv)-1
                    res[i,i]   = A.dv[i]
                    res[i+1,i] = A.ev[i]
                end
            end
            res[end,end] = A.dv[end]
        end
    end
    res
end

MemoryLayout(::Type{<:Bidiagonal{<:Any,DV,EV}}) where {DV,EV} = bidiagonallayout(MemoryLayout(DV), MemoryLayout(EV))
bandwidths(A::Bidiagonal) = A.uplo == 'U' ? (0,1) : (1,0)
bidiagonaluplo(A::Bidiagonal) = A.uplo

diagonaldata(D::Bidiagonal) = D.dv
supdiagonaldata(D::Bidiagonal) = D.uplo == 'U' ? D.ev : throw(ArgumentError("$D is lower-bidiagonal"))
subdiagonaldata(D::Bidiagonal) = D.uplo == 'L' ? D.ev : throw(ArgumentError("$D is upper-bidiagonal"))

permutedims(B::Bidiagonal) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? 'L' : 'U')

function BlockArrays.sizes_from_blocks(A::Bidiagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    (size.(A.dv, 1), size.(A.dv,2))
end


text:
# This file is a part of ChangesOfVariables.jl, licensed under the MIT License (MIT).

import Test
import ChangesOfVariables
import Documenter

Test.@testset "Package ChangesOfVariables" begin
    include("test_test.jl")
    include("test_with_ladj.jl")

    # doctests
    Documenter.DocMeta.setdocmeta!(
        ChangesOfVariables,
        :DocTestSetup,
        :(using ChangesOfVariables);
        recursive=true,
    )
    Documenter.doctest(ChangesOfVariables)
end # testset


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Lz4_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Lz4")
JLLWrappers.@generate_main_file("Lz4", UUID("5ced341a-0733-55b8-9ab6-a4889d929147"))
end  # module Lz4_jll


text:
using LocalitySensitive
using Test
using Random

@testset "MinHash" begin
    @testset "estimate jaccard distance" begin
        mh = MinHash()
        test_fp = fingerprint(mh, shingle("test", size=3))
        @test estimate_jaccard(test_fp, test_fp) == 1
        @test estimate_jaccard(test_fp, fingerprint(mh, shingle("1234", size=3))) == 0
        @test estimate_jaccard(test_fp, fingerprint(mh, shingle("testr", size=3))) > 0.5

        @test estimate_jaccard(fingerprint(mh, [1,2,3,4]), fingerprint(mh, [1,2,3,4])) == 1
        @test estimate_jaccard(fingerprint(mh, [1,2,3,4]), fingerprint(mh, [1,2,3])) > 0.5
        @test estimate_jaccard(fingerprint(mh, [1,2,3,4]), fingerprint(mh, [5,6])) == 0

        @test estimate_jaccard(fingerprint(mh, 1:4), fingerprint(mh, 1:4)) == 1
        @test estimate_jaccard(fingerprint(mh, 1:4), fingerprint(mh, 5:6)) == 0

        sets = [collect(i:j) for (i,j) in zip(1:100, 3:102)]
        rng = MersenneTwister(1)
        for i in 1:10
            s1 = sets[randperm(rng, length(sets))][1:i*10]
            @test (i-1)/10 < estimate_jaccard(fingerprint(mh, sets), fingerprint(mh, s1)) < (i+1)/10
        end
    end

    @testset "Special function for strings" begin
        mh = MinHash()
        test_fp = fingerprint(mh, "test")
        @test estimate_jaccard(test_fp, test_fp) == 1
        @test estimate_jaccard(test_fp, fingerprint(mh, "1234")) == 0
        @test estimate_jaccard(test_fp, fingerprint(mh, "testr")) > 0.5
    end

    @testset "push document to index" begin
        mh = MinHash()
        mhind = MinHashIndex(fingerprinter=mh, threshold=0.1)
        signature = fingerprint(mh, "test")
        push!(mhind, signature)
        signature = fingerprint(mh, "testsim")
        @test [1] == find_similar(mhind, signature)
        push!(mhind, fingerprint(mh,"58849"))
        signature = fingerprint(mh, "58849sim")
        @test [2] == find_similar(mhind, signature)

        sets = [collect(i:j) for (i,j) in zip(1:100, 3:102)]
        mh = MinHash(400, MersenneTwister(1))
        for i in 3:6
            mhind = MinHashIndex(fingerprinter=mh, threshold=(i-2)/10)
            push!(mhind, fingerprint(mh, sets))
            s1 = sets[1:i*10]
            @test [1] == find_similar(mhind, fingerprint(mh, s1))
            mhind = MinHashIndex(fingerprinter=mh, threshold=(i+3)/10)
            push!(mhind, fingerprint(mh, sets))
            @test [] == find_similar(mhind, fingerprint(mh, s1))
        end
    end

    @testset "Find all pair indices" begin
        mh = MinHash()
        mhind = MinHashIndex(fingerprinter= mh, threshold=0)
        push!(mhind, fingerprint(mh, "test"))
        push!(mhind, fingerprint(mh, "testSIM"))
        push!(mhind, fingerprint(mh, "SIMtest"))
        push!(mhind, fingerprint(mh, "58849"))
        push!(mhind, fingerprint(mh, "58849sim"))
        push!(mhind, fingerprint(mh, "sim58849"))
        @test sort([(1,2),(1,3),(2,3),(4,5),(4,6),(5,6)]) == sort(similar_pairs(mhind))
    end

    @testset "Fingerprint all" begin
        mh = MinHash()
        documents = ["als", "kjsqer", "al;kjslkj", "sld"]
        fingerprints = fingerprint_all(mh, [shingle(d, size=3) for d in documents])
        mhind = MinHashIndex(fingerprinter=mh, threshold = 0.9)
        for f in fingerprints
            push!(mhind, f)
        end
        @test [1] == find_similar(mhind, fingerprint(mh, "als"))
        @test [] == find_similar(mhind, fingerprint(mh, "al"))
    end
end


text:
abstract type HMCProposal end

@with_kw struct FixedStepNumber <: HMCProposal
    nsteps::Int64 = 10
end

@with_kw struct FixedTrajectoryLength <: HMCProposal
    trajectory_length::Float64 = 2.0
end

@with_kw struct NUTS <: HMCProposal
    sampling::Symbol = :MultinomialTS
    nuts::Symbol = :ClassicNoUTurn
end

function AHMCProposal(
    proposal::FixedStepNumber,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.StaticTrajectory(integrator, proposal.nsteps)
end

function AHMCProposal(
    proposal::FixedTrajectoryLength,
    integrator::AdvancedHMC.AbstractIntegrator
)
    return AdvancedHMC.HMCDA(integrator, proposal.trajectory_length)
end

function AHMCProposal(
    proposal::NUTS,
    integrator::AdvancedHMC.AbstractIntegrator
)
    sampling_type = getfield(AdvancedHMC, proposal.sampling)
    nuts_type = getfield(AdvancedHMC, proposal.nuts)

    return AdvancedHMC.NUTS{sampling_type, nuts_type}(integrator)
end


text:
module MolecularDynamicsIntegrators

import Reexport: @reexport

include("./Utils.jl")

include("./Types.jl")
@reexport using .Types

include("./Constraints.jl")
@reexport using .Constraints

include("./Verlet/Verlet.jl")
@reexport using .Verlet

include("./RATTLE/RATTLE.jl")
@reexport using .RATTLE

end


text:
module DataTables

import TableShowUtils, TableTraitsUtils, ReadOnlyArrays

using DataValues

export DataTable, NA, isna

struct DataTable{T, TCOLS} <: AbstractVector{T}
    columns::TCOLS
end

function fromNT(nt)
    nt = map(i->i isa ReadOnlyArrays.ReadOnlyArray ? i : ReadOnlyArrays.ReadOnlyArray(i), nt)
    tx = typeof(nt)
    et = NamedTuple{propertynames(nt), Tuple{(eltype(fieldtype(tx, i)) for i in 1:fieldcount(typeof(nt)))...}}
    return DataTable{et, typeof(nt)}(nt)
end

swap_dva_in(A) = A
swap_dva_in(A::Array{<:DataValue}) = DataValueArray(A)

function DataTable(;cols...)
    return fromNT(map(col -> swap_dva_in(col), values(cols)))
end

function DataTable(table)
    cols, colnames = TableTraitsUtils.create_columns_from_iterabletable(table)

    return fromNT(NamedTuple{tuple(colnames...)}(tuple(cols...)))
end

columns(dt::DataTable) = getfield(dt, :columns)

@inline Base.getproperty(dt::DataTable, name::Symbol) = getproperty(columns(dt), name)

Base.size(dt::DataTable) = size(columns(dt)[1])

Base.IndexStyle(::Type{T}) where {T<:DataTable}= Base.IndexLinear()

function Base.checkbounds(::Type{Bool}, dt::DataTable, i)
    cols = columns(dt)
    if length(cols)==0
        return true
    else
        return checkbounds(Bool, cols[1], i)
    end
end

@inline function Base.getindex(dt::DataTable{T}, i::Int) where {T}
    @boundscheck checkbounds(dt, i)
    return map(col -> @inbounds(getindex(col, i)), columns(dt))
end

function Base.show(io::IO, dt::DataTable)
    TableShowUtils.printtable(io, dt, "DataTable")
end

function Base.show(io::IO, ::MIME"text/plain", dt::DataTable)
    TableShowUtils.printtable(io, dt, "DataTable")
end

function Base.show(io::IO, ::MIME"text/html", dt::DataTable)
    TableShowUtils.printHTMLtable(io, dt)
end

Base.showable(::MIME"text/html", dt::DataTable) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", dt::DataTable)
    TableShowUtils.printdataresource(io, dt)
end

Base.showable(::MIME"application/vnd.dataresource+json", dt::DataTable) = true

end


text:
module MaterialReconstruction
using CorrelationFunctions.Directional
using AnnealingAPI
using Base.Iterators
using Statistics: mean, std
using LsqFit: curve_fit, coef
using PoissonRandom: pois_rand

include("utilities.jl")
include("modifiers.jl")
include("initialization.jl")
include("furnace.jl")
include("cost.jl")
include("cooldown.jl")
include("annealing.jl")

export
    # Cost functions
    euclid_mean,
    euclid_directional,
    euclid_mean_weighted,
    euclid_directional_weighted,
    čapek_cost, generalized_čapek_cost,
    # Initialization
    initialize_random,
    initialize_spheres,
    # Modifier types
    AbstractModifier,
    Swapper, Flipper, BatchModifier,
    modify!, reject!,
    # Sampler types
    AbstractSampler,
    UniformSampler, InterfaceSampler, DPNSampler,
    sample, update_pre!, update_post!,
    # Cooldown schedules
    exponential_cooldown,
    aarts_korst_cooldown,
    frost_heineman_cooldown,
    # Annealing
    Furnace, annealing_step

end # module


text:
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1,T2,W} = promote_type(T1,T2)

function _ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1, T2, W}
    T_canon = promote_type(T1,T2)
    ifelse(lt(pick_vector_width(T_canon), StaticInt{W}()), T1, T_canon)
end
function __ff_maybe_promote_int(::True, ::Type{T}, ::Type{T_canon}, ::Val{W}) where {T,T_canon,W}
    ifelse(eq(static_sizeof(T_canon), static_sizeof(T)), T, pick_integer(Val{W}(), T_canon))
end
__ff_maybe_promote_int(::False, ::Type{T1}, ::Type{T_canon}, ::Val{W}) where {T1,T_canon,W} = T_canon
function _ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1<:Integer, T2<:Integer, W}
    T_canon = promote_type(T1,T2)
    __ff_maybe_promote_int(lt(pick_vector_width(T_canon), StaticInt{W}()), T1, T_canon, Val{W}())
end
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: Integer, T2 <: Integer,W} = _ff_promote_rule(T1,T2,Val{W}())
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: FloatingTypes, T2<:FloatingTypes,W} = _ff_promote_rule(T1,T2,Val{W}())

Base.promote_rule(::Type{V}, ::Type{T2}) where {W,T1,T2<:NativeTypes,V<:AbstractSIMDVector{W,T1}} = Vec{W,ff_promote_rule(T1,T2,Val{W}())}
Base.promote_rule(::Type{V}, ::Type{Bool}) where {V<:AbstractMask} = V

_assemble_vec_unroll(::Val{N}, ::Type{V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = VecUnroll{N,W,T,V}
_assemble_vec_unroll(::Val{N}, ::Type{T}) where {N,T<:NativeTypes} = VecUnroll{N,1,T,T}
Base.promote_rule(::Type{VecUnroll{N,W,T1,V}}, ::Type{T2}) where {N,W,T1,V,T2<:NativeTypes} = _assemble_vec_unroll(Val{N}(), promote_type(V,T2))
Base.promote_rule(::Type{VecUnroll{N,W,T,V1}}, ::Type{V2}) where {N,W,T,V1,T2,V2<:AbstractSIMDVector{W,T2}} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))
# Base.promote_rule(::Type{VecUnroll{N,W,T,V1}}, ::Type{V2}) where {N,W,T,V1,V2<:AbstractSIMDVector{W}} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))
Base.promote_rule(::Type{VecUnroll{N,W,T,V1}}, ::Type{V2}) where {N,W,T,V1,V2<:AbstractMask{W}} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))
Base.promote_rule(::Type{VecUnroll{N,W,T1,V1}}, ::Type{VecUnroll{N,W,T2,V2}}) where {N,W,T1,T2,V1,V2} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))

Base.promote_rule(::Type{VecUnroll{N,W,T1,V1}}, ::Type{VecUnroll{N,1,T2,T2}}) where {N,W,T1,T2,V1} = _assemble_vec_unroll(Val{N}(), promote_type(V1,T2))
Base.promote_rule(::Type{VecUnroll{N,1,T1,T1}}, ::Type{VecUnroll{N,W,T2,V2}}) where {N,W,T1,T2,V2} = _assemble_vec_unroll(Val{N}(), promote_type(T1,V2))
Base.promote_rule(::Type{VecUnroll{N,1,T1,T1}}, ::Type{VecUnroll{N,1,T2,T2}}) where {N,T1,T2} = _assemble_vec_unroll(Val{N}(), promote_type(T1,T2))

Base.promote_rule(::Type{Mask{W,U}}, ::Type{EVLMask{W,U}}) where {W,U} = Mask{W,U}
Base.promote_rule(::Type{EVLMask{W,U}}, ::Type{Mask{W,U}}) where {W,U} = Mask{W,U}
Base.promote_rule(::Type{Bit}, ::Type{T}) where {T <: Number} = T

issigned(x) = issigned(typeof(x))
issigned(::Type{<:Signed}) = True()
issigned(::Type{<:Unsigned}) = False()
issigned(::Type{<:AbstractSIMD{<:Any,T}}) where {T} = issigned(T)
issigned(::Type{T}) where {T} = nothing
"""
 Promote, favoring <:Signed or <:Unsigned of first arg.
"""
@inline promote_div(x::Union{Integer,AbstractSIMD{<:Any,<:Integer}}, y::Union{Integer,AbstractSIMD{<:Any,<:Integer}}) = promote_div(x, y, issigned(x))
@inline promote_div(x, y) = promote(x, y)
@inline promote_div(x, y, ::Nothing) = promote(x, y) # for Integers that are neither Signed or Unsigned, e.g. Bool
@inline function promote_div(x::T1, y::T2, ::True) where {T1,T2}
    T = promote_type(T1, T2)
    signed(x % T), signed(y % T)
end
@inline function promote_div(x::T1, y::T2, ::False) where {T1,T2}
    T = promote_type(T1, T2)
    unsigned(x % T), unsigned(y % T)
end
itosize(i::Union{I,AbstractSIMD{<:Any,I}}, ::Type{J}) where {I,J} = signorunsign(i % J, issigned(I))
signorunsign(i, ::True) = signed(i)
signorunsign(i, ::False) = unsigned(i)

# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{T2}) where {M,N,W,T1,T2<:NativeTypes} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{Vec{W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecUnroll{M,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecTile{M,N,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}

@generated function _ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}, ::StaticInt{RS}) where {T1 <: IntegerTypes, T2 <: FloatingTypes,W, RS}
    T_canon = promote_type(T1,T2)
    (sizeof(T_canon) * W ≤ RS) && return T_canon
    @assert sizeof(T1) * W ≤ RS
    @assert sizeof(T1) == 4
    Float32
end
@inline function ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: IntegerTypes, T2 <: FloatingTypes,W}
    _ff_promote_rule(T1, T2, Val{W}(), register_size())
end
@generated function _promote_rule(
    ::Type{V1}, ::Type{V2}, ::StaticInt{RS}
) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2},RS}
    T = promote_type(T1,T2) # `T1` and `T2` should be defined in `Base`
    if RS ≥ W * sizeof(T)
        return :(Vec{$W,$T})
    end
    if T === Float64 || T === Float32
        N, r1 = (sizeof(T) * W) ÷ RS
        Wnew, r2 = divrem(W, N)
        @assert iszero(r)

        return :(VecUnroll{$(N-1),$Wnew,$T,Vec{$Wnew,$T}})
        # Should we demote `Float64` -> `Float32`?
        # return :(Vec{$W,$T})
        # don't demote to smaller than `Float32`
        # return :(Vec{$W,Float32})
    end
    # They're both of integers
    V1MM = V1 <: MM
    V2MM = V2 <: MM
    if V1MM ⊻ V2MM
        V1MM ? :(Vec{$W,$T2}) : :(Vec{$W,$T1})
    else # either both are `MM` or neither are
        B = W ÷ sizeof(T)
        if !V1MM # if neither are
            B = max(4, B)
        end
        I = integer_of_bytes_symbol(B, T <: Unsigned)
        :(Vec{$W,$I})
    end
end
@inline function Base.promote_rule(::Type{V1}, ::Type{V2}) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2}}
    _promote_rule(V1, V2, register_size(promote_type(T1,T2)))
end

maybethrow(::True) = throw(ArgumentError("The arguments were invalid."))
maybethrow(::False) = nothing

# not @generated, because calling `promote_type` on vector types
@inline function Base.promote_rule(
    ::Type{VecUnroll{Nm1,Wsplit,T,V1}}, ::Type{V2}
) where {Nm1,Wsplit,T,V1,T2,W,V2<:AbstractSIMDVector{W,T2}}
    maybethrow(ArrayInterface.ne(StaticInt{Nm1}() * StaticInt{Wsplit}() + StaticInt{Wsplit}(), StaticInt{W}()))
    V3 = Vec{Wsplit,T2}
    _assemble_vec_unroll(Val{Nm1}(), promote_type(V1,V3))
end
@inline function Base.promote_rule(
    ::Type{VecUnroll{Nm1,Wsplit,T,V1}}, ::Type{V2}
) where {Nm1,Wsplit,T,V1,W,V2<:AbstractMask{W}}
    maybethrow(ArrayInterface.ne(StaticInt{Nm1}() * StaticInt{Wsplit}() + StaticInt{Wsplit}(), StaticInt{W}()))
    V3 = Mask{Wsplit,mask_type(StaticInt{Wsplit}())}
    _assemble_vec_unroll(Val{Nm1}(), promote_type(V1,V3))
end

@inline function Base.promote_rule(
    ::Type{VecUnroll{Nm1,1,T,T}}, ::Type{V2}
) where {Nm1,T,T2,W,V2<:AbstractSIMDVector{W,T2}}
    _assemble_vec_unroll(Val{Nm1}(), promote_type(T,V2))
end
@inline function Base.promote_rule(
    ::Type{VecUnroll{Nm1,1,T,T}}, ::Type{V2}
) where {Nm1,T,W,V2<:AbstractMask{W}}
    _assemble_vec_unroll(Val{Nm1}(), promote_type(T,V2))
end


text:
"""
GLeakyIntegrator: Leaky Integration for the Gaussian task (section 2.3.1)
"""
struct GLeakyIntegrator <: GLearner
    etaleak::Float64
    N::Array{Float64, 1}
    Ny::Array{Float64, 1}
    mu_n::Array{Float64, 1}
end
function GLeakyIntegrator(;etaleak = .9)
    N = zeros(1)
    Ny = zeros(1)
    mu_n = zeros(1)
    GLeakyIntegrator(etaleak, N, Ny, mu_n)
end
export GLeakyIntegrator
""" X[t] = etaleak * X[t-1] + etaleak * I[.] """
function update!(learnerG::GLeakyIntegrator, y)
    learnerG.N[1] *= learnerG.etaleak # Discount
    learnerG.N[1] += learnerG.etaleak # Increase
    learnerG.Ny[1] *= learnerG.etaleak # Discount all
    learnerG.Ny[1] += learnerG.etaleak * y # Increase observed

    computemu_n!(learnerG)
end
export update!

function computemu_n!(learnerG::GLeakyIntegrator)
    learnerG.mu_n[1] = learnerG.Ny[1] ./ learnerG.N[1]
end


text:
## boolean conversions ##

convert(::Type{Bool}, x::Bool) = x
convert(::Type{Bool}, x::Real) = (x!=0)

# promote Bool to any other numeric type
promote_rule{T<:Number}(::Type{Bool}, ::Type{T}) = T

bool(x::Bool) = x
bool(x::Number) = convert(Bool, x)

sizeof(::Type{Bool}) = 1

typemin(::Type{Bool}) = false
typemax(::Type{Bool}) = true

## boolean operations ##

!(x::Bool) = box(Bool,not_int(unbox(Bool,x)))

(~)(x::Bool) = !x
(&)(x::Bool, y::Bool) = box(Bool,and_int(unbox(Bool,x),unbox(Bool,y)))
(|)(x::Bool, y::Bool) = box(Bool,or_int(unbox(Bool,x),unbox(Bool,y)))
($)(x::Bool, y::Bool) = (x!=y)

signbit(x::Bool) = false
sign(x::Bool) = x
abs(x::Bool) = x
abs2(x::Bool) = x

<(x::Bool, y::Bool) = y&!x
<=(x::Bool, y::Bool) = y|!x

## do arithmetic as Int ##

+(x::Bool) =  int(x)
-(x::Bool) = -int(x)

+(x::Bool, y::Bool) = int(x) + int(y)
-(x::Bool, y::Bool) = int(x) - int(y)
*(x::Bool, y::Bool) = x & y
^(x::Bool, y::Bool) = x | !y
^(x::Integer, y::Bool) = ifelse(y, x, one(x))

function +{T<:FloatingPoint}(x::Bool, y::T)
    ifelse(x, one(promote_type(Bool,T)) + convert(promote_type(Bool,T),y),
           convert(promote_type(Bool,T),y))
end
+(y::FloatingPoint, x::Bool) = x + y

function *{T<:Number}(x::Bool, y::T)
    ifelse(x, convert(promote_type(Bool,T),y),
           ifelse(signbit(y), -zero(promote_type(Bool,T)), zero(promote_type(Bool,T))))
end
function *{T<:Unsigned}(x::Bool, y::T)
    ifelse(x, convert(promote_type(Bool,T),y), zero(promote_type(Bool,T)))
end
*(y::Number, x::Bool) = x * y

div(x::Bool, y::Bool) = y ? x : throw(DivideError())
fld(x::Bool, y::Bool) = div(x,y)
rem(x::Bool, y::Bool) = y ? false : throw(DivideError())
mod(x::Bool, y::Bool) = rem(x,y)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia