text:
# Compute ϕm(A)b using the Krylov subspace Km{A,b}

############################
# Cache for expv
mutable struct ExpvCache{T}
    mem::Vector{T}
    ExpvCache{T}(maxiter::Int) where {T} = new{T}(Vector{T}(undef, maxiter^2))
end
function Base.resize!(C::ExpvCache{T}, maxiter::Int) where {T}
    C.mem = Vector{T}(undef, maxiter^2 * 2)
    return C
end
function get_cache(C::ExpvCache, m::Int)
    m^2 > length(C.mem) && resize!(C, m) # resize the cache if needed
    reshape(@view(C.mem[1:m^2]), m, m)
end

############################
# Expv
"""
    expv(t,A,b; kwargs) -> exp(tA)b

Compute the matrix-exponential-vector product using Krylov.

A Krylov subspace is constructed using `arnoldi` and `exp!` is called
on the Hessenberg matrix. Consult `arnoldi` for the values of the
keyword arguments. An alternative algorithm, where an error estimate
generated on-the-fly is used to terminate the Krylov iteration, can be
employed by setting the kwarg `mode=:error_estimate`.

    expv(t,Ks; cache) -> exp(tA)b

Compute the expv product using a pre-constructed Krylov subspace.
"""
function expv(t, A, b; mode=:happy_breakdown, kwargs...)
    # Master dispatch function
    if mode == :happy_breakdown
        _expv_hb(t, A, b; kwargs...)
    elseif mode == :error_estimate
        _expv_ee(t, A, b; kwargs...)
    else
        throw(ArgumentError("Unknown Krylov iteration termination mode, $(mode)"))
    end
end
function _expv_hb(t::Tt, A, b; cache=nothing, kwargs_arnoldi...) where Tt
    # Happy-breakdown mode: first construct Krylov subspace then expv!
    Ks = arnoldi(A, b; kwargs_arnoldi...)
    w = similar(b, promote_type(Tt, eltype(A), eltype(b)))
    expv!(w, t, Ks; cache=cache)
end
function _expv_ee(t::Tt, A, b; m=min(30, size(A, 1)), tol=1e-7, rtol=√(tol),
    ishermitian::Bool=LinearAlgebra.ishermitian(A)) where Tt
    # Error-estimate mode: construction of Krylov subspace and expv! at the same time
    n = size(A,1)
    T = promote_type(typeof(t), eltype(A), eltype(b))
    U = ishermitian ? real(T) : T
    Ks = KrylovSubspace{T,U}(n, m)
    w = similar(b, promote_type(Tt, eltype(A), eltype(b)))
    expv!(w, t, A, b, Ks, get_subspace_cache(Ks); atol=tol, rtol=rtol)
end
function expv(t::Tt, Ks::KrylovSubspace{T, U}; kwargs...) where {Tt, T, U}
    n = size(getV(Ks), 1)
    w = Vector{promote_type(Tt, T)}(undef, n)
    expv!(w, t, Ks; kwargs...)
end
"""
    expv!(w,t,Ks[;cache]) -> w

Non-allocating version of `expv` that uses precomputed Krylov subspace `Ks`.
"""
function expv!(w::AbstractVector{Tw}, t::Real, Ks::KrylovSubspace{T, U};
               cache=nothing) where {Tw, T, U}
    m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
    @assert length(w) == size(V, 1) "Dimension mismatch"
    if cache == nothing
        cache = Matrix{U}(undef, m, m)
    elseif isa(cache, ExpvCache)
        cache = get_cache(cache, m)
    else
        throw(ArgumentError("Cache must be an ExpvCache"))
    end
    copyto!(cache, @view(H[1:m, :]))
    if ishermitian(cache)
        # Optimize the case for symtridiagonal H
        F = eigen!(SymTridiagonal(cache))
        expHe = F.vectors * (exp.(lmul!(t,F.values)) .* @view(F.vectors[1, :]))
    else
        lmul!(t, cache); expH = cache
        _exp!(expH)
        expHe = @view(expH[:, 1])
    end
    lmul!(beta, mul!(w, @view(V[:, 1:m]), expHe)) # exp(A) ≈ norm(b) * V * exp(H)e
end

# NOTE: Tw can be Float64, while t is ComplexF64 and T is Float32
#       or Tw can be Float64, while t is ComplexF32 and T is Float64
#       thus they can not share the same TypeVar.
function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspace{T, U};
                cache=nothing) where {Tw, Tt, T, U}
    m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
    @assert length(w) == size(V, 1) "Dimension mismatch"
    if cache == nothing
        cache = Matrix{U}(undef, m, m)
    elseif isa(cache, ExpvCache)
        cache = get_cache(cache, m)
    else
        throw(ArgumentError("Cache must be an ExpvCache"))
    end
    if iszero(Ks.beta)
        w .= false
        return w
    end
    copyto!(cache, @view(H[1:m, :]))
    if ishermitian(cache)
        # Optimize the case for symtridiagonal H
        F = eigen!(SymTridiagonal(real(cache)))
        expHe = F.vectors * (exp.(t * F.values) .* @view(F.vectors[1, :]))
    else
        expH = _exp!(t * cache)
        expHe = @view(expH[:, 1])
    end
    lmul!(beta, mul!(w, @view(V[:, 1:m]), expHe)) # exp(A) ≈ norm(b) * V * exp(H)e
end

############################
# Cache for phiv
mutable struct PhivCache{T}
    mem::Vector{T}
    function PhivCache{T}(maxiter::Int, p::Int) where {T}
        numelems = maxiter + maxiter^2 + (maxiter + p)^2 + maxiter*(p + 1)
        new{T}(Vector{T}(undef, numelems))
    end
end
function Base.resize!(C::PhivCache{T}, maxiter::Int, p::Int) where {T}
    numelems = maxiter + maxiter^2 + (maxiter + p)^2 + maxiter*(p + 1)
    C.mem = Vector{T}(undef, numelems * 2)
    return C
end
function get_caches(C::PhivCache, m::Int, p::Int)
    numelems = m + m^2 + (m + p)^2 + m*(p + 1)
    numelems^2 > length(C.mem) && resize!(C, m, p) # resize the cache if needed
    e = @view(C.mem[1:m]); offset = m
    Hcopy = reshape(@view(C.mem[offset + 1:offset + m^2]), m, m); offset += m^2
    C1 = reshape(@view(C.mem[offset + 1:offset + (m+p)^2]), m+p, m+p); offset += (m+p)^2
    C2 = reshape(@view(C.mem[offset + 1:offset + m*(p+1)]), m, p+1)
    return e, Hcopy, C1, C2
end

############################
# Phiv
"""
    phiv(t,A,b,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using Krylov. `k` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_{k+1}(z) = \\frac{\\varphi_k(z) - 1}{z}
```

A Krylov subspace is constructed using `arnoldi` and `phiv_dense` is called
on the Hessenberg matrix. If `correct=true`, then phi_0 through phi_k-1 are
updated using the last Arnoldi vector v_m+1 [^1]. If `errest=true` then an
additional error estimate for the second-to-last phi is also returned. For
the additional keyword arguments, consult `arnoldi`.

  phiv(t,Ks,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating
the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
Formula (10).
"""
function phiv(t, A, b, k; cache=nothing, correct=false, errest=false, kwargs_arnoldi...)
    Ks = arnoldi(A, b; kwargs_arnoldi...)
    w = Matrix{eltype(b)}(undef, length(b), k+1)
    phiv!(w, t, Ks, k; cache=cache, correct=correct, errest=errest)
end
function phiv(t, Ks::KrylovSubspace{T, U}, k; kwargs...) where {T, U}
    n = size(getV(Ks), 1)
    w = Matrix{T}(undef, n, k+1)
    phiv!(w, t, Ks, k; kwargs...)
end
"""
    phiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]

Non-allocating version of 'phiv' that uses precomputed Krylov subspace `Ks`.
"""
function phiv!(w::AbstractMatrix{T}, t::Number, Ks::KrylovSubspace{T, U}, k::Integer;
               cache=nothing, correct=false, errest=false) where {T <: Number, U <: Number}
    m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
    @assert size(w, 1) == size(V, 1) "Dimension mismatch"
    @assert size(w, 2) == k + 1 "Dimension mismatch"
    if cache == nothing
        cache = PhivCache{T}(m, k)
    elseif !isa(cache, PhivCache)
        throw(ArgumentError("Cache must be a PhivCache"))
    end
    e, Hcopy, C1, C2 = get_caches(cache, m, k)
    lmul!(t, copyto!(Hcopy, @view(H[1:m, :])))
    fill!(e, zero(T)); e[1] = one(T) # e is the [1,0,...,0] basis vector
    phiv_dense!(C2, Hcopy, e, k; cache=C1) # C2 = [ϕ0(H)e ϕ1(H)e ... ϕk(H)e]
    lmul!(beta, mul!(w, @view(V[:, 1:m]), C2)) # f(A) ≈ norm(b) * V * f(H)e
    if correct
        # Use the last Arnoldi vector for correction with little additional cost
        # correct_p = beta * h_{m+1,m} * (em^T phi_p+1(H) e1) * v_m+1
        betah = beta * H[end,end] * t
        vlast = @view(V[:,end])
        @inbounds for i = 1:k
            axpy!(betah * C2[end, i+1], vlast, @view(w[:, i]))
        end
    end
    if errest
        err = abs(beta * H[end, end] * t * C2[end, end])
        return w, err
    else
        return w
    end
end


text:
module RandomStrainDistributions
using Reexport

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

include("ShearFunctions.jl")
@reexport using .ShearFunctions

include("RandomDislocations.jl")

include("DisorderConfigurations.jl")
@reexport using .DisorderConfigurations

# ===================================================================== #
#  RandomStrainDistributions module code is below
# ===================================================================== #

using Distributions
export RandomStrainDistribution, RandomDislocationDistribution, collect_dislocations

"""
Interface type for the random strain distributions
"""
abstract type RandomStrainDistribution end

"""
    collect_dislocations( rsd::RandomStrainDistribution )

Infterface function to generate dislocations from a `RandomStrainDistribution`.
"""
collect_dislocations( rsd::RandomStrainDistribution ) = error("No implementation for $(typeof(rsd)) has been defined.")

"""
    system_size( rsd::RandomStrainDistribution )

Interface function which returns the system size of a `RandomStrainDistribution`.
"""
system_size( rsd::RandomStrainDistribution ) = error("No implementation for $(typeof(rsd)) has been defined.")

# function calculate_strains( eval_r::Vector2D, rsd::RandomStrainDistribution )
#     all_dislocations = collect_dislocations( rsd )
#     axes = system_size( rsd )
#     num_dis = length( all_dislocations[:, Int(BurgersVector)] )
#     b1g = zero(eltype(eval_r.vec))
#     b2g = zero(eltype(eval_r.vec))
#     temp_eval_r = Vector2D( eval_r )
#     for dis_idx ∈ 1:num_dis
#         Vector2D!( temp_eval_r, eval_r )
#         new_strains = bxg_shears!( temp_eval_r, all_dislocations[dis_idx, Int(BurgersVector)];
#                                    source_r = all_dislocations[dis_idx, Int(DislocationOrigin) ],
#                                    diff = (A, B) -> subtract_PBC!( A, B; Lx = axes[1], Ly = axes[2] ) )
#         b1g += new_strains[1]
#         b2g += new_strains[2]
#     end
#     return ( b1g, b2g )
# end

"""
Struct containing all information required for random strains generated by edge dislocations.

# Struct Members
* `truncated::Bool`: whether one should check that the number of dislocations lie on the interval ``[1, system_size]``.
* `concentration::AbstractFloat`: the concentration of dislocations in a 2D square lattice
* `vector_diff::Function`: the difference function used to for `Vector2D` types depending on boundary conditions
* `rand_num_dis::Dis`: the `Dis <: Distribution` that returns a random number of dislocations per disorder configuration
* `burgers_vector_dist::RBVD`: the distribution `RBVD <: RandomDislocation` of random Burger's vectors and locations in the `axes`
"""
struct RandomDislocationDistribution{Dis <: Distribution, RBVD <: RandomDislocation} <: RandomStrainDistribution
    truncated::Bool
    concentration::AbstractFloat
    vector_diff::Function
    rand_num_dis::Dis
    burgers_vector_dist::RBVD
end

