text:
# The SLOPE R library is GPL so cannot be directly included in this package. Instead
# we provide an installation and build function that users need to call manually.

SLOPE_CODE_URL = "http://www-stat.stanford.edu/~candes/SortedL1/SLOPE_code.tgz"
SLOPEDestDir = joinpath(dirname(@__FILE__()), "slope_code")
#SLOPEDestDir = "/Users/feldt/dev/RegularizedRegression.jl/src/utils/slope_code"
SLOPEUnpackDir = "SLOPE_code" # Name of the dir where the tar file unpacks to

using Requests
using RCall

# When this global var has been set to a string path the SLOPE code has been successfully loaded
global SLOPECodeInPathLoaded = false

function install_SLOPE_R_code(url = SLOPE_CODE_URL;
    destdir = SLOPEDestDir, unpackdir = SLOPEUnpackDir)

    # Make destdir if not already there
    if !isdir(destdir)
        mkdir(destdir)
    end

    cd(destdir) do
        tarball = download_SLOPE_code_tarball_and_save_to_file(url)

        # Unpack tarball
        res1 = readstring(`tar xvfz $tarball`)

        # Go into dir and...
        cd(unpackdir) do
            #...build c code into R callable code
            res2 = readstring(`R CMD SHLIB cproxSortedL1.c proxSortedL1.c`)

            #...indicate that we have built the SLOPE code
            global SLOPECodeInPathLoaded
            SLOPECodeInPathLoaded = joinpath(destdir, unpackdir)
        end
    end
end

function download_SLOPE_code_tarball_and_save_to_file(url)
    filename = "SLOPE_code.tgz"
    open(filename, "w") do fh
        write(fh, Requests.readbytes(get(url)))
    end
    return filename
end

function is_installed_SLOPE()
    global SLOPECodeInPathLoaded
    if SLOPECodeInPathLoaded != false
        return isfile(joinpath(SLOPECodeInPathLoaded, "cproxSortedL1.so"))
    else
        default_path = joinpath(SLOPEDestDir, SLOPEUnpackDir)
        if isfile(joinpath(default_path, "cproxSortedL1.so"))
            SLOPECodeInPathLoaded = default_path
            return true
        else
            return false
        end
    end
end

function load_SLOPE_R_lib()
    if !is_installed_SLOPE()
        install_SLOPE_R_code()
    end

    global SLOPECodeInPathLoaded
    RCall.reval("dyn.load('$SLOPECodeInPathLoaded/cproxSortedL1.so')")
    RCall.reval("source('$SLOPECodeInPathLoaded/proxSortedL1.R')")
    RCall.reval("source('$SLOPECodeInPathLoaded/Adlas.R')")
    res = RCall.reval("Adlas")
    typeof(res) == RCall.RObject{RCall.ClosSxp}
end


text:
struct PeriodicTable
    links::Dict{Element, String}
end

function asa(::Type{DataFrame}, pt::PeriodicTable)::DataFrame
    df = DataFrame(Markdown.MD,10,18)
    for y in 1:10, x in 1:18
        df[y,x]=md""
    end
    for z in elements
        link = get(link,z,nothing)
        df[z.ypos, z.xpos] = isnothing(link) ? Markdown.parse("$(z.symbol)") : Markdown.parse("[$(z.symbol)]($link)")
    end
    return df
end

ext: jl


text:
# -------------------------------------------------
#
# An explicit tensor constructor of the coefficients
# for \prod_i f_i (x_i) given the coefficients of
# each f(x_i).
#
# -------------------------------------------------

# This is a useful way to reduce a difficult O(N log^(d-1) N)
# space of integrations down to O(N d)
# for specifying coefficients of initial conditions

# Efficiency criticality: MEDIUM

# Accuracy criticality: LOW
# Most of the accuracy dependent on functions
# called from other scripts

function tensor_construct{T<:Real}(D::Int, k::Int, n::Int,
						  coeff_array::Array{Dict{CartesianIndex{1}, Array{Array{T,1},1}}, 1};
						  scheme="sparse")
	cutoff = get_cutoff(scheme, D, n)
	coeffs = Dict{CartesianIndex{D}, Array{Array{T,D},D}}()
	# nD_coeffs = zeros(T, D)
	ls::NTuple{D, Int}    = ntuple(i-> (n+1), D)
	modes::NTuple{D, Int} = ntuple(i-> k, D)

	for level in CartesianRange(ls)
		cutoff(level) && continue

		cells::NTuple{D, Int} = ntuple(i -> 1<<max(0, level[i]-2), D)
		level_coeffs = Array{Array{T,D}}(cells)
		for cell in CartesianRange(cells)
			cell_coeffs = Array{T}(modes)
			for mode in CartesianRange(modes)
				val = one(T)
				for d in 1:D
					coeff = coeff_array[d]
					tup::NTuple{1, Int} = (level[d],) # This is somehow still slow
					l = CartesianIndex{1}(tup);
					c = cell[d];
					m = mode[d]
					val *= coeff[l][c][m]
				end
				cell_coeffs[mode] = val
			end
			level_coeffs[cell] = cell_coeffs
		end
		coeffs[level] = level_coeffs
	end
	return coeffs
end


text:
import Statistics

"""
    mean_square_error(ytrue, ypred)
"""
function mean_square_error(
        ytrue::AbstractVector{<:Real},
        ypred::AbstractVector{<:Real},
        )
    result = Statistics.mean(abs2, ytrue .- ypred)
    return result
end

"""
    root_mean_square_error(ytrue, ypred)
"""
root_mean_square_error(ytrue,ypred) = sqrt(mean_square_error(ytrue,ypred))

# convenience aliases for mean squared error:
mean_squared_error(ytrue,ypred) = mean_square_error(ytrue,ypred)
mean_square_deviation(ytrue,ypred) = mean_square_error(ytrue,ypred)
mean_squared_deviation(ytrue,ypred) = mean_square_error(ytrue,ypred)

# convenience aliases for root mean squared error:
root_mean_squared_error(ytrue,ypred) = root_mean_square_error(ytrue,ypred)
root_mean_square_deviation(ytrue,ypred) = root_mean_square_error(ytrue,ypred)
root_mean_squared_deviation(ytrue,ypred) = root_mean_square_error(ytrue,ypred)


text:
struct K2ABC{MMD,K,P,Pa} <: AbstractABC
    kernel::K
    epsilon::Float64
    prior::P
    K::Int
    weights::Vector{Float64}
    particles::Pa
    function K2ABC{MMD}(; kernel, epsilon, prior, K) where MMD
        particles = init_particles(prior, K)
        new{MMD,typeof(kernel),typeof(prior),typeof(particles)}(kernel,
                                                                epsilon,
                                                                prior, K,
                                                                ones(K)./K,
                                                                particles)
    end
end
"""
    K2ABC(; kernel, epsilon, prior, K)

Creates a `K2ABC` structure.
# Example
```julia
using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = K2ABC(kernel = Kernel(),
          epsilon = .1,
          prior = MultivariateNormal([0.], [5.]),
          K = 10^3)
result = run!(k, model, data)
mean(k)
```

[Park, M., Jitkrittum, W., and Sejdinovic, D. (2016),
K2-ABC: Approximate Bayesian Computation with Kernel Embeddings,
Proceedings of Machine Learning Research, 51:398-407
](http://proceedings.mlr.press/v51/park16.html)
"""
K2ABC(; kernel, epsilon, prior, K) = K2ABC{StandardMMD}(kernel = kernel, epsilon = epsilon, prior = prior, K = K)

function run!(rng::Random.AbstractRNG, k::K2ABC{MMD}, model, data;
              verbose = true, maxfevals = Inf, callback = () -> nothing) where MMD
    update!(k.kernel, data)
    for i in 1:k.K
        rand!(rng, k.prior, k.particles[i])
        k.weights[i] = exp(-mmd(MMD, model(k.particles[i]), data)/k.epsilon, k.kernel)
    end
    k.weights ./= sum(k.weights)
    (weights = k.weights, particles = k.particles)
end

struct KernelABC{Ke,P,Pa} <: AbstractABC
    kernel::Ke
    prior::P
    K::Int
    delta::Float64
    weights::Vector{Float64}
    particles::Pa
end
"""
    KernelABC(; kernel, prior, K, delta)

Creates a `KernelABC` structure.
# Example
```julia
using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = KernelABC(kernel = Kernel(),
              delta = 1e-8,
              prior = MultivariateNormal([0.], [5.]),
              K = 10^3)
result = run!(k, model, data)
mean(k)
```

Fukumizu, K., Song, L. and Gretton, A. (2013),
Kernel Bayes' Rule: Bayesian Inference with Positive Definite Kernels,
Journal of Machine Learning Research, 14:3753-3783,
http://jmlr.org/papers/v14/fukumizu13a.html

See also
[Muandet, K., Fukumizu, K., Sriperumbudur, B. and Schölkopf, B. (2017),
Kernel Mean Embedding of Distributions: A Review and Beyond,
Foundations and Trends® in Machine Learning, 10:1–141
](http://dx.doi.org/10.1561/2200000060)
"""
function KernelABC(; kernel, prior, K, delta)
    KernelABC(kernel, prior, K, delta, ones(K)./K, init_particles(prior, K))
end
function kernelabc!(k, model, data; updatekernel = true)
    y = model.(k.particles)
    updatekernel && update!(k.kernel, y, data)
    kystar = [k.kernel(yi, data) for yi in y]
    G = zeros(length(y), length(y))
    @simd for i in eachindex(y)
        for j in i:length(y)
            @inbounds G[i, j] = G[j, i] = k.kernel(y[i], y[j])
        end
    end
    G = Symmetric(G + k.K*k.delta*I)
    k.weights .= clamp.(G \ kystar, 0, Inf)
end
function run!(rng::Random.AbstractRNG, k::KernelABC, model, data;
              verbose = true, maxfevals = Inf, callback = () -> nothing)
    for i in 1:k.K
        rand!(rng, k.prior, k.particles[i])
    end
    kernelabc!(k, model, data)
    (weights = k.weights, particles = k.particles)
end


text:
using SimpleDifferentialOperators
using Test, LinearAlgebra, PATHSolver, Suppressor, Parameters, Random
using DualNumbers

@elapsed begin
    @time @testset "Operators with boundary conditions" begin include("operators.jl") end
    @time @testset "Operators without boundary conditions" begin include("operators-without-bc.jl") end
    @time @testset "Linear Complementarity Problems" begin include("lcp.jl") end
    @time @testset "Boundary extrapolation" begin include("utilities/extrapolatetoboundary.jl") end
end


text:
export OpenStreetMapSpace, OSMSpace, OSM

"""
    OSM
Submodule for functionality related to `OpenStreetMapSpace`.
See the docstring of the space for more info.
"""
module OSM # OpenStreetMap
using Agents
using LightOSM
using Graphs
using Pkg.Artifacts
using LinearAlgebra: dot, norm
using DataStructures
using Downloads

export test_map,
    random_road_position,
    plan_route!,
    distance,
    road_length,
    random_route!,
    latlon,
    intersection,
    road

# Stores information about an agent's path
struct OpenStreetMapPath
    route::Vector{Int}      # node IDs along path from `start` to `dest`
    start::Tuple{Int,Int,Float64} # Initial position of the agent
    dest::Tuple{Int,Int,Float64}    # Destination. `dest[1] == dest[2]` if this is an intersection
    return_route::Vector{Int}   # node IDs along path from `dest` to `start`
    has_to_return::Bool
end

# NOTE: All positions are indexed by vertex number and _not_ node id

