text:
export weight_tuple

add_vector_field(net::CachedNet, grad) = add_vector_field(net, grad)
add_vector_field(net, grad) = last(_vectorize_gradient(net, grad))

#TODO I can plug here to customize this for every network type,
# ie by passing a dict of fieldnames...
_vectorize_gradient(net, x::NamedTuple) =
    _vectorize_gradient(x, zeros(out_type(net), compute_len(x)))
function _vectorize_gradient(x::NamedTuple, vec::Vector, start=1)
    i = 0
    d=Dict{Symbol, Any}()
    for f=fieldnames(typeof(x))
        di, val = _vectorize_gradient(getfield(x,f), vec, start+i)
        i += di; push!(d, f=>val)
    end
    start == 1 && push!(d, :tuple_all_weights=>[vec])
    #@views start != 1 && push!(d, :tuple_all_weights=>vec[start:start+i-1])
    i, (;d...)
end

function _vectorize_gradient(x::Tuple, vec::Vector, start)
    i = 0
    d=Vector()
    for f=fieldnames(typeof(x))
        di, val = vectorize_tuple(getfield(x,f), vec, start+i)
        i += di; push!(d, val)
    end
    i, Tuple(d)
end

function _vectorize_gradient(x::AbstractArray{<:Number}, vec::Vector, start)
    @views data_vec = vec[start:start+length(x)-1]
    reshpd_params = reshape(data_vec, size(x))
    reshpd_params .= x
    return length(x), reshpd_params
end

#=
function compute_len(vals::Union{NamedTuple,Tuple})
    tot = 0
    for f=fieldnames(typeof(vals))
        tot += compute_len(getfield(vals, f))
    end
    tot
end =#
"""
    compute_len(vals::Union{NamedTuple, Tuple})

Computes the total length of all the arrays and matrices present in this
structure and it's children.
"""
compute_len(vals::Union{NamedTuple, Tuple}) =
    mapreduce(f->compute_len(getfield(vals, f)), +, fieldnames(typeof(vals)))
compute_len(vals::AbstractArray) = mapreduce(compute_len, +, vals)
compute_len(vals::AbstractArray{<:Number}) = length(vals)

derivative_tuple(cnet::CachedNet) = derivative_tuple(cnet.net)
derivative_tuple(net) = weight_tuple(net, fieldnames(typeof(net)),
                                     Vector{out_type(net)}())[2]

"""
    weight_tuple(net)

Returns a named tuple holding all the fields in the network, and an
extra field named `tuple_all_weights` who has the same type
"""
weight_tuple(cnet::CachedNet) = weight_tuple(cnet.net)
weight_tuple(net) = weight_tuple(net, fieldnames(typeof(net)))[2]
function weight_tuple(x, fnames, vec=Vector{weight_type(x)}(), start=1)
    i = 0
    d=Dict{Symbol, Any}()
    for f=fnames
        di, val = weight_tuple(getfield(x,f), vec, start+i)
        i += di; push!(d, f=>val)
    end
    #start == 1 && push!(d, :tuple_all_weights=>[vec])
    i, (;d...)
end

function weight_tuple(x::Tuple, vec::Vector, start)
    i = 0
    d=Vector()
    for f=fieldnames(typeof(x))
        di, val = weight_tuple(getfield(x,f), vec, start+i)
        i += di; push!(d, val)
    end
    i, Tuple(d)
end

function weight_tuple(x::AbstractArray{<:Number}, vec::Vector, start)
    length(vec) < start+length(x)-1 && resize!(vec, start+length(x)-1)
    @views data_vec = vec[start:start+length(x)-1]
    reshpd_params = reshape(data_vec, size(x))
    reshpd_params .= x
    return length(x), reshpd_params
end

##
function batched_weight_tuple(grad_tup, bsz=1)
    all_weights = grad_tup.tuple_all_weights
    all_weghts_new = [similar(w, length(w), bsz) for w=all_weights][1]
    return batched_weight_tuple(grad_tup, fieldnames(typeof(grad_tup)), all_weghts_new)[2]
end

function batched_weight_tuple(x, fnames, vec, start=1)
    i = 0
    d=Dict{Symbol, Any}()
    for f=fnames
        f==:tuple_all_weights && continue
        di, val = batched_weight_tuple(getfield(x,f), vec, start+i)
        i += di; push!(d, f=>val)
    end
    start == 1 && push!(d, :tuple_all_weights=>[vec])
    i, (;d...)
end

function batched_weight_tuple(x::Tuple, vec::Matrix, start)
    i = 0
    d=Vector()
    for f=fieldnames(typeof(x))
        di, val = batched_weight_tuple(getfield(x,f), vec, start+i)
        i += di; push!(d, val)
    end
    i, Tuple(d)
end

function batched_weight_tuple(x::AbstractArray{<:Number}, vec::Matrix, start)
    @assert length(vec) >= start+length(x)-1
    bsz = size(vec, 2)

    @views data_vec = vec[start:start+length(x)-1, :]
    reshpd_params = reshape(data_vec, size(x)..., bsz)
    reshpd_params .= x
    if reshpd_params isa Base.ReshapedArray
        reshpd_params = StridedView(reshpd_params)
    end
    return length(x), reshpd_params
end # ? stridedView?


text:
struct ALConstraintSet{T}
    constraints::Vector{ALConstraint{T}}
    c_max::Vector{T}
    μ_max::Vector{T}
    Zref::Ref{SampledTrajectory}
    Eref::Ref{CostExpansion}
end

function ALConstraintSet{T}() where T
    constraints = ALConstraint{T}[]
    c_max = T[]
    μ_max = T[]
    Z = SampledTrajectory(KnotPoint{Any,Any,Vector{T},T}[])
    Zref = Ref{RD.SampledTrajectory}(Z)
    E = CostExpansion{T}(0,0,0)
    Eref = Ref{CostExpansion}(E)

    ALConstraintSet{T}(constraints, c_max, μ_max, Zref, Eref)
end

function initialize!(conset::ALConstraintSet{T}, cons::TO.ConstraintList,
                     Z::SampledTrajectory, opts::SolverOptions,
                     costs, E=CostExpansion{T}(RD.dims(Z)[1:2]...)) where T
    N = length(Z)
    nx,nu = cons.nx, cons.nu
    @assert all(k->RD.state_dim(Z,k) == nx[k], 1:N)
    @assert all(k->RD.control_dim(Z,k) == nu[k], 1:N)
    @assert length(Z) == length(cons.p)
    ncon = length(cons)
    for i = 1:ncon
        alcon = ALConstraint{T}(Z, cons[i], cons.inds[i], costs, E, sig=cons.sigs[i],
                                diffmethod=cons.diffs[i], solveropts=opts)
        push!(conset.constraints, alcon)
    end
    resize!(conset.c_max, length(Z))
    resize!(conset.μ_max, ncon)
    conset.c_max .= 0
    conset.μ_max .= 0
    conset.Zref[] = Z
    conset.Eref[] = E
    return
end

# Indexing and Iteration
@inline Base.length(conset::ALConstraintSet) = length(conset.constraints)
@inline Base.getindex(conset::ALConstraintSet, i::Integer) = conset.constraints[i]
Base.firstindex(::ALConstraintSet) = 1
Base.lastindex(conset::ALConstraintSet) = length(conset.constraints)

function Base.iterate(conset::ALConstraintSet)
    isempty(conset.constraints) ? nothing : (conset.constraints[1], 1)
end

function Base.iterate(conset::ALConstraintSet, state::Int)
    state >= length(conset) ? nothing : (conset.constraints[state+1], state+1)
end

Base.IteratorSize(::ALConstraintSet) = Base.HasLength()
Base.IteratorEltype(::ALConstraintSet) = Base.HasEltype()
Base.eltype(::ALConstraintSet{T}) where T = ALConstraint{T}

# Methods
for method in (:evaluate_constraints!, :constraint_jacobians!)
    @eval function $method(conset::ALConstraintSet)
        for alcon in conset.constraints
            $method(alcon)
        end
    end
    @eval function $method(conset::ALConstraintSet, Z::SampledTrajectory)
        # Hack to avoid an allocation from using a function barrier w/ more than 1 arg
        # Store a pointer to the trajectory in each conval
        # Only update the trajectory if it's different than the current one
        if Z !== conset.Zref[]
            settraj!(conset, Z)  # this allocates
        end
        $method(conset)  # this does not, but something like
                         # `evaluate_constraints!(conset, Z)` does
    end
end

function alcost(conset::ALConstraintSet{T}) where T
    for alcon in conset.constraints
        alcost(alcon)
    end
end

for method in (:algrad!, :alhess!, :dualupdate!, :penaltyupdate!,
        :max_penalty, :reset_duals!, :reset_penalties!)
    @eval function $method(conset::ALConstraintSet, args...)
        for alcon in conset.constraints
            $method(alcon, args...)
        end
    end
end

function add_alcost_expansion!(conset::ALConstraintSet, E::CostExpansion)
    for alcon in conset.constraints
        add_alcost_expansion!(alcon)
    end
end

function settraj!(conset::ALConstraintSet, Z::SampledTrajectory)
    conset.Zref[] = Z
    for alcon in conset.constraints
        settraj!(alcon, Z)
    end
end

function normviolation!(conset::ALConstraintSet, p=2)
    isempty(conset) && return 0.0
    conset.c_max .= 0
    for i = 1:length(conset)
        normviolation!(conset.constraints[i], p, conset.c_max)
    end
    return norm(conset.c_max, p)
end

# Need to duplicate due to allocation
function max_violation(conset::ALConstraintSet)
    isempty(conset) && return 0.0
    conset.c_max .= 0
    for i = 1:length(conset)
        normviolation!(conset.constraints[i], Inf, conset.c_max)
    end
    return norm(conset.c_max, Inf)
end

function max_penalty(conset::ALConstraintSet)
    for i = 1:length(conset)
        conset.μ_max[i] = max_penalty(conset.constraints[i])
    end
    return maximum(conset.μ_max)
end

function reset!(conset::ALConstraintSet)
    for con in conset.constraints
        resetparams!(con)
        reset_duals!(con)
        reset_penalties!(con)
    end
end

"""
	findmax_violation(conSet)

Return details on the where the largest violation occurs. Returns a string giving the
constraint type, time step index, and index into the constraint.
"""
function findmax_violation(conSet::ALConstraintSet{T}) where T
    c_max0 = -Inf
    j_con = -Inf
    for (i,alcon) in enumerate(conSet.constraints)
        max_violation!(alcon)
        v = maximum(alcon.c_max)
        if v > c_max0
            c_max0 = v
            j_con = i
        end
    end
	if c_max0 < eps()
		return "No constraints violated"
	end
	conval = conSet[j_con]
	k_con = argmax(conval.c_max) # which index
    i_con = searchsortedfirst(conval.inds, k_con)
	c_max, i_max = findmax(abs,conval.viol[i_con])  # index into constraint
	@assert c_max == c_max0
	con_name = string(typeof(conval.con).name.name)
	return con_name * " at time step $k_con at " * TO.con_label(conval.con, i_max)
end

ext: jl


text:
module Knk

export sctenifoldknk

import ScTenifold.Net

function sctenifoldknk(A::AbstractMatrix{T}, k::Integer) where T<:Real
    A1=copy(A)
    A1[k,:].=0
    return A1
    #d,_,_=Net.manialn(A,A1)
    #fc,p,adjp=Net.drgenes(d)
end

function sctenifoldknk(A::AbstractMatrix{T},genelist::Array{String,1},targetgene::String) where T<:Real
    # genelist=vec(readdlm("genelist.txt",String))
    k=findall(genelist.==targetgene)[1]
    # fc,p,adjp=sctenifoldknk(A, k)
    return sctenifoldknk(A, k)
end

function sctenifoldknk(A::AbstractMatrix{T},genelist::Array{String,1},k::Integer) where T<:Real
    # genelist=vec(readdlm("genelist.txt",String))
    # k=findall(genelist.==targetgene)[1]
    # fc,p,adjp=sctenifoldknk(A, k)
    return sctenifoldknk(A, k)
end

end

#=
function t1(genelist::Vector{String})
    genelist[1]
end
function t2(genelist::Array{String,1})
    genelist[1]
end

using DelimitedFiles
cd(dirname(@__FILE__))
X0=readdlm("X.txt",',',Int16)
X1=copy(X0)
genelist=vec(readdlm("genelist.txt",String))
idx=findall(genelist.=="Cftr")
X1[idx[1],:].=0

lbsz=sum(X0,dims=1)
X0=(X0./lbsz)*1e4
lbsz=sum(X1,dims=1)
X1=(X1./lbsz)*1e4