"""
    RandomDislocationDistribution(; expected_num_dislocations::Int, Lx, Ly = Lx, vector_diff = subtract_PBC!, burgers_vectors )

Convenient keyword constructor for the `RandomDislocationDistribution` struct which makes sure the distributions agree with other members.

# Additional Information
* This keyword constructor assumes that the `RandomDislocation` is a `UniformBurgersVector`.
* This does not use `StatsBase.Truncated` because it's about an order of magnitude slower than my version in `collect_dislocations`.
"""
function RandomDislocationDistribution(; expected_num_dislocations::Int, Lx, Ly = Lx, vector_diff = subtract_PBC!, burgers_vectors, truncated = true )
    diff = (A, B) -> vector_diff( A, B; Lx = Lx, Ly = Ly )
    if vector_diff != subtract_PBC!
        diff = vector_diff
    end
    bv_dist = UniformBurgersVector(; Lx = Lx, Ly = Ly, burgers_vectors = burgers_vectors )

    num_plaquettes = system_size(bv_dist)
    concentration = convert(Float64, expected_num_dislocations / num_plaquettes)
    rand_num = Binomial( num_plaquettes, concentration )

    return RandomDislocationDistribution(  truncated,
                                           concentration,
                                           diff,
                                           rand_num,
                                           bv_dist )
end

"""
    system_size( rdd::RandomDislocationDistribution )

Alias call to the `RandomDislocation` `system_size` function.
"""
system_size( rdd::RandomDislocationDistribution ) = system_size( rdd.burgers_vector_dist )

"""
    collect_dislocations( rdd::RandomDislocationDistribution )

Alias call to the `RandomDislocation` `collect_dislocations` function.

# Additional Information
* Note that this distribution will not allow for configurations with zero dislocations or for a number of dislocations larger than the `system_size`.
* This keyword constructor assumes that the `rand_num_dis` is a truncated `Binomial` `Distribution` on the interval ``[1, Lx \\times Ly]`` when `truncated == true`.
"""
function collect_dislocations( rdd::RandomDislocationDistribution )
    num_dis = rand(rdd.rand_num_dis)
    if rdd.truncated
        while num_dis == zero(typeof(num_dis)) || num_dis > system_size(rdd)
            num_dis = rand(rdd.rand_num_dis)
        end
    end
    return collect_dislocations( rdd.burgers_vector_dist, num_dis )
end

end # module RandomStrainDistributions

ext: jl


text:
function SeisProcessHeaders(in,out,functions,parameters;group="some",key=[],ntrace=1000000,update_tracenum=true)

	if (group=="all")
		h1 = SeisReadHeaders(in,group=group,key=key,itrace=1,ntrace=ntrace)
		for ifunc = 1 : length(functions)
			f = functions[ifunc]
			p = parameters[ifunc]
			h2 = f(h1;p...)
			h1 = copy(h2)
		end
		SeisWriteHeaders(out,h1)
	else
		itrace_in = 1
		itrace_out = 1
		NX = GetNumTraces(in)
		while itrace_in <= NX
			nx = NX - itrace_in + 1
			ntrace = nx > ntrace ? ntrace : nx
			h1 = SeisReadHeaders(in,group=group,key=key,itrace=itrace_in,ntrace=ntrace)
			num_traces_in = length(h1)
			for ifunc = 1 : length(functions)
				f = functions[ifunc]
				p = parameters[ifunc]
				h2 = f(h1;p...)
				h1 = copy(h2)
			end
			num_traces_out = length(h1)
			SeisWriteHeaders(out,h1,itrace=itrace_out,update_tracenum=update_tracenum)
			itrace_in += num_traces_in
			itrace_out += num_traces_out
		end
	end
end


text:
using Test, ImageCore
function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") in ("true", "True") ||
    get(ENV, "CI", "") in ("true", "True")
end

if is_ci()
    @info "CI testing is disabled due to the need for a Spinnaker-compatible camera during tests."
else
    using Spinnaker

    camlist = CameraList()
    if length(camlist) < 1
        error("""Spinnaker could not find a compatible camera.
        Tests require that a Spinnaker-compatible camera is connected and discoverable.
        If such a camera is connected, check that it is discoverable in SpinView""")
    else
        @testset "Camera Interaction" begin
            cam = camlist[0]
            @test typeof(cam) == Spinnaker.Camera

            @testset "Set continuous mode" begin
                triggermode!(cam, "Off")
                @test triggermode(cam) == "Off"

                acquisitionmode!(cam, "Continuous")
                @test acquisitionmode(cam) == "Continuous"

                framerate!(cam, 60.0)
                @test isapprox(framerate(cam), 60.0)
            end

            @testset "Set exposure" begin
                exposure!(cam, 4000)
                exp, mode = exposure(cam)
                @test isapprox(exp, 4000, rtol=0.1)
                @test mode == "Off"

                exposure!(cam)
                exp, mode = exposure(cam)
                @test mode == "Continuous"
            end

            @testset "Read image" begin
                @testset "Mono8" begin
                    pixelformat!(cam, "Mono8")
                    @assert pixelformat(cam) == "Mono8"
                    start!(cam)
                    img = getimage(cam, Gray{N0f8}, normalize=true)
                    @test all(size(img) .> 0)
                    stop!(cam)
                end
            end
        end
    end
end


text:
cd(@__DIR__)
using DiscreteEvents, DataStructures, Distributions, Random, BenchmarkTools
using Plots
const _bench = [false]

mutable struct Server
    capacity::Int
    ids::Vector{Int}
end

Server(capacity) = Server(capacity, Int[])

is_full(s) = length(s.ids) ≥ s.capacity

remove!(s, id) = filter!(x-> x != id, s.ids)

add!(s, id) = push!(s.ids, id)

Random.seed!(8710) # set random number seed for reproducibility
num_customers = 10 # total number of customers generated
num_servers = 2 # number of servers
mu = 1.0 / 2 # service rate
lam = 0.9 # arrival rate
arrival_dist = Exponential(1 / lam) # interarrival time distriubtion
service_dist = Exponential(1 / mu); # service time distribution

function enter_line(clock, server, id, service_dist, arrival_time)
    delay!(clock, arrival_time)
    _bench[end] ||  now!(clock, fun(println, "Customer $id arrived: ", clock.time))
    if is_full(server)
        _bench[end] || now!(clock, fun(println, "Customer $id is waiting: ", clock.time))
        wait!(clock, ()->!is_full(server))
    end
    _bench[end] ||  now!(clock, fun(println,"Customer $id starting service: ", clock.time))
    add!(server, id)
    tΔ = rand(service_dist)
    delay!(clock, tΔ)
    leave(clock, server, id)
    # now!(clock, fun(println, "servers available ", server.capacity-length(server.ids), " ", clock.time))
    # now!(clock, fun(println, "servers full ", is_full(server), " ", clock.time))
end

function leave(clock, server, id)
    _bench[end] || now!(clock, fun(println, "Customer $id finishing service: ", clock.time))
    remove!(server, id)
end

function initialize!(clock, arrival_dist, service_dist, num_customers, server)
    arrival_time = 0.0
    for i in 1:num_customers
        arrival_time += rand(arrival_dist)
        process!(clock, Prc(i, enter_line, server, i, service_dist, arrival_time), 1)
    end
end

function run_model(arrival_dist, service_dist, num_customers, num_servers, t)
    clock = Clock()
    server = Server(num_servers)
    initialize!(clock, arrival_dist, service_dist, num_customers, server)
    run!(clock, t)
end

function run_model(num_customers)
    result = @benchmark run_model($arrival_dist, $service_dist, $num_customers, $num_servers, 2000000)
    return mean(result).time*1e-9
end

_bench[end] = true
N = [10,1000,2000,4000]

times = run_model.(N)
plot(N, times, xlabel="Customers", ylabel="Time (seconds)", leg=false, grid=false)
savefig("DiscreteEvents_Queue.png")


text:
"""
    BalanceLaw

An abstract type representing a PDE balance law of the form

elements for balance laws of the form

```math
q_{,t} + Σ_{i=1,...d} F_{i,i} = s
```

Subtypes `L` should define the following methods:
- `vars_aux(::L)`: a tuple of symbols containing the auxiliary variables
- `vars_state(::L)`: a tuple of symbols containing the state variables
- `vars_state_for_transform(::L)`: a tuple of symbols containing the state variables which are passed to the `transform!` function.
- `vars_gradient(::L)`: a tuple of symbols containing the transformed variables of which gradients are computed
- `vars_diffusive(::L)`: a tuple of symbols containing the diffusive variables
- `flux_nondiffusive!(::L, flux::Grad, state::State, auxstate::State, t::Real)`
- `flux_diffusive!(::L, flux::Grad, state::State, diffstate::State, auxstate::State, t::Real)`
- `gradvariables!(::L, transformstate::State, state::State, auxstate::State, t::Real)`
- `diffusive!(::L, diffstate::State, ∇transform::Grad, auxstate::State, t::Real)`
- `source!(::L, source::State, state::State, diffusive::Vars, auxstate::State, t::Real)`
- `wavespeed(::L, nM, state::State, aux::State, t::Real)`
- `boundary_state!(::NumericalFluxGradient, ::L, stateP::State, auxP::State, normalM, stateM::State, auxM::State, bctype, t)`
- `boundary_state!(::NumericalFluxNonDiffusive, ::L, stateP::State, auxP::State, normalM, stateM::State, auxM::State, bctype, t)`
- `boundary_state!(::NumericalFluxDiffusive, ::L, stateP::State, diffP::State, auxP::State, normalM, stateM::State, diffM::State, auxM::State, bctype, t)`
- `init_aux!(::L, aux::State, coords, args...)`
- `init_state!(::L, state::State, aux::State, coords, args...)`

"""
abstract type BalanceLaw end # PDE part

# function stubs
function vars_state end
function vars_aux end
function vars_gradient end
function vars_diffusive end
vars_integrals(::BalanceLaw, FT) = @vars()
vars_reverse_integrals(::BalanceLaw, FT) = @vars()

num_aux(m::BalanceLaw, FT) = varsize(vars_aux(m,FT))
num_state(m::BalanceLaw, FT) = varsize(vars_state(m,FT)) # nstate
num_gradient(m::BalanceLaw, FT) = varsize(vars_gradient(m,FT))  # number_gradient_states
num_diffusive(m::BalanceLaw, FT) = varsize(vars_diffusive(m,FT)) # number_viscous_states
num_integrals(m::BalanceLaw, FT) = varsize(vars_integrals(m,FT))
num_reverse_integrals(m::BalanceLaw, FT) = varsize(vars_reverse_integrals(m,FT))

function update_aux! end
function integral_load_aux! end
function integral_set_aux! end
function reverse_integral_load_aux! end
function reverse_integral_set_aux! end
function flux_nondiffusive! end
function flux_diffusive! end
function gradvariables! end
function diffusive! end
function source! end
function wavespeed end
function boundary_state! end
function init_aux! end
function init_state! end

function calculate_dt end

function create_state(bl::BalanceLaw, grid, commtag)
  topology = grid.topology
  # FIXME: Remove after updating CUDA
  h_vgeo = Array(grid.vgeo)
  FT = eltype(h_vgeo)
  Np = dofs_per_element(grid)
  DA = arraytype(grid)

  weights = view(h_vgeo, :, grid.Mid, :)
  weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

  state = MPIStateArray{FT}(topology.mpicomm, DA, Np, num_state(bl,FT),
                            length(topology.elems),
                            realelems=topology.realelems,
                            ghostelems=topology.ghostelems,
                            vmaprecv=grid.vmaprecv, vmapsend=grid.vmapsend,
                            nabrtorank=topology.nabrtorank,
                            nabrtovmaprecv=grid.nabrtovmaprecv,
                            nabrtovmapsend=grid.nabrtovmapsend,
                            weights=weights, commtag=commtag)
  return state
end