"""
    OpenStreetMapSpace(path::AbstractString; kwargs...)
Create a space residing on the Open Street Map (OSM) file provided via `path`.
A sample file is provided using [`OSM.test_map`](@ref). Additional maps can
be downloaded using the [functions provided by LightOSM.jl](https://deloittedigitalapac.github.io/LightOSM.jl/docs/download_network/).
The functionality related to Open Street Map spaces is in the submodule `OSM`.

This space represents the underlying map as a *continuous* entity choosing accuracy over
performance. The map is represented as a graph, consisting of nodes connected by edges. Nodes
are not necessarily intersections, and there may be multiple nodes on a road joining two
intersections. The length of an edge between two nodes is specified in the units of the
map's `weight_type` as listed in the documentation for
[`LightOSM.OSMGraph`](https://deloittedigitalapac.github.io/LightOSM.jl/docs/types/#LightOSM.OSMGraph).
The possible `weight_type`s are:
- `:distance`: The distance in kilometers of an edge
- `:time`: The time in hours to travel along an edge
- `:lane_efficiency`: Time scaled by number of lanes

An example of its usage can be found in [Zombie Outbreak](@ref).

Much of the functionality of this space is provided by interfacing with
[LightOSM.jl](https://github.com/DeloitteDigitalAPAC/LightOSM.jl).

For details on how to obtain an OSM file for your use case, consult the LightOSM.jl documentation.
We provide a function `OSM.test_map` to use for testing.

All keywords are passed to
[`LightOSM.graph_from_file`](https://deloittedigitalapac.github.io/LightOSM.jl/docs/create_graph/#LightOSM.graph_from_file).

## The OSMAgent

The base properties for an agent residing on an `OSMSpace` are as follows:
```julia
mutable struct OSMAgent <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int,Float64}
end
```

Current `pos`ition tuple is represented as
`(start intersection index, finish intersection index, distance travelled)`.
The distance travelled is in the units of `weight_type`.

Further details can be found in [`OSMAgent`](@ref).

## Routing

There are two ways to generate a route, depending on the situation.
1. Use [`plan_route!`](@ref) to plan a route from an agent's current position to a target
   destination. This also has the option of planning a return trip.
2. [`random_route!`](@ref), choses a new random destination and plans a path to it.

Both of these functions override any pre-existing route that may exist for an agent.
"""
struct OpenStreetMapSpace <: Agents.AbstractSpace
    map::OSMGraph
    s::Vector{Vector{Int}}
    routes::Dict{Int,OpenStreetMapPath} # maps agent ID to corresponding path
end

function OpenStreetMapSpace(
    path::AbstractString;
    kwargs...
)
    m = graph_from_file(path; kwargs...)
    agent_positions = [Int[] for _ in 1:Agents.nv(m.graph)]
    return OpenStreetMapSpace(m, agent_positions, Dict())
end

function Base.show(io::IO, s::OpenStreetMapSpace)
    print(
        io,
        "OpenStreetMapSpace with $(length(s.map.highways)) ways " *
        "and $(length(s.map.nodes)) nodes",
    )
end

"""
    OSM.test_map()

Download a small test map of [Göttingen](https://nominatim.openstreetmap.org/ui/details.html?osmtype=R&osmid=191361&class=boundary)
as an artifact. Return a path to the downloaded file.
"""
function test_map()
    artifact_toml = joinpath(@__DIR__, "../../Artifacts.toml")
    map_hash = artifact_hash("osm_map_gottingen", artifact_toml)
    if isnothing(map_hash) || !artifact_exists(map_hash)
        map_hash = create_artifact() do artifact_dir
            Downloads.download(
                "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/artifacts/agents/osm_map_gottingen.json",
                joinpath(artifact_dir, "osm_map_gottingen.json")
            )
        end

        bind_artifact!(artifact_toml, "osm_map_gottingen", map_hash; force = true)
    end

    return joinpath(artifact_path(map_hash), "osm_map_gottingen.json")
end

#######################################################################################
# Custom functions for OSMSpace
#######################################################################################

# EXPORTED
"""
    OSM.random_road_position(model::ABM{<:OpenStreetMapSpace})

Similar to [`random_position`](@ref), but rather than providing only intersections, this method
returns a location somewhere on a road heading in a random direction.
"""
function random_road_position(model::ABM{<:OpenStreetMapSpace})
    # pick a random source and destination, and then a random distance on that edge
    s = Int(rand(model.rng, 1:Agents.nv(model)))
    if isempty(all_neighbors(model.space.map.graph, s))
        return (s, s, 0.0)
    end
    d = Int(rand(model.rng, all_neighbors(model.space.map.graph, s)))
    dist = rand(model.rng) * road_length(s, d, model)
    return (s, d, dist)
end

"""
    OSM.random_route!(agent, model::ABM{<:OpenStreetMapSpace}; kwargs...)

Plan a new random route for the agent, by selecting a random destination and
planning a route from the agent's current position. Overwrite any existing route.

The keyword `limit = 10` specifies the limit on the number of attempts at planning
a random route. Returns `true` if a route was successfully planned, `false` otherwise.
All other keywords are passed to [`plan_route!`](@ref)
"""
function random_route!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A};
    return_trip = false,
    limit = 10,
    kwargs...
) where {A<:AbstractAgent}
    tries = 0
    while tries < limit && !plan_route!(
        agent,
        random_road_position(model),
        model;
        return_trip,
        kwargs...
    )
        tries += 1
    end

    return tries < limit
end

"""
    plan_route!(agent, dest, model::ABM{<:OpenStreetMapSpace};
                   return_trip = false, kwargs...)

Plan a route from the current position of `agent` to the location specified in `dest`, which
can be an intersection or a point on a road.

If `return_trip = true`, a route will be planned from start ⟶ finish ⟶ start. All other
keywords are passed to
[`LightOSM.shortest_path`](https://deloittedigitalapac.github.io/LightOSM.jl/docs/shortest_path/#LightOSM.shortest_path).

Returns `true` if a path to `dest` exists, and `false` if it doesn't. Specifying
`return_trip = true` also requires the existence of a return path for a route to be
planned.
"""
function Agents.plan_route!(
    agent::A,
    dest::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace,A};
    return_trip = false,
    kwargs...
) where {A<:AbstractAgent}
    if agent.pos[1] == agent.pos[2] == dest[1] == dest[2] ||    # identical start and end
       agent.pos == dest ||
       agent.pos == get_reverse_direction(dest, model)

        return true
    end

    if agent.pos[1:2] == dest[1:2] || agent.pos[1:2] == dest[2:-1:1]  # start and end on same road
        if agent.pos[1:2] == dest[2:-1:1]
            dest = get_reverse_direction(dest, model)
        end
        if agent.pos[3] < dest[3]   # same direction
            model.space.routes[agent.id] = OpenStreetMapPath(
                Int[],
                agent.pos,
                dest,
                Int[],
                return_trip
            )
        else    # opposite direction
            move_agent!(agent, get_reverse_direction(agent.pos, model), model)
            model.space.routes[agent.id] = OpenStreetMapPath(
                Int[],
                agent.pos,
                get_reverse_direction(dest, model),
                Int[],
                return_trip,
            )
        end
        return true
    end

    start_node = if agent.pos[1] == agent.pos[2] || 2.0 * agent.pos[3] < road_length(agent.pos, model)
        agent.pos[1]
    else
        agent.pos[2]
    end
    end_node = if dest[1] == dest[2] || 2.0 * dest[3] < road_length(dest, model)
        dest[1]
    else
        dest[2]
    end

    if start_node == end_node   # LightOSM.shortest_path fails in this case
        # either one of start or end is a node and the other an edge incident on it
        if agent.pos[1] == agent.pos[2]
            if dest[2] == agent.pos[1]
                dest = get_reverse_direction(dest, model)
            end
            move_agent!(agent, (dest[1], dest[2], 0.0), model)
            model.space.routes[agent.id] = OpenStreetMapPath(
                Int[],
                agent.pos,
                dest,
                return_trip ? Int[dest[1]] : Int[],
                return_trip,
            )
        elseif dest[1] == dest[2]
            if agent.pos[1] == dest[1]
                move_agent!(agent, get_reverse_direction(agent.pos, model), model)
            end
            model.space.routes[agent.id] = OpenStreetMapPath(
                Int[dest[1]],
                agent.pos,
                dest,
                Int[],
                return_trip,
            )
            # or both are edges incident on a common node
        else
            # swap around directions so that agent is moving toward common node
            # and destination is in the direction from common node to other node
            if agent.pos[1] == dest[2] || agent.pos[2] == dest[2]
                dest = get_reverse_direction(dest, model)
            end
            if agent.pos[1] == dest[1]
                move_agent!(agent, get_reverse_direction(agent.pos, model), model)
            end
            model.space.routes[agent.id] = OpenStreetMapPath(
                Int[dest[1]],
                agent.pos,
                dest,
                return_trip ? Int[dest[1]] : Int[],
                return_trip,
            )
        end
        return true
    end
    route = Int[]

    try
        route = shortest_path(
            model.space.map,
            model.space.map.index_to_node[start_node],
            model.space.map.index_to_node[end_node];
            kwargs...
        )
    catch
        return false
    end

    # convert back to graph indices
    for i in 1:length(route)
        route[i] = Int(model.space.map.node_to_index[route[i]])
    end

    # route won't be empty, those cases are already handled
    reverse!(route)

    # starting from this intersection, so remove it from route
    if agent.pos[1] == agent.pos[2]
        pop!(route)
        # move in reverse direction
    elseif length(route) > 1 && agent.pos[1] == route[end] && agent.pos[2] != route[end-1]
        move_agent!(agent, get_reverse_direction(agent.pos, model), model)
    end

    return_route = Int[]
    if return_trip
        try
            return_route =
                shortest_path(
                    model.space.map,
                    model.space.map.index_to_node[end_node],
                    model.space.map.index_to_node[start_node];
                    kwargs...
                )
        catch
            return false
        end
    end

    if return_trip
        # convert back to graph indices
        for i in 1:length(return_route)
            return_route[i] = Int(model.space.map.node_to_index[return_route[i]])
        end

        reverse!(return_route)
        # analogous case to forward route
        dest[1] == dest[2] && pop!(return_route)
        # will not check other case since dest shouldn't be flipped until
        # actually doing return trip
    end

    model.space.routes[agent.id] = OpenStreetMapPath(
        route,
        agent.pos,
        dest,
        return_route,
        return_trip,
    )
    return true
end

# Allows passing destination as an index
Agents.plan_route!(agent::A, dest::Int, model; kwargs...) where {A<:AbstractAgent} =
    plan_route!(agent, (dest, dest, 0.0), model; kwargs...)