include("scTenifoldNet.jl")
using .scTenifoldNet
#X0=rand(100,1000);
#X1=rand(100,1000);
@show Threads.nthreads()
@time Z0=scTenifoldNet.pcnet(X0')
@time Z1=scTenifoldNet.pcnet(X1')

@time d,aln0,aln1=scTenifoldNet.manialn(Z0,Z1)
fc,p,adjp=scTenifoldNet.drgenes(d)

writedlm("qvalues.txt",[fc p adjp])

using StatsPlots, Distributions
x=rand(Chisq(1), length(fc))
qqplot(x, fc)
end
=#

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using Plots
gr()
x = LinRange(0,2π,100)
plot(x, sin.(x), c = :black, legend =:bottomleft, label ="sin(x)",
    fg_legend = :red, bg_legend = nothing,xlab ="x",ylab ="f(x)",
    title = "Sine", guidefont=(18,"Times"),tickfont=(14, "Times"),
    legendfont=(13, "Times"), titlefont=(18, "Times"))
savefig("figSin1.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
# HALF-PLANES

"""Half plane of equation ax + by + c <= 0."""
struct HalfPlane{T} <: Surface
    a::T
    b::T
    c::T
end

HalfPlane(a, b, c) = HalfPlane(promote(a, b, c)...)

function HalfPlane{T}(p1::Point{T}, p2::Point{T}, side) where T
    a = - p2[2] + p1[2]
    b = + p2[1] - p1[1]
    c = - a*p1[1] - b*p1[2]
    if side==:right
        HalfPlane{T}(a, b, c)
    else
        HalfPlane{T}(-a, -b, -c)
    end
end
HalfPlane(p1, p2, side::Symbol=:right) = (T = eltype(p1); HalfPlane{T}(Point{T}(p1...), Point{T}(p2...), side))

# POLAR COORDINATES
# The line is defined in polar coordinates around the center (xc, yc).
# The angle θ is the angle of normal vector.
#
#   ┏━━━━━━┓          ┏━━━━━━┓           ┏━━━━━━┓           ┏━━━━━━┓
#   ┃▒▒▒   ┃          ┃▒     ┃           ┃      ┃           ┃     ▒┃
#   ┃▒▒▒   ┃ θ = 0    ┃▒▒    ┃ θ = π/4   ┃      ┃ θ = π/2   ┃    ▒▒┃ θ = 3π/4
#   ┃▒▒▒   ┃          ┃▒▒▒▒  ┃           ┃▒▒▒▒▒▒┃           ┃  ▒▒▒▒┃
#   ┃▒▒▒   ┃          ┃▒▒▒▒▒ ┃           ┃▒▒▒▒▒▒┃           ┃ ▒▒▒▒▒┃
#   ┗━━━━━━┛          ┗━━━━━━┛           ┗━━━━━━┛           ┗━━━━━━┛
#   ┏━━━━━━┓          ┏━━━━━━┓           ┏━━━━━━┓           ┏━━━━━━┓
#   ┃   ▒▒▒┃          ┃ ▒▒▒▒▒┃           ┃▒▒▒▒▒▒┃           ┃▒▒▒▒▒ ┃
#   ┃   ▒▒▒┃ θ = π    ┃  ▒▒▒▒┃ θ = 5π/4  ┃▒▒▒▒▒▒┃ θ = 3π/2  ┃▒▒▒▒  ┃ θ = 7π/4 = -π/4
#   ┃   ▒▒▒┃          ┃    ▒▒┃           ┃      ┃           ┃▒▒    ┃
#   ┃   ▒▒▒┃          ┃     ▒┃           ┃      ┃           ┃▒     ┃
#   ┗━━━━━━┛          ┗━━━━━━┛           ┗━━━━━━┛           ┗━━━━━━┛
PolarHalfPlane(r, θ; center=Point(0.0, 0.0)) = HalfPlane(cos(θ), sin(θ), -cos(θ)*center[1] - sin(θ)*center[2] + r)
PolarHalfPlane(θ; center=Point(0.0, 0.0)) = PolarHalfPlane(0.0, θ; center=center)

convert(::Type{HalfPlane{T}}, h::HalfPlane{U}) where {T, U} = HalfPlane{T}(convert(T, h.a), convert(T, h.b), convert(T, h.c))
promote_rule(::Type{HalfPlane{T}}, ::Type{HalfPlane{U}}) where {T, U} = HalfPlane{promote_type(T, U)}

isapprox(h1::HalfPlane, h2::HalfPlane; kw...) = all((isapprox(h1.a, h2.a; kw...), isapprox(h1.b, h2.b; kw...), isapprox(h1.c, h2.c; kw...)))

equation(h::HalfPlane) = (x, y) -> h.a*x + h.b*y + h.c
in(p, h::HalfPlane) = equation(h)(p[1], p[2]) <= 0.0
signed_distance(p, h::HalfPlane) = equation(h)(p[1], p[2])/hypot(h.a, h.b)
distance(p, h::HalfPlane) = max(signed_distance(p, h), 0.0)

outward_normal(h::HalfPlane{T}) where T = SVector{2, T}(h.a, h.b)/hypot(h.a, h.b)
angle(h::HalfPlane{T}) where T = mod(atan(h.b, h.a), 2*T(π))

translate(h::HalfPlane, v) = (n = outward_normal(h); HalfPlane(h.a, h.b, h.c - v[1]*n[1] - v[2]*n[2]))
rotate(h::HalfPlane, ϕ; center=Point(0.0, 0.0)) = PolarHalfPlane(signed_distance(center, h), angle(h) + ϕ; center=center)
complement(h::HalfPlane) = HalfPlane(-h.a, -h.b, -h.c)
exchange_x_and_y(h::HalfPlane) = HalfPlane(h.b, h.a, h.c)

function corner_point(h1::HalfPlane{T}, h2::HalfPlane{U}) where {T, U}
    S = promote_type(U, T)
    idet = one(S) / (h1.a * h2.b - h2.a * h1.b)
    return Point{S}((-h2.b * h1.c + h1.b * h2.c) * idet, (-h1.a * h2.c + h2.a * h1.c) * idet)
end

isempty(h::HalfPlane) = false
area(h::HalfPlane) = Inf


text:
using LinearAlgebra
include("poisson_linear_system.jl")
include("simpson_integrate.jl")

function poisson_error(numelements)
    xL = -1.0
    xR = +1.0
    wavenumber = pi / 2
    bcleft = bcright = 0.0

    numpts = numelements + 1
    stepsize = (xR - xL) / numelements
    xrange = range(xL, stop = xR, length = numpts)

    rhs = cos.(wavenumber * xrange)
    op = -poisson_linear_system(numpts)
    rhs[1] = bcleft / stepsize^2
    rhs[end] = bcright / stepsize^2
    sol = stepsize^2 * (op \ rhs)
    exactsolution =
        1.0 / wavenumber^2 * (cos.(wavenumber * xrange) .- cos(wavenumber))
    err = simpson_integrate(sol - exactsolution, stepsize)
    return err
end

function convergence_rate(err,dx)
    return diff(log.(err)) ./ diff(log.(dx))
end

powers = 1:10
numelements = 2 .^ powers
stepsize = 2.0 ./ numelements
err = poisson_error.(numelements)
rate = convergence_rate(err,stepsize)

# using Plots
# plot(stepsize,err,xscale = :log, yscale = :log)
# scatter!(stepsize,err)

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
import Reactive
using BenchmarkTools
using StaticArrays
using Statistics
Reactive.async_mode.x = false

@testset "benchmark" begin
    println("")
    Signals.async_mode(false)
    n = 200
    signals = Vector{Signal}()
    push!(signals, Signal(1))
    for i = 1:n
        A = Signal(signals[i]) do x
            x + 1
        end
        push!(signals, A)
    end
    A = signals[1]
    bench = @benchmark begin; $A[] = 0; $signals[end](); end;
    ptime = median(bench).time
    println("Signals function call overhead on pull = $(ptime/n)ns") == nothing

    bench = @benchmark begin $A(1); end;
    ptime = median(bench).time
    println("Signals function call overhead on push = $(ptime/n)ns") == nothing

    signals = Vector{Reactive.Signal}()
    push!(signals, Reactive.Signal(1))
    for i = 1:n
        A = Reactive.foreach(signals[i]) do x
            x + 1
        end
        push!(signals, A)
    end
    A = signals[1]
    sleep(2)

    bench = @benchmark Reactive.push!($A, Reactive.value($A))
    rtime = median(bench).time

    println("Reactive function call overhead = $(rtime/n)ns") == nothing
    # @test ptime < rtime

    println("")

    typ = SMatrix{4, 4, Float64, 16}
    A = Signal(typ(rand(4, 4)))
    B = Signal(typ(rand(4, 4)))
    C = Signal(typ(rand(4, 4)))
    D = Signal(typ(rand(4, 4)))
    E = Signal(A, B) do a, b
        a*b
    end
    F = Signal(C, D) do a, b
        a*b
    end
    G = Signal(E, F) do e, f
        e*f
    end
    Z = typ(zeros(4, 4))
    bench = @benchmark begin; $A[] = $Z; $G(); end
    ptime = median(bench).time
    println("Signals function call time on pull (4x4 SArray multiply)= $(ptime/2)ns") == nothing
    bench = @benchmark begin $A($Z); end
    ptime = median(bench).time
    println("Signals function call time on push (4x4 SArray multiply) = $(ptime/2)ns") == nothing

    A = Reactive.Signal(typ(rand(4, 4)))
    B = Reactive.Signal(typ(rand(4, 4)))
    C = Reactive.Signal(typ(rand(4, 4)))
    D = Reactive.Signal(typ(rand(4, 4)))
    E = Reactive.foreach(A, B) do a, b
        a*b
    end
    F = Reactive.foreach(C, D) do a, b
        a*b
    end
    G = Reactive.foreach(E, F) do e, f
        e*f
    end
    Z = typ(zeros(4, 4))
    Reactive.async_mode.x = false

    bench = @benchmark Reactive.push!($A, Reactive.value($A))
    rtime = median(bench).time
    println("Reactive function call time on push (4x4 SArray multiply) = $(rtime/2)ns") == nothing
    # @test ptime < rtime

    println("")

    typ = Matrix
    A = Signal(typ(rand(4, 4)))
    B = Signal(typ(rand(4, 4)))
    C = Signal(typ(rand(4, 4)))
    D = Signal(typ(rand(4, 4)))
    E = Signal(A, B) do a, b
        a*b
    end
    F = Signal(C, D) do a, b
        a*b
    end
    G = Signal(E, F) do e, f
        e*f
    end
    Z = typ(zeros(4, 4))
    bench = @benchmark begin; $A[] = $Z; $G(); end
    ptime = median(bench).time
    println("Signals function call time on pull (4x4 Matrix)= $(ptime/2)ns") == nothing
    bench = @benchmark begin $A($Z); end
    ptime = median(bench).time
    println("Signals function call time on push (4x4 Matrix) = $(ptime/2)ns") == nothing

    A = Reactive.Signal(typ(rand(4, 4)))
    B = Reactive.Signal(typ(rand(4, 4)))
    C = Reactive.Signal(typ(rand(4, 4)))
    D = Reactive.Signal(typ(rand(4, 4)))
    E = Reactive.foreach(A, B) do a, b
        a*b
    end
    F = Reactive.foreach(C, D) do a, b
        a*b
    end
    G = Reactive.foreach(E, F) do e, f
        e*f
    end
    Z = typ(zeros(4, 4))
    bench = @benchmark Reactive.push!($A, Reactive.value($A))
    rtime = median(bench).time
    println("Reactive function call time on push (4x4 Matrix) = $(rtime/2)ns") == nothing
    # @test ptime < rtime

end


text:
module ReferenceFEInterfacesTests

using Test
using Gridap.Fields
using Gridap.Polynomials
using Gridap.ReferenceFEs

D = 2
T = Float64
order = 1
prebasis = MonomialBasis{D}(T,order)

polytope = QUAD
x = get_vertex_coordinates(polytope)
dofs = LagrangianDofBasis(T,x)

ndofs = 4

face_own_dofs = [[1],[2],[3],[4],Int[],Int[],Int[],Int[],Int[]]
a = [1]
b = Int[]
face_own_dofs_permutations = [[a],[a],[a],[a],[b,b],[b,b],[b,b],[b,b],fill(b,8)]
face_dofs = [[1],[2],[3],[4],[1,2],[3,4],[1,3],[2,4],[1,2,3,4]]
metadata = nothing

struct MockConformity <: Conformity end

conf = MockConformity()
reffe = GenericRefFE{typeof(conf)}(
  ndofs, polytope, prebasis, dofs, conf, metadata ,face_dofs)

function ReferenceFEs.get_face_own_dofs(reffe::GenericRefFE{MockConformity},::MockConformity)
  face_own_dofs
end

function ReferenceFEs.get_face_own_dofs_permutations(reffe::GenericRefFE{MockConformity},::MockConformity)
  face_own_dofs_permutations
end

test_reference_fe(reffe)

@test get_face_own_dofs(reffe,L2Conformity()) == [[], [], [], [], [], [], [], [], [1, 2, 3, 4]]
@test get_face_own_dofs_permutations(reffe,L2Conformity()) == [[[]], [[]], [[]], [[]], [[]], [[]], [[]], [[]], [[1, 2, 3, 4]]]

shapefuns = get_shapefuns(reffe)

@test evaluate(shapefuns,x) == [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]

@test get_own_dofs_permutations(reffe) == face_own_dofs_permutations[end]

end # module


text:
"""
    distance_lat_lon_euclidean(coords1::LatLonCoords, coords2::LatLonCoords)

Compute the Euclidean distance in kilometresbetween two lat-long coordinates close to each other.
"""
function distance_lat_lon_euclidean(coords1::LatLonCoords, coords2::LatLonCoords)
    deglen = 110.25
    x = coords1.lat - coords2.lat
    y = (coords1.lon - coords2.lon)*cos(coords2.lat)
    return deglen*sqrt(x^2 + y^2)
end

function is_in_bounds(params::CityParams, coords::LatLonCoords)

    return (coords.lat >= params.lat_start && coords.lat <= params.lat_end) &&
            (coords.lon >= params.lon_start && coords.lon <= params.lon_end)
end

function hhmmss_to_seconds(time_str::String)

    str_arr = split(time_str, ":")
    int_arr = [parse(Int64, s) for s in str_arr]
    time_arr = [3600., 60., 1.]

    return dot(int_arr, time_arr)
end

"""
Sample locations from a Uniform distribution within the bounding box.
"""
function generate_sites(coords_bb_lo::LatLonCoords, coords_bb_hi::LatLonCoords,
                        n_sites::Int64, rng::RNG) where {RNG <: AbstractRNG}

    sites = Vector{LatLonCoords}(undef, n_sites)

    for i = 1:n_sites

        lat_dist = Uniform(coords_bb_lo.lat, coords_bb_hi.lat)
        lon_dist = Uniform(coords_bb_lo.lon, coords_bb_hi.lon)

        sites[i] = LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist)))

    end

    return sites

end # function