function create_auxstate(bl, grid, commtag=222)
  topology = grid.topology
  Np = dofs_per_element(grid)

  h_vgeo = Array(grid.vgeo)
  FT = eltype(h_vgeo)
  DA = arraytype(grid)

  weights = view(h_vgeo, :, grid.Mid, :)
  weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

  auxstate = MPIStateArray{FT}(topology.mpicomm, DA, Np, num_aux(bl,FT),
                               length(topology.elems),
                               realelems=topology.realelems,
                               ghostelems=topology.ghostelems,
                               vmaprecv=grid.vmaprecv,
                               vmapsend=grid.vmapsend,
                               nabrtorank=topology.nabrtorank,
                               nabrtovmaprecv=grid.nabrtovmaprecv,
                               nabrtovmapsend=grid.nabrtovmapsend,
                               weights=weights, commtag=commtag)

  dim = dimensionality(grid)
  polyorder = polynomialorder(grid)
  vgeo = grid.vgeo
  device = typeof(auxstate.data) <: Array ? CPU() : CUDA()
  nrealelem = length(topology.realelems)
  @launch(device, threads=(Np,), blocks=nrealelem,
          initauxstate!(bl, Val(dim), Val(polyorder), auxstate.data, vgeo, topology.realelems))
  MPIStateArrays.start_ghost_exchange!(auxstate)
  MPIStateArrays.finish_ghost_exchange!(auxstate)

  return auxstate
end

function create_diffstate(bl, grid, commtag=111)
  topology = grid.topology
  Np = dofs_per_element(grid)

  h_vgeo = Array(grid.vgeo)
  FT = eltype(h_vgeo)
  DA = arraytype(grid)

  weights = view(h_vgeo, :, grid.Mid, :)
  weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

  # TODO: Clean up this MPIStateArray interface...
  diffstate = MPIStateArray{FT}(topology.mpicomm, DA, Np, num_diffusive(bl,FT),
                                length(topology.elems),
                                realelems=topology.realelems,
                                ghostelems=topology.ghostelems,
                                vmaprecv=grid.vmaprecv,
                                vmapsend=grid.vmapsend,
                                nabrtorank=topology.nabrtorank,
                                nabrtovmaprecv=grid.nabrtovmaprecv,
                                nabrtovmapsend=grid.nabrtovmapsend,
                                weights=weights, commtag=commtag)

  return diffstate
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""EndpointSlice represents a subset of the endpoints that implement a service. For a given service there may be multiple EndpointSlice objects, selected by labels, which must be joined to produce the full set of endpoints.

    IoK8sApiDiscoveryV1beta1EndpointSlice(;
        addressType=nothing,
        apiVersion=nothing,
        endpoints=nothing,
        kind=nothing,
        metadata=nothing,
        ports=nothing,
    )

    - addressType::String : addressType specifies the type of address carried by this EndpointSlice. All addresses in this slice must be the same type. This field is immutable after creation. The following address types are currently supported: * IPv4: Represents an IPv4 Address. * IPv6: Represents an IPv6 Address. * FQDN: Represents a Fully Qualified Domain Name.
    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - endpoints::Vector{IoK8sApiDiscoveryV1beta1Endpoint} : endpoints is a list of unique endpoints in this slice. Each slice may include a maximum of 1000 endpoints.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata.
    - ports::Vector{IoK8sApiDiscoveryV1beta1EndpointPort} : ports specifies the list of network ports exposed by each endpoint in this slice. Each port must have a unique name. When ports is empty, it indicates that there are no defined ports. When a port is defined with a nil port value, it indicates \"all ports\". Each slice may include a maximum of 100 ports.
"""
mutable struct IoK8sApiDiscoveryV1beta1EndpointSlice <: SwaggerModel
    addressType::Any # spec type: Union{ Nothing, String } # spec name: addressType
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    endpoints::Any # spec type: Union{ Nothing, Vector{IoK8sApiDiscoveryV1beta1Endpoint} } # spec name: endpoints
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiDiscoveryV1beta1EndpointPort} } # spec name: ports

    function IoK8sApiDiscoveryV1beta1EndpointSlice(;addressType=nothing, apiVersion=nothing, endpoints=nothing, kind=nothing, metadata=nothing, ports=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("addressType"), addressType)
        setfield!(o, Symbol("addressType"), addressType)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("endpoints"), endpoints)
        setfield!(o, Symbol("endpoints"), endpoints)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        o
    end
end # type IoK8sApiDiscoveryV1beta1EndpointSlice

const _property_map_IoK8sApiDiscoveryV1beta1EndpointSlice = Dict{Symbol,Symbol}(Symbol("addressType")=>Symbol("addressType"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("endpoints")=>Symbol("endpoints"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("ports")=>Symbol("ports"))
const _property_types_IoK8sApiDiscoveryV1beta1EndpointSlice = Dict{Symbol,String}(Symbol("addressType")=>"String", Symbol("apiVersion")=>"String", Symbol("endpoints")=>"Vector{IoK8sApiDiscoveryV1beta1Endpoint}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("ports")=>"Vector{IoK8sApiDiscoveryV1beta1EndpointPort}")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1EndpointSlice))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1EndpointSlice[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1EndpointSlice[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1EndpointSlice)
    (getproperty(o, Symbol("addressType")) === nothing) && (return false)
    (getproperty(o, Symbol("endpoints")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, name::Symbol, val)
end


text:
module SkeletonGrids

using Gridap

export SkeletonGrid
import Gridap: points
import Gridap: cells
import Gridap: celltypes
import Gridap: cellorders
import Gridap: Triangulation
import Gridap: SkeletonTriangulation

struct SkeletonGrid{D,Z} <: Grid{D,Z}
  grid::Grid{D,Z}
  descriptor1::BoundaryDescriptor
  descriptor2::BoundaryDescriptor

  function SkeletonGrid(
    grid::Grid{D,Z},
    descriptor1::BoundaryDescriptor,
    descriptor2::BoundaryDescriptor) where {D,Z}
    @assert D == Z + 1
    new{D,Z}(grid,descriptor1,descriptor2)
  end

end

for op in (:points,:cells,:celltypes,:cellorders)
  @eval begin
    $op(g::SkeletonGrid) = $op(g.grid)
  end
end

function Triangulation(grid::SkeletonGrid)
  trian = Triangulation(grid.grid)
  SkeletonTriangulation(trian,grid.descriptor1,grid.descriptor2)
end

function SkeletonGrid(model::DiscreteModel)
  d = celldim(model)
  graph = GridGraph(model)
  facet_to_cells = connections(graph,d-1,d)
  nfacets = length(facet_to_cells)
  mask = fill(false,nfacets)
  _find_on_skeleton!(mask,facet_to_cells)
  SkeletonGrid(model,mask)
end

function _find_on_skeleton!(mask,facet_to_cells)
  for (facet,cells) in enumerate(facet_to_cells)
    if length(cells) == 2
      mask[facet] = true
    end
  end
end

function SkeletonGrid(model::DiscreteModel, tags)
  cell1 = 1
  cell2 = 2
  bgrid1 = BoundaryGrid(model,tags,cell1)
  bgrid2 = BoundaryGrid(model,tags,cell2)
  grid = bgrid1.grid
  descriptor1 = bgrid1.descriptor
  descriptor2 = bgrid2.descriptor
  SkeletonGrid(grid,descriptor1,descriptor2)
end

function SkeletonTriangulation(model::DiscreteModel,tags)
  grid = SkeletonGrid(model,tags)
  Triangulation(grid)
end

function SkeletonTriangulation(model::DiscreteModel)
  grid = SkeletonGrid(model)
  Triangulation(grid)
end

end # module


text:
"""
    HookList{TNext, THandler, TPlugin}

Provides fast, inlinable call to the implementations of a specific hook.

You can get a HookList by calling `hooklist()` directly, or using `hooks()`.

The `HookList` can be called with arbitrary number of extra arguments. If any of the
plugins referenced in the list fails to handle the extra arguments, the call will raise a `MethodError`
"""
struct HookList{TNext, THandler, TPlugin}
    next::TNext
    handler::THandler
    plugin::TPlugin
end

struct HookListTerminal end

@inline function (hook::HookList)(args...)::Bool
    if hook.handler(hook.plugin, args...) !== true
        return hook.next(args...)
    end
    return true
end

(hook::HookListTerminal)(_...) = false

@inline function call_optional(hook::Plugins.HookList, args...)
    if !isnothing(hook.handler)
        val = hook.handler(hook.plugin, args...)
        if !isnothing(val)
            return val
        end
    end
    return call_optional(hook.next, args...)
end
call_optional(hook::HookListTerminal, _...) = nothing

length(l::HookListTerminal) = 0
length(l::HookList) = 1 + length(l.next)

iterate(l::HookList) = (l, l.next)
iterate(::HookList, state) = isnothing(state) ? nothing : (state, state.next)
iterate(::HookList, ::HookListTerminal) = nothing

"""
    function hooklist(plugins, hookfn)

Create a HookList which allows fast, inlinable call to the merged implementations of `hookfn`
by the given plugins.

A plugin of type `TPlugin` found in plugins will be referenced in the resulting HookList if there is a method
that matches the following signature: `hookfn(::TPlugin, ...)`
"""
function hooklist(plugins, hookfn)
    if length(plugins) == 0
        return HookListTerminal()
    end
    plugin = plugins[1]
    if length(methods(hookfn, (typeof(plugin), Vararg{Any}))) > 0
        return HookList(hooklist(plugins[2:end], hookfn), hookfn, plugin)
    end
    return hooklist(plugins[2:end], hookfn)
end
hooklist(stack::PluginStack, hookfn) = hooklist(stack.plugins, hookfn)

"""
    hooks(app)
    hooks(pluginstack::PluginStack)

Create or get a hook cache for `stack`.

The first form can be used when `pluginstack` is stored in `app.plugins` (the recommended pattern).

When this function is called first time on a `PluginStack`, the hooks cache will be created by calling
`hook_cache()`, and stored in `pluginstack` for quick access later.
"""
function hooks end

@inline hooks(app) = hooks(app.plugins)
@inline hooks(stack::PluginStack) = stack.hookcache

function rebuild_cache!(stack::PluginStack)
    stack.symbolcache = symbolcache(stack.plugins)
    stack.hookcache = hook_cache(stack)
end

"""
    hook_cache(stack::PluginStack, hookfns)
    hook_cache(plugins, hookfns)

Create a cache of `HookList`s for a PluginStack or from lists of plugins and hook functions.

Returns a NamedTuple with an entry for every handler.

# Examples