"""
    OSM.distance(pos_1, pos_2, model::ABM{<:OpenStreetMapSpace}; kwargs...)

Return the distance between the two positions along the shortest path joining them in the given
model. Returns `Inf` if no such path exists.

All keywords are passed to
[`LightOSM.shortest_path`](https://deloittedigitalapac.github.io/LightOSM.jl/docs/shortest_path/#LightOSM.shortest_path).
"""
function distance(
    pos_1::Tuple{Int,Int,Float64},
    pos_2::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace};
    kwargs...
)
    # positions are identical
    if pos_1[1] == pos_1[2] == pos_2[1] == pos_2[2] ||
       pos_1 == pos_2 ||
       pos_1 == get_reverse_direction(pos_2, model)
        return 0.0
    end

    # positions on same road
    if pos_1[1:2] == pos_2[1:2]
        return abs(pos_1[3] - pos_2[3])
    elseif pos_1[1:2] == pos_2[2:-1:1]
        return abs(pos_1[3] - road_length(pos_1, model) + pos_2[3])
    end

    # starting vertex
    st_node = if pos_1[1] == pos_1[2] || pos_1[3] < road_length(pos_1, model) / 2
        pos_1[1]
    else
        pos_1[2]
    end

    # ending vertex
    en_node = if pos_2[1] == pos_2[2] || pos_2[3] < road_length(pos_2, model) / 2
        pos_2[1]
    else
        pos_2[2]
    end

    # Case where they are same
    if st_node == en_node
        r1 = (pos_1[1] == pos_1[2] ? 0.0 : pos_1[3])
        r2 = (pos_2[1] == pos_2[2] ? 0.0 : pos_2[3])
        if pos_1[1] != pos_1[2] && st_node == pos_1[2]
            r1 = road_length(pos_1, model) - r1
        end
        if pos_2[1] != pos_2[2] && en_node == pos_2[2]
            r2 = road_length(pos_2, model) - r2
        end
        return r1 + r2
    end

    route = Int[]

    # get route
    try
        route = shortest_path(
            model.space.map,
            model.space.map.index_to_node[st_node],
            model.space.map.index_to_node[en_node];
            kwargs...
        )
    catch
        return Inf
    end

    # distance along route
    dist = sum(weights_from_path(model.space.map, route))

    # cases where starting or ending position is partway along a road
    # route may or may not pass through that road, so all cases need to be handled
    if pos_1[1] != pos_1[2]
        if route[1] == pos_1[1]
            if route[2] == pos_1[2]
                dist -= pos_1[3]
            else
                dist += pos_1[3]
            end
        else
            if route[2] == pos_1[1]
                dist -= road_length(pos_1, model) - pos_1[3]
            else
                dist += road_length(pos_1, model) - pos_1[3]
            end
        end
    end

    if pos_2[1] != pos_2[2]
        if route[end] == pos_2[1]
            if route[end-1] == pos_2[2]
                dist -= pos_2[3]
            else
                dist += pos_2[3]
            end
        else
            if route[end-1] == pos_2[1]
                dist -= road_length(pos_2, model) - pos_2[3]
            else
                dist += road_length(pos_2, model) - pos_2[3]
            end
        end
    end

    return dist
end

function distance(
    pos_1::Int,
    pos_2::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace};
)
    distance((pos_1, pos_1, 0.0), pos_2, model)
end

function distance(
    pos_1,
    pos_2::Int,
    model::ABM{<:OpenStreetMapSpace}
)
    distance(pos_1, (pos_2, pos_2, 0.0), model)
end

"""
    OSM.latlon(pos, model)
    OSM.latlon(agent, model)

Return `(latitude, longitude)` of current road or intersection position.
"""
latlon(pos::Int, model::ABM{<:OpenStreetMapSpace}) =
    Tuple(model.space.map.node_coordinates[pos])

function latlon(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace})
    # extra checks to ensure consistency between both versions of `latlon`
    if pos[3] == 0.0 || pos[1] == pos[2]
        return latlon(pos[1], model)
    elseif pos[3] == road_length(pos, model)
        return latlon(pos[2], model)
    else
        gloc1 = get_geoloc(pos[1], model)
        gloc2 = get_geoloc(pos[2], model)
        dist = norm(LightOSM.to_cartesian(gloc1) .- LightOSM.to_cartesian(gloc2))
        dir = heading(gloc1, gloc2)
        geoloc = calculate_location(gloc1, dir, pos[3] / road_length(pos, model) * dist)
        return (geoloc.lat, geoloc.lon)
    end
end

latlon(agent::A, model::ABM{<:OpenStreetMapSpace,A}) where {A<:AbstractAgent} =
    latlon(agent.pos, model)

"""
    OSM.intersection(latlon::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})

Return the nearest intersection position to (latitude, longitude).
Quicker, but less precise than [`OSM.road`](@ref).
"""
function intersection(ll::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})
    vert = Int(model.space.map.node_to_index[nearest_node(model.space.map, [GeoLocation(ll..., 0.0)])[1][1][1]])
    return (vert, vert, 0.0)
end

"""
    OSM.road(latlon::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})

Return a location on a road nearest to (latitude, longitude). Significantly slower, but more
precise than [`OSM.intersection`](@ref).
"""
function road(ll::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})
    best_sq_dist = Inf
    best = (-1, -1, -1.0)
    pt = LightOSM.to_cartesian(GeoLocation(ll..., 0.0))
    for e in edges(model.space.map.graph)
        s = LightOSM.to_cartesian(GeoLocation(model.space.map.node_coordinates[src(e)]..., 0.0))
        d = LightOSM.to_cartesian(GeoLocation(model.space.map.node_coordinates[dst(e)]..., 0.0))
        road_vec = d .- s

        # closest point on line segment requires checking if perpendicular from point lies on line
        # segment. If not, use the closest end of the line segment
        if dot(pt .- s, road_vec) < 0.0
            int_pt = s
        elseif dot(pt .- d, road_vec) > 0.0
            int_pt = d
        else
            int_pt = s .+ (dot(pt .- s, road_vec) / dot(road_vec, road_vec)) .* road_vec
        end

        sq_dist = dot(int_pt .- pt, int_pt .- pt)

        if sq_dist < best_sq_dist
            best_sq_dist = sq_dist
            rd_dist = norm(int_pt .- s) / norm(road_vec) * road_length(Int(src(e)), Int(dst(e)), model)

            best = (
                Int(src(e)),
                Int(dst(e)),
                rd_dist,
            )
        end
    end

    return best
end

"""
    OSM.road_length(start::Int, finish::Int, model)
    OSM.road_length(pos::Tuple{Int,Int,Float64}, model)

Return the road length between two intersections. This takes into account the
direction of the road, so `OSM.road_length(pos_1, pos_2, model)` may not be the
same as `OSM.road_length(pos_2, pos_1, mode)`. Units of the returned quantity
are as specified by the underlying graph's `weight_type`.
"""
road_length(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace}) =
    road_length(pos[1], pos[2], model)
function road_length(p1::Int, p2::Int, model::ABM{<:OpenStreetMapSpace})
    len = model.space.map.weights[p1, p2]
    if len == 0.0 || len == Inf
        len = model.space.map.weights[p2, p1]
    end
    return len
end

function Agents.is_stationary(agent, model::ABM{<:OpenStreetMapSpace})
    return !haskey(model.space.routes, agent.id)
end

"""
    OSM.get_geoloc(pos::Int, model::ABM{<:OpenStreetMapSpace})

Return `GeoLocation` corresponding to node `pos`
"""
get_geoloc(pos::Int, model::ABM{<:OpenStreetMapSpace}) = GeoLocation(model.space.map.node_coordinates[pos]..., 0.0)

"""
    OSM.get_reverse_direction(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace})

Returns the same position, but with `pos[1]` and `pos[2]` swapped and `pos[3]` updated accordingly
"""
get_reverse_direction(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace}) =
    (pos[2], pos[1], road_length(pos, model) - pos[3])

#######################################################################################
# Agents.jl space API
#######################################################################################

function Agents.random_position(model::ABM{<:OpenStreetMapSpace})
    vert = Int(rand(model.rng, 1:Agents.nv(model)))
    return (vert, vert, 0.0)
end

function Agents.add_agent_to_space!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
) where {A<:AbstractAgent}
    push!(model.space.s[agent.pos[1]], agent.id)
    return agent
end

function Agents.remove_agent_from_space!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
) where {A<:AbstractAgent}
    prev = model.space.s[agent.pos[1]]
    ai = findfirst(i -> i == agent.id, prev)
    deleteat!(prev, ai)
    return agent
end

function Agents.move_agent!(
    agent::A,
    pos::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace,A},
) where {A<:AbstractAgent}
    if pos[1] == agent.pos[1]
        agent.pos = pos
        return agent
    end
    Agents.remove_agent_from_space!(agent, model)
    agent.pos = pos
    Agents.add_agent_to_space!(agent, model)
end

"""
    move_along_route!(agent, model::ABM{<:OpenStreetMapSpace}, distance::Real)

Move an agent by `distance` along its planned route. Units of distance are as specified
by the underlying graph's weight_type.
"""
function Agents.move_along_route!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
    distance::Real,
) where {A<:AbstractAgent}
    if is_stationary(agent, model)
        return nothing
    end

    # branching here corresponds to nesting of the following cases:
    # - Is the agent moving to the end of the road or somewhere in the middle (isempty(osmpath.route))
    # - Is `distance` more than the distance to the next waypoint / end of route (will it overshoot)?
    #   - Did we reach the end of the route? If so, does the agent need to return?
    #     - If the agent returns, what direction does it go in?
    #   - If we have another waypoint after this, move there
    # It might be easier to formulate this as a recursive structure, so recursive calls are annotated where necessary
    # instead of the loop this currently runs in. These annotations are marked with `##` just to make it clear.

    osmpath = model.space.routes[agent.id]
    while distance > 0
        if isempty(osmpath.route)   # last leg of route, to the middle of a road
            # distance left to reach dest
            dist_left = osmpath.dest[3] - agent.pos[3]
            if dist_left < distance # can overshoot destination
                distance -= dist_left
                move_agent!(agent, osmpath.dest, model) # reach destination
                if osmpath.has_to_return    # if we have to return
                    # empty return route implies we just have to go reverse on this edge
                    if isempty(osmpath.return_route)
                        osmpath = OpenStreetMapPath(
                            osmpath.return_route,
                            get_reverse_direction(osmpath.dest, model),
                            get_reverse_direction(osmpath.start, model),
                            Int[],
                            false,
                        )
                        move_agent!(agent, get_reverse_direction(agent.pos, model), model)
                        break
                        ## return
                    end

                    # non-empty return route
                    osmpath = OpenStreetMapPath(    # construct return path
                        osmpath.return_route,
                        osmpath.dest,
                        osmpath.start,
                        Int[],
                        false
                    )
                    # get next waypoint on return path
                    # this will either one of the endpoints of this road
                    next_wp = osmpath.route[end]
                    if next_wp == agent.pos[1]
                        # need to go back along same road, so reverse direction
                        move_agent!(agent, get_reverse_direction(agent.pos, model), model)
                    end
                    # move remaining distance along return path
                    continue
                    ## return Agents.move_along_route!(agent, model, distance)
                end
                # delete path data so this agent is recognized as stationary
                delete!(model.space.routes, agent.id)
                break
                ## return
            end
            # can't directly reach destination, so move distance along road
            # ensure we don't overshoot the destination
            result_pos = min(agent.pos[3] + distance, osmpath.dest[3])
            move_agent!(agent, (agent.pos[1:2]..., result_pos), model)
            break
            ## return
        end
        # don't need an else clause, since everything going inside the if clause will break
        # or continue

        # distance left till next node
        dist_left = road_length(agent.pos, model) - agent.pos[3]
        if dist_left < distance # can overshoot
            distance -= dist_left   # leftover distance
            node_a = pop!(osmpath.route)    # remove the node we just reached
            if isempty(osmpath.route)
                # this was the last node, so either we reached the end or dest is on an outgoing road
                if osmpath.dest[1] == osmpath.dest[2]   # we reached the end
                    if osmpath.has_to_return    # need to return from here
                        # empty return route, so reverse along same edge
                        if isempty(osmpath.return_route)
                            move_agent!(agent, (dest[1], agent.pos[1], 0.0), model)
                            osmpath = OpenStreetMapPath(
                                osmpath.return_route,
                                osmpath.dest,
                                get_reverse_direction(osmpath.start, model),
                                Int[],
                                false,
                            )
                            break
                            ## return
                        end
                        osmpath = OpenStreetMapPath(
                            osmpath.return_route,
                            osmpath.dest,
                            osmpath.start,
                            Int[],
                            false,
                        )
                        # next node on return path
                        node_b = osmpath.route[end]
                        # set to move along reverse path
                        move_agent!(agent, (node_a, node_b, 0.0), model)
                        # move rest of distance along return route
                        ## return Agents.move_along_route!(agent, model, distance)
                        continue
                    end
                    # move to end
                    move_agent!(agent, osmpath.dest, model)
                    # remove route so agent is marked as stationary
                    delete!(model.space.routes, agent.id)
                    ## return
                    break
                end

                # destination is on an outgoing road from this last waypoint
                # move to beginning of this road
                move_agent!(agent, (osmpath.dest[1:2]..., 0.0), model)
                # move rest of distance to destination
                ## return Agents.move_along_route!(agent, model, distance)
                continue
            end

            # there is a further waypoint
            node_b = osmpath.route[end]
            move_agent!(agent, (node_a, node_b, 0.0), model)
            # move rest of distance to next waypoint
            ## return Agents.move_along_route!(agent, model, distance)
            continue
        end

        # will not overshoot
        result_pos = min(agent.pos[3] + distance, road_length(agent.pos, model))
        move_agent!(agent, (agent.pos[1:2]..., result_pos), model)
        ## return
        break
    end