"""
    get_non_dominated_trip_points(env::MAPFTransitEnv, avoid_vertex_idxs::Set{Int64},
                                  vtx::MAPFTransitVertexState, tid::Int64)

Use Kung's algorithm for d = 2 (very simple) to get the stops on a trip that
are non-dominated by other stops on the same trip for the (elapsed time, flight distance) metric
 - See Appendix for details.
"""
function get_non_dominated_trip_points(env::MAPFTransitEnv, avoid_vertex_idxs::Set{Int64},
                                       vtx::MAPFTransitVertexState, tid::Int64)

    trip_vtx_range = env.trip_to_vtx_range[tid]

    time_dist_set = Vector{Tuple{Int64,Float64,Float64}}(undef, 0)

    # First run through and put in list if reachable by agent
    for seq = trip_vtx_range[1]:trip_vtx_range[2]
        dist = env.dist_fn(vtx.state.location, env.state_graph.vertices[seq].state.location)
        tdiff = env.state_graph.vertices[seq].state.time - vtx.state.time

        if dist < env.drone_params.max_distance && env.drone_params.avg_speed * tdiff > dist &&
            ~(seq in avoid_vertex_idxs)
            push!(time_dist_set, (seq, tdiff, dist))
        end
    end

    if isempty(time_dist_set)
        return []
    end

    # YOU KNOW THIS IS ALREADY SORTED IN INCREASING ETA AND REACHABLE
    non_dom_idxs = Set{Int64}(1)

    # Start from second element
    for (i, elem) in enumerate(time_dist_set[2:end])
        dom = false
        for ndi in non_dom_idxs
            # If dist is worse, then dominated
            if elem[3] >= time_dist_set[ndi][3]
                dom = true
                break
            end
        end

        if dom == false
            push!(non_dom_idxs, i)
        end
    end

    # Now return all seq points from ndi
    seqs = [time_dist_set[ndi][1] for ndi in non_dom_idxs]

end # function

"""
    load_depots(depot_file::String)

Reads in a depot file and returns a vector of depot locations.
"""
function load_depots(depot_file::String)

    open(depot_file, "r") do df
        depot_dict = JSON.parse(df)
    end

    n_depots = length(collect(keys(depot_dict)))

    depots = Vector{LatLonCoords}(undef, n_depots)

    # Enforcing keys to be integers
    for (idx_str, latlon) in depot_dict
        idx = parse(Int64, id_str)
        depots[idx] = LatLonCoords((lat = latlon["lat"], lon = latlon["lon"]))
    end

    return depots
end

"""
    load_stop_to_location(::Type{LatLonCoords}, stop_coords_file::String)

Loads up the "stop_to_coords.json" dictionary file from GTFS parsing.
"""
function load_stop_to_location(::Type{LatLonCoords}, stop_coords_file::String)

    stop_coord_dict = Dict()
    open(stop_coords_file, "r") do f
        stop_coord_dict = JSON.parse(f)
    end

    stop_to_location = Dict{Int64,LatLonCoords}()

    for (id_str, latlon) in stop_coord_dict
        stop_to_location[parse(Int64, id_str)] = LatLonCoords((lat = latlon["lat"], lon = latlon["lon"]))
    end

    return stop_to_location
end

"""
    load_transit_routes_single_time(trip_file::String)

Loads the transit routes from the trips.json file from GTFS parsing.
"""
function load_transit_routes_single_time(trip_file::String)

    trips_dict = Dict()
    open(trip_file, "r") do f
        trips_dict = JSON.parse(f)
    end

    n_trips = length(collect(keys(trips_dict)))
    transit_trips = Vector{Vector{RouteWaypoint}}(undef, n_trips)

    for trip_num = 1:n_trips

        trip = trips_dict[string(trip_num)]
        n_stops = length(collect(keys(trip)))

        this_trip = Vector{RouteWaypoint}(undef, n_stops)

        for stop_seq = 1:n_stops

            # Cross-reference stop sequence from trip dictionary
            id_time = trip[string(stop_seq)]

            # Get info for route waypoint
            waypt = RouteWaypoint(stop_id = id_time["stop_id"],
                                  arrival_time = id_time["arrival_time"])

            this_trip[stop_seq] = waypt
        end

        transit_trips[trip_num] = this_trip
    end

    return transit_trips
end

## For nearestneighbor metric
struct EuclideanLatLong <: Metric
end

function Distances.evaluate(::EuclideanLatLong,
                            x::AbstractVector{Float64},
                            y::AbstractVector{Float64})
    coords1 = (lat=x[1], lon=x[2])
    coords2 = (lat=y[1], lon=y[2])
    return distance_lat_lon_euclidean(coords1, coords2)
end

"""
    get_halton_value(index::Int64, base::Int64)

Returns the halton sequence value for a given index and base.
"""
function get_halton_value(index::Int64, base::Int64)

    res = 0
    f = 1

    while index > 0
        f = f*1.0/base
        res += f*(index%base)
        index = div(index, base)
    end

    return res
end

"""
    get_halton_sequence(n::Int64, bases::Vector{Int64}, lower::Vector{Float64}, upper::Vector{Float64},
                        discard::Int64)

Sample d-dimensional Halton points between lower and upper limits. Discard the first x points,
where x is the value of the discard::Int64 parameter.
"""
function get_halton_sequence(n::Int64, bases::Vector{Int64}, lower::Vector{Float64}, upper::Vector{Float64},
                             discard::Int64)

    @assert length(bases) == length(lower) && length(lower) == length(upper)

    ndims = length(bases)
    diff = upper - lower
    grid_points = [SVector{ndims,Float64}([get_halton_value(idx, b) for b in bases])  for idx = 1:n+discard]

    grid_scaled_points = [SVector{ndims,Float64}(lower + gp.*diff) for gp in grid_points]

    return grid_scaled_points
end

"""
    plot_depots_sites!(background::Plots.Plot,side_x::Int64, side_y::Int64, bb_params::CityParams, depots::Vector{LatLonCoords}, sites::Vector{LatLonCoords},
                            depot_size::Int64=15, site_size::Int64=5)

Render the depot and package delivery locations on the map background.
"""
function plot_depots_sites!(background::Plots.Plot,side_x::Int64, side_y::Int64, bb_params::CityParams, depots::Vector{LatLonCoords}, sites::Vector{LatLonCoords},
                            depot_size::Int64=15, site_size::Int64=5)

    londiff = bb_params.lon_end - bb_params.lon_start
    latdiff = bb_params.lat_end - bb_params.lat_start

    for dep in depots

        xpos = convert(Int64, round((dep.lon - bb_params.lon_start)*side_x/londiff))
        ypos = convert(Int64, round((dep.lat - bb_params.lat_start)*side_y/latdiff))

        scatter!(background, [xpos], [ypos], markershape=:pentagon, markersize=depot_size, markercolor=:black, markerstrokewidth=0)
    end

    for site in sites

        xpos = convert(Int64, round((site.lon - bb_params.lon_start)*side_x/londiff))
        ypos = convert(Int64, round((site.lat - bb_params.lat_start)*side_y/latdiff))

        scatter!(background, [xpos], [ypos], markershape=:rect, markersize=site_size, markercolor=:grey, markerstrokecolor=:black)
    end
end

"""
    render_drones(background::Plots.Plot, side_x::Int64, side_y::Int64, bb_params::CityParams, time_val::Float64,
                       solution::Vector{PR}, drone_size::Int64=10) where {PR <: PlanResult}

Takes as input the static background obtained from plot_depots_sites! and plots the drone states.
"""
function render_drones(background::Plots.Plot, side_x::Int64, side_y::Int64, bb_params::CityParams, time_val::Float64,
                       solution::Vector{PR}, drone_size::Int64=10) where {PR <: PlanResult}

    p = deepcopy(background)

    # Loop over agents and add if still flying
    n_agents = length(solution)

    londiff = bb_params.lon_end - bb_params.lon_start
    latdiff = bb_params.lat_end - bb_params.lat_start

    for idx = 1:n_agents

        agt_soln_states = solution[idx].states

        pre_state_idx = 0
        post_state_idx = 0

        for (si, (state, t)) in enumerate(agt_soln_states)
            if t >= time_val
                post_state_idx = si
                break
            end
        end

        if post_state_idx > 0 && post_state_idx <= length(agt_soln_states)

            @assert post_state_idx > 1 "$(time_val) is too short for agent $(idx)!"

            pre_state_idx = post_state_idx - 1

            pre_time = agt_soln_states[pre_state_idx][2]
            post_time = agt_soln_states[post_state_idx][2]

            pre_loc = agt_soln_states[pre_state_idx][1].state.location
            post_loc = agt_soln_states[post_state_idx][1].state.location

            # Get interpolated position
            interp_factor = (time_val - pre_time)/(post_time - pre_time)
            new_lat = pre_loc.lat + interp_factor*(post_loc.lat - pre_loc.lat)
            new_lon = pre_loc.lon + interp_factor*(post_loc.lon - pre_loc.lon)

            # Only on car if previous and current are both r
            pre_str = agt_soln_states[pre_state_idx][1].vertex_str
            post_str = agt_soln_states[post_state_idx][1].vertex_str

            if startswith(pre_str, "r") && startswith(post_str, "r")
                dr_color = :red
            else
                dr_color = :blue
            end

            # Now obtain normalized position
            # Longitude corresponds to x
            xpos = convert(Int64, round((new_lon - bb_params.lon_start)*side_x/londiff))
            ypos = convert(Int64, round((new_lat - bb_params.lat_start)*side_y/latdiff))

            scatter!(p, [xpos], [ypos], markershape=:circle, markersize=drone_size, markercolor=dr_color, markerstrokecolor=:black)

        end
    end

    return p
end


text:
module QuasiGLM

using DataFrames, Distributions, GLM, PrettyTables

include("PeripheralFunctions.jl")
include("AdjustQuasiGLM.jl")

# Exports

export PearsonResiduals
export EstimateDispersionParameter
export coefarray
export AdjustQuasiGLM

end


text:
# # Geostrophic adjustment in the hydrostatic Boussinesq equations
#
# This example simulates a one-dimensional geostrophic adjustement problem
# using the `ClimateMachine.Ocean` subcomponent to solve the hydrostatic
# Boussinesq equations.
#
# First we `ClimateMachine.init()`.

using ClimateMachine

ClimateMachine.init()

# # Domain setup
#
# We formulate our problem in a Cartesian domain 100 km in ``x, y`` and 400 m
# deep, and discretized on a grid with 100 fourth-order elements in ``x``, and 1
# fourth-order element in ``y, z``,

using ClimateMachine.Ocean.Domains

domain = RectangularDomain(
    Ne = (25, 1, 1),
    Np = 4,
    x = (0, 1e6),
    y = (0, 1e6),
    z = (-400, 0),
    periodicity = (false, true, false),
)

# # Physical parameters
#
# We use a Coriolis parameter appropriate for mid-latitudes,

f = 1e-4 # s⁻¹, Coriolis parameter
nothing # hide

# and Earth's gravitational acceleration,

using CLIMAParameters: AbstractEarthParameterSet, Planet
struct EarthParameters <: AbstractEarthParameterSet end

g = Planet.grav(EarthParameters()) # m s⁻²

# # An unbalanced initial state
#
# We use a Gaussian, partially-balanced initial condition with parameters

U = 0.1              # geostrophic velocity (m s⁻¹)
L = domain.L.x / 40  # Gaussian width (m)
a = f * U * L / g    # amplitude of the geostrophic surface displacement (m)
x₀ = domain.L.x / 4  # Gaussian origin (m, recall that x ∈ [0, Lx])

# and functional form

Gaussian(x, L) = exp(-x^2 / (2 * L^2))

## Geostrophic ``y``-velocity: f V = g ∂_x η
vᵍ(x, y, z) = -U * (x - x₀) / L * Gaussian(x - x₀, L)

## Geostrophic surface displacement
ηᵍ(x, y, z) = a * Gaussian(x - x₀, L)

# We double the initial surface displacement so that the surface is half-balanced,
# half unbalanced,

ηⁱ(x, y, z) = 2 * ηᵍ(x, y, z)

# In summary,

using ClimateMachine.Ocean.OceanProblems: InitialConditions

initial_conditions = InitialConditions(v = vᵍ, η = ηⁱ)

@info """ Parameters for the Geostrophic adjustment problem are...

    Coriolis parameter:                            $f s⁻¹
    Gravitational acceleration:                    $g m s⁻²
    Geostrophic velocity:                          $U m s⁻¹
    Width of the initial geostrophic perturbation: $L m
    Amplitude of the initial surface perturbation: $a m
    Rossby number (U / f L):                       $(U / (f * L))

"""

# # Boundary conditions, Driver configuration, and Solver configuration
#
# Next, we configure the `HydrostaticBoussinesqModel` and build the `DriverConfiguration`.
# We configure our model in a domain which is bounded in the ``x`` direction.
# Both the boundary conditions in ``x`` and in ``z`` require boundary conditions,
# which we define:

using ClimateMachine.Ocean:
    Impenetrable, Penetrable, FreeSlip, Insulating, OceanBC

solid_surface_boundary_conditions = OceanBC(
    Impenetrable(FreeSlip()), # Velocity boundary conditions
    Insulating(),             # Temperature boundary conditions
)

free_surface_boundary_conditions = OceanBC(
    Penetrable(FreeSlip()),   # Velocity boundary conditions
    Insulating(),             # Temperature boundary conditions
)

boundary_conditions =
    (solid_surface_boundary_conditions, free_surface_boundary_conditions)

# We refer to these boundary conditions by their indices in the `boundary_tags` tuple
# when specifying the boundary conditions for the `state`; in other words, "1" corresponds to
# `solid_surface_boundary_conditions`, while `2` corresponds to `free_surface_boundary_conditions`,

boundary_tags = (
    (1, 1), # (west, east) boundary conditions
    (0, 0), # (south, north) boundary conditions
    (1, 2), # (bottom, top) boundary conditions
)

# We're now ready to build the model.