```julia
cache = hook_cache([Plugin1(), Plugin2()], [hook1, hook2])
cache.hook1()
```
"""
function hook_cache(plugins, hookfns)
    return (;(nameof(hook) => hooklist(plugins, hook) for hook in hookfns)...)
end
hook_cache(stack::PluginStack) = hook_cache(stack.plugins, stack.hookfns)

struct LifecycleHook
    op::Function
end

Base.string(lfhook::LifecycleHook) = string(lfhook.op)

function (lfhook::LifecycleHook)(stack::PluginStack, data...)
    allok = true
    results = []
    for plugin in reverse(stack.plugins)
        try
            @debug "Calling $(string(lfhook)) lifecycle hook of $(typeof(plugin))"
            pushfirst!(results, lfhook.op(plugin, data...))
        catch e
            allok = false
            @debug "Got error $(e) from $(string(lfhook)) lifecycle hook of $(typeof(plugin))"
            pushfirst!(results, (e, catch_backtrace()))
        end
    end
    return (allok = allok, results = results)
end

function create_lifecyclehook(op::Function)
    return LifecycleHook(op)
end


text:
include("training_script.jl")

init_training(; save_tag="31-12-20_resnet152", learning_rate_by_epochs = [(0.005, 1), (0.02, 17), (0.002, 12), (0.001, 12)], should_use_pmap = true)


text:
using Colors, GLVisualize
using Gadfly, GLVisualize.ComposeBackend

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end
const static_example = true

p = plot([sin, cos], 0, 25)

draw(composebackend, p)
if !isdefined(:runtests)
renderloop(window)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation <: SwaggerModel
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    url::Any # spec type: Union{ Nothing, String } # spec name: url

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation(;description=nothing, url=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation, Symbol("url"), url)
        setfield!(o, Symbol("url"), url)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation = Dict{Symbol,Symbol}(Symbol("description")=>Symbol("description"), Symbol("url")=>Symbol("url"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation = Dict{Symbol,String}(Symbol("description")=>"String", Symbol("url")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1ExternalDocumentation }, name::Symbol, val)
end


text:
"""
Truncation in a box between `a` and `b`.
"""
struct BoxTruncationRegion <: TruncationRegion
    a::Vector{Float64}
    b::Vector{Float64}
end

intruncationregion(r::BoxTruncationRegion, x::AbstractArray) = all(r.a .<= x) && all(x .<= r.b)

"""
An Elliptical Truncation region.
"""
struct EllipticalTruncationRegion <: TruncationRegion
    H::PDMat
    h::Vector{Float64}
    c::Float64
end

intruncationregion(r::EllipticalTruncationRegion, x::AbstractArray) = (x-r.h)'*r.H*(x-r.h) <= r.c

abstract type PolytopeTruncationRegion <: TruncationRegion end


text:
function terms(orb::RelativisticOrbital, w::Int=one(Int))
    j = kappa_to_j(orb.κ)
    w <= 2j+1 || throw(ArgumentError("w=$w too large for $orb orbital"))

    2w ≥ 2j+1 && (w = convert(Int, 2j) + 1 - w)
    w == 0 && return [zero(HalfInteger)]
    w == 1 && return [j]

    # Forms full Cartesian product of all mⱼ, not necessarily the most
    # performant method.
    mⱼs = filter(allunique, collect(allchoices([-j:j for i = 1:w])))
    MJs = map(x -> reduce(+, x), mⱼs) # TODO: make map(sum, mⱼs) work

    Js = HalfInteger[]

    while !isempty(MJs)
        # Identify the maximum MJ and associate it with J.
        MJmax = maximum(MJs)
        N = count(isequal(MJmax), MJs)
        append!(Js, repeat([MJmax], N))
        # Remove the -MJ:MJ series, N times.
        for MJ = -MJmax:MJmax
            deleteat!(MJs, findall(isequal(MJ), MJs)[1:N])
        end
    end

    # Do we really want unique here?
    sort(unique(Js))
end

# This is a workaround until seniority number are implemented for
# jj-coupled subshells.
intermediate_terms(orb::RelativisticOrbital, w::Int=one(Int)) =
    terms(orb, w)

export terms, intermediate_terms


text:
export update_graph, is_conducting, swap_index

using SparseArrays

"""
```julia
swap_index(a)
```
Takes an iterable collection of iterable collections of the
form `(a,b,c,...)` and returns an object with "reversed" nesting
of the form `((a[1],b[1],c[1],...),(a[2],b[2],c[2],...),...)`.
"""
function swap_index(a)
    # if a is empty then set n=0 otherwise n=length(a[1])
    n = length(a)
    n==0 || (n=length(a[1]))
    collect(broadcast(x->getindex(x,i),a) for i ∈ 1:n)
end

"""
```julia
update_graph(nwn::NWN_circuit{:JDA,T,N})
```
Generates the graph of the electrical connections under the Junction Dominated Assumption model, `:JDA`, in the nanowire network `nwn` by mutating the `graph`, `ls`, and `ns` fields of nwn.
"""
function update_graph(nwn::NWN_circuit{:JDA,T,N}) where {T,N}
    ww_src, ww_dst = connections(nwn.lines, nwn.dims)
    we_src, we_dst = connections(nwn.lines, nwn.elecs, nwn.dims)
    wg_src, wg_dst = connections(nwn.lines, nwn.grnds, nwn.dims)
    nₗ, nₑ, n₀ = length(nwn.lines), length(nwn.elecs), length(nwn.grnds)
    n = nₗ + nₑ + n₀
    mₗ, mₑ, m₀ = length(ww_src), length(we_src), length(wg_src)
    src = vcat(ww_src, we_src, wg_src)
    dst = vcat(ww_dst, we_dst.+nₗ, wg_dst.+(nₗ+nₑ))
    wgt = vcat(ones(mₗ)./nwn.Rⱼ,ones(mₑ+m₀)./nwn.Rₑ)
    sp = sparse(vcat(src,dst), vcat(dst,src), vcat(wgt,wgt),n,n,+)
    nwn.graph=SimpleWeightedGraph(sp)
    nwn.ls = ones(Int,n)
    nwn.ns = cumsum(nwn.ls)
    nwn
end

"""
```julia
update_graph(nwn::NWN_circuit{:MNR,T,N})
```
Generates the graph of the electrical connections under the MNR model, `:MNR`, in the nanowire network `nwn` by mutating the `graph`, `ls`, and `ns` fields of nwn.
"""
function update_graph(nwn::NWN_circuit{:MNR,T,N}) where {T,N}
    nₗ, nₑ, n₀ = length(nwn.lines), length(nwn.elecs), length(nwn.grnds)
    n = nₗ + nₑ + n₀
    line_ip = SVector{n}([T[] for i ∈ 1:n])
    juncts = NamedTuple{(:R,:c₁,:c₂),Tuple{T,Tuple{Int,Int},Tuple{Int,Int}}}[]
    each_intersect(nwn.lines, nwn.dims, []) do i₁,i₂,arr,dims,ps,p,arg
        push!(line_ip[i₁],ps[1])
        push!(line_ip[i₂],ps[2])
        push!(juncts,(R=nwn.Rⱼ,c₁=(i₁,length(line_ip[i₁])),c₂=(i₂,length(line_ip[i₂]))))
    end
    each_intersect(nwn.lines, vcat(nwn.elecs,nwn.grnds), nwn.dims, []) do i₁,i₂,arr₁,arr₂,dims,ps,p,arg
        i₂ += nₗ
        push!(line_ip[i₁],ps[1])
        push!(line_ip[i₂],ps[2])
        push!(juncts,(R=nwn.Rₑ,c₁=(i₁,length(line_ip[i₁])),c₂=(i₂,length(line_ip[i₂]))))
    end
    nwn.ls = length.(line_ip) # number of intersections on each wire
    nwn.ns = cumsum(nwn.ls) # number of nodes up to wire i
    # to count number of nodes involved in inner-wire resistances:
    # if ls[i]==0 then line i does not add to the node count or add a resistance
    # if ls[1]==1 then line i has no inner wire resistances; -1 to nodes in inner-wires
    # if ls[1]>1 then line i has i-1 inner wire resistances; -1 to nodes in inner-wires
    # ⟹ if ls[i]>0 then -1 to nodes in inner-wires
    ms = vcat([0],cumsum(nwn.ls[1:nₗ].>0))
    nᵢ = nwn.ns[nₗ]-ms[end]
    nⱼ = length(juncts)
    perms = sortperm.(line_ip)
    invperms = invperm.(perms)
    permute!.(line_ip,perms)

    src, dst, wgt = zeros(Int,nᵢ+nⱼ), zeros(Int,nᵢ+nⱼ), zeros(T,nᵢ+nⱼ)
    for i ∈ 1:nₗ
        js = 1:(nwn.ls[i]-1)
        slice = js.+(nwn.ns[i]-nwn.ls[i]-ms[i])
        Rₗ = 4*len(nwn.lines[i])*nwn.props[i].ρ/(π*nwn.props[i].D*nwn.props[i].D)
        src[slice] = js.+(nwn.ns[i]-nwn.ls[i])
        dst[slice] = js.+(1+nwn.ns[i]-nwn.ls[i])
        wgt[slice] = 1 ./(Rₗ.*(diff(line_ip[i])))
    end
    for (i,ju) ∈ enumerate(juncts)
        src[i+nᵢ] = nwn.ns[ju.c₁[1]]-nwn.ls[ju.c₁[1]]+invperms[ju.c₁[1]][ju.c₁[2]]
        dst[i+nᵢ] = nwn.ns[ju.c₂[1]]-nwn.ls[ju.c₂[1]]+invperms[ju.c₂[1]][ju.c₂[2]]
        wgt[i+nᵢ] = 1/ju.R
    end
    sp = sparse(vcat(src,dst), vcat(dst,src), vcat(wgt,wgt),nwn.ns[end],nwn.ns[end],+)
    nwn.graph = SimpleWeightedGraph(sp)
    nwn
end

"""
```julia
is_conducting(nwn::NWN_circuit)
```
Returns a boolean value for if the given network has a path between at least one electrode node and one ground node.
"""
function is_conducting(nwn::NWN_circuit)
    nₗ, nₑ = length(nwn.lines), length(nwn.elecs)
    for i ∈ (nwn.ns[nₗ]+1):nwn.ns[nₗ+nₑ]
        for j ∈ (nwn.ns[nₗ+nₑ]+1):nwn.ns[end]
            if has_path(nwn.graph,i,j)
                return true
            end
        end
    end
    return false
end


text:
#=

# 201: 2D Laplace equation
([source code](SOURCE_URL))

=#

module Example201_Laplace2D

using VoronoiFVM,ExtendableGrids
using GridVisualize

## Flux function which describes the flux
## between neigboring control volumes
function g!(f,u,edge)
    f[1]=u[1,1]-u[1,2]
end

function main(;Plotter=nothing,n=5)
    nspecies=1
    ispec=1
    X=collect(0:1.0/n:1)
    grid=VoronoiFVM.Grid(X,X)

    physics=VoronoiFVM.Physics(flux=g!)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,ispec,[1])
    boundary_dirichlet!(sys,ispec,1,0.0)
    boundary_dirichlet!(sys,ispec,3,1.0)
    inival=unknowns(sys,inival=0)
    solution=unknowns(sys)
    solve!(solution,inival,sys)
    nf=nodeflux(sys,solution)
    vis=GridVisualizer(Plotter=Plotter)
    scalarplot!(vis,grid,solution[1,:],clear=true,colormap=:summer)
    vectorplot!(vis,grid,nf[:,1,:],clear=false,spacing=0.1,vscale=0.5)
    reveal(vis)
    return solution[7]
end

## Called by unit test

function test()
    main() ≈ 0.2
end

end


text:
mutable struct EFITEquilibrium{T<:Real, S<:AbstractRange{Float64},
                               R<:AbstractMatrix{Float64},
                               Q<:AbstractVector{Float64}} <: AbstractEquilibrium
    cocos::COCOS       # COCOS
    r::S               # Radius/R range
    z::S               # Elevation/Z range
    psi::S             # Polodial Flux range (polodial flux from magnetic axis)
    psi_rz::R          # Polodial Flux on RZ grid (polodial flux from magnetic axis)
    g::Q               # Polodial Current
    p::Q               # Plasma pressure
    q::Q               # Q profile
    phi::Q             # Electric Potential
    axis::NTuple{2, T} # Magnetic Axis (raxis,zaxis)
    sigma::Int         # sign(dot(J,B))
end

function efit(cc::COCOS, r::S, z::S, psi::S, psi_rz, g, p,
              q, phi, axis::NTuple{2,T}, sigma::Int) where {T,S<:AbstractRange}

    psi_rz_itp = CubicSplineInterpolation((r,z), psi_rz, extrapolation_bc=Flat())
    if step(psi) > 0
        g_itp = CubicSplineInterpolation(psi, g, extrapolation_bc=Flat())
        p_itp = CubicSplineInterpolation(psi, p, extrapolation_bc=Flat())
        q_itp = CubicSplineInterpolation(psi, q, extrapolation_bc=Flat())
        phi_itp = CubicSplineInterpolation(psi, phi, extrapolation_bc=Flat())
    else # CubicSplineInterpolation doesn't like decreasing psi so reverse them
        g_itp = CubicSplineInterpolation(reverse(psi), reverse(g), extrapolation_bc=Flat())
        p_itp = CubicSplineInterpolation(reverse(psi), reverse(p), extrapolation_bc=Flat())
        q_itp = CubicSplineInterpolation(reverse(psi), reverse(q), extrapolation_bc=Flat())
        phi_itp = CubicSplineInterpolation(reverse(psi), reverse(phi), extrapolation_bc=Flat())
    end
    EFITEquilibrium(cc, r, z, psi, psi_rz_itp, g_itp, p_itp, q_itp, phi_itp, axis, Int(sigma))
end

function Base.show(io::IO, N::EFITEquilibrium)
    print(io, "EFITEquilibrium{$(eltype(N.psi))}")
end

Base.broadcastable(N::EFITEquilibrium) = (N,)

cocos(N::EFITEquilibrium) = N.cocos

function (N::EFITEquilibrium)(r,z)
    return N.psi_rz(r,z)
end

function magnetic_axis(N::EFITEquilibrium)
    return N.axis
end

function B0Ip_sign(N::EFITEquilibrium)
    return N.sigma
end

function limits(N::EFITEquilibrium)
    return extrema(N.r), extrema(N.z)
end

function psi_limits(N::EFITEquilibrium)
    return N.psi[1], N.psi[end]
end

function psi_gradient(N::EFITEquilibrium,r,z)
    return SVector{2}(Interpolations.gradient(N.psi_rz, r, z))
end

function poloidal_current(N::EFITEquilibrium, psi)
    return N.g(psi)
end

function poloidal_current_gradient(n::EFITEquilibrium, psi)
    return Interpolations.gradient(n.g, psi)[1]
end

function pressure(N::EFITEquilibrium, psi)
    return N.p(psi)
end

function pressure_gradient(n::EFITEquilibrium, psi)
    return Interpolations.gradient(n.p, psi)[1]
end

function safety_factor(N::EFITEquilibrium, psi)
    return N.q(psi)
end

function electric_potential(N::EFITEquilibrium, psi)
    return N.phi(psi)
end

function electric_potential_gradient(N::EFITEquilibrium, psi)
    return Interpolations.gradient(N.phi, psi)[1]
end


text:
# Copyright 2015 Conrad Scott

immutable Corpus{S <: AbstractString, D <: AbstractVector, I <: Integer, C <: AbstractMatrix}
    v :: Dictionary{S, D}

    # A term-document matrix: A row per term; a column per document; values are document word counts.
    c :: C

    function Corpus(v :: Dictionary{S, D}, c :: AbstractMatrix{I})
        new(v, c)
    end
end

function Corpus(v :: Dictionary, c :: AbstractMatrix)
  Corpus{eltype(v.d), typeof(v.d), eltype(c), typeof(c)}(v, c)
end

function words(c :: Corpus)
  size(c.c, 1)
end

function documents(c :: Corpus)
  size(c.c, 2)
end


text:
module RegressionTools

using Distances: euclidean
using StatsFuns: logistic, softplus, logit
using PLINK

import Base.vecnorm

export threshold!
export threshold
export update_residuals!
export update_indices!
export update_col!
export update_weights!
export update_xk!
export fill_perm!
export update_xb!
export count_partialnz
export fill_partial!
export update_partial_residuals!
export difference!
export ypatzmw!
export project_k!
export selectpermk!
export selectpermk
export mce
export logistic_loglik
export update_x!
export df_norm
export logistic_grad!
export logistic!
export log2xb!
export fit_logistic
export mask!
export cv_get_folds
export vecnorm
export issymmetric
export vec!

const Float = Union{Float32, Float64}

depwarn("NOTE: As of Julia v0.6, this package is deprecated and no longer maintained.\n")

# do not load any functions; deprecated package will now be empty
#include("regtool.jl")
#include("logistic.jl")

end # end module RegressionTools


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using GLMakie, Colors, ColorSchemes, Random
GLMakie.activate!() #HIDE
let
    Random.seed!(38)
    x = y = 1:100
    z = exp10.(2*rand(100,100))
    fig = Figure(resolution = (850, 800), fontsize = 20)
    ax = Axis(fig, xscale = log10, yscale = log10, aspect = 1, xlabel = "x",
        ylabel = "y", yminorticksvisible = true, xminorticksvisible = true,
        yminorticks = IntervalsBetween(9),xminorticks = IntervalsBetween(9))
    cmap = cgrad(:CMRmap, scale = :log)
    hmap = heatmap!(x,y,z, colormap = cmap, colorrange = (1,100))
    cbar = Colorbar(fig, hmap, label = "z values", width = 25, ticksize=25,
        minorticks = IntervalsBetween(9), minorticksvisible=true,labelsize=18,
        ticklabelsize = 18,tickalign = 1, minortickalign =1,
        height = Relative(0.87),scale = log10) # log10 colorbar ticks
    #ax.xticks = [1,10,100]
    #ax.yticks = [1,10,100]
    #cbar.ticks = [1,10,100]
    fig[1, 1] = ax
    fig[1, 2] = cbar
    #colgap!(fig.layout, 7) # HIDE
    fig
    save(joinpath(@__DIR__, "output", "heatmapLogIrregular.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie", "Colors", "ColorSchemes", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#############################################################################
# sumlargesteigs.jl
# Handles top k eigenvalues of a symmetric or Hermitian matrix
# (and imposes the constraint that its argument be symmetric or Hermitian)
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
import LinearAlgebra.eigvals

### sumlargesteigs

struct SumLargestEigs <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,AbstractExpr}
    size::Tuple{Int,Int}

    function SumLargestEigs(x::AbstractExpr, k::AbstractExpr)
        children = (x, k)
        m, n = size(x)
        if m == n
            return new(:sumlargesteigs, hash(children), children, (1, 1))
        else
            error("sumlargesteigs can only be applied to a square matrix.")
        end
    end
end

function sign(x::SumLargestEigs)
    return NoSign()
end

function monotonicity(x::SumLargestEigs)
    return (Nondecreasing(), NoMonotonicity())
end

function curvature(x::SumLargestEigs)
    return ConvexVexity()
end

function evaluate(x::SumLargestEigs)
    return eigvals(evaluate(x.children[1]))[end-x.children[2]:end]
end

function sumlargesteigs(x::AbstractExpr, k::Int)
    return k == 0 ? Constant(0) : SumLargestEigs(x, Constant(k))
end

# Create the equivalent conic problem:
#   minimize sk + Tr(Z)
#   subject to
#            Z ⪰ 0
#            A ⪰ 0
#            Z + sI ⪰ A
# See Ben-Tal and Nemirovski, "Lectures on Modern Convex Optimization"
# Example 18.c

function conic_form!(x::SumLargestEigs, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        k = x.children[2]
        m, n = size(A)
        if sign(A) == ComplexSign()
            Z = ComplexVariable(n, n)
        else
            Z = Variable(n, n)
        end
        s = Variable()
        # The two inequality constraints have the side effect of constraining A to be symmetric,
        # since only symmetric matrices can be positive semidefinite.
        p = minimize(
            s * k + real(tr(Z)),
            Z + s * Matrix(1.0I, n, n) - A ⪰ 0,
            Z ⪰ 0,
        )
        cache_conic_form!(unique_conic_forms, x, p)
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
module A
export foo, bar
export to_function

f(x) = 0
g(x) = x + 1

function to_function(expr)
    @eval f(x) = $expr
    f
end
function foo(expr)
    f = to_function(expr)
    for x=1:1000
       f(x)
    end
    f
end
function bar(expr)
    f = to_function(expr)
    for x=1:1000
        g(x)
    end
    f
end
end #module

using A

@time expr = :(x+1)
@time to_function(expr)
A.f(1)
foo(expr)
bar(expr)

for i=1:1000000
    println("$i: foo-bar")
    i2 = rand(1:i)
    i3 = rand(1:i)
    i4 = rand(1:i)
    i5 = rand(1:i)
    i6 = rand(1:i)
    i7 = rand(1:i)
    i8 = rand(1:i)
    expr = :(x+$i+$i2+$i3+$i4+$i5+$i6+$i7+$i8)
    @time foo(expr)
    @time bar(expr)
end


text:
#=
Creates a basic neuron with 3 inputs.

Associated YT NNFS tutorial: https://www.youtube.com/watch?v=Wo5dMEP_BbI
=#

inputs = [1.2, 5.1, 2.1]
weights = [3.1, 2.1, 8.7]
bias = 3.0

output = inputs[1]*weights[1] + inputs[2]*weights[2] + inputs[3]*weights[3] + bias
println(output)

max_stars_repo_name: yarikbratashchuk/NNfSiX

max_issues_repo_name: yarikbratashchuk/NNfSiX

max_forks_repo_name: yarikbratashchuk/NNfSiX


text:
using SafeTestsets

@safetestset "Day09" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Computer
    using AdventOfCode.Day09

    @testset "quine" begin
        quine = [109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99]
        outputs = run_program!(Program(quine))
        @test quine == outputs
    end

    @testset "large numbers" begin
        @test run_program!([1102,34915192,34915192,7,4,7,99,0]).output == 1219070632396864
        @test run_program!([104,1125899906842624,99]).output == 1125899906842624
    end

    btest("9.1") do
        part1() == 2494485073
    end

    btest("9.2") do
        part2() == 44997
    end
end


text:
## Methods giving basic information about graphs ##

"""
    memsize(g::Graph)