end

# Nearby positions must be intersections, since edges imply a direction.
# nearby agents/ids can be on an intersection or on a road X m away.
# We cannot simply use nearby_positions to collect nearby ids then.

# Default is searching both backwards and forwards.
# I assume it would be useful to turn off one or the other at some point.
function Agents.nearby_ids(
    pos::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace},
    distance::Real,
    args...;
    kwargs...
)
    distances = Dict{Int,Float64}()
    queue = Queue{Int}()
    nearby = Int[]

    if pos[1] == pos[2]
        distances[pos[1]] = 0.0
        enqueue!(queue, pos[1])
    else
        dist_1 = pos[3]
        dist_2 = road_length(pos, model) - pos[3]

        if dist_1 <= distance && dist_2 <= distance
            # just add to queue, all IDs on this road will be covered
            distances[pos[1]] = dist_1
            enqueue!(queue, pos[1])
            distances[pos[2]] = dist_2
            enqueue!(queue, pos[2])
        elseif dist_1 <= distance   # && dist_2 > distance
            # BFS covers IDs `distance` away from `pos[1]`, but not those in range
            # `(distance, dist_1 + distance)` away
            distances[pos[2]] = dist_2
            enqueue!(queue, pos[2])
            # push IDs that won't be covered
            for id in forward_ids_on_road(pos[1], pos[2], model)
                dist_1 < model[id].pos[3] <= dist_1 + distance && push!(nearby, id)
            end
            for id in reverse_ids_on_road(pos[1], pos[2], model)
                dist = road_length(pos, model) - model[id].pos[3]
                dist_1 < dist <= dist_1 + distance && push!(nearby, id)
            end
        elseif dist_2 <= distance # && dist_1 > distance
            # BFS covers IDs `distance` away from `pos[2]`, but not those in range
            # `(distance, dist_2 + distance)` away
            distances[pos[2]] = dist_2
            enqueue!(queue, pos[2])
            # push IDs that won't be covered
            for id in forward_ids_on_road(pos[1], pos[2], model)
                dist = road_length(pos, model) - model[id].pos[3]
                dist_2 < dist <= dist_2 + distance && push!(nearby, id)
            end
            for id in reverse_ids_on_road(pos[1], pos[2], model)
                dist_2 < model[id].pos[3] < dist_2 + distance && push!(nearby, id)
            end
        else    # neither node is within `distance` of `pos`, so simply filter IDs on this road
            for id in forward_ids_on_road(pos[1], pos[2], model)
                abs(model[id].pos[3] - dist_1) <= distance && push!(nearby, id)
            end
            for id in reverse_ids_on_road(pos[1], pos[2], model)
                abs(model[id].pos[3] - dist_2) <= distance && push!(nearby, id)
            end
        end
    end

    # NOTE: During BFS, each node is only explored once. From each node, every outgoing and incoming edge is
    # considered and the node on the other end is added to the queue if it is close enough. The edge is explored
    # for IDs by the node with the lower index (to prevent the same edge from being explored twice).
    # If the node on the other end can't be reached, then the current node explores however far it can on
    # this road. If the other node is reached through another path, it takes this into account and only explores
    # the unexplored part of this edge.
    while !isempty(queue)
        node = dequeue!(queue)
        for nb in nearby_positions(node, model; neighbor_type = :all)
            rd_len = road_length(node, nb, model)
            if rd_len == 0.0 || rd_len == Inf
                rd_len = road_length(nb, node, model)
            end

            if rd_len <= distance - distances[node]  # can reach pos_2 from this node
                if !haskey(distances, nb) # mark for exploration if not visited
                    distances[nb] = distances[node] + rd_len
                    enqueue!(queue, nb)
                end
                if node < nb # road is explored by smaller indexed node
                    append!(nearby, ids_on_road(node, nb, model))
                end
            else
                # cannot reach pos_2 from this node
                dist_to_explore = distance - distances[node]
                if haskey(distances, nb)
                    # part of it is already explored, so only explore remaining part
                    dist_to_explore = max(min(dist_to_explore, rd_len - distances[nb]), 0.0)
                end
                for id in forward_ids_on_road(node, nb, model)
                    model[id].pos[3] <= dist_to_explore && push!(nearby, id)
                end
                for id in reverse_ids_on_road(node, nb, model)
                    dist = rd_len - model[id].pos[3]
                    dist <= dist_to_explore && push!(nearby, id)
                end
            end
        end
    end

    return nearby
end

forward_ids_on_road(pos_1::Int, pos_2::Int, model::ABM{<:OpenStreetMapSpace}) =
    Iterators.filter(i -> model[i].pos[2] == pos_2, model.space.s[pos_1])

reverse_ids_on_road(pos_1::Int, pos_2::Int, model::ABM{<:OpenStreetMapSpace}) =
    forward_ids_on_road(pos_2, pos_1, model)

ids_on_road(pos_1::Int, pos_2::Int, model::ABM{<:OpenStreetMapSpace}) =
    Iterators.flatten((
        forward_ids_on_road(pos_1, pos_2, model),
        reverse_ids_on_road(pos_1, pos_2, model),
    ))

function Agents.nearby_ids(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
    args...;
    kwargs...
) where {A<:AbstractAgent}
    all = nearby_ids(agent.pos, model, args...; kwargs...)
    filter!(i -> i ≠ agent.id, all)
end

Agents.nearby_positions(pos::Tuple{Int,Int,Float64}, model, args...; kwargs...) =
    nearby_positions(pos[1], model, args...; kwargs...)

function Agents.nearby_positions(
    position::Int,
    model::ABM{<:OpenStreetMapSpace};
    neighbor_type::Symbol = :default
)
    @assert neighbor_type ∈ (:default, :all, :in, :out)
    neighborfn = if neighbor_type == :default
        Graphs.neighbors
    elseif neighbor_type == :in
        Graphs.inneighbors
    elseif neighbor_type == :out
        Graphs.outneighbors
    else
        Graphs.all_neighbors
    end
    Int.(neighborfn(model.space.map.graph, position))
end

end # module OSM

const OpenStreetMapSpace = OSM.OpenStreetMapSpace

const OSMSpace = OSM.OpenStreetMapSpace


text:
function parse_geom(string::AbstractString)
    if isfile(string)
        return read_xyz_file(string)
    else
        return read_xyz_string(string)
    end
end

function read_xyz_line(line::AbstractString)
    s, x, y, z = split(line)[1:4]
    position = parse.(Float64, [x, y, z]) / CODATA08_BOHR_TO_A
    number = findfirst(==(titlecase(s)), ELEMENTS)
    return position, number
end

function read_xyz_file(file::AbstractString)
    io = open(file, "r")
    n = parse(Int, readline(io))
    readline(io)
    positions = Vector{Vec3{Float64}}(undef, n)
    numbers = Vector{Int}(undef, n)
    for i = 1:n
        line = readline(io)
        positions[i], numbers[i] = read_xyz_line(line)
    end
    close(io)
    return positions, numbers
end

function read_xyz_string(string::AbstractString)
    lines = split(strip(string), "\n")
    n = length(lines)
    positions = Vector{Vec3{Float64}}(undef, n)
    numbers = Vector{Int}(undef, n)
    for (i, line) in enumerate(lines)
        positions[i], numbers[i] = read_xyz_line(line)
    end
    return positions, numbers
end

function parse_pc(string::AbstractString)
    if isfile(string)
        return read_pc_file(string)
    else
        return read_pc_string(string)
    end
end

function read_pc_line(line::AbstractString)
    q, x, y, z = split(line)[1:4]
    position = parse.(Float64, [x, y, z]) / CODATA08_BOHR_TO_A
    charge = parse(Float64, q)
    return position, charge
end

function read_pc_file(file::AbstractString)
    io = open(file, "r")
    n = parse(Int, readline(io))
    positions = Vector{Vec3{Float64}}(undef, n)
    charges = Vector{Float64}(undef, n)
    for i = 1:n
        line = readline(io)
        positions[i], charges[i] = read_pc_line(line)
    end
    close(io)
    return positions, charges
end

function read_pc_string(string::AbstractString)
    lines = split(strip(string), "\n")
    n = length(lines)
    positions = Vector{Vec3{Float64}}(undef, n)
    charges = Vector{Float64}(undef, n)
    for (i, line) in enumerate(lines)
        positions[i], charges[i] = read_pc_line(line)
    end
    return positions, charges
end


text:
function Base.read!(f, block::ConstantBlockHeader{T}) where T
    block.val
end

function Base.read!(f, block::PlainMeshBlockHeader{T,D}) where {T,D}
    offset = block.base_header.data_location

    raw_data = ntuple(Val(D)) do i
        Array{T, 1}(undef, block.dims[i])
    end
    @inbounds for i in eachindex(raw_data)
        seek(f, offset)
        offset += sizeof(T) * block.dims[i]
        read!(f, raw_data[i])
    end

    return raw_data
end

function Base.read!(f, block::PointMeshBlockHeader{T,D}) where {T,D}
    offset = block.base_header.data_location

    raw_data = ntuple(Val(D)) do _
        Array{T, 1}(undef, block.np)
    end
    @inbounds for i in eachindex(raw_data)
        seek(f, offset)
        offset += sizeof(T) * block.np
        read!(f, raw_data[i])
    end

    return raw_data
end

function Base.read!(f, block::PlainVariableBlockHeader{T}) where T
    dim = prod(block.dims)
    raw_data = Array{T, 1}(undef, dim)

    seek(f, block.base_header.data_location)
    read!(f, raw_data)

    reshape(raw_data, block.dims)
end

function Base.read!(f, block::PointVariableBlockHeader{T}) where T
    raw_data = Array{T, 1}(undef, block.np)
    seek(f, block.base_header.data_location)

    read!(f, raw_data)
end

function Base.read!(f, block::RunInfoBlockHeader{T,D}) where {T,D}
end


text:
# Copyright (c) 2019 Harshavardhan P K
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

function getcoords(arr::Array{Tuple{Char,Int64}})::Array{Tuple{Int64,Int64}}
    coords::Array{Tuple{Int64,Int64}} = Array([])
    x, y = 0, 0
    for (dir, steps) ∈ arr
        for _ ∈ 1:steps
            if dir == 'R'
                x += 1
            elseif dir == 'L'
                x -= 1
            elseif dir == 'U'
                y += 1
            elseif dir == 'D'
                y -= 1
            end
            push!(coords, (x, y))
        end
    end
    return coords
end

function mindelay(arr1::Array{Tuple{Int64,Int64}}, arr2::Array{Tuple{Int64,Int64}})::Int64
    common_coords = Set(arr1) ∩ Set(arr2)
    delay = indexin(common_coords, arr1) + indexin(common_coords, arr2)
    return minimum(delay)
end

function getData(fname::String)
    f = open(fname)
    line = readline(f)
    line = split(strip(line), ',')
    w1 = map(x -> (x[1], parse(Int64, x[2:end])), line)
    line = readline(f)
    line = split(strip(line), ',')
    w2 = map(x -> (x[1], parse(Int64, x[2:end])), line)
    return w1, w2
end