using ClimateMachine.Ocean

model = Ocean.HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = 2.0,
    initial_conditions = initial_conditions,
    parameters = EarthParameters(),
    turbulence_closure = (νʰ = 0, κʰ = 0, νᶻ = 0, κᶻ = 0),
    coriolis = (f₀ = f, β = 0),
    boundary_tags = boundary_tags,
    boundary_conditions = boundary_conditions,
)

# !!! info "Horizontallly-periodic boundary conditions"
#     To set horizontally-periodic boundary conditions with
#     `(solid_surface_boundary_conditions, free_surface_boundary_conditions)`
#     in the vertical direction use `periodicity = (true, true, false)` in
#     the `domain` constructor and `boundary_tags = ((0, 0), (0, 0), (1, 2))`
#     in the constructor for `HydrostaticBoussinesqSuperModel`.

# # Animating the solution
#
# To animate the `ClimateMachine.Ocean` solution, we'll create a callback
# that draws a plot and stores it in an array. When the simulation is finished,
# we'll string together the plotted frames into an animation.

using Printf
using Plots
using ClimateMachine.GenericCallbacks: EveryXSimulationSteps
using ClimateMachine.Ocean.Fields: assemble
using ClimateMachine.Ocean: current_step, current_time

u, v, η, θ = model.fields

## Container to hold the plotted frames
movie_plots = []

plot_every = 200 # iterations

plot_maker = EveryXSimulationSteps(plot_every) do

    @info "Steps: $(current_step(model)), time: $(current_time(model))"

    assembled_u = assemble(u.elements)
    assembled_v = assemble(v.elements)
    assembled_η = assemble(η.elements)

    umax = 0.5 * max(maximum(abs, u), maximum(abs, v))
    ulim = (-umax, umax)

    u_plot = plot(
        assembled_u.x[:, 1, 1],
        [assembled_u.data[:, 1, 1] assembled_v.data[:, 1, 1]],
        xlim = domain.x,
        ylim = (-0.7U, 0.7U),
        label = ["u" "v"],
        linewidth = 2,
        xlabel = "x (m)",
        ylabel = "Velocities (m s⁻¹)",
    )

    η_plot = plot(
        assembled_η.x[:, 1, 1],
        assembled_η.data[:, 1, 1],
        xlim = domain.x,
        ylim = (-0.01a, 1.2a),
        linewidth = 2,
        label = nothing,
        xlabel = "x (m)",
        ylabel = "η (m)",
    )

    push!(movie_plots, (u = u_plot, η = η_plot, time = current_time(model)))

    return nothing
end

# # Running the simulation and animating the results
#
# Finally, we run the simulation,

hours = 3600.0

model.solver_configuration.timeend = 2hours

result = ClimateMachine.invoke!(
    model.solver_configuration;
    user_callbacks = [plot_maker],
)

# and animate the results,

animation = @animate for p in movie_plots
    title = @sprintf("Geostrophic adjustment at t = %.2f hours", p.time / hours)
    frame = plot(
        p.u,
        p.η,
        layout = (2, 1),
        size = (800, 600),
        title = [title ""],
    )
end

gif(animation, "geostrophic_adjustment.gif", fps = 8) # hide

max_forks_repo_name: jm-c/CLIMA


text:
using OVERTVerify

state_vars = tora_input_vars
control_vars = tora_control_vars
# input_set_hyperrect = Hyperrectangle(low=[0.6, -0.7, -0.4, 0.5], high=[0.7, -0.6, -0.3, 0.6])
ranges = [[z...] for z in zip([0.6, -0.7, -0.4, 0.5], [0.7, -0.6, -0.3, 0.6])]
input_dict = Dict(zip(state_vars, ranges))
controller_file = "nnet_files/jmlr/tora_smallest_controller.nnet"
dynamics_map = Dict(state_vars[1]=>state_vars[2],
                    state_vars[2]=>tora_dim2,
                    state_vars[3]=>state_vars[4],
                    state_vars[4]=>tora_control_vars[1]
                    )
dt = 0.1
output_constraints = [:(x1 <= -2.0), :(x1 >= 2.0)] # (avoid set)
N_steps=15
experiment_name = "tora_smallest_controller"
dirname="src/examples/jmlr/comparisons/"

ΔT = compare_to_dreal(state_vars, control_vars, input_dict, controller_file, dynamics_map, dt, output_constraints, dirname, experiment_name, N_steps; jobs=28, dreal_path="/opt/dreal/4.20.12.1/bin/dreal")


text:
function get_normalization(block::T) where T
    hasproperty(block, :mult) ? block.mult : block.mults
end

Base.nameof(block::AbstractBlockHeader) = nameof(block.base_header)

labels(block::PointMeshBlockHeader) = block.labels
labels(block::PlainMeshBlockHeader) = block.labels

function add_units(raw_data::NTuple, block)
    𝟙 = one(eltype(block))
    T = typeof(𝟙*get_units(block.units)[1])
    map(data->reinterpret(T, data), raw_data)
end

function add_units(raw_data, block)
    𝟙 = one(eltype(block))
    reinterpret(typeof(𝟙*get_units(block.units)), raw_data)
end

function apply_normalization!(raw_data, block)
    𝟙 = one(eltype(block))

    n = get_normalization(block)
    if n isa Number
        default_factor = 𝟙
    else
        default_factor = ntuple(one, ndims(block))
    end

    if n ≠ default_factor
        raw_data .*= get_normalization(block)
    end

    return nothing
end


text:
module SeisOptimization

    using Printf, LinearAlgebra, SparseArrays

    import Base.convert, Base.show

    include("linear/linear.jl")

end # module


text:
module FluxBot

using GitHub, Glob
using GitHub.JSON, GitHub.HTTP
using Pkg.TOML

include("buildkite.jl")
include("trial.jl")
include("make_reply.jl")
include("zoo_ci.jl")

end # module


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays

@testset "Kronecker products and sums" begin
    @testset "Kronecker product" begin
        A = rand(ComplexF64, 3, 3)
        B = rand(ComplexF64, 2, 2)
        K = kron(A, B)
        LA = LinearMap(A)
        LB = LinearMap(B)
        LK = @inferred kron(LA, LB)
        @test kron(LA, 2LB) isa LinearMaps.ScaledMap
        @test kron(3LA, LB) isa LinearMaps.ScaledMap
        @test kron(3LA, 2LB) isa LinearMaps.ScaledMap
        @test kron(3LA, 2LB).λ == 6
        @test kron(kron(LA, LA), 2LB) isa LinearMaps.ScaledMap
        @test kron(3LA, kron(LB, LB)) isa LinearMaps.ScaledMap
        @test_throws ErrorException LinearMaps.KroneckerMap{Float64}((LA, LB))
        @test occursin("6×6 LinearMaps.KroneckerMap{$(eltype(LK))}", sprint((t, s) -> show(t, "text/plain", s), LK))
        @test @inferred size(LK) == size(K)
        @test LinearMaps.MulStyle(LK) === LinearMaps.ThreeArg()
        for i in (1, 2)
            @test @inferred size(LK, i) == size(K, i)
        end
        @test LK isa LinearMaps.KroneckerMap{ComplexF64}
        L = ones(3) ⊗ ones(ComplexF64, 4)'
        v = rand(4)
        @test Matrix(L) == ones(3,4)
        @test L*v ≈ fill(sum(v), 3)

        for transform in (identity, transpose, adjoint)
            @test Matrix(transform(LK)) ≈ transform(Matrix(LK)) ≈ transform(K)
            @test Matrix(kron(transform(LA), transform(LB))) ≈ transform(K)
            @test Matrix(transform(LinearMap(LK))) ≈ transform(Matrix(LK)) ≈ transform(K)
        end
        @test kron(LA, kron(LA, B)) == kron(LA, LA, LB)
        @test kron(kron(LA, LB), kron(LA, LB)) == kron(LA, LB, LA, LB) == ⊗(LA, LB, LA, LB)
        @test kron(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) ≈ Matrix(@inferred LA^⊗(3)) ≈ Matrix(@inferred A^⊗(3))
        LAs = LinearMap(sparse(A))
        K = @inferred kron(A, A, A, LAs)
        @test K isa LinearMaps.KroneckerMap
        @test Matrix(K) ≈ kron(A, A, A, A)
        @test convert(AbstractMatrix, K) isa SparseMatrixCSC
        @test convert(AbstractMatrix, K) ≈ kron(A, A, A, A)
        @test sparse(K) ≈ kron(A, A, A, A)
        @test Matrix(@inferred K*K) ≈ kron(A, A, A, A)*kron(A, A, A, A)
        K4 = @inferred kron(A, B, B, LB)
        # check that matrices don't get Kronecker-multiplied, but that everything is lazy
        @test K4.maps[1].lmap === A
        @test @inferred kron(LA, LB)' == @inferred kron(LA', LB')
        @test (@inferred kron(LA, B)) == (@inferred kron(LA, LB)) == (@inferred kron(A, LB))
        @test @inferred ishermitian(kron(LA'LA, LB'LB))
        A = rand(2, 5); B = rand(4, 2)
        K = @inferred kron(A, LinearMap(B))
        v = rand(size(K, 2))
        @test K*v ≈ kron(A, B)*v
        @test Matrix(K) ≈ kron(A, B)
        K = @inferred kron(LinearMap(B), A)
        v = rand(size(K, 2))
        @test K*v ≈ kron(B, A)*v
        @test Matrix(K) ≈ kron(B, A)
        A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B)
        @test @inferred issymmetric(kron(LA'LA, LB'LB))
        @test @inferred ishermitian(kron(LA'LA, LB'LB))
        # use mixed-product rule
        K = kron(LA, LB) * kron(LA, LB) * kron(LA, LB)
        @test Matrix(K) ≈ kron(A, B)^3
        # example that doesn't use mixed-product rule
        A = rand(3, 2); B = rand(2, 3)
        K = @inferred kron(A, LinearMap(B))
        @test Matrix(@inferred K*K) ≈ kron(A, B)*kron(A, B)
        A = rand(3, 2); B = rand(4, 3)
        @test Matrix(kron(LinearMap(A), B, [A A])*kron(LinearMap(A), B, A')) ≈ kron(A, B, [A A])*kron(A, B, A')

        m = 3
        A = rand(m, m)
        F = LinearMap(x -> A*x, m, m)
        S = sparse(I, m, m)
        J = LinearMap(I, m)
        @test kron(J, J) == LinearMap(I, m*m)
        v = rand(m^3)
        for (K, M) in ((⊗(A, J, J), kron(A, S, S)),
                       (⊗(J, A, J), kron(S, A, S)),
                       (⊗(J, J, A), kron(S, S, A)),
                       (⊗(F, J, J), kron(A, S, S)),
                       (⊗(J, F, J), kron(S, A, S)),
                       (⊗(J, J, F), kron(S, S, A)))
            @test K * v ≈ M * v
            @test Matrix(K) ≈ M
        end
    end

    @testset "Kronecker sum" begin
        for elty in (Float64, ComplexF64)
            A = rand(elty, 3, 3)
            B = rand(elty, 2, 2)
            LA = LinearMap(A)
            LB = LinearMap(B)
            KS = @inferred kronsum(LA, B)
            @test occursin("6×6 LinearMaps.KroneckerSumMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), KS))
            @test_throws ArgumentError kronsum(LA, [B B]) # non-square map
            KSmat = kron(A, Matrix(I, 2, 2)) + kron(Matrix(I, 3, 3), B)
            @test Matrix(KS) ≈ Matrix(kron(A, LinearMap(I, 2)) + kron(LinearMap(I, 3), B))
            @test KS * ones(size(KS, 2)) ≈ KSmat * ones(size(KS, 2))
            @test size(KS) == size(kron(A, Matrix(I, 2, 2)))
            for transform in (identity, transpose, adjoint)
                @test Matrix(transform(KS)) ≈ transform(Matrix(KS)) ≈ transform(KSmat)
                @test Matrix(kronsum(transform(LA), transform(LB))) ≈ transform(KSmat)
                @test Matrix(transform(LinearMap(kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat)
            end
            @test @inferred(kronsum(A, A, LB)) == @inferred(⊕(A, A, B))
            @test Matrix(@inferred LA^⊕(3)) == Matrix(@inferred A^⊕(3)) ≈ Matrix(kronsum(LA, A, A))
            @test @inferred(kronsum(LA, LA, LB)) == @inferred(kronsum(LA, kronsum(LA, LB))) == @inferred(kronsum(A, A, B))
            @test Matrix(@inferred kronsum(A, B, A, B, A, B)) ≈ Matrix(@inferred kronsum(LA, LB, LA, LB, LA, LB))
            T = typeof(kron(Diagonal(rand(elty, 3)), sprand(3, 3, 0.3)))
            @test convert(AbstractMatrix, kronsum(sparse(A), sparse(B), sparse(A))) isa T
            @test convert(AbstractMatrix, kronsum(A, B, A)) == Matrix(kronsum(A, B, A))
            @test sparse(kronsum(sparse(A), B, A)) == Matrix(kronsum(sparse(A), B, A))
        end
    end
end


text:
immutable Constant <: Degree{0} end

"""
Constant b-splines are *nearest-neighbor* interpolations, and effectively
return `A[round(Int,x)]` when interpolating
"""
Constant

"""
`define_indices_d` for a constant b-spline calculates `ix_d = round(Int,x_d)`
"""
function define_indices_d(::Type{BSpline{Constant}}, d, pad)
    symix, symx = Symbol("ix_",d), Symbol("x_",d)
    :($symix = clamp(round(Int, $symx), 1, size(itp, $d)))
end

"""
`coefficients` for a constant b-spline simply sets `c_d = 1` for compatibility
with the general b-spline framework
"""
function coefficients(::Type{BSpline{Constant}}, N, d)
    sym, symx = Symbol("c_",d), Symbol("x_",d)
    :($sym = 1)
end

"""
`gradient_coefficients` for a constant b-spline simply sets `c_d = 0` for
compatibility with the general b-spline framework
"""
function gradient_coefficients(::Type{BSpline{Constant}}, d)
    sym, symx = Symbol("c_",d), Symbol("x_",d)
    :($sym = 0)
end

"""
`hessian_coefficients` for a constant b-spline simply sets `c_d = 0` for
compatibility with the general b-spline framework
"""
function hessian_coefficients(::Type{BSpline{Constant}}, d)
    sym = Symbol("c_",d)
    :($sym = 0)
end

function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Constant}}, ::Type{IT}, N::Integer, offsets...)
    if (length(offsets) < N)
        d = length(offsets)+1
        sym = Symbol("c_", d)
        return :($sym * $(index_gen(IT, N, offsets..., 0)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end


text:
isna(X::DataValueArray, I::Int...) = X.isna[I...]
Base.values(X::DataValueArray, I::Int...) = X.values[I...]

"""
    size(X::DataValueArray, [d::Real])

Return a tuple containing the lengths of each dimension of `X`, or if `d` is
specific, the length of `X` along dimension `d`.
"""
Base.size(X::DataValueArray) = size(X.values)

function Base.similar(x::AbstractArray, ::Type{DataValue{T}}, dims::Dims) where {T}
    return DataValueArray{T}(dims)
end

function Base.similar(x::Array, ::Type{DataValue{T}}, dims::Dims) where {T}
    return DataValueArray{T}(dims)
end

function Base.similar(x::SubArray, ::Type{DataValue{T}}, dims::Dims) where {T}
    return DataValueArray{T}(dims)
end

"""
    copy(X::DataValueArray)

Return a shallow copy of `X`; the outer structure of `X` will be copied, but
all elements will be identical to those of `X`.
"""
function Base.copy(X::DataValueArray{T}) where {T}
    return Base.copyto!(similar(X, DataValue{T}), X)
end

# DA TODO This was my version, not clear which one is better
# function Base.copy!{T}(dest::DataValueArray{T},
#     src::DataValueArray{T})
#     length(dest) >= length(src) || throw(BoundsError())

#     n = length(src)

#     if isbitstype(T)
#         unsafe_copy!(pointer(dest.values, 1), pointer(src.values, 1), n)
#     else
#         ccall(:jl_array_ptr_copy, Nothing, (Any, Ptr{Nothing}, Any, Ptr{Nothing}, Int),
#             dest.values, pointer(dest.values, 1), src.values, pointer(src.values, 1), n)
#     end
#     unsafe_copy!(pointer(dest.isna, 1), pointer(src.isna, 1), n)
#     return dest
# end

"""
    copyto!(dest::DataValueArray, src::DataValueArray)

Copy the initialized values of a source DataValueArray into the respective
indices of the destination DataValueArray. If an entry in `src` is null, then
this method nullifies the respective entry in `dest`.
"""
function Base.copyto!(dest::DataValueArray, src::DataValueArray)
    if isbitstype(eltype(dest)) && isbitstype(eltype(src))
        copyto!(dest.values, src.values)
    else
        dest_values = dest.values
        src_values = src.values
        length(dest_values) >= length(src_values) || throw(BoundsError())
        # copy only initilialized values from src into dest
        for i in 1:length(src_values)
            @inbounds !(src.isna[i]) && (dest.values[i] = src.values[i])
        end
    end
    copyto!(dest.isna, src.isna)
    return dest
end

"""
    fill!(X::DataValueArray, x::DataValue)

Fill `X` with the value `x`. If `x` is empty, then `fill!(X, x)` nullifies each
entry of `X`. Otherwise, `fill!(X, x)` fills `X.values` with the value of `x`
and designates each entry of `X` as present.
"""
function Base.fill!(X::DataValueArray, x::DataValue)
    if isna(x)
        fill!(X.isna, true)
    else
        fill!(X.values, get(x))
        fill!(X.isna, false)
    end
    return X
end

"""
    fill!(X::DataValueArray, x::DataValue)

Fill `X` with the value `x` and designate each entry as present. If `x` is an
object reference, all elements will refer to the same object. Note that
`fill!(X, Foo())` will return `X` filled with the result of evaluating `Foo()`
once.
"""
function Base.fill!(X::DataValueArray, x::Any)
    fill!(X.values, x)
    fill!(X.isna, false)
    return X
end

"""
    Base.deepcopy(X::DataValueArray)

Return a `DataValueArray` object whose internal `values` and `isna` fields are
deep copies of `X.values` and `X.isna` respectively.
"""
function Base.deepcopy(X::DataValueArray)
    return DataValueArray(deepcopy(X.values), deepcopy(X.isna))
end

"""
    resize!(X::DataValueVector, n::Int)

Resize a one-dimensional `DataValueArray` `X` to contain precisely `n` elements.
If `n` is greater than the current length of `X`, then each new entry will be
designated as null.
"""
function Base.resize!(X::DataValueArray{T,1}, n::Int) where {T}
    resize!(X.values, n)
    oldn = length(X.isna)
    resize!(X.isna, n)
    X.isna[oldn+1:n] .= true
    return X
end

function Base.reshape(X::DataValueArray, dims::Dims)
    return DataValueArray(reshape(X.values, dims), reshape(X.isna, dims))
end

"""
    ndims(X::DataValueArray)

Returns the number of dimensions of `X`.
"""
Base.ndims(X::DataValueArray) = ndims(X.values)

"""
    length(X::DataValueArray)

Returns the maximum index `i` for which `getindex(X, i)` is valid.
"""
Base.length(X::DataValueArray) = length(X.values)

"""
    lastindex(X::DataValueArray)

Returns the last entry of `X`.
"""
Base.lastindex(X::DataValueArray) = lastindex(X.values)

# DA TODO Unclear whether I want that
# function Base.find(X::DataValueArray{Bool})
#     ntrue = 0
#     @inbounds for (i, isna) in enumerate(X.isna)
#         ntrue += !isna && X.values[i]
#     end
#     res = Array{Int}(ntrue)
#     ind = 1
#     @inbounds for (i, isna) in enumerate(X.isna)
#         if !isna && X.values[i]
#             res[ind] = i
#             ind += 1
#         end
#     end
#     return res
# end

"""
    dropna(X::AbstractVector)

Return a vector containing only the non-missing entries of `X`,
unwrapping `DataValue` entries. A copy is always returned, even when
`X` does not contain any missing values.
"""
function dropna(X::AbstractVector{T}) where {T}
    if !(DataValue <: T) && !(T <: DataValue)
        return copy(X)
    else
        Y = filter(x->!isna(x), X)
        res = similar(Y, eltype(T))
        for i in eachindex(Y, res)
            @inbounds res[i] = isa(Y[i], DataValue) ? Y[i].value : Y[i]
        end
        return res
    end
end
dropna(X::DataValueVector) = X.values[(!).(X.isna)]

"""
    dropna!(X::AbstractVector)

Remove missing entries of `X` in-place and return a `Vector` view of the
unwrapped `DataValue` entries. If no missing values are present, this is a no-op
and `X` is returned.
"""
function dropna!(X::AbstractVector{T}) where {T}                 # -> AbstractVector
    if !(DataValue <: T) && !(T <: DataValue)
        return X
    else
        deleteat!(X, findall(isna, X))
        res = similar(X, eltype(T))
        for i in eachindex(X, res)
            @inbounds res[i] = isa(X[i], DataValue) ? X[i].value : X[i]
        end
        return res
    end
end

# TODO: replace `find(X.isna)` with `X.isna` when
# https://github.com/JuliaLang/julia/pull/20465 is merged and part of
# current release (either v0.6 or v1.0)
"""
    dropna!(X::DataValueVector)

Remove missing entries of `X` in-place and return a `Vector` view of the
unwrapped `DataValue` entries.
"""
dropna!(X::DataValueVector) = deleteat!(X, (LinearIndices(X.isna))[findall(X.isna)]).values # -> Vector

# DA TODO I don't think we want this
# """
#     isnan(X::DataValueArray)

# Test whether each entry of `X` is null and if not, test whether the entry is
# not a number (`NaN`). Return the results as `DataValueArray{Bool}`. Note that
# null entries of `X` will be reflected by null entries of the resultant
# `DataValueArray`.
# """
# function Base.isnan(X::DataValueArray) # -> DataValueArray{Bool}
#     return DataValueArray(isnan.(X.values), copy(X.isna))
# end

# DA TODO I don't think we want this
# """
#     isfinite(X::DataValueArray)

# Test whether each entry of `X` is null and if not, test whether the entry is
# finite. Return the results as `DataValueArray{Bool}`. Note that
# null entries of `X` will be reflected by null entries of the resultant
# `DataValueArray`.
# """
# function Base.isfinite(X::DataValueArray) # -> DataValueArray{Bool}
#     res = Array{Bool}(size(X))
#     for i in eachindex(X)
#         if !X.isna[i]
#             res[i] = isfinite(X.values[i])
#         end
#     end
#     return DataValueArray(res, copy(X.isna))
# end

"""
    convert(T, X::DataValueArray)

Convert `X` to an `AbstractArray` of type `T`. Note that if `X` contains any
null entries then calling `convert` without supplying a replacement value for
null entries will result in an error.

Currently supported return type arguments include: `Array`, `Array{T}`,
`Vector`, `Matrix`.

    convert(T, X::DataValueArray, replacement)

Convert `X` to an `AbstractArray` of type `T` and replace all null entries of
`X` with `replacement` in the result.
"""
function Base.convert(::Type{Array{S, N}}, X::DataValueArray{T, N}) where {S,T,N}
    if any(isna, X)
        throw(DataValueException())
    else
        return convert(Array{S, N}, X.values)
    end
end

function Base.convert(::Type{Array{S}}, X::DataValueArray{T, N}) where {S,T,N} # -> Array{S, N}
    return convert(Array{S, N}, X)
end

function Base.convert(::Type{Vector}, X::DataValueVector{T}) where {T} # -> Vector{T}
    return convert(Array{T, 1}, X)
end

function Base.convert(::Type{Matrix}, X::DataValueMatrix{T}) where {T} # -> Matrix{T}
    return convert(Array{T, 2}, X)
end

function Base.convert(::Type{Array},
                            X::DataValueArray{T, N}) where {T,N} # -> Array{T, N}
    return convert(Array{T, N}, X)
end

# Conversions with replacements for handling null values

function Base.convert(::Type{Array{S, N}},
                               X::DataValueArray{T, N},
                               replacement::Any) where {S,T,N} # -> Array{S, N}
    replacementS = convert(S, replacement)
    res = Array{S}(undef, size(X))
    for i in 1:length(X)
        if X.isna[i]
            res[i] = replacementS
        else
            res[i] = X.values[i]
        end
    end
    return res
end

function Base.convert(::Type{Vector},
                         X::DataValueVector{T},
                         replacement::Any) where {T} # -> Vector{T}
    return convert(Array{T, 1}, X, replacement)
end

function Base.convert(::Type{Matrix},
                         X::DataValueMatrix{T},
                         replacement::Any) where {T} # -> Matrix{T}
    return convert(Array{T, 2}, X, replacement)
end

function Base.convert(::Type{Array},
                            X::DataValueArray{T, N},
                            replacement::Any) where {T,N} # -> Array{T, N}
    return convert(Array{T, N}, X, replacement)
end

Base.any(::typeof(isna), X::DataValueArray) = Base.any(X.isna)
Base.all(::typeof(isna), X::DataValueArray) = Base.all(X.isna)


text:
# ---
# title: 1178. Number of Valid Words for Each Puzzle
# id: problem1178
# author: Qling
# date: 2021-02-26
# difficulty: Hard
# categories: Hash Table, Bit Manipulation
# link: <https://leetcode.com/problems/number-of-valid-words-for-each-puzzle/description/>
# hidden: true
# ---
#
# With respect to a given `puzzle` string, a `word` is _valid_  if both the
# following conditions are satisfied:
#
#   * `word` contains the first letter of `puzzle`.
#   * For each letter in `word`, that letter is in `puzzle`.
# For example, if the puzzle is "abcdefg", then valid words are "faced",
# "cabbage", and "baggage"; while invalid words are "beefed" (doesn't include
# "a") and "based" (includes "s" which isn't in the puzzle).
#
# Return an array `answer`, where `answer[i]` is the number of words in the
# given word list `words` that are valid with respect to the puzzle
# `puzzles[i]`.
#
#
#
# **Example :**
#
#
#
#     Input:
#     words = ["aaaa","asas","able","ability","actt","actor","access"],
#     puzzles = ["aboveyz","abrodyz","abslute","absoryz","actresz","gaswxyz"]
#     Output: [1,1,3,2,4,0]
#     Explanation:
#     1 valid word for "aboveyz" : "aaaa"
#     1 valid word for "abrodyz" : "aaaa"
#     3 valid words for "abslute" : "aaaa", "asas", "able"
#     2 valid words for "absoryz" : "aaaa", "asas"
#     4 valid words for "actresz" : "aaaa", "asas", "actt", "access"
#     There're no valid words for "gaswxyz" cause none of the words in the list contains letter 'g'.
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 10^5`
#   * `4 <= words[i].length <= 50`
#   * `1 <= puzzles.length <= 10^4`
#   * `puzzles[i].length == 7`
#   * `words[i][j]`, `puzzles[i][j]` are English lowercase letters.
#   * Each `puzzles[i] `doesn't contain repeated characters.
#
#
## @lc code=start
using LeetCode

function get_bit_mask(word::String)::Int
    mask = 0

    for c in word
        i = Int(c) - Int('a')
        mask |= 1 << i
    end

    return mask
end

function find_num_of_valid_words(
    words::Vector{String},
    puzzles::Vector{String},
)::Vector{Int}
    letter_frequencies = Dict{Int,Int}()

    for word in words
        mask = get_bit_mask(word)
        letter_frequencies[mask] = get(letter_frequencies, mask, 0) + 1
    end

    solution = fill(0, length(puzzles))

    for (i, puzzle) in enumerate(puzzles)
        mask = get_bit_mask(puzzle)
        sub_mask = mask
        first_bit_index = Int(puzzle[1]) - Int('a')

        while true
            if (sub_mask >> first_bit_index & 1) == 1
                solution[i] += get(letter_frequencies, sub_mask, 0)
            end
            (sub_mask == 0) && break

            sub_mask = (sub_mask - 1) & mask
        end
    end

    return solution
end
## @lc code=end


text:
function erase_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.erase_mode = getproperty(han.b["erase_button"],:active,Bool)

    nothing
end

function erase_start(han::Tracker_Handles,x,y)

    plot_image(han,han.current_frame')
    r = Gtk.getgc(han.c)
    Cairo.save(r)
    ctxcopy = copy(r)

    plot_whiskers(han)

    push!((han.c.mouse, :button1motion),  (c, event) -> erase_move(han, event.x, event.y, ctxcopy))
    push!((han.c.mouse, :motion), Gtk.default_mouse_cb)
    push!((han.c.mouse, :button1release), (c, event) -> erase_stop(han, event.x, event.y, ctxcopy))

    nothing
end

function erase_move(han::Tracker_Handles, x,y,ctxcopy)

    r=Gtk.getgc(han.c)

    #check for whisker overlap and erase points that are overlapping
    keep=trues(han.wt.whiskers[han.woi_id].len)
    for i=1:han.wt.whiskers[han.woi_id].len
        if (x+5.0>han.wt.whiskers[han.woi_id].x[i])&(x-5.0<han.wt.whiskers[han.woi_id].x[i])
            if (y+5.0>han.wt.whiskers[han.woi_id].y[i])&(y-5.0<han.wt.whiskers[han.woi_id].y[i])
                keep[i]=false
            end
        end
    end

    if length(find(keep.==false))>0
        keep[1:findfirst(keep.==false)] .= false
    end

    han.wt.whiskers[han.woi_id].x=han.wt.whiskers[han.woi_id].x[keep]
    han.wt.whiskers[han.woi_id].y=han.wt.whiskers[han.woi_id].y[keep]
    han.wt.whiskers[han.woi_id].thick=han.wt.whiskers[han.woi_id].thick[keep]
    han.wt.whiskers[han.woi_id].scores=han.wt.whiskers[han.woi_id].scores[keep]

    han.wt.whiskers[han.woi_id].len=length(han.wt.whiskers[han.woi_id].x)

    #redraw whisker
    set_source(r,ctxcopy)
    paint(r)

    plot_whiskers(han)

    nothing
end

function erase_stop(han::Tracker_Handles,x,y,ctxcopy)

    assign_woi(han)

    pop!((han.c.mouse, :button1motion))
    pop!((han.c.mouse, :motion))
    pop!((han.c.mouse, :button1release))

    nothing
end


text:
#=
Author : Ng Yeesian
Modified : Guy Farmer / Michael Hilborn / Zander Tedjo / Will Worthington
Defines the DeploymentProblem struct and constructor. Defines params and naive_solution
=#
using AmbulanceDeployment
abstract type BM <: AbstractMatrix{Bool} end
abstract type IM <: AbstractMatrix{Int} end

mutable struct DeploymentProblem
        nambulances::Int
        nlocations::Int
        nregions::Int
        demand::Array{Int64,2}      # nperiod x nregion !! can potentially break !!
        train::Vector{Int}
        test::Vector{Int}
        coverage::Array{Bool,2}    # nregion x nlocation !! can potentially break !!
        adjacency::Array{Bool,2}   # nregion x nregion !! can potentially break !!
end

function DeploymentProblem(
        nregions::Int,
        hourly_calls::DataFrame,
        adjacent_nbhd::DataFrame,
        coverage::Array{Bool,2}; #!! changed from original !!
        namb = 30,
        train_filter = (hourly_calls[:year] .== 2020) .* (hourly_calls[:month] .<= 3)
    )
    #regions = Int[parse(Int,string(x)) for x in names(hourly_calls[!,6:end - 8])]
    regions = collect(1:nregions) #recall julia indexes from 1
    locations = collect(1:size(coverage,2))
    adjacent = Matrix(adjacent_nbhd[!,2:end]) #adjacent = convert(Array, adjacent_nbhd[!,2:end])[regions,regions] .> 0.5
    adjacent = adjacent[regions,regions] .> 0.5
    #adjacent = adjacent[1,:,1,:] #the previous bit operation changes the shape to 1x200x1x200

    demand = Matrix(hourly_calls[:,6:7]) #demand = convert(Array,hourly_calls[:,6:end]) #convert is deprecated

    indices = 1:nrow(hourly_calls)
    train_indices = indices[train_filter]
    test_indices = indices[.!train_filter]

    # println(namb)
    # println(length(locations))
    # println(length(regions))
    # println(size(demand))
    # println(size(train_indices))
    # println(size(test_indices))
    # println(size(coverage[:,:]))
    # println(typeof(coverage))
    # println(size(Array{Bool,2}(adjacent)))
    # println(typeof(adjacent))

    DeploymentProblem(
        namb,
        length(locations),
        length(regions),
        demand,
        train_indices,
        test_indices,
        coverage[:,:],
        Array{Bool,2}(adjacent)
    )
end
function naive_solution(p::DeploymentProblem)
    # evenly distribute the ambulances over all the locations
    x = zeros(Int, p.nlocations)
    for i in 0:p.nambulances-1
        x[1 + (i % p.nlocations)] += 1
    end
    x
end

mutable struct DispatchProblem
    emergency_calls::DataFrame
    hospitals::DataFrame
    stations::DataFrame
    coverage::Matrix{Bool} # (nbhd x stns)
    turnaround::LogNormal
    deployment::Vector{Int}

    shortfalls::Int
    wait_queue::Vector{Vector{Int}} # length nbhd
    available::Vector{Int} # length stns
    redeploy_events::Vector{Tuple{Int,Int,Int,Int}} # amb,from,to,time

    DispatchProblem(emergency_data::DataFrame,
                        hospitals::DataFrame,
                        stations::DataFrame,
                        coverage::Array{Bool,2};
                        turnaround::LogNormal = LogNormal(3.65, 0.3),
                        ) =
        new(emergency_data, hospitals, stations, coverage, turnaround)
end

function initialize!(problem::DispatchProblem, deployment::Vector{Int})
    problem.wait_queue = [Int[] for i in 1:size(problem.coverage,2)]
    problem.available = copy(deployment)
    problem.deployment = deepcopy(deployment)
    problem.redeploy_events = Tuple{Int,Int,Int,Int}[]
    problem.shortfalls = 0

    problem.emergency_calls[!,:arrival_seconds] =
        cumsum(problem.emergency_calls[!,:interarrival_seconds])

    problem
end

function DispatchProblem(
        emergency_data::DataFrame,
        hospitals::DataFrame,
        stations::DataFrame,
        coverage::Array{Bool,2},
        deployment::Vector{Int};
        turnaround::LogNormal = LogNormal(3.65, 0.3),
    )
    problem = DispatchProblem(
        emergency_data,
        hospitals,
        stations,
        coverage,
        turnaround=turnaround,
    )
    initialize!(problem, deployment)
end

function returned_to!(problem::DispatchProblem, location::Int, t::Int)
    @assert problem.available[location] >= 0
    problem.available[location] += 1
end

struct Params
    α::Float64 # Probabilistic Guarantee
    ε::Float64 # Convergence
    δ::Float64 # Solver
    nperiods::Int # for
    maxiter::Int # for RobustDeployment
end
params = Params(0.01, 0.5, 1e-6, 500, 50)

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
using Polyline

# Load gpx file and parse.
gpxFile = readGPX("route.gpx")
gpxRoute = parseGPX(gpxFile)

# Encode polyline
polyline = encodePolyline(gpxRoute)

# Decode polyline
route = decodePolyline(polyline)

# Plot the route in Goole maps. You would need to obtain a token from the static maps API
# which you can find here:
# https://developers.google.com/maps/documentation/maps-static/get-api-key
# Then either pass it as an argument or set it as the environment variable GOOGLE_MAPS_API
#url = mapsURL(polyline)

# Or set the API token as an argument.
url = mapsURL(polyline; token="...")

# With the custom defaults you can plot the route.
getMapImage(url; pathFig="/tmp/fig1.png")

# Full customization.
url = mapsURL(polyline; token="...",
              type = "roadmap", # https://developers.google.com/maps/documentation/maps-static/dev-guide#MapTypes
              size = 1000,
              scale = 2, # https://developers.google.com/maps/documentation/maps-static/dev-guide#scale_values
              MarkersStart = (gpxRoute[1, 1], gpxRoute[1, 2]),
              MarkersEnd = (gpxRoute[end, 1], gpxRoute[end, 2]))

getMapImage(url; pathFig="/tmp/fig2.png")


text:
function addRecords!(dest::FASTX.FASTA.Record, src::FASTX.FASTA.Record)
  append!(dest.data, src.data[src.sequence])
  dest.filled = UnitRange{Int}(dest.filled.start, dest.filled.stop + src.filled.stop - src.sequence.start + 1)
  dest.sequence = UnitRange{Int}(dest.sequence.start, dest.sequence.stop + src.sequence.stop - src.sequence.start + 1)
end

function concatenate(infiles::Vector{String}, outfile::String)
  concat = Dict{String, FASTX.FASTA.Record}()
  allLabels = Set(String[])
  alignLen = 0
  alphabet=""
  firstrecord=true
  for (i, file) in enumerate(infiles)
    theseLabels = String[]
    if i == 1 # first file will have only new ids
      reader = open(FASTX.FASTA.Reader, file)
      for record in reader
        if firstrecord==true
          alignLen += length(record.sequence)
          firstrecord=false
          alphabet = typeof(FASTX.FASTA.sequence(record))
        end
        id = FASTX.FASTA.identifier(record)
        push!(allLabels, id)
        concat[id] = record
      end
      close(reader)
    else # later records need to check for id existance
      record = FASTX.FASTA.Record()
      open(FASTX.FASTA.Reader, file) do reader
        while !eof(reader)
          read!(reader, record)
          id = FASTX.FASTA.identifier(record)
          push!(theseLabels, id)
          if id in allLabels
            addRecords!(concat[id], record)
          else # need to create and back-propogate gaps
            push!(allLabels, id)
            newRecord = FASTX.FASTA.Record(id,repeat("-",alignLen))
            pop!(newRecord.data)
            concat[id] = newRecord
            addRecords!(concat[id], record)
          end
        end
        alignLen += length(record.sequence)
        # assume last sequence will be representative of alphabet
        @assert alphabet == typeof(FASTX.FASTA.sequence(record))
      end
      missingLabels = setdiff(allLabels, theseLabels)
      for label in missingLabels
        gapRecord = copy(record)
        gapRecord.data[gapRecord.sequence] = repeat(UInt8[0x2d], length(gapRecord.sequence))
        addRecords!(concat[label], gapRecord)
      end
    end
  end
  open(FASTX.FASTA.Writer, outfile) do w
    for (key, value) in concat
      write(w, value)
    end
  end
  return concat
end

#using BenchmarkTools
#tmpdir = tempdir()
#infiles = simAlignments(30,100,tmpdir)
#@benchmark concatenate($simAlignments(30,100,tmpdir),$mktemp(tmpdir)[1])
#rm(tmpdir, recursive=true)

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
module MLJBase

# ===================================================================
# IMPORTS

import Base: ==, precision, getindex, setindex!
import Base.+, Base.*

# Scitype
import MLJScientificTypes.ScientificTypes: TRAIT_FUNCTION_GIVEN_NAME
using MLJScientificTypes
using MLJModelInterface

# Interface
import MLJModelInterface: fit, update, update_data, transform,
    inverse_transform, fitted_params, predict,
    predict_mode, predict_mean, predict_median, predict_joint,
    evaluate, clean!, is_same_except,
    save, restore, is_same_except, istransparent,
    params

# Macros
using Parameters

# Containers & data manipulation
using Tables
import PrettyTables
using DelimitedFiles
using OrderedCollections
using CategoricalArrays
import InvertedIndices: Not
import JLSO
import Dates

# Distributed computing
using Distributed
using ComputationalResources
using ComputationalResources: CPUProcesses
using ProgressMeter
import .Threads

# Operations & extensions
import LossFunctions
import LossFunctions: DistanceLoss, MarginLoss, SupervisedLoss
import StatsBase
import StatsBase: fit!, mode, countmap
import Missings: levels
import Distributions
import Distributions: pdf, logpdf, sampler
const Dist = Distributions

# from Standard Library:
using Statistics, LinearAlgebra, Random, InteractiveUtils

# ===================================================================
## EXPORTS

# -------------------------------------------------------------------
# re-exports from MLJModelInterface, MLJScientificTypes
# NOTE: MLJBase does **not** re-export UnivariateFinite to avoid
# ambiguities between the raw constructor (MLJBase.UnivariateFinite)
# and the general method (MLJModelInterface.UnivariateFinite)

# MLJ model hierarchy
export MLJType, Model, Supervised, Unsupervised,
    Probabilistic, JointProbabilistic, Deterministic, Interval, Static,
    ProbabilisticComposite, JointProbabilisticComposite, DeterministicComposite,
    IntervalComposite, UnsupervisedComposite, StaticComposite,
    ProbabilisticSurrogate, JointProbabilisticSurrogate, DeterministicSurrogate,
    IntervalSurrogate, UnsupervisedSurrogate, StaticSurrogate,
    Surrogate, Composite

export UnivariateFinite

# MLJType equality
export is_same_except

# model constructor + metadata
export @mlj_model, metadata_pkg, metadata_model

# model api
export fit, update, update_data, transform, inverse_transform,
       fitted_params, predict, predict_mode, predict_mean, predict_median, predict_joint,
       evaluate, clean!

# model/measure matching:
export Checker, matching

# model traits
export input_scitype, output_scitype, target_scitype,
       is_pure_julia, package_name, package_license,
       load_path, package_uuid, package_url,
       is_wrapper, supports_weights, supports_online,
       docstring, name, is_supervised,
       prediction_type, implemented_methods, hyperparameters,
       hyperparameter_types, hyperparameter_ranges

# data operations
export matrix, int, classes, decoder, table,
       nrows, selectrows, selectcols, select

# re-exports from MLJScientificTypes
export Unknown, Known, Finite, Infinite,
       OrderedFactor, Multiclass, Count, Continuous, Textual,
       Binary, ColorImage, GrayImage, Image, Table
export scitype, scitype_union, elscitype, nonmissing, trait
export coerce, coerce!, autotype, schema, info

# -------------------------------------------------------------------
# exports from this module, MLJBase

export DeterministicComposite,
    ProbabilisticComposite,
    UnsupervisedComposite, @load

# computational_resources.jl:
export default_resource

# one_dimensional_ranges.jl:
export ParamRange, NumericRange, NominalRange, iterator, scale

# parameter_inspection.jl:
export params # note this is *not* an extension of StatsBase.params

# data.jl:
export partition, unpack, complement, restrict, corestrict

# utilities.jl:
export @set_defaults, flat_values, recursive_setproperty!,
       recursive_getproperty, pretty, unwind

# show.jl
export HANDLE_GIVEN_ID, @more, @constant, @bind, color_on, color_off

# univariate_finite/
export average, UnivariateFiniteArray, UnivariateFiniteVector

# info_dict.jl:
export info_dict

# datasets.jl:
export load_boston, load_ames, load_iris,
       load_reduced_ames, load_crabs, load_smarket,
       @load_boston, @load_ames, @load_iris,
       @load_reduced_ames, @load_crabs, @load_smarket

# sources.jl:
export source, Source, CallableReturning

# machines.jl:
export machine, Machine, fit!, report, fit_only!

# datasets_synthetics.jl
export make_blobs, make_moons, make_circles, make_regression

# composition:
export machines, sources, anonymize!, @from_network, @pipeline,
    glb, @tuple, node, @node, sources, origins, return!,
    nrows_at_source, machine,
    rebind!, nodes, freeze!, thaw!, models, Node, AbstractNode,
    DeterministicSurrogate, ProbabilisticSurrogate, UnsupervisedSurrogate,
    DeterministicComposite, ProbabilisticComposite, UnsupervisedComposite

# aliases to the above,  kept for backwards compatibility:
export  DeterministicNetwork, ProbabilisticNetwork, UnsupervisedNetwork

# resampling.jl:
export ResamplingStrategy, Holdout, CV, StratifiedCV,
       evaluate!, Resampler, PerformanceEvaluation

# openml.jl:
export OpenML

# -------------------------------------------------------------------
# exports from MLJBase specific to Measure (these may go in their
# specific MLJMeasureInterface package in some future)

# measures/registry.jl:
export measures, metadata_measure

# measure/measures.jl:
export orientation, reports_each_observation,
    is_feature_dependent, aggregation,
    aggregate, default_measure, value,
    spports_weights, prediction_type

# measures/continuous.jl:
export mav, mae, mape, rms, rmsl, rmslp1, rmsp, l1, l2, log_cosh

# measures/confusion_matrix.jl:
export confusion_matrix, confmat

# measures/finite.jl
export cross_entropy, BrierScore, brier_score,
       misclassification_rate, mcr, accuracy,
       balanced_accuracy, bacc, bac,
       matthews_correlation, mcc

# measures/finite.jl -- binary order independent:
export auc, area_under_curve, roc_curve, roc

# measures/finite.jl -- binary order dependent:
export TruePositive, TrueNegative, FalsePositive, FalseNegative,
       TruePositiveRate, TrueNegativeRate, FalsePositiveRate,
       FalseNegativeRate, FalseDiscoveryRate, Precision, NPV, FScore,
       # standard synonyms
       TPR, TNR, FPR, FNR, FDR, PPV,
       Recall, Specificity, BACC,
       # instances and their synonyms
       truepositive, truenegative, falsepositive, falsenegative,
       true_positive, true_negative, false_positive, false_negative,
       truepositive_rate, truenegative_rate, falsepositive_rate,
       true_positive_rate, true_negative_rate, false_positive_rate,
       falsenegative_rate, negativepredictive_value,
       false_negative_rate, negative_predictive_value,
       positivepredictive_value, positive_predictive_value,
       tpr, tnr, fpr, fnr,
       falsediscovery_rate, false_discovery_rate, fdr, npv, ppv,
       recall, sensitivity, hit_rate, miss_rate,
       specificity, selectivity, f1score, fallout

# -------------------------------------------------------------------
# re-export from Random, StatsBase, Statistics, Distributions,
# CategoricalArrays, InvertedIndices:
export pdf, sampler, mode, median, mean, shuffle!, categorical, shuffle,
       levels, levels!, std, Not, support, logpdf

# ===================================================================
## CONSTANTS

# the directory containing this file: (.../src/)
const MODULE_DIR = dirname(@__FILE__)

# horizontal space for field names in `MLJType` object display:
const COLUMN_WIDTH = 24
# how deep to display fields of `MLJType` objects:
const DEFAULT_SHOW_DEPTH = 0
const DEFAULT_AS_CONSTRUCTED_SHOW_DEPTH = 2
const INDENT = 4

const Arr = AbstractArray
const Vec = AbstractVector

const MMI = MLJModelInterface
const FI  = MLJModelInterface.FullInterface
# ===================================================================
# Computational Resource
# default_resource allows to switch the mode of parallelization

default_resource()    = DEFAULT_RESOURCE[]
default_resource(res) = (DEFAULT_RESOURCE[] = res)

# stub for @load (extended by MLJModels)
macro load end

# ===================================================================
# Includes

include("init.jl")
include("utilities.jl")
include("show.jl")
include("info_dict.jl")
include("interface/data_utils.jl")
include("interface/model_api.jl")

include("univariate_finite/types.jl")
include("univariate_finite/methods.jl")
include("univariate_finite/arrays.jl")

include("sources.jl")
include("machines.jl")

include("composition/abstract_types.jl")
include("composition/learning_networks/nodes.jl")
include("composition/learning_networks/inspection.jl")
include("composition/learning_networks/machines.jl")
@static if VERSION ≥ v"1.3.0-"
    include("composition/learning_networks/arrows.jl")
end

include("composition/models/methods.jl")
include("composition/models/from_network.jl")
include("composition/models/inspection.jl")
include("composition/models/pipelines.jl")
include("composition/models/_wrapped_function.jl")

include("operations.jl")
include("resampling.jl")

include("hyperparam/one_dimensional_ranges.jl")
include("hyperparam/one_dimensional_range_methods.jl")

include("data/data.jl")
include("data/datasets.jl")
include("data/datasets_synthetic.jl")

include("matching.jl")
include("measures/measures.jl")
include("measures/measure_search.jl")

include("openml.jl")

end # module


text:
struct ClebschGordan{A,B,C,j₁T,j₂T}
    j₁::A
    m₁::A
    j₂::B
    m₂::B
    J::C
    M::C
end

Base.zero(::Type{ClebschGordan{A,B,C,j₁T,j₂T}}) where {A,B,C,j₁T,j₂T} =
    ClebschGordan{A,B,C,j₁T,j₂T}(zero(A),zero(A),zero(B),zero(B),zero(C),zero(C))

const ClebschGordanℓs{I<:Integer} = ClebschGordan{I,Half{I},Half{I},:ℓ,:s}
ClebschGordanℓs(j₁::I, m₁::I, j₂::R, m₂::R, J::R, M::R) where {I<:Integer,R<:Half{I}} =
    ClebschGordanℓs{I}(j₁, m₁, j₂, m₂, J, M)

Base.convert(::Type{T}, cg::ClebschGordan) where {T<:Real} =
    clebschgordan(cg.j₁,cg.m₁,cg.j₂,cg.m₂,cg.J,cg.M)

# We have chosen the Clebsch–Gordan coefficients to be real.
Base.adjoint(cg::ClebschGordan) = cg

Base.show(io::IO, cg::ClebschGordan) =
    write(io, "⟨$(cg.j₁),$(cg.j₂);$(cg.m₁),$(cg.m₂)|$(cg.J),$(cg.M)⟩")

function Base.show(io::IO, cg::ClebschGordanℓs)
    spin = cg.m₂ > 0 ? "↑" : "↓"
    write(io, "⟨$(spectroscopic_label(cg.j₁));$(cg.m₁),$(spin)|$(cg.J),$(cg.M)⟩")
end

#=
Relevant equations taken from

- Sakurai, J. J., & Napolitano, J. J. (2010). Modern quantum mechanics
(2nd edition). : Addison-Wesley.
=#

#=

Repeated action [Eq. (3.8.45)] with the ladder operators on the jmⱼ
basis

J±|ℓs;jmⱼ⟩

gives the following recursion relations, Eq. (3.8.49):

√((j∓mⱼ)(j±mⱼ+1))⟨ℓs;mₗ,mₛ|ℓs;j,mⱼ±1⟩
= √((ℓ∓mₗ+1)(ℓ±mₗ))⟨ℓs;mₗ∓1,mₛ|ℓs;j,mⱼ⟩
+ √((s∓mₛ+1)(s±mₛ))⟨ℓs;mₗ,mₛ∓1|ℓs;j,mⱼ⟩.

We know [Eq. (3.8.58)] that

⟨ℓs;±ℓ,±1/2|ℓs;ℓ±1/2,ℓ±1/2⟩ = 1.

=#
function rotate!(blocks::Vector{M}, orbs::RelativisticOrbital...) where {T,M<:AbstractMatrix{T}}
    2length(blocks) == sum(degeneracy.(orbs))-2 ||
        throw(ArgumentError("Invalid block size for $(orbs...)"))
    ℓ = κ2ℓ(first(orbs).κ)
    # We sort by mⱼ and remove the first and last elements since they
    # are pure and trivially unity.
    ℓms = sort(vcat([[(ℓ,m,s) for m ∈ -ℓ:ℓ] for s = -half(1):half(1)]...)[2:end-1], by=((ℓms)) -> (ℓms[2],+(ℓms[2:3]...)))
    jmⱼ = sort(vcat([[(j,mⱼ) for mⱼ ∈ -j:j] for j ∈ [κ2j(o.κ) for o in orbs]]...), by=last)[2:end-1]
    for (a,(ℓ,m,s)) in enumerate(ℓms)
        bi = cld(a, 2)
        o = 2*(bi-1)
        for (b,(j,mⱼ)) in enumerate(jmⱼ)
            b-o ∉ 1:2 && continue
            blocks[bi][a-o,b-o] = convert(T, ClebschGordanℓs(ℓ,m,half(1),s,j,mⱼ))
        end
    end
    blocks
end

"""
    jj2ℓsj([T=Float64, ]orbs...)

Generates the block-diagonal matrix that transforms jj-coupled
configurations to lsj-coupled ones.

The blocks correspond to invariant subspaces, which rotate among
themselves to form new linear combinations.

They are sorted according to n,ℓ and within the blocks, the columns
are sorted according to mⱼ,j (increasing) and the rows according to
s,ℓ,m (increasing).

E.g. the p-block will have the following structure:

```
 ⟨p;-1,↓|3/2,-3/2⟩  │       ⋅                  ⋅             │       ⋅                ⋅           │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │  ⟨p;0,↓|1/2,-1/2⟩   ⟨p;0,↓|3/2,-1/2⟩   │       ⋅                ⋅           │       ⋅
      ⋅             │  ⟨p;-1,↑|1/2,-1/2⟩  ⟨p;-1,↑|3/2,-1/2⟩  │       ⋅                ⋅           │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │       ⋅                  ⋅             │  ⟨p;1,↓|1/2,1/2⟩  ⟨p;1,↓|3/2,1/2⟩  │       ⋅
      ⋅             │       ⋅                  ⋅             │  ⟨p;0,↑|1/2,1/2⟩  ⟨p;0,↑|3/2,1/2⟩  │       ⋅
 ───────────────────┼────────────────────────────────────────┼────────────────────────────────────┼─────────────────
      ⋅             │       ⋅                  ⋅             │       ⋅                ⋅           │  ⟨p;1,↑|3/2,3/2⟩
```

"""
function jj2ℓsj(::Type{T}, orbs::RelativisticOrbital...) where T
    nℓs = map(o -> (o.n, κ2ℓ(o.κ)), sort([orbs...]))
    blocks = map(unique(nℓs)) do (n,ℓ)
        i = findall(isequal((n,ℓ)), nℓs)
        subspace = orbs[i]
        mⱼ = map(subspace) do orb
            j = convert(Rational, κ2j(orb.κ))
            -j:j
        end

        jₘₐₓ = maximum([κ2j(o.κ) for o in subspace])
        pure = [Matrix{T}(undef,1,1),Matrix{T}(undef,1,1)]
        pure[1][1] = convert(T, ClebschGordanℓs(ℓ,-ℓ,half(1),-half(1),jₘₐₓ,-jₘₐₓ))
        pure[2][1] = convert(T, ClebschGordanℓs(ℓ,ℓ,half(1),half(1),jₘₐₓ,jₘₐₓ))

        n = sum(length.(mⱼ))-2
        if n > 0
            nblocks = div(n,2)
            b = [zeros(T,2,2) for i = 1:nblocks]
            rotate!(b, subspace...)
            [pure[1],b...,pure[2]]
        else
            pure
        end
    end |> b -> vcat(b...)
    rows = size.(blocks,1)
    N = sum(rows)
    R = BlockBandedMatrix(Zeros{T}(N,N), rows,rows, (0,0))
    for (i,b) in enumerate(blocks)
        R[Block(i,i)] .= b
    end
    R
end
jj2ℓsj(orbs::RelativisticOrbital...) = jj2ℓsj(Float64, orbs...)

angular_integral(a::SpinOrbital{<:Orbital}, b::SpinOrbital{<:Orbital}) =
    a.orb.ℓ == b.orb.ℓ && a.m == b.m

angular_integral(a::SpinOrbital{<:RelativisticOrbital}, b::SpinOrbital{<:RelativisticOrbital}) =
    a.orb.ℓ == b.orb.ℓ && a.m == b.m

angular_integral(a::SpinOrbital{<:Orbital}, b::SpinOrbital{<:RelativisticOrbital}) =
    Int(a.orb.ℓ == b.orb.ℓ) * clebschgordan(a.orb.ℓ, a.m[1], half(1), a.m[2], b.orb.j, b.m[1])

angular_integral(a::SpinOrbital{<:RelativisticOrbital}, b::SpinOrbital{<:Orbital}) =
    conj(angular_integral(b, a))

function jj2ℓsj(sorb::SpinOrbital{<:Orbital})
    orb,(mℓ,ms) = sorb.orb,sorb.m
    @unpack n,ℓ = orb
    mj = mℓ + ms
    map(max(abs(mj),ℓ-half(1)):(ℓ+half(1))) do j
        ro = SpinOrbital(RelativisticOrbital(n,ℓ,j),mj)
        ro => angular_integral(ro, sorb)
    end
end

function jj2ℓsj(sorb::SpinOrbital{<:RelativisticOrbital})
    orb,(mj,) = sorb.orb,sorb.m
    @unpack n,ℓ,j = orb
    map(max(-ℓ, mj-half(1)):min(ℓ, mj+half(1))) do mℓ
        ms = mj - mℓ
        o = SpinOrbital(Orbital(n,ℓ),(Int(mℓ),ms))
        o => angular_integral(o, sorb)
    end
end

function jj2ℓsj(sorbs::AbstractVector{<:SpinOrbital{<:RelativisticOrbital}})
    @assert issorted(sorbs)
    # For each jj spin-orbital, find all corresponding ls
    # spin-orbitals with their CG coeffs.
    couplings = map(jj2ℓsj, sorbs)

    CGT = typeof(couplings[1][1][2])
    LSOT = eltype(reduce(vcat, map(c -> map(first, c), couplings)))
    blocks = Vector{Tuple{Int,Int,Matrix{CGT}}}()
    tmp_blocks = Vector{Matrix{CGT}}(undef, length(sorbs))

    # Sort the ls spin-orbitals such that appear in the same place as
    # those jj orbitals with the same mⱼ; this is essential to allow
    # applying the JJ ↔ LS transform in-place via orbital
    # rotations. Also generate the rotation matrices which are either
    # 1×1 (for pure states) or 2×2.
    orb_map = Dict{LSOT,Int}()
    ls_orbitals = Vector{LSOT}(undef, length(sorbs))
    for (i,cs) in enumerate(couplings)
        if length(cs) == 1
            (o,coeff) = first(cs)
            ls_orbitals[i] = o
            M = zeros(CGT, 1, 1)
            M[1] = coeff
            push!(blocks, (i,i,M))
        else
            (a,ca),(b,cb) = cs
            mi,ma = minmax(a,b)
            j = get!(orb_map, mi, i)
            ls_orbitals[i] = if i == j # First time we see this block
                tmp_blocks[i] = zeros(CGT, 2, 2)
                mi
            else # Block finished
                push!(blocks, (j, i, tmp_blocks[j]))
                ma
            end
            tmp_blocks[j][(i == j ? 1 : 2),:] .= a < b ? (ca,cb) : (cb,ca)
        end
    end

    ls_orbitals, blocks
end

export jj2ℓsj, ClebschGordan, ClebschGordanℓs


text:
include(joinpath("dataProcessing", "anomalies.jl"))
include(joinpath("dataProcessing", "causeSelection.jl"))
include(joinpath("dataProcessing", "dataSelection.jl"))
include(joinpath("dataProcessing", "processOverflow.jl"))
include(joinpath("dataProcessing", "siteSelection.jl"))
include(joinpath("dataProcessing", "RainYday.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
# Test of transmission and distribution decoupling

# Note 1
# Only distribution part is implemented at the moment: given a distribution network having some
# candidate storage / flexible loads, a set of flexibility candidates is returned.
# The entire distribution network is seen as a flexibility candidate by transmission. Multiple
# candidates related to the same distribution network represent alternative planning options, which
# vary in amount of flexibility provided and in cost.

# Note 2
# Among flexibility sources in distribution, only storage is considered at the moment.

## Import packages and choose a solver

using Memento
_LOGGER = Logger(basename(@__FILE__)[1:end-3]) # A logger for this script, also used by included files.

import FlexPlan; const _FP = FlexPlan
include("../io/load_case.jl")
include("../io/td_decoupling.jl")
import Cbc
optimizer = _FP.optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0)

## Input parameters

# Test case parameters
flex_load         = false # Toggles flexibility of loads
ne_storage        =  true # Toggles candidate storage
scale_gen         =   1.0 # Scaling factor of all generators, wind included
scale_wind        =   6.0 # Further scaling factor of wind generator
scale_load        =   1.0 # Scaling factor of loads
energy_cost       =  50.0 # Cost of energy exchanged with transmission network [€/MWh]
year_scale_factor =    10 # How many years a representative year should represent
number_of_hours   =    24 # Number of hourly optimization periods
start_period      =     1 # First period of profile data to use

# Procedure parameters
number_of_candidates = 4 # Number of flexibility candidates for each distribution network to be returned

# Script parameters
out_dir = "./test/data/output_files/td_decoupling/" # Directory of output files
plot = true # Toggles plotting of results

## Load distribution network instance

d_mn_data = load_cigre_mv_eu(; flex_load, ne_storage, scale_gen, scale_wind, scale_load, energy_cost, year_scale_factor, number_of_hours, start_period)

## Solve problem

dist_candidates = _FP.solve_td_decoupling_distribution(d_mn_data; optimizer, number_of_candidates)

## Result analysis and output

# Candidates to be processed and plotted. Passing a subset of candidates can be useful in making plots less crowded. Order of candidates is respected in plots.
candidate_ids = sort(collect(keys(dist_candidates)))

mkpath(out_dir)

# Default Plots backend can be changed here
#import Plots
#Plots.plotlyjs()

# Kwargs: `plot_ext` can be used to set plot file extension; also all Plots kwargs are accepted. Example: `plot_ext="png", dpi=300`
report_dist_candidates_pcc_power(dist_candidates, out_dir; plot, candidate_ids)
report_dist_candidates_branch(dist_candidates, out_dir, d_mn_data; plot, candidate_ids)
report_dist_candidates_storage(dist_candidates, out_dir; plot, candidate_ids)

report_dist_candidates_investment(dist_candidates, out_dir; candidate_ids)
report_dist_candidates_cost(dist_candidates, out_dir; candidate_ids)
#report_dist_candidates_nw_summary(dist_candidates, out_dir)

println("Test completed. Results saved in $out_dir")


text:
using Test
using DelimitedFiles
using SeuratRDS
using Conda

dn = joinpath(@__DIR__,"data")
testfn = joinpath(dn,"testSeur.rds")
checkfn = joinpath(dn,"dataSeur.csv")
@testset "SeuratRDS.jl" begin
    modality = "RNA"
    assay = "data"
    metadata = "nCount_RNA"

    env = initR() # make the env

    @test Conda.channels(env) == ["r","defaults"]

    dat = loadSeur(testfn,env,modality,assay,metadata)
    check,ccols = readdlm(checkfn,header=true)

    @test check == dat.dat

    closeR(env)
    @test !isdir(env)
end


text:
include("../MDMKP/MDMKP.jl")

function MDMKP_id_dataframe()
	tight_problems = MDMKP.load_folder()
	loose_problems = MDMKP.loosen.(tight_problems)
	all_problems = vcat(tight_problems..., loose_problems...)
	labeled_problems = MDMKP.set_decision.(all_problems)

	ids = [p.id for p in labeled_problems]

	df = DataFrame(
		id = Int[],
		dataset = Int[],
		instance = Int[],
		case = Int[],
		tightness = Number[],
		n_vars = Int[],
		n_demands = Int[],
		n_dimensions = Int[],
		mixed_obj = Bool[],
		loosened = Bool[],
		category = String[]
	)

	get_row(id) = [id.id, id.dataset, id.instance, id.case, id.tightness,
		id.n_vars, id.n_demands, id.n_dimensions, id.mixed_obj, id.loosened,
		"$(id.category)"]

	[push!(df, get_row(id)) for id in ids]
	df
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
# Functions to replace BLAS kernels in gmres.jl and orthogonalize.jl for
# counting flops with GFlops.jl

import LinearAlgebra: dot, mul!, norm
import Base: *

# gmres228, orthogonalize55, orthogonalize70
function dot(x::StridedVector{Float64}, y::StridedVector{Float64})
	z = zero(Float64)
	for i=1:length(x)
		z += x[i] * y[i]
	end
	z
end

# orthogonalize14, orthogonalize42
function mul!(y::StridedVector{Float64}, A::Adjoint{<:Any, <:StridedMatrix}, x::StridedVector{Float64})
	for i in 1:length(y)
		tmp = zero(Float64)
		for j in 1:length(x)
			tmp += A[i,j] * x[j]
		end
		y[i] = tmp
	end
end

# orthogonalize26
function *(A::Adjoint{<:Any, <:StridedMatrix{Float64}}, x::StridedVector{Float64})
	y = zeros(size(A)[1])
	for i in 1:length(y)
		tmp = zero(Float64)
		for j in 1:length(x)
			tmp += A[i,j] * x[j]
		end
		y[i] = tmp
	end
	y
end

# gmres267, gmres272, orthogonalize15, orthogonalize29, orthogonalize43
function mul!(y::StridedVector{Float64}, A::StridedMatrix{Float64}, x::StridedVector{Float64}, alpha::Float64, beta::Float64)
	for i in 1:length(y)
		tmp = zero(Float64)
		for j in 1:length(x)
			tmp += alpha * A[i,j] * x[j]
		end
		y[i] = tmp + beta * y[i]
	end
end

# gmres252, orthogonalize17, orthogonalize22, orthogonalize28, orthogonalize32,
# orthogonalize45, orthogonalize61, orthogonalize75
function norm(x::StridedVector{Float64})
	result = 0.0
	for i = 1:length(x)
		result += x[i] * x[i]
	end
	sqrt(result)
end

max_stars_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_issues_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_forks_repo_name: jkrch/JuliaParGMRES-BenchCLI


text: println("Hello, Wold!")

max_stars_repo_name: Susensio/programming-languages

max_issues_repo_name: Susensio/programming-languages

max_forks_repo_name: Susensio/programming-languages


text:
module Debugger

using Markdown
using Base.Meta: isexpr
using REPL
using REPL.LineEdit
using REPL.REPLCompletions

using JuliaInterpreter: JuliaInterpreter, JuliaStackFrame, @lookup, Compiled, JuliaProgramCounter, JuliaFrameCode,
      finish!, enter_call_expr, step_expr!

# TODO: Work on better API in JuliaInterpreter and rewrite Debugger.jl to use it
# These are undocumented functions from JuliaInterpreter.jl used by Debugger.jl`
using JuliaInterpreter: _make_stack, pc_expr,isassign, getlhs, do_assignment!, maybe_next_call!, is_call, _step_expr!, next_call!,  moduleof,
                        iswrappercall, next_line!, linenumber

const SEARCH_PATH = []
function __init__()
    append!(SEARCH_PATH,[joinpath(Sys.BINDIR,"../share/julia/base/"),
            joinpath(Sys.BINDIR,"../include/")])
    return nothing
end

export @enter

include("LineNumbers.jl")
using .LineNumbers: SourceFile, compute_line

include("operations.jl")
include("repl.jl")
include("printing.jl")
include("commands.jl")

macro enter(arg)
    quote
        let stack = $(_make_stack(__module__,arg))
            RunDebugger(stack)
        end
    end
end

end # module


text:
# TODO: Add HigherPrecision.jl tests here?
# TODO: Test Interior Point algorithm as well
@testset "Arbitrary Precision" begin
    prob = MVP.UnconstrainedProblems.examples["Rosenbrock"]

    f = MVP.objective(prob)
    g! = MVP.gradient(prob)
    h! = MVP.hessian(prob)

    @testset "BigFloat" begin
        x0 = big.(prob.initial_x)
        res = optimize(f, x0)
        debug_printing && @show res
        @test Optim.converged(res) == true
        @test Optim.minimum(res) < 1e-8
        @test Optim.minimizer(res) ≈ [1.0, 1.0] atol=1e-4 rtol=0

        res = optimize(f, g!, x0)
        debug_printing && @show res
        @test Optim.converged(res) == true
        @test Optim.minimum(res) < 1e-16
        @test Optim.minimizer(res) ≈ [1.0, 1.0] atol=1e-10 rtol=0

        res = optimize(f, g!, h!, x0)
        debug_printing && @show res
        @test Optim.converged(res) == true
        @test Optim.minimum(res) < 1e-16
        @test Optim.minimizer(res) ≈ [1.0, 1.0] atol=1e-10 rtol=0

        lower = big.([-Inf, -Inf])
        upper = big.([0.5, 1.5])
        res = optimize(f, g!, lower, upper, x0)
        debug_printing && @show res
        @test Optim.converged(res) == true
        @test Optim.minimum(res) ≈ 0.25 atol=1e-10 rtol=0
        @test Optim.minimizer(res) ≈ [0.5, 0.25] atol=1e-10 rtol=0

        res = optimize(f, lower, upper, x0)
        debug_printing && @show res
        @test Optim.converged(res) == true
        @test Optim.minimum(res) ≈ 0.25 atol=1e-10 rtol=0
        @test Optim.minimizer(res) ≈ [0.5, 0.25] atol=1e-10 rtol=0
    end
end