Approximate size of `g` in memory, in bytes.
"""
memsize(g::Graph) =
    sizeof(g) + sizeof(g.nodes) +
    sum(Int[sizeof(n.forward) + sizeof(n.reverse) for n = g.nodes]) +
    sizeof(g.nodeattrs) + sizeof(g.edgeattrs)

"""
    nodes(g::Graph)

Indices of all nodes in the graph.
"""
nodes(g::Graph) = 1:length(g.nodes)

"""
    nodecount(g::Graph)

Count nodes in the graph.
"""
nodecount(g::Graph) = length(g.nodes)

# Matrix-stype size of Graph
size(g::Graph) = (nodecount(g), nodecount(g))
size(g::Graph, i::Integer) = i <= 2 ? nodecount(g) : 1

"""
    edgecount(g::Graph)

Count edges in the graph.
"""
edgecount(g::Graph) = g.edgecount

"""
    maxedgecount(g::Graph)

Maximum number of edges possible in a graph having no multiple edges and self-loops.
"""
maxedgecount(g::DirectedGraph)   = (n = nodecount(g); n * (n - 1))
maxedgecount(g::UndirectedGraph) = (n = nodecount(g); div(n * (n - 1), 2))


text:
using MultipleScattering
using Plots; pyplot()
using Random
using LinearAlgebra

function box_size_convergence(  m = 4,
                                volfrac = 0.05,
                                radius = 1.0,
                                times = [20.0,30.0,40.0,50.0,60.0,70.0],
                                ω = collect(LinRange(0.01,1.0,100))
                                )

    host_medium = Acoustic(1.0, 1.0, 2)
    particle_medium = Acoustic(1.0, 0.0, 2)
    particle_shape = Circle(radius)

    listener_position = [-10.0; 0.0]
    bigshape = TimeOfFlightPlaneWaveToPoint(listener_position, maximum(times))

    seed = MersenneTwister(1).seed
    allparticles = random_particles(particle_medium, particle_shape;
                                            region_shape = bigshape,
                                            volume_fraction = volfrac,
                                            seed = seed
                                    )

    source =  plane_source( host_medium; position = [0.0, 0.0], direction = [1.0, 0.0] )

    simulations = map(times) do t
        println("Calculating response with particles at a maximum of $t away")
        shape = TimeOfFlightPlaneWaveToPoint(listener_position, t)
        particles = filter(p -> p ⊆ shape, allparticles)
        return run(particles, source, [listener_position], ω; basis_order = m)
    end

    return map(x -> frequency_to_time(x; t_vec = LinRange(0.0,100.0,201)), simulations)
end

function plot_box_size_convergence(simulations; times = [20,30,40,50,60,70])

    plot(xguide = "Time (t)", yguide ="Response")
    for s in eachindex(simulations)
        plot!(simulations[s], label = "Box cut at t=$(times[s])")
        plot!(size=(800,600))
    end
    plot!(title="Time response from random particles from increasingly large boxes")
end

simulations = box_size_convergence()
plot_box_size_convergence(simulations)


text:
using Knet, KnetModules
import KnetModules.forward

type Fire <: KnetModule
    squeeze::Conv
    expand1x1::Conv
    expand3x3::Conv
end

function Fire(inp::Int, sq::Int, e1x1::Int, e3x3::Int)
    squeeze   = Conv(1, 1, inp, sq)
    expand1x1 = Conv(1, 1, sq, e1x1)
    expand3x3 = Conv(3, 3, sq, e3x3; padding=1)
    return Fire(squeeze, expand1x1, expand3x3)
end

function forward(ctx, f::Fire, x)
    o = relu.(@mc f.squeeze(x))
    return ccat(relu.(@mc f.expand1x1(o)),
                relu.(@mc f.expand3x3(o)))
end

# FIXME: too many memory allocations
# FIXME: model is too slow in gpu, possibly because of this
function ccat(a1, a2)
    @assert size(a1, 1,2,4) == size(a2, 1,2,4)
    a1_ = mat(permutedims(a1, (1,2,4,3)))
    a2_ = mat(permutedims(a2, (1,2,4,3)))
    out = hcat(a1_, a2_)
    out = reshape(out, (size(a1, 1,2,4)...,
                        size(out, 2)))
    return permutedims(out, (1,2,4,3))
end

"""
# Usage
    s = SqueezeNet(;num_classes=1000, pdrop=0.5, trained=true)
    @run s(x)