w1, w2 = getData("input")
c1, c2 = getcoords(w1), getcoords(w2)
println(mindelay(c1, c2))


text:
module FixedTypes
include("create_struct.jl")
@define_struct Param range value unit
@define_struct Out range value unit
@define_struct Compression compr_method detail_th compr_aux
@define_struct Signal y dt t compression
@define_struct MyData param aux comment out wfm
export Param, Out, Signal, Compression, MyData, @define_struct
end


text:
using Distributions
using PyPlot

function rsemicts(n:UInt64, pzero:Real, cts_density:String, cts_params)

end

max_stars_repo_name: psaiku/semictsdistributions

max_issues_repo_name: psaiku/semictsdistributions

max_forks_repo_name: psaiku/semictsdistributions


text:
module ScottishBenefits
#=
# This module is a holder for all the Scotland-specific
# benefits - currently SCP and bedroom tax
# mitigation,
=#
using ScottishTaxBenefitModel
using .ModelHousehold:
    BenefitUnit,
    Household,
    Person,
    count,
    get_head,
    get_spouse,
    le_age

using .STBParameters:
    ScottishChildPayment

using .STBIncomes
using .Definitions

using .Intermediate
using .Results: BenefitUnitResult, HouseholdResult, has_any

export calc_scottish_child_payment!,
    calc_bedroom_tax_mitigation!

function calc_scottish_child_payment!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit        :: BenefitUnit,
    intermed            :: MTIntermediate,
    scpsys              :: ScottishChildPayment )
    scp = 0.0
    bu = benefit_unit
    bur = benefit_unit_result # shortcuts
    nkids = count( bu, le_age, scpsys.maximum_age )
    if( nkids > 0 ) && has_any( bur, scpsys.qualifying_benefits... )
        scp = nkids * scpsys.amount
        spouse = get_spouse( bu )
        target_pid = BigInt(-1)
        if spouse === nothing
            target_pid = get_head( bu ).pid
        else
            target_pid = spouse.pid
        end
        bur.pers[target_pid].income[SCOTTISH_CHILD_PAYMENT] = scp
    end
end

"""

This is the Scottish Government's bedroom tax
removal for social renters.
DISCRESIONARY_HOUSING_PAYMENT can be used for other things, too, but not so far
in this model.

I'm basically taking a wild guess at how this works for UC
since I can't find any documentation on what 'housing element'
means in this context - it can't mean all of it
regardless of total UC entitlement.

Assigns a DISCRESIONARY_HOUSING_PAYMENT of the
minimum of hb/uc hosts and the amount of
rooms reduction to whoever recieves uc/hb in the 1st benefit unit.
'But for weekly or monthly payments, the amount of a discretionary housing payment cannot be more than the amount of universal credit or housing benefit that you get to help with your rent.'

See: https://cpag.org.uk/scotland/welfare-rights/scottish-benefits/discretionary-housing-payments-scotland
TODO parameterise this somehow - mitigate fully or whatever.
"""
function calc_bedroom_tax_mitigation!(
    hr    :: HouseholdResult,
    hh    :: Household )
    # DISCRESIONARY_HOUSING_PAYMENT
    if ! is_social_renter( hh.tenure )
        return
    end
    hrep = hr.bus[1].legacy_mtbens.hb_recipient
    urep = hr.bus[1].uc.recipient
    uche = hr.bus[1].uc.housing_element
    benred = hr.bus[1].bencap.reduction
    hb = 0.0
    if hrep > 0
        hb = hr.bus[1].pers[hrep].income[HOUSING_BENEFIT]
    end
    uc = 0.0
    if urep > 0
        uc = hr.bus[1].pers[urep].income[UNIVERSAL_CREDIT]
    end
    rrd = hr.housing.rooms_rent_reduction
    if rrd > 0
        # the > housing tests here are kinda redundant, but still ..
        if uc > 0
            # !!! FIXME!! Assume uc housing element is the 1st part to be withdrawn
            # with income so eligibiliity for DHP ceases
            # when income > housing element.
            # FIXME I think that's how it must work but I don't know from
            # any documentation I have that this is right
            # uc itself is here in case this hhls is caught
            # by the benefit cap.
            uchousing = max(
                0.0,
                uche - hr.bus[1].uc.total_income - benred )
            hr.bus[1].pers[urep].income[DISCRESIONARY_HOUSING_PAYMENT] =
                min( uchousing, rrd )
        elseif hrep > 0 && hb > 0
            # can't exceed housing benefit
            hr.bus[1].pers[hrep].income[DISCRESIONARY_HOUSING_PAYMENT] =
                min( hb , rrd )
        end
    end
end # calc_ ..

end # module

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using StaticArrays
using Binary_BB_Trees
using Binary_BB_Trees: findMinSVSV, findMaxSVSV

p1 = SVector(1.0, 5.0, 9.0)
p2 = 2 * p1
p3 = 3 * p1
p4 = 4 * p1

@testset "findMinSVSV/findMaxSVSV" begin
    @test p1 == findMinSVSV(p1, p2) == findMinSVSV(p2, p1)
    @test p1 == findMinSVSV(SVector(p1, p2, p3)) == findMinSVSV(SVector(p3, p2, p1)) == findMinSVSV(SVector(p2, p1, p3))
    @test p1 == findMinSVSV(SVector(p1, p2, p3, p4)) == findMinSVSV(SVector(p4, p3, p2, p1)) == findMinSVSV(SVector(p2, p4, p1, p3))

    @test p4 == findMaxSVSV(p4, p3) == findMaxSVSV(p3, p4)
    @test p4 == findMaxSVSV(SVector(p4, p3, p2)) == findMaxSVSV(SVector(p2, p3, p4)) == findMaxSVSV(SVector(p1, p4, p3))
    @test p4 == findMaxSVSV(SVector(p1, p2, p3, p4)) == findMaxSVSV(SVector(p4, p3, p2, p1)) == findMaxSVSV(SVector(p2, p4, p1, p3))
end

@testset "minMaxToCenterExtent" begin
    center, extent = Binary_BB_Trees.minMaxToCenterExtent(p1, p2)
    @test center == (p1 + p2) * 0.5
    @test extent == (p2 - p1) * 0.5
end

@testset "calc_min_max" begin
    min_, max_ = calc_min_max(p1, p2)
    @test min_ == p1
    @test max_ == p2
    min_, max_ = calc_min_max(p2, p1)
    @test min_ == p1
    @test max_ == p2
end


text:
"num_obs <- Internal function used to determine the number of observations in the input dataset"
(num_obs(data::AbstractVector{T})::Int) where {T} = length(data)
(num_obs(data::AbstractMatrix{T})::Int) where {T} = size(data, 1)
(num_obs(data::Vector{Vector{T}})::Int) where {T} = (isempty(data) || any(length.(data) .!= length(data[1]))) ? error("Input dataset is empty, or inner vectors of input dataset do not have matching length: $(length.(data))") : length(data[1])
(num_obs(data::DataFrames.DataFrame)::Int) = size(data, 1)
(num_obs(data::TimeSeries.TimeArray{T,1})::Int) where {T} = size(data, 1)
(num_obs(data::TimeSeries.TimeArray{T,2})::Int) where {T} = size(data, 1)

"num_var <- Internal function used to determine the number of variables in the input dataset"
(num_var(data::AbstractVector{T})::Int) where {T} = 1
(num_var(data::AbstractMatrix{T})::Int) where {T} = size(data, 2)
(num_var(data::Vector{Vector{T}})::Int) where {T} = isempty(data) ? error("Input dataset is empty") : length(data)
(num_var(data::DataFrames.DataFrame)::Int) = size(data, 2)
(num_var(data::TimeSeries.TimeArray{T,1})::Int) where {T} = 1
(num_var(data::TimeSeries.TimeArray{T,2})::Int) where {T} = size(data, 2)

"local_get_var <- Internal function used to get the ith variable in dataset"
(local_get_var(data::AbstractVector{T}, i::Int)::Vector{T}) where {T} = i == 1 ? data[:] : error("Invalid index $(i) given data $(typeof(data))")
(local_get_var(data::AbstractMatrix{T}, i::Int)::Vector{T}) where {T} = (1 <= i <= size(data,2)) ? data[:, i] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(size(data, 2))")
(local_get_var(data::Vector{Vector{T}}, i::Int)::Vector{T}) where {T} = (1 <= i <= length(data)) ? data[i] : error("Invalid index $(i) given data $(typeof(data)) with outer length: $(length(data))")
(local_get_var(data::DataFrames.DataFrame, i::Int)) = (1 <= i <= num_var(data)) ? data[:, i] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")
(local_get_var(data::TimeSeries.TimeArray{T,1}, i::Int)::Vector{T}) where {T} = i == 1 ? values(data) : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")
(local_get_var(data::TimeSeries.TimeArray{T,2}, i::Int)::Vector{T}) where {T} = (1 <= i <= num_var(data)) ? values(data)[:,i] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")

"local_get_index <- Internal function used to resample the dataset data using the input resampling index inds"
(local_get_index(data::AbstractVector{T}, inds::Vector{Int})::Vector{T}) where {T} = data[inds]
(local_get_index(data::AbstractMatrix{T}, inds::Vector{Int})::Matrix{T}) where {T} = data[inds, :]
(local_get_index(data::Vector{Vector{T}}, inds::Vector{Int})::Vector{Vector{T}}) where {T} = [ y[inds] for y in data ]
(local_get_index(data::DataFrames.DataFrame, inds::Vector{Int})::DataFrames.DataFrame) = data[inds, :]
(local_get_index(data::TimeSeries.TimeArray{T,1}, inds::Vector{Int})::TimeSeries.TimeArray{T,1}) where {T} = TimeSeries.TimeArray(TimeSeries.timestamp(data), values(data)[inds], TimeSeries.colnames(data) ; unchecked=true)
(local_get_index(data::TimeSeries.TimeArray{T,2}, inds::Vector{Int})::TimeSeries.TimeArray{T,2}) where {T} = TimeSeries.TimeArray(TimeSeries.timestamp(data), values(data)[inds, :], TimeSeries.colnames(data) ; unchecked=true)

"""
    dbootdata_one(data::T, bi::BootInput)::T
    dbootdata_one(data::T; kwargs...)::T

Get a single resampled dataset of the input data using the dependent boostrap
methodology defined in BootInput.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.

Note, the output type will always be the same as the type of the input data.
"""
(dbootdata_one(data::Td, bi::BootInput)::Td) where {Td} = local_get_index(data, dbootinds_one(bi))
(dbootdata_one(data::Td ; kwargs...)::Td) where {Td} = dbootdata_one(data, BootInput(data ; kwargs...))
# function dbootdata_one_infl(data_infl::Td, bi::BootInput{BootTapered})::TD where {Td}
#     dataout = apply_inds_to_data(data_infl, dbootinds_one(bi))
#     dboot_weight!(dataout, bi)
#     return dataout
# end
# function (dbootdata_one(data::Td, bi::BootInput{BootTapered})::Td) where {Td}
#     return dbootdata_one_infl(apply_influence_function(data), bi, bm)
# end

"""
    dbootdata(data::T , bi::BootInput)::Vector{T}
    dbootdata(data::T ; kwargs...)::Vector{T}

Get the resampled datasets of the input data using the dependent bootstrap
methodology defined in BootInput.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.

Note, this function should always have output type Vector{T}.
"""
(dbootdata(data::Td, bi::BootInput)::Vector{Td}) where {Td} = [ local_get_index(data, dbootinds_one(bi)) for j = 1:bi.numresample ]
(dbootdata(data::Td ; kwargs...)::Vector{Td}) where {Td} = dbootdata(data, BootInput(data ; kwargs...))
# function dbootdata(data::Td, bi::BootInput{BootTapered})::Vector{Td} where {Td}
#     data_infl = apply_influence_function(data)
#     return [ dbootdata_one_infl(data_infl, bi, bm) for j = 1:bi.numresample ]
# end