"""
type SqueezeNet <: KnetModule
    features::Sequential
    last_conv::Conv
    pdrop
end

function SqueezeNet(;num_classes=1000, pdrop=0.5, trained=true)
    features = Sequential(
        Conv(3, 3, 3, 64; stride=2),
        ReLU(),
        MaxPool(;window=(3,3), stride=(2,2)),
        Fire(64, 16, 64, 64),
        Fire(128, 16, 64, 64),
        #Fire(128, 32, 128, 128),
        MaxPool(;window=(3,3), stride=(2,2)),
        Fire(128, 32, 128, 128),
        Fire(256, 32, 128, 128),
        MaxPool(;window=(3,3), stride=(2,2)),
        Fire(256, 48, 192, 192),
        Fire(384, 48, 192, 192),
        Fire(384, 64, 256, 256),
        Fire(512, 64, 256, 256)
    )
    last_conv = Conv(1, 1, 512, num_classes;
                     winit=(et, a...)->et(0.01) .* rand(et, a...))
    return SqueezeNet(features, last_conv, pdrop)
end

function forward(ctx, sn::SqueezeNet, x)
    o = @mc sn.features(x)
    o = dropout(o, sn.pdrop)
    o = @mc sn.last_conv(o)
    o = pool(o; window=size(o,1,2), mode=2)
    return mat(o)
end


text:
export rule

@rule AR(:θ, Marginalisation) (q_y_x::MultivariateNormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta) = begin
    order = getorder(meta)
    F     = getvform(meta)

    myx, Vyx = mean(q_y_x), cov(q_y_x)
    my, Vy   = ar_slice(F, myx, 1:order), ar_slice(F, Vyx, 1:order, 1:order)
    mx, Vx   = ar_slice(F, myx, order+1:2order), ar_slice(F, Vyx, order+1:2order, order+1:2order)
    Vyx      = ar_slice(F, Vyx, order+1:2order, 1:order)

    mγ = mean(q_γ)

    D = mγ * (Vx + mx * mx')
    c = ar_unit(getvform(meta), order)

    mθ, Vθ = cholinv(D) * (Vyx + mx * my') * mγ * c, cholinv(D)

    return convert(promote_variate_type(F, NormalMeanVariance), mθ, Vθ)
end


text:
@testset "Mortal walks" begin
    walk = WalkF(status_type = Mortal())

    lattice_rates = LatticeVar{1}(Exponential(1e6))
    waiting_times = LatticeVar{1}(Pareto(0.5))

    ql = Lattice(waiting_times, lattice_rates)

    latwalk = LatticeWalk(ql, walk)

    times = logrange(1:.3:6)
    step_limit_actor = StepLimitActor(10^4)
    storing_actor = storing_nsteps_position_actor(times)

    actors = ActorSet(step_limit_actor, storing_actor)
    walk!(latwalk, actors)
    @test  typeof(Actors.get_values(storing_actor)) == Tuple{Vector{Int64}, Vector{Float64}}

    storing_actor = storing_nsteps_actor(times)
    actors = ActorSet(step_limit_actor, storing_actor)
    walk!(latwalk, actors)
    @test  typeof(Actors.get_values(storing_actor)) == Tuple{Vector{Int64}}

    storing_actor = storing_position_actor(times)
    actors = ActorSet(step_limit_actor, storing_actor)
    walk!(latwalk, actors)
    @test  typeof(Actors.get_values(storing_actor)) == Vector{Float64}

    lwp = WalkPlan(latwalk, ActorSet(step_limit_actor, storing_actor, NullActor()))
    walk!(lwp)
    @test  typeof(Actors.get_values(storing_actor)) == Vector{Float64}
end


text:
# # Simple example of a homogeneous ensemble using learning networks

# In this simple example, no bagging is used, so every atomic model
# gets the same learned parameters, unless the atomic model training
# algorithm has randomness, eg, DecisionTree with random subsampling
# of features at nodes.

# Note that MLJ has a built in model wrapper called `EnsembleModel`
# for creating bagged ensembles with a few lines of code.

# ## Definition of composite model type

using MLJ, PyPlot
import Statistics

# learning network (composite model spec):

Xs = source()
ys = source(kind=:target)

atom = @load DecisionTreeRegressor
atom.n_subfeatures = 4 # to ensure diversity among trained atomic models

machines = (machine(atom, Xs, ys) for i in 1:100)

# overload `mean` for nodes:
Statistics.mean(v...) = mean(v)
Statistics.mean(v::AbstractVector{<:AbstractNode}) = node(mean, v...)

yhat = mean([predict(m, Xs) for  m in machines]);

# new composite model type and instance:

one_hundred_models = @from_network OneHundredModels(atom=atom) <= yhat

# ## Application to data

X, y = @load_boston;

# tune regularization parameter for a *single* tree:

r = range(atom,
          :min_samples_split,
          lower=2,
          upper=100, scale=:log)

mach = machine(atom, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e1.svg")) # hide

# \fig{e1.svg}

# tune regularization parameter for all trees in ensemble simultaneously:

r = range(one_hundred_models,
          :(atom.min_samples_split),
          lower=2,
          upper=100, scale=:log)

mach = machine(one_hundred_models, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e2.svg")) # hide

# \fig{e2}

#-

PyPlot.close_figs() # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
#! /usr/bin/env julia

HERE = @__FILE__
println(@__FILE__)
println(splitdir(HERE))

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
using Krylov, Stopping

#https://github.com/JuliaSmoothOptimizers/Krylov.jl

include("random_coordinate_descent_method.jl")
include("stop_random_coordinate_descent_method.jl")
include("cheap_stop_random_coordinate_descent_method.jl")
include("instate_coordinate_descent_method.jl")

#using ProfileView,
using BenchmarkTools

n=5000;
A=rand(n,n)
sol=rand(n)
b=A*sol

x0   = zeros(size(A,2))
pb   = LinearSystem(A,b)
s0   = GenericState(x0, similar(b))
mcnt = Stopping._init_max_counters_linear_operators()

###############################################################################
#
# The original algorithm
#
###############################################################################
@time x, OK, k = RandomizedCD(A, b, max_iter = 100)

###############################################################################
#
# Stopping with memory optimization using the State and with cheap_stop!
#
###############################################################################
@time meta3 = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 100,
                        retol = false,
                        tol_check = (atol, rtol, opt0) -> (atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s3   = GenericState(zeros(size(A,2)), similar(b))
@time stp3 = LAStopping(pb, meta3, s3)
@time StopRandomizedCD2(stp3)

###############################################################################
#
# Stopping version of the algorithm
#
###############################################################################
@time meta2 = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 100,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s2   = GenericState(zeros(size(A,2)), similar(b))
@time stp2 = LAStopping(pb, meta2, s2)
@time StopRandomizedCD(stp2)

###############################################################################
#
# Stopping version of the algorithm with cheap remote control
#
###############################################################################
@time meta1 = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 100,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s1   = GenericState(zeros(size(A,2)), similar(b))
@time stp1 = LAStopping(pb, meta1, cheap_stop_remote_control(), s1)
@time StopRandomizedCD(stp1)

###############################################################################
#
# Stopping version of the algorithm with cheap remote control
#
###############################################################################
@time meta = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 100,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s0   = GenericState(zeros(size(A,2)), similar(b))
@time stp = LAStopping(pb, meta, StopRemoteControl(), s0)
@time StopRandomizedCD(stp)

###############################################################################
#
# Stopping with memory optimization using the State and with cheap_stop!
# uses the @ macro instate.
#
#DOESN'T WORK ??
#
###############################################################################
#=
@time meta4 = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 100,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s4   = GenericState(zeros(size(A,2)), similar(b))
@time stp4 = LAStopping(pb, meta4, s4)
@time stp4 = StopRandomizedCD3(stp4)
=#

Lnrm = [norm(stp.current_state.current_score),
        norm(stp1.current_state.current_score),
        norm(stp2.current_state.current_score),
        norm(stp3.current_state.current_score), #norm(stp4.current_state.current_score),
        norm(b-A*x)
       ]

using Test
@test Lnrm ≈ minimum(Lnrm) * ones(length(Lnrm)) atol = 1e-7

nothing;


text:
module Judycon

"""
    QuickFind
"""
struct QuickFind
    id :: Vector{Int}
end

function QuickFind(n::Int)
    id = collect(1:n)
    return QuickFind(id)
end

function find(qf::QuickFind, p::Int)
    id = qf.id
    return id[p]
end

function isconnected(qf::QuickFind, p::Int, q::Int)
    i = find(qf, p)
    j = find(qf, q)
    return isequal(i, j)
end

function connect!(qf::QuickFind, p::Int, q::Int)
    id = qf.id
    pid = id[p]
    qid = id[q]
    n = length(id)
    for i=1:n
        if isequal(getindex(id, i), pid)
            setindex!(id, qid, i)
        end
    end
    return
end

"""
    QuickUnion
"""
struct QuickUnion
    id :: Vector{Int}
    sz :: Vector{Int}
    weighted :: Bool
    compress :: Bool
end

function QuickUnion(n::Int, weighted=true, compress=true)
    id = collect(1:n)
    sz = zeros(Int, n)
    return QuickUnion(id, sz, weighted, compress)
end

function find(qu::QuickUnion, p::Int)
    id = qu.id
    while p != id[p]
        if qu.compress
            id[p] = id[id[p]]
        end
        p = id[p]
    end
    return p
end

function isconnected(qu::QuickUnion, p::Int, q::Int)
    i = find(qu, p)
    j = find(qu, q)
    return isequal(i, j)
end

"""
    connect(G, p, q)

Connect p and q in G.

If weighting is used, connect such a way that the depth of the tree is minimized.
"""
function connect!(qu::QuickUnion, p::Int, q::Int)
    id = qu.id
    i = find(qu, p)
    j = find(qu, q)
    if qu.weighted
        sz = qu.sz
        if (sz[i] < sz[j])
            id[i] = j
            sz[j] += sz[i]
        else
            id[j] = i
            sz[i] += sz[j]
        end
    else
        id[i] = j
    end
end

export QuickFind, QuickUnion, connect!, isconnected, find

end


text:
"Supertype of tensor product dictionaries."
abstract type TensorProductDict{S,T} <: Dictionary{S,T} end

# component(dict::TensorProductDict, range::AbstractRange) = tensorproduct(dict.dicts[range]...)
function component(dict::TensorProductDict, range::AbstractRange)
    # TODO: remove this method in the future
    @warn "Range selection of product dictionaries now simply returns the components"
    [component(dict,i) for i in range]
end

factors(dict::TensorProductDict) = components(dict)

product_domaintype(dicts::Dictionary...) = Tuple{map(domaintype, dicts)...}
function product_domaintype(dicts::Vararg{Dictionary{<:Number},N}) where {N}
    T = promote_type(map(domaintype, dicts)...)
    SVector{N,T}
end

function TensorProductDict(dicts::Dictionary...)
    N = length(dicts)
    S = product_domaintype(dicts...)
    T = promote_type(map(codomaintype, dicts)...)
    C = promote_type(map(coefficienttype, dicts)...)
    dicts2 = map(dict->ensure_coefficienttype(C,dict),dicts)
    DT = typeof(dicts2)
    TupleProductDict{N,DT,S,T}(dicts2)
end

tolerance(dict::TensorProductDict) = minimum(map(tolerance,components(dict)))

dimensions(d::TensorProductDict) = map(dimensions, components(d))

coefficienttype(d::TensorProductDict) = promote_type(map(coefficienttype,components(d))...)

IndexStyle(d::TensorProductDict) = IndexCartesian()

^(d::Dictionary, n::Int) = tensorproduct(d, n)

similar(dict::TensorProductDict, ::Type{T}, size::Int) where {T} =
    similar(dict, T, approx_length(dict, size))

## Properties

for op in (:isreal, :isbasis, :isframe)
    @eval $op(s::TensorProductDict) = mapreduce($op, &, components(s))
end

for op in (:isorthogonal, :isorthonormal)
    @eval $op(s::TensorProductDict, m::ProductWeight) = mapreduce($op, &, components(s), components(m))
end

for op in (:isorthogonal, :isorthonormal)
    @eval BasisFunctions.$op(s::TensorProductDict, m::BasisFunctions.DiscreteProductWeight) = mapreduce($op, &, components(s), components(m))
end

for op in (:isorthogonal, :iscompatible)
    @eval BasisFunctions.$op(s::TensorProductDict, m::BasisFunctions.ProductGrid) = mapreduce($op, &, components(s), components(m))
end

## Feature methods

for op in (:hasinterpolationgrid, :hasextension, :hasderivative, :hasantiderivative)
    @eval $op(s::TensorProductDict) = mapreduce($op, &, components(s))
end
hasderivative(Φ::TensorProductDict, order) =
    mapreduce(hasderivative, &, components(Φ), order)
hasantiderivative(Φ::TensorProductDict, order) =
    mapreduce(hasantiderivative, &, components(Φ), order)

hasgrid_transform(s::TensorProductDict, gb, grid::ProductGrid) =
    mapreduce(hastransform, &, components(s), components(grid))

hasgrid_transform(s::TensorProductDict, gb, grid::AbstractGrid) = false

for op in (:derivative_dict, :antiderivative_dict)
    @eval $op(s::TensorProductDict, order; options...) =
        tensorproduct( map( (el,ord) -> $op(el, ord; options...), components(s), order)... )
end

for op in (:differentiation, :antidifferentiation)
    @eval function $op(::Type{T}, src::TensorProductDict, dest::TensorProductDict, order; options...) where {T}
        @assert length(order) == dimension(src)
        @assert length(order) == dimension(dest)
        tensorproduct(map( (el_s,el_d,ord) -> $op(T, el_s, el_d, ord; options...), components(src), components(dest), order)...)
    end
end

diff(dict::TensorProductDict, order; options...) =
    TensorProductDict( (diff(e,o) for (e,o) in zip(components(dict),order))...)

resize(d::TensorProductDict, n::Int) = resize(d, approx_length(d, n))
resize(d::TensorProductDict, dims) = TensorProductDict(map(resize, components(d), dims)...)

# Delegate dict_in_support to _dict_in_support with the composing dicts as extra arguments,
# in order to avoid extra memory allocation.
dict_in_support(dict::TensorProductDict, idx, x) =
    _dict_in_support(dict, components(dict), idx, x)
# catch CartesianIndex, convert to tuple, so that iteration works
dict_in_support(dict::TensorProductDict, idx::CartesianIndex, x) =
    _dict_in_support(dict, components(dict), Tuple(idx), x)

_dict_in_support(::TensorProductDict, dicts, idx, x) = mapreduce(in_support, &, dicts, idx, x)

function approx_length(s::TensorProductDict, n::Int)
    # Rough approximation: distribute n among all dimensions evenly, rounded upwards
    N = dimension(s)
    m = ceil(Int, n^(1/N))
    tuple([approx_length(component(s, j), m^dimension(s, j)) for j in 1:ncomponents(s)]...)
end

extensionsize(s::TensorProductDict) = map(extensionsize, components(s))

getindex(s::TensorProductDict, ::Colon, i::Int) = (@assert ncomponents(s)==2; component(s,1))
getindex(s::TensorProductDict, i::Int, ::Colon) = (@assert ncomponents(s)==2; component(s,2))
getindex(s::TensorProductDict, ::Colon, ::Colon) = (@assert ncomponents(s)==2; s)

getindex(s::TensorProductDict, ::Colon, i::Int, j::Int) = (@assert ncomponents(s)==3; component(s,1))
getindex(s::TensorProductDict, i::Int, ::Colon, j::Int) = (@assert ncomponents(s)==3; component(s,2))
getindex(s::TensorProductDict, i::Int, j::Int, ::Colon) = (@assert ncomponents(s)==3; component(s,3))
getindex(s::TensorProductDict, ::Colon, ::Colon, i::Int) =
    (@assert ncomponents(s)==3; TensorProductDict(component(s,1),component(s,2)))
getindex(s::TensorProductDict, ::Colon, i::Int, ::Colon) =
    (@assert ncomponents(s)==3; TensorProductDict(component(s,1),component(s,3)))
getindex(s::TensorProductDict, i::Int, ::Colon, ::Colon) =
    (@assert ncomponents(s)==3; TensorProductDict(component(s,2),component(s,3)))
getindex(s::TensorProductDict, ::Colon, ::Colon, ::Colon) = (@assert ncomponents(s)==3; s)

interpolation_grid(s::TensorProductDict) =
    ProductGrid(map(interpolation_grid, components(s))...)
gauss_rule(s::TensorProductDict) =
    productmeasure(map(gauss_rule, components(s))...)

support(s::TensorProductDict) = productdomain(map(support, components(s))...)
support(s::TensorProductDict, idx::LinearIndex) = support(s, native_index(s,idx))
support(s::TensorProductDict, idx::ProductIndex) = productdomain(map(support, components(s), indextuple(idx))...)

# We pass on the elements of s as an extra argument in order to avoid
# memory allocations in the lines below
# unsafe_eval_element(set::TensorProductDict, idx, x) = _unsafe_eval_element(set, components(set), indexable_index(set, idx), x)
unsafe_eval_element(dict::TensorProductDict, idx::ProductIndex, x) =
    _unsafe_eval_element(dict, components(dict), idx, x)

_unsafe_eval_element(dict::TensorProductDict, dicts, i, x) =
    mapreduce(unsafe_eval_element, *, dicts, i, x)
_unsafe_eval_element(dict::TensorProductDict, dicts, i::CartesianIndex, x) =
    mapreduce(unsafe_eval_element, *, dicts, Tuple(i), x)

unsafe_eval_element_derivative(dict::TensorProductDict, idx::ProductIndex, x, order) =
    _unsafe_eval_element_derivative(dict, components(dict), idx, x, order)

_unsafe_eval_element_derivative(dict::TensorProductDict, dicts, i, x, order) =
    mapreduce(unsafe_eval_element_derivative, *, dicts, i, x, order)
_unsafe_eval_element_derivative(dict::TensorProductDict, dicts, i::CartesianIndex, x, order) =
    mapreduce(unsafe_eval_element_derivative, *, dicts, Tuple(i), x, order)

hasmeasure(dict::TensorProductDict) = mapreduce(hasmeasure, &, components(dict))
measure(dict::TensorProductDict) = productmeasure(map(measure, components(dict))...)

innerproduct_native(Φ1::TensorProductDict, i, Φ2::TensorProductDict, j, measure::ProductWeight; options...) =
    mapreduce(innerproduct, *, components(Φ1), Tuple(i), components(Φ2), Tuple(j), components(measure))

evaluation(::Type{T}, dict::TensorProductDict, gb::GridBasis, grid::ProductGrid; options...) where {T} =
    tensorproduct(map( (d,g) -> evaluation(T, d, g; options...), components(dict), components(grid))...)

dual(dict::TensorProductDict, measure::Union{ProductWeight,DiscreteProductWeight}=measure(dict); options...) =
    TensorProductDict([dual(dicti, measurei; options...) for (dicti, measurei) in zip(components(dict),components(measure))]...)

gram(::Type{T}, dict::TensorProductDict, measure::Union{ProductWeight,DiscreteProductWeight}; options...) where {T} =
    TensorProductOperator(map((x,y)->gram(T, x,y; options...), components(dict), components(measure))...)
mixedgram(::Type{T}, dict1::TensorProductDict, dict2::TensorProductDict, measure::Union{ProductWeight,DiscreteProductWeight}; options...) where {T} =
    TensorProductOperator(map((x,y,z)->mixedgram(T, x,y,z; options...), components(dict1), components(dict2), components(measure))...)

Display.combinationsymbol(d::TensorProductDict) = Display.Symbol('⊗')
Display.displaystencil(d::TensorProductDict) = composite_displaystencil(d)
show(io::IO, mime::MIME"text/plain", d::TensorProductDict) = composite_show(io, mime, d)
show(io::IO, d::TensorProductDict) = composite_show_compact(io, d)

"A flat tensor product dict has `N` scalar components."
abstract type FlatTensorProductDict{N,S,T} <: TensorProductDict{S,T}
end

const TensorProductDict1{DT,S,T} = FlatTensorProductDict{1,S,T}
const TensorProductDict2{DT,S,T} = FlatTensorProductDict{2,S,T}
const TensorProductDict3{DT,S,T} = FlatTensorProductDict{3,S,T}
const TensorProductDict4{DT,S,T} = FlatTensorProductDict{4,S,T}

similar(dict::FlatTensorProductDict{N}, ::Type{T}, size::Vararg{Int,N}) where {T,N} =
    TensorProductDict(map(similar, components(dict), T.parameters, size)...)

## Native indices are of type ProductIndex

# The native indices of a tensor product dict are of type ProductIndex
ordering(d::FlatTensorProductDict{N}) where {N} = ProductIndexList{N}(size(d))

native_index(d::TensorProductDict, idx) = product_native_index(size(d), idx)

# We have to amend the boundscheck ecosystem to catch some cases:
# - This line will catch indexing with tuples of integers, and we assume
#   the user wanted to use a CartesianIndex
checkbounds(::Type{Bool}, d::FlatTensorProductDict{N}, idx::NTuple{N,Int}) where {N} =
    checkbounds(Bool, d, CartesianIndex(idx))
checkbounds(::Type{Bool}, d::FlatTensorProductDict{2}, idx::NTuple{2,Int}) =
    checkbounds(Bool, d, CartesianIndex(idx))
# - Any other tuple we assume is a recursive native index, which we convert
#   elementwise to a tuple of linear indices
checkbounds(::Type{Bool}, d::TensorProductDict, idx::Tuple) =
    checkbounds(Bool, d, map(linear_index, components(d), idx))

"""
    struct TupleProductDict{N,DT,S,T} <: Dictionary{S,T}