"""
    dbootlevel1(data::T1, bi::BootInput)
    dbootlevel1(data::T1; kwargs...)

Get the level 1 bootstrapped statistics associated with dataset in data, and bootstrap methodology in BootInput.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.

Note, the return type is determined by bi.flevel1, which must be a function that accepts T1,
ie typeof(data), as input. It may return any output type T2, as long as bi.flevel2 will
accept Vector{T2} as input.

For example, if data is a Vector{<:Number} then bi.flevel1 might be the function mean,
which in this case will return Float64, so bi.flevel2 must be some function that can
accept Vector{Float64} as input.

A more complicated example: if data is Matrix{<:Number} then bi.flevel1 might be the anonymous
function x->mean(x,dims=1), which in this case will return a single row Matrix{Float64}, and
so bi.flevel2 must be some function that can accept Vector{Matrix{Float64}} as input.
"""
dbootlevel1(data, bi::BootInput) = [ bi.flevel1(dbootdata_one(data, bi)) for j = 1:bi.numresample ]
dbootlevel1(data ; kwargs...) = dbootlevel1(data, BootInput(data ; kwargs...))

"""
    dboot(data, bi::BootInput)
    dboot(data ; kwargs...)

Get the level 2 bootstrapped statistics associated with dataset in data, and bootstrap methodology in BootInput.

A keyword method that calls the keyword constructor for BootInput is also provided. Please
use ?BootInput at the REPL for more detail on feasible keywords.

Note, the return type of the output will be determined by bi.flevel2, which must be a function that accepts
Vector{T}, where T is the output type of bi.flevel1.

For example, if data is a Vector{<:Number} and bi.flevel1 is mean, then in this case, bi.flevel1 will return
Float64, and so bi.flevel2 must be some function that accepts Vector{Float64} as input (and can have any output
type.)

Alternatively, bi.flevel2 could be the anonymous function (x -> quantile(x, [0.025, 0.975])), in which case
the input should be Vector{Float64}, and so bi.flevel1 should return Float64. Note, the output of bi.flevel2
in this case will be a 2-element Vector{Float64} with elements corresponding bootstrapped 95% confidence interval
for the level1 statistic of the input dataset
"""
dboot(data, bi::BootInput) = bi.flevel2(dbootlevel1(data, bi))
dboot(data ; kwargs...) = dboot(data, BootInput(data ; kwargs...))

"dbootlevel2 <- Identical to the dboot function. This function is only included for naming consistency with dbootlevel1"
dbootlevel2(data, bi::BootInput)= dboot(data, bi)
dbootlevel2(data ; kwargs...)= dbootlevel2(data, BootInput(data ; kwargs...))

"dbootvar <- Identical to dboot but with the level 2 statistic set to variance"
dbootvar(data ; kwargs...) = dboot(data, BootInput(data ; flevel2=var, kwargs...))

"dbootconf <- Identical to dboot but with the level 2 statistic set to a confidence interval with width determined by keyword alpha. Default alpha=0.05 corresponds to a 95% confidence interval."
dbootconf(data ; alpha::Float64=0.05, kwargs...) = (0.0 < alpha < 0.5) ? dboot(data, BootInput(data ; flevel2=(x -> quantile(x, [alpha/2, 1-(alpha/2)])), kwargs...)) : error("Invalid alpha of $(alpha) for confidence interval. alpha must lie on the interval (0.0, 0.5)")


text:
const HBAR = 1.0
const IOHBAR = (0.0 + 1.0im) / HBAR


text:
# Plotting functions, for illustrating plume rise evolution under different situations
#
# Scritto da: Patrizia Favaron
#
# Questo codice è coperto dalla licenza MIT.

using Plots

function show_max_x_vs(
    us::Float64,            # Velocità del vento alla quota del vertice della ciminiera (m/s)
    vs_min::Float64,        # Velocità verticale minima dei fumi all'uscita dalla ciminiera (m/s)
    vs_max::Float64,        # Velocità verticale massima dei fumi all'uscita dalla ciminiera (m/s)
    ds::Float64,            # Diametro della ciminiera al vertice (m)
    Ts::Float64,            # Temperatura dei fumi all'uscita dalla ciminiera (K)
    Ta::Float64             # Temperatura ambiente (K)
)

    gr()

    vs = range(vs_min, vs_max; length=100)

    max_x_vals = Array{Float64}(undef, length(vs), 2)
    for i in 1:length(vs)
        max_x_vals[i,1] = max_x(1, us, vs[i], ds, Ts, Ta) # Categoria A: convettività
        max_x_vals[i,2] = max_x(6, us, vs[i], ds, Ts, Ta) # Categoria F: stabilità
    end

    p = plot(vs, max_x_vals, label = ["Cat. A" "Cat. F"]);
    xlabel!("Vel. fumi (m/s)");
    ylabel!("x alla max altezza (m)");
    display(p)

    println(max_x_vals[1,1])
    println(max_x_vals[2,2])

    println(max_x_vals[length(vs),1])
    println(max_x_vals[length(vs),2])

end

max_stars_repo_name: mafavaron/books-pbl-software

max_issues_repo_name: mafavaron/books-pbl-software

max_forks_repo_name: mafavaron/books-pbl-software


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# create a restart file

trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"))

###############################################################################
# adapt the parameters that have changed compared to "elixir_advection_extended.jl"

# Note: If you get a restart file from somewhere else, you need to provide
# appropriate setups in the elixir loading a restart file

restart_filename = joinpath("out", "restart_000040.h5")
mesh = load_mesh(restart_filename)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

tspan = (load_time(restart_filename), 2.0)
ode = semidiscretize(semi, tspan, restart_filename);

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using PkgBenchmark
using Random
using MIToS.Utils
using MIToS.MSA

@benchgroup "Utils" begin
    include("Utils/GeneralUtils.jl")
end

@benchgroup "MSA" begin
    include("MSA/Residues.jl")
    include("MSA/Annotations.jl")
end


text:
function solve()
    n, q = [parse(Int, x) for x in split(readline())]
    h = w = n
    ret = (n-2) * (n-2)
    a = ones(Int, n) * n
    b = ones(Int, n) * n
    for _=1:q
        t, x = [parse(Int, x) for x in split(readline())]
        if t == 1
            if x < w
                b[x:w-1] .= h
                w = x
            end
            ret -= b[x] - 2
        else
            if x < h
                a[x:h-1] .= w
                h = x
            end
            ret -= a[x] - 2
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module OldSK

using StaticArrays, JuLIP
using JuLIP: AbstractAtoms

"""
`abstract type SKHamiltonian{NORB}` : supertype for Hamiltonian
operators of Slater-Koster type
"""
abstract type SKHamiltonian{NORB} end

# TODO: general implementation
"""
`nbonds(H)` : number of hopping functions needed to evaluate the
hamiltonian.
"""
nbonds(::SKHamiltonian{1}) = 1
nbonds(::SKHamiltonian{4}) = 4
nbonds(::SKHamiltonian{9}) = 10

"""
`isorthogonal(H) -> Bool` : specify whether the tb model is orthogonal
or non-orthogonal (has non-trivial overlap). Default is `false`.
"""
isorthogonal(H::SKHamiltonian) = false

"""
`hop(H::SKHamiltonian, r, i)`: where `r` is real, `i` integer, this
should return the
"""
function hop end

# TODO: we don't know yet what this means!
# hop_d(H::SKHamiltonian, r, i) = ForwardDiff.derivative(s -> hop(H,s,i), r)

"""
`hop!(H::SKHamiltonian, r, bonds)` : fill `bonds` with `hop(H, r, i)`
"""
function hop!(H::SKHamiltonian, r, bonds)
   for i = 1:nbonds(H)
      bonds[i] = hop(H, r, i)
   end
   return bonds
end

# # TODO: hop_d should return hop and hop_d, since we need both!!!!
# function hop_d!(H::SKHamiltonian, r, b, db)
#    for i = 1:nbonds(H)
#       b[i] = hop(H, r, i)
#       db[i] = hop_d(H, r, i)
#    end
#    return b, db
# end

"""
`overlap(H::SKHamiltonian, r, i)`: where `r` is real, `i` integer, this
should return the overlap SK parameter
"""
function overlap end

# overlap_d(H::SKHamiltonian, r::Real, i) = ForwardDiff.derivative(s->overlap(H,s,i), r)

"""
`overlap!(H::SKHamiltonian, r, bonds)` : fill `bonds` with `overlap(H, r, i)`
"""
function overlap!(H::SKHamiltonian, r, bonds)
   for i = 1:nbonds(H)
      bonds[i] = overlap(H, r, i)
   end
   return bonds
end

# function overlap_d!(H::SKHamiltonian, r, b, db)
#    for i = 1:nbonds(H)
#       b[i] = overlap(H, r, i)
#       db[i] = overlap_d(H, r, i)
#    end
#    return b, db
# end

# TODO: discuss the onsite terms
# rewrite onsite! as diagonals???
# but we don't provide AD for these since they depend on too many variables,
# so AD will necessarily be inefficient.
function onsite! end
# function onsite_grad! end

# Matrix Assembly

const SKBlock{NORB} = SMatrix{NORB, NORB, Float64}

"""
`struct SparseSKH` :
a triplet sparse matrix kind of thing that stores pre-computed
hamiltonian blocks, and can efficiently generate k-dependent Hamiltonians,
either sparse of full.

### Methods:
* `full(A::SparseSKH, k)`: returns full `H`, `M` for given `k`-vector.
* `full!(out, A, k)`: same as `full`, but in-place

### Methods to be implemented:
* `sparse(A::SparseSKH, k)`: same but sparse
* `collect(A::SparseSKH, k)`: chooses full (small systems) or sparse (large
                           systems) based on some simple heuristic

### Fields:
* `H` : the hamiltonian used to construct it
* `at` : the associated atoms object
* `i, j` : row and column indices
* `first` : `first[n]` is the index in `i, j` for which `i[idx] == n`
* `vH` : Hamiltonian blocks
* `vM` : overlap blocks
* `Rcell` : each hamiltonian block is associated with an e^{i k ⋅ S} multiplier
            from the Bloch transform; this S is stored in Rcell.
"""
struct SparseSKH{HT, TV}  # v0.6: require that TV <: SKBlock{NORB}
   H::HT
   at::AbstractAtoms
   i::Vector{Int32}
   j::Vector{Int32}
   first::Vector{Int32}
   vH::Vector{TV}
   vM::Vector{TV}
   Rcell::Vector{JVecF}
end

Base.length(skh::SparseSKH) = length(skh.i)

"""
Assemble a hamiltonian block for a Slater-Koster type
hamiltonian, with 4 orbitals (s, p).

### Parameters:
* `U` : R / |R|, orientation of the bond, must be a 3-vector
* `hop` : vector of hopping function values, must be a 4-vector
* `mat` : output matrix, must be at least 4 x 4
"""
function sk4!(U, hop, mat)
   l, m, n = U[1], U[2], U[3]
   # 4 orbitals are s, px, py, pz, these are the mat-indices
   # 4 bond types are : ssσ,spσ,ppσ,ppπ, these are the hop-indices
   #                     1   2   3   4
   mat[1,1] = hop[1]                            # E_ss = V_ssσ
   mat[2,2] = l*l * hop[3] + (1-l*l) * hop[4]   # E_xx = l² V_ppσ + (1-l²) V_ppπ
   mat[3,3] = m*m * hop[3] + (1-m*m) * hop[4]   # E_yy = m² V_ppσ + (1-m²) V_ppπ
   mat[4,4] = n*n * hop[3] + (1-n*n) * hop[4]   # E_zz = n² V_ppσ + (1-n²) V_ppπ
   mat[1,2] = l * hop[2]                        # E_sx = l V_spσ
   mat[1,3] = m * hop[2]                        # E_sy = m V_spσ
   mat[1,4] = n * hop[2]                        # E_sz = n V_spσ
   mat[2,1] = - mat[1,2]                        # E_xs = - E_sx
   mat[3,1] = - mat[1,3]                        # E_ys = - E_sy
   mat[4,1] = - mat[1,4]                        # E_zs = - E_sz
   mat[2,3] = l * m * (hop[3] - hop[4])         # E_xy = l m (V_ppσ - V_ppπ)
   mat[2,4] = l * n * (hop[3] - hop[4])         # E_xz = l n (V_ppσ - V_ppπ)
   mat[3,4] = m * n * (hop[3] - hop[4])         # E_yz = n m (V_ppσ - V_ppπ)
   mat[3,2] =  mat[2,3]                         # E_yx = E_xy
   mat[4,2] =  mat[2,4]                         # E_zx = E_xz
   mat[4,3] =  mat[3,4]                         # E_zy = E_yz
   return mat
end

"""
Assemble a hamiltonian block for a slater-koster type
hamiltonian, with 9 orbitals (s, p, d).

* `U` : orientation of the bond
* `hop` : hopping functions for bond type, 10-vector
* `mat` : output matrix 9 x 9
"""
function sk9!(U, hop, mat)
   # fill the [1:4, 1:4] block
   sk4!(U, hop, mat)
   # and then all the rest
   l, m, n = U[1], U[2], U[3]
   # ll, mm, nn = l*l, m*m, n*n    # TODO: for faster evaluation???
   # lm, ln, mn = l*m, l*n, m*n
   # rt3 = √3

   # sd
   mat[1,5] = √3 * l * m * hop[5]
   mat[1,6] = √3 * m * n * hop[5]
   mat[1,7] = √3 * l * n * hop[5]
   mat[1,8] = √3/2 * (l^2 - m^2) * hop[5]
   mat[1,9] = ( n^2 - (l^2 + m^2)/2 ) * hop[5]
   mat[5,1] = mat[1,5]
   mat[6,1] = mat[1,6]
   mat[7,1] = mat[1,7]
   mat[8,1] = mat[1,8]
   mat[9,1] = mat[1,9]

   # pd
   mat[2,5] = √3 * l * l * m * hop[6] + m * (1.0 - 2.0 * l^2) * hop[7]
   mat[2,6] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[2,7] = √3 * l * l * n * hop[6] + n * (1.0 - 2.0 * l^2) * hop[7]
   mat[2,8] = √3/2 * l * (l^2 - m^2) * hop[6] + l * (1.0 - l^2 + m^2) * hop[7]
   mat[2,9] = l * (n^2 - (l^2 + m^2)/2) * hop[6] - √3 * l * n^2 * hop[7]
   mat[5,2] = - mat[2,5]
   mat[6,2] = - mat[2,6]
   mat[7,2] = - mat[2,7]
   mat[8,2] = - mat[2,8]
   mat[9,2] = - mat[2,9]
   mat[3,5] = √3 * l * m * m * hop[6] + l * (1.0 - 2.0 * m^2) * hop[7]
   mat[3,6] = √3 * m * m * n * hop[6] + n * (1.0 - 2.0 * m^2) * hop[7]
   mat[3,7] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[3,8] = √3/2 * m * (l^2 - m^2) * hop[6] - m * (1.0 + l^2 - m^2) * hop[7]
   mat[3,9] = m * (n^2 - (l^2 + m^2)/2) * hop[6] - √3 * m * n^2 * hop[7]
   mat[5,3] = - mat[3,5]
   mat[6,3] = - mat[3,6]
   mat[7,3] = - mat[3,7]
   mat[8,3] = - mat[3,8]
   mat[9,3] = - mat[3,9]
   mat[4,5] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[4,6] = √3 * m * n * n * hop[6] + m * (1.0 - 2.0 * n^2) * hop[7]
   mat[4,7] = √3 * l * n * n * hop[6] + l * (1.0 - 2.0 * n^2) * hop[7]
   mat[4,8] = √3/2 * n * (l^2 - m^2) * hop[6] - n * (l^2 - m^2) * hop[7]
   mat[4,9] = n * (n^2 - (l^2 + m^2)/2) * hop[6] + √3 * n * (l^2 + m^2) * hop[7]
   mat[5,4] = - mat[4,5]
   mat[6,4] = - mat[4,6]
   mat[7,4] = - mat[4,7]
   mat[8,4] = - mat[4,8]
   mat[9,4] = - mat[4,9]

   # dd
   mat[5,5] = 3.0 * l^2 * m^2 * hop[8] + (l^2 + m^2 - 4.0 * l^2 * m^2) * hop[9] +
               (n^2 + l^2 * m^2) * hop[10]
   mat[6,6] = 3.0 * m^2 * n^2 * hop[8] + (m^2 + n^2 - 4.0 * m^2 * n^2) * hop[9] +
               (l^2 + m^2 * n^2) * hop[10]
   mat[7,7] = 3.0 * l^2 * n^2 * hop[8] + (l^2 + n^2 - 4.0 * l^2 * n^2) * hop[9] +
               (m^2 + l^2 * n^2) * hop[10]
   mat[8,8] = 3.0/4 * (l^2 - m^2)^2 * hop[8] + (l^2 + m^2 - (l^2 - m^2)^2) * hop[9] +
               (n^2 + (l^2 - m^2)^2 /4 ) * hop[10]
   mat[9,9] = (n^2 - (l^2 + m^2) /2)^2 * hop[8] + 3.0 * n^2 * (l^2 + m^2) * hop[9] +
               3.0/4 * (l^2 + m^2)^2 * hop[10]
   mat[5,6] = 3.0 * l * m^2 * n * hop[8] + l * n * (1.0 - 4.0 * m^2) * hop[9] +
               l * n * (m^2 - 1.0) * hop[10]
   mat[5,7] = 3.0 * l^2 * m * n * hop[8] + m * n * (1.0 - 4.0 * l^2) * hop[9] +
               m * n * (l^2 - 1.0) * hop[10]
   mat[5,8] = 3.0/2 * l * m * (l^2 - m^2) * hop[8] + 2.0 * l * m * (m^2 - l^2) * hop[9] +
               1.0/2 * l * m * (l^2 - m^2) * hop[10]
   mat[5,9] = √3 * l * m * (n^2 - (l^2 + m^2)/2) * hop[8] - 2.0*√3 * l * m * n^2 * hop[9] +
               √3/2 * l * m * (1.0 + n^2) * hop[10]
   mat[6,7] = 3.0 * l * m * n^2 * hop[8] + l * m * (1.0 - 4.0 * n^2) * hop[9] +
               l * m * (n^2 - 1.0) * hop[10]
   mat[6,8] = 3.0/2 * m * n * (l^2 - m^2) * hop[8] -
               m * n * (1.0 + 2.0 * (l^2 - m^2)) * hop[9] +
               m * n * (1.0 + (l^2 - m^2) /2) * hop[10]
   mat[6,9] = √3 * m * n * (n^2 - (l^2 + m^2)/2) * hop[8] +
               √3 * m * n * (l^2 + m^2 - n^2) * hop[9] -
               √3/2 * m * n * (l^2 + m^2) * hop[10]
   mat[7,8] = 3.0/2 * l * n * (l^2 - m^2) * hop[8] +
               l * n * (1.0 - 2.0 * (l^2 - m^2)) * hop[9] -
               l * n * (1.0 - (l^2 - m^2) /2) * hop[10]
   mat[7,9] = √3 * l * n * (n^2 - (l^2 + m^2) /2) * hop[8] +
               √3 * l * n * (l^2 + m^2 - n^2) * hop[9] -
               √3/2 * l * n * (l^2 + m^2) * hop[10]
   mat[8,9] = √3/2 * (l^2 - m^2) * (n^2 - (l^2 + m^2) /2) * hop[8] +
               √3 * n^2 * (m^2 - l^2) * hop[9] +
               √3/4 * (1.0 + n^2) * (l^2 - m^2) * hop[10]
   mat[6,5] = mat[5,6]
   mat[7,5] = mat[5,7]
   mat[8,5] = mat[5,8]
   mat[9,5] = mat[5,9]
   mat[7,6] = mat[6,7]
   mat[8,6] = mat[6,8]
   mat[9,6] = mat[6,9]
   mat[8,7] = mat[7,8]
   mat[9,7] = mat[7,9]
   mat[9,8] = mat[8,9]
   return mat
end

function sk4_d!(U, r, hh, dhh, dmat)
   # U = ∇r, so it appears here both as the orientation and as ∇r
   l, m, n = U[1], U[2], U[3]
   ll, lm, ln, mm, mn, nn = l*l, l*m, l*n, m*m, m*n, n*n
   dl = ( (1.0-ll)/r ,     - lm/r ,     - ln/r )
   dm = (     - lm/r , (1.0-mm)/r ,     - mn/r )
   dn = (     - ln/r ,     - mn/r , (1.0-nn)/r )

   for d = 1:3
      dmat[d,1,1] = dhh[1] * U[d]
      dmat[d,2,2] = (ll * dhh[3] + (1 .- ll) * dhh[4]) * U[d] + 2*l * hh[3] * dl[d] - 2*l * hh[4] * dl[d]
      dmat[d,3,3] = (mm * dhh[3] + (1 .- mm) * dhh[4]) * U[d] + 2*m * hh[3] * dm[d] - 2*m * hh[4] * dm[d]
      dmat[d,4,4] = (nn * dhh[3] + (1 .- nn) * dhh[4]) * U[d] + 2*n * hh[3] * dn[d] - 2*n * hh[4] * dn[d]
      dmat[d,1,2] = l * dhh[2] * U[d] + hh[2] * dl[d]
      dmat[d,1,3] = m * dhh[2] * U[d] + hh[2] * dm[d]
      dmat[d,1,4] = n * dhh[2] * U[d] + hh[2] * dn[d]
      dmat[d,2,1] = - dmat[d,1,2]
      dmat[d,3,1] = - dmat[d,1,3]
      dmat[d,4,1] = - dmat[d,1,4]
      dmat[d,2,3] = lm * (dhh[3] - dhh[4]) * U[d] + (dl[d] * m + l * dm[d]) * (hh[3] - hh[4])
      dmat[d,2,4] = ln * (dhh[3] - dhh[4]) * U[d] + (dl[d] * n + l * dn[d]) * (hh[3] - hh[4])
      dmat[d,3,4] = mn * (dhh[3] - dhh[4]) * U[d] + (dm[d] * n + m * dn[d]) * (hh[3] - hh[4])
      dmat[d,3,2] =  dmat[d,2,3]
      dmat[d,4,2] =  dmat[d,2,4]
      dmat[d,4,3] =  dmat[d,3,4]
   end
   return dmat
end

sk!(out, H::SKHamiltonian{4}, U, bonds) = sk4!(U, bonds, out)
sk!(out, H::SKHamiltonian{9}, U, bonds) = sk4!(U, bonds, out)

using StaticArrays

using JuLIP: AbstractAtoms, positions, cutoff, neighbourlist
using NeighbourLists: sites, npairs
using LinearAlgebra: dot, diagind, Hermitian

"""
`indexblock`: auxiliary function to compute indices of Slater Koster orbitals,
this is returned as an SVector, i.e. it is generated on the stack so that no
heap memory is allocated.
"""
indexblock(n::Integer, H::SKHamiltonian{NORB}) where {NORB} =
   SVector{NORB, Int}( ((n-1)*NORB+1):(n*NORB) )

ndofs(H::SKHamiltonian{NORB}, at::AbstractAtoms) where {NORB} =
      ndofs(NORB, length(at))
ndofs(norb::Integer, nat::Integer) =
      norb * nat