Parameters:
- DT is a tuple of types, representing the (possibly different) types of the dicts.
- N is the dimension of the product (equal to the length of the DT tuple)
- S is the domain type
- T is the codomain type.
"""
struct TupleProductDict{N,DT,S,T} <: FlatTensorProductDict{N,S,T}
    dicts   ::  DT
    size    ::  NTuple{N,Int}

    TupleProductDict{N,DT,S,T}(dicts::DT) where {N,DT,S,T} = new(dicts, map(length, dicts))
end

# Generic functions for composite types:
components(d::TupleProductDict) = d.dicts

size(d::TupleProductDict) = d.size

"Return a list of all tensor product indices (1:s+1)^n."
index_set_tensorproduct(s,n) = CartesianIndices(ntuple(k->1,n), tuple(ntuple(k->s+1,n)))

"Return a list of all indices of total degree at most s, in n dimensions."
function index_set_total_degree(s, n)
    # We make a list of arrays because that is easier
    I = _index_set_total_degree(s, n)
    # and then we convert to tuples
    [tuple((1+i)...) for i in I]
end

function _index_set_total_degree(s, n)
    if n == 1
        I = [[i] for i in 0:s]
    else
        I = Array{Array{Int,1}}(0)
        I_rec = _index_set_total_degree(s, n-1)
        for idx in I_rec
            for m in 0:s-sum(abs.(idx))
                push!(I, [idx...; m])
            end
        end
        I
    end
end


text:
import Base: @__doc__
using PowerDynamics: @DynamicNode
import PowerDynamics: construct_vertex, dimension, symbolsof, showdefinition, AbstractNode
using LinearAlgebra: Diagonal
using NetworkDynamics: ODEVertex

begin
    @__doc__ struct DGUnitPLL <: AbstractNode
            I_r
            K_pll
            Y_n
        end
    DGUnitPLL(; I_r, K_pll, Y_n) = DGUnitPLL(I_r, K_pll, Y_n)
    function construct_vertex(par::DGUnitPLL)
        I_r = par.I_r
        K_pll = par.K_pll
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt
            i = total_current(e_s, e_d) + Y_n * u
            θ = x[3]

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = - K_pll * abs(u) * sin(θ - angle(u))

            du = i - I_r*exp(θ*1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=3, mass_matrix=Diagonal([0, 0, 1]), sym=Symbol[:u_r, :u_i, :θ])
    end
    symbolsof(::DGUnitPLL) = begin
            [:u_r, :u_i, :θ]
        end
    dimension(::DGUnitPLL) = begin
            3
        end
end

begin
    @__doc__ struct DGUnitPLLPQ <: AbstractNode
            I_r
            K_pll
            S_pq
            Y_n
        end
    DGUnitPLLPQ(; I_r, K_pll, S_pq, Y_n) = DGUnitPLLPQ(I_r, K_pll, S_pq, Y_n)
    function construct_vertex(par::DGUnitPLLPQ)
        I_r = par.I_r
        K_pll = par.K_pll
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            θ = x[3]

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = - K_pll * abs(u) * sin(θ - angle(u))

            du = i - I_r*exp(θ*1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=3, mass_matrix=Diagonal([0, 0, 1]), sym=Symbol[:u_r, :u_i, :θ])
    end
    symbolsof(::DGUnitPLLPQ) = begin
            [:u_r, :u_i, :θ]
        end
    dimension(::DGUnitPLLPQ) = begin
            3
        end
end

begin
    @__doc__ struct DGUnitPLLPQTracking <: AbstractNode
            K_pll
            K_PT1
            T_PT1
            S_pq
            Y_n
        end
    DGUnitPLLPQTracking(; K_pll, K_PT1, T_PT1, S_pq, Y_n) = DGUnitPLLPQTracking(K_pll, K_PT1, T_PT1, S_pq, Y_n)
    function construct_vertex(par::DGUnitPLLPQTracking)
        K_pll = par.K_pll
        K_PT1 = par.K_PT1
        T_PT1 = par.T_PT1
        S_pq = par.S_pq
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            # current from incoming/outgoing edges, nodal shunt, PQ load background
            i = total_current(e_s, e_d) + Y_n * u - conj(S_pq) / conj(u)
            Vamp = abs(u)

            θ = x[3]
            P_err = x[4]
            Q_err = x[5]
            P_g = x[6]
            Q_g = x[7]

            # power flow tracking
            ΔP, ΔQ = p

            dP_err = - ΔP
            dQ_err = - ΔQ
            dP_g = (K_PT1 * P_err / 2. - P_g) / T_PT1
            dQ_g = (K_PT1 * Q_err / 2. - Q_g) / T_PT1

            # current set point adjustment
            I_r = complex(P_g, -Q_g) / Vamp

            # PLL
            # (x[2] * cos(θ) - x[1] * sin(θ))
            dθ = - K_pll * Vamp * sin(θ - angle(u))

            du = i - I_r*exp(θ*1im) # = 0 constraint

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dθ
                dx[4] = dP_err
                dx[5] = dQ_err
                dx[6] = dP_g
                dx[7] = dQ_g
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=7, mass_matrix=Diagonal([0, 0, 1, 1, 1, 1, 1]), sym=Symbol[:u_r, :u_i, :θ, :P_err, :Q_err, :P_g, :Q_g])
    end
    symbolsof(::DGUnitPLLPQTracking) = begin
            [:u_r, :u_i, :θ, :P_err, :Q_err, :P_g, :Q_g]
        end
    dimension(::DGUnitPLLPQTracking) = begin
            7
        end
end

using OrdinaryDiffEq: ODEProblem, DiscreteCallback
import PowerDynamics: NodeShortCircuit

function (nsc::NodeShortCircuit)(powergrid::ControlledPowerGrid)
    # Currently this assumes that the lines are PiModels...
    lines = copy(powergrid.lines)
    l_idx = findfirst(l -> (l.from == nsc.node || l.to == nsc.node) && l isa PiModelLine, lines)

    if l_idx == nothing
        @warn "Node needs to be connected to a PiModelLine to implement NodeShortCircuit"
        return nothing
    end

    if lines[l_idx].from == nsc.node
        lines[l_idx] = PiModelLine(;from=lines[l_idx].from, to=lines[l_idx].to, y = lines[l_idx].y, y_shunt_km = 1/nsc.R,  y_shunt_mk = lines[l_idx].y_shunt_mk)
    elseif lines[l_idx].to == nsc.node
        lines[l_idx] = PiModelLine(;from=lines[l_idx].from, to=lines[l_idx].to, y = lines[l_idx].y, y_shunt_km = lines[l_idx].y_shunt_km,  y_shunt_mk = 1/nsc.R)
    end
    ControlledPowerGrid(powergrid.graph, powergrid.nodes, lines, powergrid.controller)
end

function simulate_nsc(nsc::NodeShortCircuit, powergrid::ControlledPowerGrid, x1, timespan)
    @assert first(timespan) <= nsc.sc_timespan[1] "fault cannot begin in the past"
    @assert nsc.sc_timespan[2] <= last(timespan) "fault cannot end in the future"
    nsc_powergrid = nsc(powergrid)

    problem = ODEProblem{true}(rhs(powergrid), x1, timespan)

    function errorState(integrator)
        sol1 = integrator.sol
        x2 = find_valid_initial_condition(nsc_powergrid, sol1[end]) # Jump the state to be valid for the new system.
        println(x2)
        integrator.f = rhs(nsc_powergrid)
        integrator.u = x2
    end

    function regularState(integrator)
        sol2 = integrator.sol
        x3 = find_valid_initial_condition(nsc_powergrid, sol2[end]) # Jump the state to be valid for the new system.
        println(x3)
        integrator.f = rhs(powergrid)
        integrator.u = x3
    end

    cb1 = DiscreteCallback(((u,t,integrator) -> t in nsc.sc_timespan[1]), errorState)
    cb2 = DiscreteCallback(((u,t,integrator) -> t in nsc.sc_timespan[2]), regularState)
    sol = solve(problem, Rodas4(), callback = CallbackSet(cb1, cb2), tstops=[nsc.sc_timespan[1];nsc.sc_timespan[2]])
    return sol
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
function update_obj!(problem, s, y, ∇fx, z, ∇fz, B, scheme, scale = nothing)
    fz, ∇fz = upto_gradient(problem, ∇fz, z)
    # add Project gradient

    # Update y
    @. y = ∇fz - ∇fx

    # Update B
    if scale == nothing
        if isa(scheme.approx, Direct)
            yBy = dot(y, B * y)
            if !iszero(yBy)
                Badj = dot(y, s) / yBy .* B
            end
        else
            ys = dot(y, s)
            if !iszero(ys)
                Badj = dot(y, B * y) / ys .* B
            else
                return fz, ∇fz, B, s, y
            end
        end
    else
        Badj = B
    end
    # Quasi-Newton update
    B = update!(scheme, Badj, s, y)

    return fz, ∇fz, B, s, y
end

function update_obj!(problem, s, y, ∇fx, z, ∇fz, B, scheme::Newton, scale = nothing)
    fz, ∇fz, B = upto_hessian(problem, ∇fz, B, z)

    return fz, ∇fz, B, s, s
end

function update_obj(problem, s, ∇fx, z, ∇fz, B, scheme, scale = nothing)
    fz, ∇fz = upto_gradient(problem, ∇fz, z)
    # add Project gradient

    # Update y
    y = ∇fz - ∇fx

    # Update B
    if scale == nothing
        if isa(scheme.approx, Direct)
            yBy = dot(y, B * y)
            if !iszero(yBy) && isfinite(yBy)
                Badj = dot(y, s) / yBy * B # this is different than above?
            else
                Badj = B
            end
        else
            ys = dot(y, s)
            if !iszero(ys) && isfinite(ys)
                Badj = dot(y, B * y) / ys * B
            else
                Badj = B
            end
        end
    else
        Badj = B
    end

    # Quasi-Newton update
    B = update(scheme, Badj, s, y)

    return fz, ∇fz, B, s, y
end

function update_obj(problem, s, ∇fx, z, ∇fz, B, scheme::Newton, is_first = nothing)
    fz, ∇fz, B = upto_hessian(problem, ∇fx, B, z)

    return fz, ∇fz, B, s, nothing
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
abstract type Game{B,M} end

"""
    boardtype(game)

Return the type of the board for game
"""
boardtype(::Game{B,M}) where {B,M} = B
boardtype(g::Type{G}) where G <: Game = boardtype(g())

"""
    movetype(game)

Return the type of the move for game
"""
movetype(::Game{B,M}) where {B,M} = M
movetype(g::Type{G}) where G <: Game = movetype(g())

"""
    play(board, move)

Play move on the board and return the resulting board
"""
function play end

"""
    play!(board, move)

Play move on the board, changing the argument
"""
function play! end

"""
    getmoves(board)

Return the available moves. If the array is empty, the game is assumed over
"""
function getmoves end

"""
    hasmoves(board) -> Bool

Return if there are any moves left or if the game is over
"""
function hasmoves end
hasmoves(board) = !isempty(getmoves(board))

"""
    initialboard(game)
"""
function initialboard end

initialboard(game::Game{B,M}) where {B,M} = B()

"""
    playerturn(board) -> Int

Return the player who's to move next
"""
function playerturn end

"""
    winner(board) -> Int

Return the player who just won the game
"""
function winner end

"""
    nplayers(game)

Return the number of players in the game
"""
function nplayers end
nplayers(game::G) where {G<:Game} = nplayers(G)


text:
macro sayhello()
    return :( println("Hello, world!") )
end

@sayhello

macro sayhello(name)
    return :( println("Hello, ", $name) )
end

@sayhello "MAC110"

ex = macroexpand(Main, :(@sayhello("human")) )

@macroexpand @sayhello "human"

@macroexpand @macroexpand @sayhello "human"

macro twostep(arg)
    println("I execute at parse time. The argument is: ", arg)
    return :(println("I execute at runtime. The argument is: ", $arg))
end

ex = macroexpand(Main, :(@twostep :(1, 2, 3)) );

typeof(ex)

ex

eval(ex)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


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

using CellwiseEditingDifferentiation
using Distributions
using Gadfly
using Lora

hyperpars = Dict{Symbol, Any}(:λ=>fill(10.0, nsites))

prior = Dict{Symbol, Any}(:v=>Function[(m::Float64, v::Float64, a::Float64, b::Float64)->
  vpcprior(m, v, a, b, hyperpars[:λ][i]) for i in 1:nsites])

prior_target = Dict{Symbol, Any}(:w=>Array(Function, nsites))
for i in 1:nsites
  prior_target[:w][i] = function (w::Vector{Float64})
    v = data[:m][i]*(1-data[:m][i])*inv_logit(w[1])
    a, b = beta_pars_from_mv(data[:m][i], v)
    prior[:v][i](data[:m][i], v, a, b)+logdvdw(w[1], data[:m][i])
  end
end

prior_init = Dict{Symbol, Any}(:w=>fill(0.0, nsites))

for i in 1:nsites
  np = 100000

  vposterior = vec(readdlm(joinpath(OUTDIR, @sprintf("vchain_%s_site%02d.txt", string(simulationid), i))))

  pposterior = Array(Float64, np)

  for j in 1:np
    a, b = beta_pars_from_mv(data[:m][i], rand(vposterior))
    pposterior[j] = rand(Beta(a, b))
  end

  println("Mean of pposterior for site $i = $(mean(pposterior))")
  println("Var of pposterior for site $i = $(var(pposterior))")

  mcmodel = model(prior_target[:w][i], init=[prior_init[:w][i]])
  mcsampler = RAM()
  mcrunner = SerialMC(nsteps=int(1.01*np), burnin=int(0.1*np))
  mcchain = run(mcmodel, mcsampler, mcrunner)
  mcchain = map(w->data[:m][i]*(1-data[:m][i])*inv_logit(w), vec(mcchain.samples))

  pprior = Array(Float64, np)

  for j in 1:np
    a, b = beta_pars_from_mv(data[:m][i], rand(mcchain))
    pprior[j] = rand(Beta(a, b))
  end

  println("Mean of pprior for site $i = $(mean(pprior))")
  println("Var of pprior for site $i = $(var(pprior))")

  colors = Dict{Symbol, Color.RGB{Float64}}(:posterior=>color("red"), :prior=>color("blue"))

  layers = Layer[]

  push!(layers, layer(
    x=pposterior,
    Stat.histogram(bincount=50, density=true),
    Geom.line,
    Theme(default_color=colors[:posterior])
  )[1])

  push!(layers, layer(
    x=pprior,
    Stat.histogram(bincount=50, density=true),
    Geom.line,
    Theme(default_color=colors[:prior])
  )[1])

  vplot = plot(
    layers,
    Guide.xlabel("p$i"),
    Guide.title("Marginal histogram of p$i")
    # Guide.manual_color_key("Distribution", [string(k) for k in keys(colors)], [c for c in values(colors)])
  )

  draw(PDF(joinpath(OUTDIR, @sprintf("marginal_phist_%s_site%02d.pdf", string(simulationid), i)), 4inch, 3inch), vplot)
end


text:
include("types.jl")
include("score.jl")
for line in eachline()
    hands,results = split(line, ':')
    results = reshape(map(hex2int, collect(reverse(results))), (4,5))
    (par,bid,decl) = parscore(results)
    println(join((line, par), ' '))
end

max_stars_repo_name: JuliaTagBot/play

max_issues_repo_name: JuliaTagBot/play

max_forks_repo_name: JuliaTagBot/play