function SparseSKH(H::SKHamiltonian{NORB}, at::AbstractAtoms) where {NORB}

   # here is a little hack that will turn an abstract type into a concrete type
   SKB = typeof(zero(SKBlock{NORB}))

   # get back positions since we still can't do without to get the cell-shifts
   # TODO: >>> maybe implement an alternative in JuLIP?
   X = positions(at)

   # ALLOCATIONS: TODO: move to a different method?
   # neighbourlist
   nlist = neighbourlist(at, cutoff(H))
   #      off-diagonal + diagonal
   nnz = npairs(nlist) + length(at)
   # allocate space for ordered triplet format
   i = zeros(Int32, nnz)
   j = zeros(Int32, nnz)
   first = zeros(Int32, length(at))
   vH = zeros(SKB, nnz)
   vM = zeros(SKB, nnz)
   Rcell = zeros(JVecF, nnz)
   # index into the triplet format
   idx = 0

   # allocate space to assemble the hamiltonian blocks, we use MMatrix
   # here, but is this really necessary? Matrix should do fine?
   H_nm = zero(MMatrix{NORB, NORB, Float64})
   M_nm = zero(MMatrix{NORB, NORB, Float64})
   bonds = zeros(nbonds(H))     # temporary array for storing the potentials

   # loop through sites
   # for (n, neigs, r, R) in sites(nlist)
   for (n, neigs, R) in sites(nlist)
      r = norm.(R)
      first[n] = idx+1          # where do the triplet entries for atom n start?

      # --------- on-site terms ----------
      # add the diagonal/on-site entries
      onsite!(H, r, R, H_nm, M_nm)
      # onsite_overlap!(H, r, R, M_nm)   # TODO: should this be just (H, M)
      idx += 1
      i[idx], j[idx], vH[idx], vM[idx] = n, n, SKB(H_nm), SKB(M_nm)
      # Rvector corresponding to this block is just X[n]-X[n] = 0
      Rcell[idx] = zero(JVecF)

      # loop through the neighbours of the current atom (i.e. the bonds)
      for m = 1:length(neigs)
         U = R[m]/r[m]
         # hamiltonian block; TODO: new framework plugs into here!!!! <<<<<<<<<< !!!!!!
         sk!(H_nm, H, U, hop!(H, r[m], bonds))
         sk!(M_nm, H, U, overlap!(H, r[m], bonds))
         idx += 1
         i[idx], j[idx], vH[idx], vM[idx] = n, neigs[m], SKB(H_nm), SKB(M_nm)
         # compute the Rcell vector for these blocks
         Rcell[idx] = R[m] - (X[neigs[m]] - X[n])
      end
   end

   return SparseSKH(H, at, i, j, first, vH, vM, Rcell)
end

_alloc_full(skh::SparseSKH) = _alloc_full(skh.H, skh.at)

_alloc_full(H::SKHamiltonian, at::AbstractAtoms) =
         ( Matrix{ComplexF64}(undef, ndofs(H, at), ndofs(H, at)),
           Matrix{ComplexF64}(undef, ndofs(H, at), ndofs(H, at)) )

Base.Matrix(H::SparseSKH, k::AbstractVector = zero(JVecF)) =
         _full!(_alloc_full(H), H, k)

_full!(out::Tuple, H::SparseSKH, k::AbstractVector = zero(JVecF)) =
         _full!(out[1], out[2], H, k, H.H)

function _full!(Hout, Mout, skh, k, H::SKHamiltonian{NORB}) where {NORB}
   fill!(Hout, 0.0)
   fill!(Mout, 0.0)
   k = JVecF(k)
   for (i, j, H_ij, M_ij, S) in zip(skh.i, skh.j, skh.vH, skh.vM, skh.Rcell)
      eikR = exp( im * dot(k, S) )
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      @inbounds for a = 1:NORB, b = 1:NORB
         Hout[Ii[a], Ij[b]] += H_ij[a, b] * eikR
         Mout[Ii[a], Ij[b]] += M_ij[a, b] * eikR
      end
   end
   # TODO: enforce numerical self-adjointness here???
   #       e.g. by wrapping into Hermition???
   Hout[diagind(Hout)] = real(Hout[diagind(Hout)])
   Mout[diagind(Mout)] = real(Mout[diagind(Mout)])
   return Hout, Mout
end

hamiltonian(H::SKHamiltonian, at::AbstractAtoms, k = zero(JVecF)) =
            Hermitian.(Matrix(SparseSKH(H, at), k))

using Parameters

using JuLIP.Potentials: fcut, fcut_d, SplineCutoff,
   PairPotential, ZeroPairPotential, EAM, @analytic

import JuLIP: cutoff

# export KwonHamiltonian

"""
`KwonHamiltonian <: SKHamiltonian{4}`

Hamiltonian for an orthogonal sp TB model of Si developed by Kwon et al [1].

This implementation deviates  from [1] in how the cut-off is applied:
instead of "patching" a cubic spline between r1 and rcut, we simply multiply
with a quintic spline on the interval [0.5 (rcut + r0), rcut].

[1] I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho and C. M. Soukoulis.
Transferable tight-binding models for silicon.
Phys Rev B 49 (11), 1994.
"""
@with_kw struct KwonHamiltonian <: SKHamiltonian{4}
   r0::Float64 = 2.360352   # Å
   Es::Float64 = -5.25       # eV
   Ep::Float64 = 1.2        # eV
   E0::Float64 = 8.7393204  # eV   8.7393204  is the original value; but it is just a constant
   # ------------------------------- Electronic Parameters
   # α   1     2    3     4
   #    ssσ   spσ  ppσ   ppπ
   hr0::NTuple{4, Float64} = (-2.038, 1.745, 2.75, -1.075)
   nc::NTuple{4, Float64} = (9.5, 8.5, 7.5, 7.5)
   rc::NTuple{4, Float64} = (3.4, 3.55, 3.7, 3.7)
   # ------------------------------- 2-body parameters
   m::Float64 = 6.8755
   mc::Float64 = 13.017
   dc::Float64 = 3.66995     # Å
   C::NTuple{4, Float64} = (2.1604385, -0.1384393, 5.8398423e-3, -8.0263577e-5)
   # --------------------------------
   r1::Float64 = 3.260176     # Å  (start of cut-off) >>> different r1 from Kwon paper
   rcut::Float64 = 4.16       # Å  (end of cut-off)
end

cutoff(H::KwonHamiltonian) = H.rcut

isorthogonal(H::KwonHamiltonian) = true

kwon_hop(H::KwonHamiltonian, r::Real, α) = ( H.hr0[α] * (H.r0 / r)^2 *
               exp( - 2 * (r/H.rc[α])^H.nc[α] + 2 * (H.r0/H.rc[α])^H.nc[α] ) )

hop(H::KwonHamiltonian, r::Real, α) = kwon_hop(H, r, α) * fcut(r, H.r1, H.rcut)

overlap(H::KwonHamiltonian, r, i) = 0.0

function onsite!(H::KwonHamiltonian, _r, _R, H_nn, M_nn)
   fill!(H_nn, 0.0)
   fill!(M_nn, 0.0)
   H_nn[1,1] = H.Es
   H_nn[2,2] = H_nn[3,3] = H_nn[4,4] = H.Ep
   M_nn[1,1] = M_nn[2,2] = M_nn[3,3] = M_nn[4,4] = 1.0
   return H_nn, M_nn
end

end


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

function datadown(thepath, thefile, theurl)
    filepath = joinpath(thepath, thefile)
    try
        download(theurl, filepath)
    catch
        if isfile(filepath)
            @warn("Failed to update $thefile, but you already have it. Things might continue to work, but if you would like to make sure you have the latest version of $thefile, use may use your web-browser to download it from $theurl and place it in $_pkg_deps.")
        else
            error("Failed to download $thefile from $theurl. You may use your web-browser to download it from $theurl and place it in $_pkg_deps.")
        end
    end
end

datadown(_pkg_deps, "plotschema.json",
         "https://api.plot.ly/v2/plot-schema?sha1")
datadown(_pkg_assets, "plotly-latest.min.js",
         "https://cdn.plot.ly/plotly-latest.min.js")

include("make_schema_docs.jl")
include("find_attr_groups.jl")
AttrGroups.main()


text:
using Test, Random, Statistics
using Flux, DistributedFlux, MPI

@testset "DistributedFlux.jl" begin
    Random.seed!(1234)

    x = randn(Float32, 100, 16)
    y = mean(x, dims = 1)

    model = Chain(Dense(100, 100), Dense(100, 1))

    function loss(x, y)
        ŷ = model(x)
        Flux.mse(y, ŷ)
    end

    ps = Flux.params(model)
    opt = ADAM(1e-3, (0.9, 0.999))

    MPI.Init()

    data = Flux.Data.DataLoader((x, y), batchsize = 4)
    Flux.train!(loss, ps, data, opt, verbose = true, sync = true)
    Flux.train!(loss, ps, data, opt, verbose = true, sync = false)

    MPI.Finalize()
end


text:
module SpringCollab2020MadelineSilhouetteTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/MadelineSilhouetteTrigger" MadelineSilhouetteTrigger(x::Integer, y::Integer,
    width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, enable::Bool=true)

const placements = Ahorn.PlacementDict(
    "Madeline Silhouette (Spring Collab 2020)" => Ahorn.EntityPlacement(
        MadelineSilhouetteTrigger,
        "rectangle"
    )
)

end


text:
module SimpleHypergraphs

using LightGraphs
using StatsBase
using DataStructures
using PyCall
using Conda
using PyPlot
using JSON3
using Random
using LinearAlgebra

export Hypergraph, getvertices, gethyperedges
export add_vertex!, add_hyperedge!, remove_vertex!, remove_hyperedge!
export prune_hypergraph!, prune_hypergraph
export set_vertex_meta!, get_vertex_meta
export set_hyperedge_meta!, get_hyperedge_meta
export adjacency_matrix, edge_adjacency_matrix

export BipartiteView, TwoSectionView
export shortest_path
export get_twosection_adjacency_mx, get_twosection_weighted_adjacency_mx
export dual
export random_model, random_kuniform_model, random_dregular_model, random_preferential_model

export Abstract_HG_format, HGF_Format, JSON_Format
export hg_load, hg_save

export modularity
export randompartition
export AbstractCommunityFinder, CFModularityRandom, CFModularityCNMLike, CFLabelPropagationFinder
export findcommunities

export nhv, nhe
export random_walk
export get_connected_components
export conductance
export AbstractDistance, SnodeDistanceDijkstra, SedgeDistanceDijkstra
export distance

export HyperNetX, GraphBased
export draw

const hnx  = PyNULL()
const nx = PyNULL()
const has_plotting = Ref(false)

function __init__()
    has_networkx = false
    has_hypernetx = false
    try
        copy!(nx, pyimport("networkx"))
        has_networkx = true
    catch e; end
    try
        copy!(hnx, pyimport("hypernetx"))
        has_hypernetx = true
    catch e; end
    has_plotting[] = has_networkx && has_hypernetx
    if !has_plotting[]
        @warn "The plotting functionality of HyperNetX will not work!\n"*
		(has_networkx ? "" : "Conda Python networkx not found.\n")*
        (has_hypernetx ? "" : "Conda Python HyperNetX not found.\n")*
        "To test your installation try running `using PyCall;pyimport(\"networkx\");pyimport(\"hypernetx\")`"
    end
end

function support_hypernetx()
    return has_plotting[]
end

include("hypergraph.jl")
include("io.jl")

include("models/bipartite.jl")
include("models/twosection.jl")
include("models/random-models.jl")
include("models/dual.jl")

include("algorithms/conductance.jl")
include("algorithms/distance.jl")

include("algorithms/community/modularity.jl")
include("algorithms/community/label-propagation.jl")

include("viz/drawing.jl")
include("viz/widget.jl")

end # module