text:
# This is the startup.jl file that should reside in
# ~/.julia/config/startup.jl
# A simple way to install it is to run the ./install.sh script
# in this directory.

"""
    SimplePackageLoader(pkg_name,verbose=false)

This is a function form of the `using` command in Julia.
Instead of `using SomePackage` we can instead call
`SimplePackageLoader("SomePackage")`. With `verbose` set to
`true` each `using` command is printed.

This can also be called with a list of package names.
"""
function SimplePackageLoader(pkg_name::String, verbose::Bool = false)
    command = "using $pkg_name"
    if verbose
        @info command
    end
    eval(Meta.parse(command))
    nothing
end

function SimplePackageLoader(pkgs, verbose::Bool = false)
    for p in pkgs
        SimplePackageLoader(p, verbose)
    end
    nothing
end

"""
`list_of_favorites()` returns a list of the
packages to be loaded when `SimpleFavorites()` is called.
"""
function list_of_favorites()
    packages = [
        "Permutations"
        "SimplePartitions"
        "SimpleGraphs"
        "BigCombinatorics"
        "SimpleTools"
        "LinearAlgebra"
        "ShowSet"
        "Primes"
        "SimpleRandom"
        "SimplePolynomials"
    ]
    return packages
end

"""
`list_of_graph_theory()` are the packages to use when focused on
graph computing.
"""
function list_of_graph_theory()
    packages = [
        "SimpleGraphs"
        "SimpleGraphAlgorithms"
        "ChooseOptimizer"
        "Gurobi"
        "ShowSet"
        "SimplePolynomials"
    ]
    return packages
end

"""
`list_of_graph_draw()` returns a list of the packages I use
for graph visualization.
"""
function list_of_graph_draw()
    packages = [
        "Plots"
        "SimpleDrawing"
        "DrawSimpleGraphs"
    ]
end

"""
    GraphTheory(draw_pkgs::Bool=true, verbose::Bool=true)

loads my packages for working on graphs
and sets the optmizer to be `Gurobi` (non-verbose). The list is generated by `list_of_graph_theory()`.

If `draw_pkgs` is set to false, don't load the drawing packages returned by  `list_of_graph_draw()`
"""
function GraphTheory(draw_pkgs::Bool = true, verbose::Bool = true)
    SimplePackageLoader(list_of_graph_theory(), verbose)
    eval(Meta.parse("use_Gurobi()"))

    if draw_pkgs
        SimplePackageLoader(list_of_graph_draw(), verbose)
    else
        if verbose
            @info "Graph drawing packages not loaded"
        end
    end
end
"""
    box_my_text(line)

prints out a single line of text embedded in a
box consisting of corner and connector characters.
```
julia> box_my_text("Hello")
┌───────┐
│ Hello │
└───────┘
```
"""
function box_my_text(line::String)
    UL = "$(Char(9484))"
    UR = "$(Char(9488))"
    LL = "$(Char(9492))"
    LR = "$(Char(9496))"
    HOR = "$(Char(9472))"
    VER = "$(Char(9474))"
    SPA = " "

    n = length(line)
    println(UL * HOR^(n + 2) * UR)
    println(VER * SPA * line * SPA * VER)
    println(LL * HOR^(n + 2) * LR)
    nothing
end

# When the "east" servers are not updating right, do this:

"""
    change_package_server(server::String = "us-west.pkg.julialang.org")

Sometimes the default package server doesn't get updated, so this is a way to
specify a different one that seems to be more reliable.
"""
function change_package_server(server::String = "us-west.pkg.julialang.org")
    Base.ENV["JULIA_PKG_SERVER"] = server
    println("Julia package server set to $server")
    nothing
end

# box_my_text("Use GraphTheory() to load graph theory packages")


text:
@testset "MaskedDenseArray/MaskedHashArray" begin
    for T in (Float32, Float64)
        for ArrayType in (TreeArrays.MaskedDenseArray, TreeArrays.MaskedHashArray)
            x = (@inferred ArrayType{T}(undef, 2, 3))::ArrayType{T, 2}
            @test size(x) == (2, 3)
            @test size(x.mask) == (2, 3)
            @test any(x.mask) == false

            # getindex, checkactive
            @test_throws Exception x[i]

            # setindex
            x[1] = 2
            @test TreeArrays.isactive(x, 1) == true
            @test x[1] == 2

            # setindex nothing
            x[1] = nothing
            @test TreeArrays.isactive(x, 1) == false
            if x isa TreeArrays.MaskedDenseArray
                @test x.data[1] == 2 # still value is not changed
            else
                @test x.data[TreeArrays.FastHashInt(1)] == 2 # still value is not changed
            end
        end
    end
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 7a0128c0-a13f-11eb-1098-ef9c00274a0d
using Plots

# ╔═╡ 7c5e40c4-f84e-4527-a87e-c6ac395f5d9e
function ansatz(x,α)
	return ifelse(abs(x) > 1,1 - abs(x)^α,0)
end;

# ╔═╡ 007bb513-b173-4723-9d8b-2e858ae54dc6
function probability_density(x,α)
	return (ansatz(x,α)^2)
end;

# ╔═╡ dcc6bdb1-daa8-44f1-a6ec-285ac1916179
function E_local(x,α)
	return ifelse(abs(x)<1,(0.5*α*(α - 1)*abs(x)^(α - 2))/(1 - abs(x)^α),0)
end;

# ╔═╡ 99eff10f-5b1e-42d9-9009-56e5aeb821c2
function metropolis(N,α)
	L = 2
	x = rand(-1:0.001:1)
	E = 0
	E2 = 0
	Eln_avg = 0
	ln_avg = 0
	rejection_ratio = 0
	#start the loop
	for i in 1:N
		x_t = x + 0.01*rand(-L:0.000001:L)
		if probability_density(x_t,α) >= probability_density(x,α)
			x = x_t
		else
			dummy = rand()
			if dummy < probability_density(x_t,α)/probability_density(x,α)
				x = x_t
			else
				rejection_ratio += 1/N
			end
		end

		E += E_local(x,α)/N
		E2 += E_local(x,α)^2/N
		Eln_avg += (E_local(x,α)*((-1*log(abs(x))*abs(x)^α)/(1 - abs(x)^α)))/(N)
		ln_avg += ((-1*log(abs(x))*abs(x)^α)/(1 - abs(x)^α))/N
	end
	return E,E2,Eln_avg,ln_avg,rejection_ratio
end;

# ╔═╡ cbfdf51f-721e-4172-b7d7-1911bf666163
function iterate_alpha()
	α = 1.1
	alpha_iter = 50
	N = 500
	random_walkers = 200
	γ = 0.01
	energy = Vector{Float64}(undef,alpha_iter)
	alpha = Vector{Float64}(undef,alpha_iter)
	variance = Vector{Float64}(undef,alpha_iter)
	for i in 1:alpha_iter
		E = 0
		E2 = 0
		dE_dalpha = 0
		Eln = 0
		ln = 0
		rejection_ratio = 0
		for j in 1:random_walkers
			E_met, E2_met, Eln_met, ln_met, _ = metropolis(N,α)
			E += E_met/random_walkers
			E2 += E2_met/random_walkers
			Eln += Eln_met/random_walkers
			ln += ln_met/random_walkers
			# rejection_ratio += rejections_met/random_walkers
		end

		# Define the next α
		dE_dalpha = 2*(Eln-E*ln)
		α = α + 0.05#γ*dE_dalpha

		# Update the arrays for plotting
		energy[i] = E
		alpha[i] = α
		variance[i] = E2 - E^2
	end
	return energy,alpha,variance
end;

# ╔═╡ 20151696-5be9-4b68-afb7-8a7217c16402
energy,alpha,variance = iterate_alpha();

# ╔═╡ e5f66da9-ef79-4f82-9e1b-6ad7b27ec75d
energy

# ╔═╡ 64f872bd-cb8b-42de-8fa3-3270430e9a7f
alpha

# ╔═╡ 5ca1aa2b-099e-479c-88ad-7c9764fa6f85
begin
	plot(energy,
		label = false,
		title = "Energy expenctation Plot",
		linecolor = :blue,
		ylabel = "<E>",
		xlabel = "timestep"
	)

	hline!([π^2/8],
		label = false
	)
end

# ╔═╡ b21c39ff-68f7-4725-9bdc-4c180b5b368b
plot(variance,
	title = "Variance in energy",
	label = false,
	linecolor = :blue,
	ylabel = "Var(E)",
	xlabel = "timestep"
)

# ╔═╡ bb6cfb6d-325e-4913-9f77-3efcfcfc4368
plot(alpha,
	title = "Variational Parameter",
	label = false,
	linecolor = :blue,
	ylabel = "α",
	xlabel = "timestep"
)

# ╔═╡ 11dd75f3-e142-4712-914d-d8d594acada3
begin
	scatter(alpha,
		energy,
		title = "Energy vs α",
		xlabel = "α",
		ylabel = "<E>",
		label = false,
		markercolor = :blue
	)

	hline!([π^2/8],
		label = false
	)
end

# ╔═╡ 43a28afe-622c-4cfc-9679-9554994be6bf
variance

# ╔═╡ Cell order:
# ╠═7a0128c0-a13f-11eb-1098-ef9c00274a0d
# ╠═7c5e40c4-f84e-4527-a87e-c6ac395f5d9e
# ╠═007bb513-b173-4723-9d8b-2e858ae54dc6
# ╠═dcc6bdb1-daa8-44f1-a6ec-285ac1916179
# ╠═99eff10f-5b1e-42d9-9009-56e5aeb821c2
# ╠═cbfdf51f-721e-4172-b7d7-1911bf666163
# ╠═20151696-5be9-4b68-afb7-8a7217c16402
# ╠═e5f66da9-ef79-4f82-9e1b-6ad7b27ec75d
# ╠═64f872bd-cb8b-42de-8fa3-3270430e9a7f
# ╠═5ca1aa2b-099e-479c-88ad-7c9764fa6f85
# ╠═b21c39ff-68f7-4725-9bdc-4c180b5b368b
# ╠═bb6cfb6d-325e-4913-9f77-3efcfcfc4368
# ╠═11dd75f3-e142-4712-914d-d8d594acada3
# ╠═43a28afe-622c-4cfc-9679-9554994be6bf

max_stars_repo_name: Ved-Mahajan/VMC

max_issues_repo_name: Ved-Mahajan/VMC

max_forks_repo_name: Ved-Mahajan/VMC


text:
@testset "Orders in absolute number fields" begin
  include("NfOrd/NfOrd.jl")
  include("NfOrd/Elem.jl")
  include("NfOrd/Ideal.jl")
  include("NfOrd/FracIdl.jl")
  include("NfOrd/ResidueRing.jl")
  include("NfOrd/Clgp.jl")
  include("NfOrd/RayClassGroup.jl")
  include("NfOrd/ResidueRingMultGrp.jl")
  include("NfOrd/LinearAlgebra.jl")
end


text:
# Basic definitions of openbuildnet.jl

abstract OBNPort    # Abstract port class

# A (physical) input port
abstract OBNInputAbstract <: OBNPort

# A (physical) output port
abstract OBNOutputAbstract <: OBNPort


text:
import Distributions

export Dists
const Dists = Distributions

function Measure(dist::Dists.Distribution{F,S}) where {F,S}
    X = Base.eltype(dist)
    DistributionMeasure{F,S,X}(dist)
end

struct DistributionMeasure{F, S, X} <: AbstractMeasure
    dist :: Dists.Distribution{F, S}
end

function basemeasure(μ::DistributionMeasure{F,S,X}) where {F, S <: Dists.Continuous, X}
    return Lebesgue(X)
end

export logdensity

function logdensity(μ::DistributionMeasure{F,S,X}, x::X) where {F, S, X}
    return Dists.logpdf(μ.dist, x)
end

function Base.rand(rng::AbstractRNG, μ::DistributionMeasure{F,S,X}) where {F, S, X}
    return rand(rng, μ.dist)
end

basemeasure(d::Dists.AbstractMvNormal) = Lebesgue(ℝ)^size(d)

function basemeasure(μ::Dists.Distribution{Dists.Univariate,Dists.Continuous})
    return Lebesgue(ℝ)
end

function basemeasure(μ::Dists.Distribution{Dists.Univariate,Dists.Discrete})
    return CountingMeasure(ℤ)
end

∫(::typeof(identity), ::Dists.Distribution) = 1.0

logdensity(μ::Dists.Distribution, x) = Dists.logpdf(μ,x)

density(μ::Dists.Distribution, x) = Dists.pdf(μ,x)

testvalue(d::Dists.Distribution) = rand(d)


text:
"""
Follow a given route, longitudinal acceleration is controlled by the IntelligentDriverModel
"""
@with_kw mutable struct RouteFollowingIDM <: LaneFollowingDriver
    route::Vector{Lane{Float64}} = Lane{Float64}[]
    dir::Int64 = 1
    a::Float64 = NaN # predicted acceleration
    σ::Float64 = NaN # optional stdev on top of the model, set to zero or NaN for deterministic behavior

    k_spd::Float64 = 1.0 # proportional constant for speed tracking when in freeflow [s⁻¹]

    δ::Float64 = 4.0 # acceleration exponent [-]
    T::Float64  = 1.5 # desired time headway [s]
    v_des::Float64 = 8.0 # desired speed [m/s]
    s_min::Float64 = 1.0 # minimum acceptable gap [m]
    a_max::Float64 = 2.0 # maximum acceleration ability [m/s²]
    d_cmf::Float64 = 2.0 # comfortable deceleration [m/s²] (positive)
    d_max::Float64 = 9.0 # maximum decelleration [m/s²] (positive)
end

AutomotiveDrivingModels.get_name(::RouteFollowingIDM) = "RouteFollowingIDM"

function AutomotiveDrivingModels.set_desired_speed!(model::RouteFollowingIDM, v_des::Float64)
    model.v_des = v_des
    model
end

AutomotiveDrivingModels.reset_hidden_state!(model::RouteFollowingIDM) = model

"""
Compute acceleration according to IDM
"""
function AutomotiveDrivingModels.track_longitudinal!(model::RouteFollowingIDM, v_ego::Float64, v_oth::Float64, headway::Float64)

    if !isnan(v_oth) && headway > 0.
        # @assert !isnan(headway) && headway > 0

        Δv = v_oth - v_ego
        s_des = model.s_min + v_ego*model.T - v_ego*Δv / (2*sqrt(model.a_max*model.d_cmf))
        v_ratio = model.v_des > 0.0 ? (v_ego/model.v_des) : 1.0
        model.a = model.a_max * (1.0 - v_ratio^model.δ - (s_des/headway)^2)
    elseif headway < 0.
        model.a = -model.d_max
    else
        # no lead vehicle, just drive to match desired speed
        Δv = model.v_des - v_ego
        model.a = Δv*model.k_spd # predicted accel to match target speed
    end

    @assert !isnan(model.a)

    model.a = clamp(model.a, -model.d_max, model.a_max)

    model
end

"""
set the lane exit to take at the next junction to reach the goal
"""
function set_direction!(model::RouteFollowingIDM, scene::Scene, roadway::Roadway, egoid::Int64)
    ego = scene[findfirst(egoid, scene)]
    cur_lane = get_lane(roadway, ego)
    set_direction!(model, cur_lane, roadway)
end
function set_direction!(model::RouteFollowingIDM, lane::Lane, roadway::Roadway)
    model.dir = get_direction(lane, model.route)
end

function get_direction(lane::Lane{T}, route::Vector{Lane{T}}) where T<:Real
    dir = 1
    ind = findfirst(isequal(lane), route)
    if ind == nothing
        ind = 0
    end
    if ind == length(route)
        return dir
    end
    next_exit = route[ind+1]
    for i=1:length(lane.exits)
        if lane.exits[i].target.tag == next_exit.tag
            dir=i
            break
        end
    end
    return dir
end

function AutomotiveDrivingModels.observe!(model::RouteFollowingIDM, scene::Scene, roadway::Roadway, egoid::Int)

    vehicle_index = findfirst(egoid, scene)

    fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway,
                                        VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    v_ego = scene[vehicle_index].state.v

    if fore.ind != nothing
        headway, v_oth = fore.Δs, scene[fore.ind].state.v
    else
        headway, v_oth = NaN, NaN
    end
    track_longitudinal!(model, v_ego, v_oth, headway)
    set_direction!(model, scene, roadway, egoid)

    return model
end

function Base.rand(model::RouteFollowingIDM)
    if isnan(model.σ) || model.σ ≤ 0.0
        LonAccelDirection(model.a, model.dir)
    else
        LonAccelDirection(rand(Normal(model.a, model.σ)), model.dir)
    end
end

## Discrete distribution
# P(a | s)

function action_space(model::RouteFollowingIDM)
    return (LonAccelDirection(-model.σ, model.dir), LonAccelDirection(model.σ, model.dir), LonAccelDirection(model.σ, model.dir))
end

function get_distribution(model::RouteFollowingIDM)
    actions = action_space(model)
    probs = (1/3, 1/3, 1/3)
    return actions, probs
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Documenter, URIs

makedocs(;
    modules=[URIs],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaWeb/URIs.jl/blob/{commit}{path}#L{line}",
    sitename="URIs.jl",
    authors = "Jacob Quinn, Sam O'Connor and contributors: https://github.com/JuliaWeb/URIs.jl/graphs/contributors"
)

deploydocs(;
    repo="github.com/JuliaWeb/URIs.jl",
    push_preview=true
)


text:
import RowEchelon

#=
Homework 1
Section 1.2
Problem 3
=#

# Create matrix to be evalutaed
A = [1 2 4 8 ;
     2 4 6 8 ;
     3 6 9 12]

# Print matrix for reference
println("Input Matrix:")
println(A)
println()

# Evaluate with RowEchelon
A_rref = RowEchelon.rref(A)

# Print result
println("Reduced Row Echelon Form:")
println(A_rref)


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

mutable struct ExpressRouteServiceProviderBandwidthsOffered <: SwaggerModel
    offerName::Any # spec type: Union{ Nothing, String } # spec name: offerName
    valueInMbps::Any # spec type: Union{ Nothing, Int32 } # spec name: valueInMbps

    function ExpressRouteServiceProviderBandwidthsOffered(;offerName=nothing, valueInMbps=nothing)
        o = new()
        validate_property(ExpressRouteServiceProviderBandwidthsOffered, Symbol("offerName"), offerName)
        setfield!(o, Symbol("offerName"), offerName)
        validate_property(ExpressRouteServiceProviderBandwidthsOffered, Symbol("valueInMbps"), valueInMbps)
        setfield!(o, Symbol("valueInMbps"), valueInMbps)
        o
    end
end # type ExpressRouteServiceProviderBandwidthsOffered

const _property_map_ExpressRouteServiceProviderBandwidthsOffered = Dict{Symbol,Symbol}(Symbol("offerName")=>Symbol("offerName"), Symbol("valueInMbps")=>Symbol("valueInMbps"))
const _property_types_ExpressRouteServiceProviderBandwidthsOffered = Dict{Symbol,String}(Symbol("offerName")=>"String", Symbol("valueInMbps")=>"Int32")
Base.propertynames(::Type{ ExpressRouteServiceProviderBandwidthsOffered }) = collect(keys(_property_map_ExpressRouteServiceProviderBandwidthsOffered))
Swagger.property_type(::Type{ ExpressRouteServiceProviderBandwidthsOffered }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ExpressRouteServiceProviderBandwidthsOffered[name]))}
Swagger.field_name(::Type{ ExpressRouteServiceProviderBandwidthsOffered }, property_name::Symbol) =  _property_map_ExpressRouteServiceProviderBandwidthsOffered[property_name]

function check_required(o::ExpressRouteServiceProviderBandwidthsOffered)
    true
end

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


text:
#########
# Phase #
#########

abstract type Phase end

struct Execute <: Phase end

struct Transform <: Phase end

proceed(::Execute) = Transform()
proceed(::Transform) = Execute()

###############
# TraceConfig #
###############

get_world_age() = ccall(:jl_get_tls_world_age, UInt, ()) # ccall(:jl_get_world_counter, UInt, ())

abstract type AbstractPass end

struct TraceConfig{C<:Context,M,w,b,P<:Union{AbstractPass,Unused}}
    context::C
    metadata::M
    world::Val{w}
    boxes::Val{b}
    pass::P
end

const AnyTraceConfig{world} = TraceConfig{<:Context,<:Any,world}

function TraceConfig(context::Context;
                     metadata = unused,
                     world::Val = Val(get_world_age()),
                     boxes::Val = Val(false),
                     pass::Union{AbstractPass,Unused} = unused)
    return TraceConfig(context, metadata, world, boxes, pass)
end

@inline prehook(::AnyTraceConfig{w}, ::Vararg{Any}) where {w} = nothing
@inline posthook(::AnyTraceConfig{w}, ::Vararg{Any}) where {w} = nothing
@inline is_user_primitive(cfg::AnyTraceConfig{w}, f, args...) where {w} = false
@inline is_core_primitive(cfg::AnyTraceConfig{w}, f, args...) where {w} = _is_core_primitive(cfg, Core.Typeof(f), args...)
@inline execution(cfg::AnyTraceConfig{w}, f, args...) where {w} = unboxcall(cfg.context, f, args...)

@generated function _is_core_primitive(cfg::TraceConfig{C,M,w}, ::Type{F}, args...) where {C,M,w,F}
    ftype = unbox(C, F)
    atypes = Tuple(unbox(C, T) for T in args)
    signature = Tuple{ftype,atypes...}
    # TODO: this is slow, we should try to check whether CodeInfo is retrievable
    # rather than going through the whole process of actually retrieving it
    if isa(lookup_method_body(signature; world = w), CodeInfo)
        result = :(false)
    else
        result = :(true)
    end
    return quote
        $(Expr(:meta, :inline))
        $(result)
    end
end

###########
# Overdub #
###########

struct Overdub{P<:Phase,F,C<:TraceConfig}
    phase::P
    func::F
    config::C
    function Overdub(phase::P, func, config::C) where {P,C}
        F = Core.Typeof(func) # this yields `Type{T}` instead of `UnionAll` for constructors
        return new{P,F,C}(phase, func, config)
    end
end

@inline overdub(::Type{C}, f; phase = Execute(), kwargs...) where {C<:Context} = overdub(C(f), f; phase = phase, kwargs...)
@inline overdub(ctx::Context, f; phase = Execute(), kwargs...) = Overdub(phase, f, TraceConfig(ctx; kwargs...))

@inline proceed(o::Overdub, f = o.func) = Overdub(proceed(o.phase), f, o.config)

@inline context(o::Overdub) = o.config.context

@inline func(o::Overdub) = o.func
@inline func(f) = f

Base.show(io::IO, o::Overdub{P}) where {P} = print("Overdub{$(P.name.name)}($(typeof(context(o)).name), $(func(o)))")

####################
# Overdub{Execute} #
####################

@inline function execute(o::Overdub{Execute}, args...)
    if is_user_primitive(o.config, o.func, args...)
        return execution(o.config, o.func, args...)
    else
        return proceed(o)(args...)
    end
end

@inline function (o::Overdub{Execute})(args...)
    prehook(o.config, o.func, args...)
    output = execute(o, args...)
    posthook(o.config, o.func, output, args...)
    return output
end

######################
# Overdub{Transform} #
######################

# Note that this approach emits code in which LHS SSAValues are not
# monotonically increasing. This currently isn't a problem, but in
# the future, valid IR might require monotonically increasing LHS
# SSAValues, in which case we'll have to add an extra SSA-remapping
# pass to this function.
function overdub_pass!(method_body::CodeInfo, boxes_enabled::Bool)
    # set up new SSAValues
    self_ssa = SSAValue(method_body.ssavaluetypes)
    method_body.ssavaluetypes += 1
    ctx_ssa = SSAValue(method_body.ssavaluetypes)
    greatest_ssa_value = method_body.ssavaluetypes

    # set up replacement code
    new_code = copy_prelude_code(method_body.code)
    prelude_end = length(new_code)
    push!(new_code, :($self_ssa = $(GlobalRef(Cassette, :func))($(SlotNumber(1)))))
    push!(new_code, :($ctx_ssa = $(GlobalRef(Cassette, :context))($(SlotNumber(1)))))
    in_overdub_region = false

    # fill in replacement code
    for i in (prelude_end + 1):length(method_body.code)
        stmnt = method_body.code[i]
        if stmnt === BEGIN_OVERDUB_REGION
            in_overdub_region = true
        else
            # replace `SlotNumber(1)` references with `self_ssa`
            replace_match!(s -> self_ssa, s -> isa(s, SlotNumber) && s.id == 1, stmnt)
            if in_overdub_region
                # replace calls with overdubbed calls
                replace_match!(s -> is_call(s), stmnt) do call
                    greatest_ssa_value += 1
                    new_ssa_value = SSAValue(greatest_ssa_value)
                    new_ssa_stmnt = Expr(:(=), new_ssa_value, Expr(:call, GlobalRef(Cassette, :proceed), SlotNumber(1), call.args[1]))
                    push!(new_code, new_ssa_stmnt)
                    call.args[1] = new_ssa_value
                    return call
                end
            end
            push!(new_code, stmnt)
        end
    end

    if boxes_enabled
        # replace all `new` expressions with calls to `Cassette._newbox`
        replace_match!(x -> Base.Meta.isexpr(x, :new), new_code) do x
            return Expr(:call, GlobalRef(Cassette, :_newbox), ctx_ssa, x.args...)
        end
    end

    method_body.code = fix_labels_and_gotos!(new_code)
    method_body.ssavaluetypes = greatest_ssa_value + 1
    return method_body
end

function overdub_transform_call_generator(::Type{F}, ::Type{C}, ::Type{M}, world, boxes, pass, f, args) where {F,C,M}
    ftype = unbox(C, F)
    atypes = Tuple(unbox(C, T) for T in args)
    signature = Tuple{ftype,atypes...}
    try
        method_body = lookup_method_body(signature; world = world, pass = pass)
        if isa(method_body, CodeInfo)
            method_body = overdub_pass!(method_body, boxes)
            method_body.inlineable = true
            method_body.signature_for_inference_heuristics = Core.svec(ftype, atypes, world)
            @safe_debug "returning overdubbed codeinfo" method_body
        else
            method_body = quote
                $(Expr(:meta, :inline))
                $Cassette.execution(f.config, f.func, $(OVERDUB_ARGS_SYMBOL)...)
            end
            @safe_debug "no codeinfo found; executing as primitive"
        end
        return method_body
    catch err
        @safe_error "error compiling" signature context=C
        errmsg = "ERROR COMPILING $signature IN CONTEXT $C: " * sprint(showerror, err)
        return quote
            error($errmsg)
        end
    end
end

function overdub_transform_call_definition(pass, line, file)
    return quote
        function (f::$Cassette.Overdub{$Cassette.Transform,F,$Cassette.TraceConfig{C,M,world,boxes,pass}})($(OVERDUB_ARGS_SYMBOL)...) where {F,C,M,world,boxes,pass<:$pass}
            $(Expr(:meta,
                   :generated,
                   Expr(:new,
                        Core.GeneratedFunctionStub,
                        :overdub_transform_call_generator,
                        Any[:f, OVERDUB_ARGS_SYMBOL],
                        Any[:F, :C, :M, :world, :boxes, :pass],
                        @__LINE__,
                        QuoteNode(Symbol(file)),
                        true)))
        end
    end
end

eval(overdub_transform_call_definition(:Unused, @__LINE__, @__FILE__))


text:
module NimrodFiles

export readnimrod_hdr, readnimrod, Nimrod
export toasciigrid

"""The horizontal grid type for the data. National Grid=0, latitude/longitude,
space view, polar stereographic, xy or other."""
@enum NimrodGridType nationalgrid=0 latlon=1 space=2 polar=3 xy=4 other=5

"""
Representation of the contents of a Nimrod file. The type parameter `T` depends
on the data type and size declared in fields 12 and 13 of the header.

Fields
======

- `hdr_dump::Array{Any}` - Contents of header in the order defined in the Nimrod
  specification.
- `data::Matrix{T}` - Data contents of a Nimrod file.
- `gridtype::NimrodGridType` - Grid type of data coordinates (e.g., lat/lon, NG).
- `xpix_size::Real` - Size of pixels in the x dimension.
- `ypix_size::Real` - Size of pixels in the y dimension.
- `llcorner::Tuple{Real,Real}` - (Easting,Northing) of lower left corner of the image.
"""
type Nimrod{T<:Real}
    hdr_dump::Array{Any}
    data::Matrix{T}

    gridtype::NimrodGridType
    xpix_size::Real
    ypix_size::Real

    """(Easting, Northing) of lower left corner of image."""
    llcorner::Tuple{Real,Real}
    """(Easting, Northing) of the centre of the lower left pixel."""
    llcenter::Tuple{Real,Real}
    """(Easting, Northing) of the centre of the top right pixel."""
    trcenter::Tuple{Real,Real}
end

include("./io.jl")
include("./asciigrid.jl")
end # module


text:
using LinearAlgebra, Statistics

export embed, remove_embedding, has_embedding
export set_vertex_color, set_line_color, get_vertex_color, get_line_color

export set_fill_color

export GraphEmbedding, show
# export draw, draw_labels, unclip

export transform, scale, rotate, translate, recenter
export getxy
export set_vertex_size, get_vertex_size
export edge_length

const DEFAULT_MARKER_SIZE = 6
# MARKER_SIZE = DEFAULT_MARKER_SIZE

"""
`GraphEmbedding(G)` creates a drawing of the `SimpleGraph` `G`. It
is given a circular embedding.

**Note**: Direct usage of `GraphEmbedding` is **deprecated**!

`GraphEmbedding(G,d)` creates a drawing using the `Dict` `d` to
specify the vertex locations. `d` should map vertices of `G` to
`Vector`s that specify x,y-coordinates.
"""
mutable struct GraphEmbedding
    G::SimpleGraph
    xy::Dict{Any,Vector{Float64}}
    line_color::String
    fill_color::String
    vertex_size::Int

    function GraphEmbedding(GG::SimpleGraph)
        d = circular_embedding(GG)
        new(GG,d,"black","white",DEFAULT_MARKER_SIZE)
    end

    function GraphEmbedding(GG::SimpleGraph,d::Dict)
        T = vertex_type(GG)
        dd = Dict{T,Vector{Float64}}()

        try
            for v in GG.V
                dd[v] = d[v]
            end
        catch
            error("Dictionary doesn't include all vertices")
        end
        new(GG,dd,"black","white", DEFAULT_MARKER_SIZE)
    end
end

"""
`remove_embedding(G)` erases the graph's embedding saved in its cache.
"""
function remove_embedding(G::SimpleGraph)
  cache_clear(G,:GraphEmbedding)
  nothing
end

"""
`has_embedding(G)` returns `true` if an embedding has been created
for `G` in its cache.
"""
function has_embedding(G::SimpleGraph)
  return cache_check(G,:GraphEmbedding)
end

"""
`get_embedding_direct(G)` returns a direct reference to a
graph's embedding. This function is not exposed. Users should
use `cache_recall(G,:GraphEmbedding)` instead to retrieve
a copy of the drawing. If there is not drawing, a default
(circular) drawing is created.
"""
function get_embedding_direct(G::SimpleGraph)
  if !cache_check(G,:GraphEmbedding)
    X = GraphEmbedding(G)
    G.cache[:GraphEmbedding] = X
  end
  return G.cache[:GraphEmbedding]
end

"""
`set_embedding_direct(G,X)` overwrites (or creates) the drawing
`X` in the graph's cache. This non-exposed function should not be
called by the user.
"""
function set_embedding_direct(G::SimpleGraph, X::GraphEmbedding)
  G.cache[:GraphEmbedding] = X
  nothing
end

"""
`set_vertex_color(G,color)` sets the color that gets drawn in the
interior of the vertices. (All vertices get the same color.)
"""
function set_vertex_color(G::SimpleGraph, color::String="white")
  if !has_embedding(G)
    embed(G)
  end
  G.cache[:GraphEmbedding].fill_color=color
  nothing
end

"""
`get_vertex_color(G)` returns the color used to fill in vertices.
See `set_vertex_color`.
"""
function get_vertex_color(G::SimpleGraph)
  if !has_embedding(G)
    embed(G)
  end
  return G.cache[:GraphEmbedding].fill_color
end

"""
`set_line_color(G,color)` sets the color used to draw edges and
the circles around vertices.
"""
function set_line_color(G::SimpleGraph, color::String="black")
  if !has_embedding(G)
    embed(G)
  end
  G.cache[:GraphEmbedding].line_color = color
  nothing
end

"""
`get_line_color(G)` returns the color used to draw edges and the
circles around vertices. See `set_line_color`.
"""
function get_line_color(G::SimpleGraph)
  if !has_embedding(G)
    embed(G)
  end
  return G.cache[:GraphEmbedding].line_color
end

"""
`list2dict(list)` takes a list of `(Symbol,Any)` pairs and
converts them to a dictionary mapping the symbols to their
associated values.
"""
function list2dict(list::Vector)
  d = Dict{Symbol,Any}()
  for it in list
    x,y = it
    d[x]=y
  end
  return d
end

"""
`embed(G)` creates a new embedding for `G`. The full call is
```
embed(G,algorithm;args...)
```
The `symbol` algorithm indicates the embedding algorithm.
The `args` are a collection of possible arguments to be sent
to the algorithm.

The `algorithm` defaults to `:circular` and may be one of the following:

* `:circular`: arrange the vertices evenly in a circle.
* `:random`: arrange the vertices randomly.
* `:spring`: use the `spring` layout from `GraphLayout`. Optional argument:
  * `iterations`.
* `:stress`: use the `stress` layout from `GraphLayout`.
* `:spectral`: use the `spectral` embedding. Optional arguments:
  * `xcol` -- which eigenvalue to use for the `x` coordinate.
  * `ycol` -- which eigenvalue to use for the `y` coordinate

Note that if the graph already has (say) an embedding, that embedding may
be used as the starting point for one of the algorithms.
"""
function embed(G::SimpleGraph, algorithm::Symbol=:circular; args...)
  X = get_embedding_direct(G)
  arg_dict = list2dict(collect(args))

  ## FULL LIST OF POSSIBLE ARGUMENTS PRE-PARSED HERE

  verbose::Bool = true
  if haskey(arg_dict,:verbose)
    verbose = arg_dict[:verbose]
  end

  iterations::Int = 0
  if haskey(arg_dict,:iterations) && arg_dict[:iterations]>0
    iterations = arg_dict[:iterations]
  end

  tolerance = 0.001
  if haskey(arg_dict,:tolerance) && arg_dict[:tolerance]>0
    tolerance = arg_dict[:tolerance]
  end

  xcol::Int = 2
  if haskey(arg_dict,:xcol) && arg_dict[:xcol] > 0
    xcol = arg_dict[:xcol]
  end

  ycol::Int = 3
  if haskey(arg_dict,:ycol) && arg_dict[:ycol] > 0
    xcol = arg_dict[:ycol]
  end

  if algorithm == :circular
    circular!(X)
    return nothing
  end

  if algorithm == :random
    random!(X)
    return nothing
  end

  if algorithm == :spring
    if iterations <= 0
      spring!(X)
    else
      spring!(X,iterations)
    end
    return nothing
  end

  if algorithm == :stress
    stress!(X)
    return nothing
  end
  #
  # if algorithm == :distxy
  #   distxy!(X,tolerance,verbose)
  #   return nothing
  # end

  if algorithm == :combined
    embed(G,:spring,iterations = iterations)
    scale(G)
    embed(G,:stress)
    return nothing
  end

  if algorithm == :spectral
    spectral!(X,xcol, ycol)
    return nothing
  end

  @warn "Unknown embedding algorithm; no action taken"
  return nothing
end

"""
`embed(G,d)` specifies an embedding of the graph `G` with
a dictionary `d` that maps vertices to coordinates (as two
dimensional vectors `[x,y]`).
"""
function embed(G::SimpleGraph,d::Dict)
  X = GraphEmbedding(G,d)
  set_embedding_direct(G,X)
  nothing
end

function circular_embedding(G::SimpleGraph{T}) where T
    d = Dict{T, Vector{Float64}}()
    n = NV(G)

    vv = vlist(G)  # vertices as a list

    s = 2*sin(pi/n)

    t = collect(0:n-1)*2*pi/n
    x = map(sin,t)/s
    y = map(cos,t)/s

    for i = 1:n
        v = vv[i]
        d[v] = [x[i],y[i]]
    end
    return d
end

function show(io::IO, X::GraphEmbedding)
    print(io,"Embedding of $(X.G)")
end

"""
`set_vertex_size(G,sz)` sets the size of the circle used
when drawing vertices. The default is `$DEFAULT_MARKER_SIZE`.
See also `get_vertex_size`.
"""
function set_vertex_size(G::SimpleGraph, m::Int = DEFAULT_MARKER_SIZE)
  if m < 1
    error("Vertex size must be nonnegative")
  end
  if !has_embedding(G)
    embed(G)
  end
  G.cache[:GraphEmbedding].vertex_size = m
  nothing
end

"""
`get_vertex_size(G)` returns the size of the circle used when
drawing vertices. See also `set_vertex_size`.
"""
function get_vertex_size(G::SimpleGraph)
  if !has_embedding(G)
    embed(G)
  end
  return G.cache[:GraphEmbedding].vertex_size
end

"""
`circular!(X)` arranges the vertices of the graph held in the
`GraphEmbedding` around a circle.
"""
function circular!(X::GraphEmbedding)
    X.xy = circular_embedding(X.G)
    rescale!(X)
    return
end

function private_adj(X::GraphEmbedding)
    n = NV(X.G)
    A = zeros(Int,n,n)
    vv = vlist(X.G)
    for i=1:n-1
        u = vv[i]
        for j=(i+1):n
            w = vv[j]
            if has(X.G,u,w)
                A[i,j] = 1
                A[j,i] = 1
            end
        end
    end
    return A,vv
end

function private_dist(X::GraphEmbedding)
    n = NV(X.G)
    d = dist(X.G)

    A = zeros(n,n)
    vv = vlist(X.G)

    for i=1:n
        u = vv[i]
        for j=1:n
            w = vv[j]
            A[i,j] = d[u,w]
            if A[i,j] < 0
                A[i,j] = n/2  # should be enough to separate the comps
            end
        end
    end

    return A,vv
end

# TEMPORARILY OFF LINE DURING TRANSITION TO JULIA 0.7

include("my_spring.jl")

"""
`spring!(X)` gives the graph held in `X` with a spring embedding
(based on code in the `GraphLayout` module). If runs a default number of
iterations (100) of that algorithm. To change the number of
iterations, use `spring!(X,nits)`.
"""
function spring!(X::GraphEmbedding, nits::Int=100)
    n = NV(X.G)
    A,vv = private_adj(X)

    x,y = layout_spring_adj(A,MAXITER=nits)

    d = Dict{vertex_type(X.G), Vector{Float64}}()
    for i = 1:n
        v = vv[i]
        d[v] = [x[i], y[i]]
    end
    X.xy = d
    return
end

include("my_stress.jl")

"""
`stress!(X)` computes a stress major layout using code taken from the
`GraphLayout` package.
"""
function stress!(X::GraphEmbedding)
    n = NV(X.G)
    A,vv = private_dist(X)

    currentxy = zeros(n,2)
    for i=1:n
        v = vv[i]
        currentxy[i,:] = X.xy[v]
    end

    xy = my_layout_stressmajorize_adj(A,2,nothing,currentxy)

    d = Dict{vertex_type(X.G), Vector{Float64}}()
    for i = 1:n
        v = vv[i]
        d[v] = collect(xy[i,:])
    end
    X.xy = d
    return
end

"""
`random!(X)` gives the graph held in `X` a random embedding.
"""
function random!(X::GraphEmbedding)
    rootn = sqrt(NV(X.G))
    for v in X.G.V
        X.xy[v] = [rand(), rand()]*rootn
    end
    recenter!(X)
end

function transform!(X::GraphEmbedding,
                    A::Array{S,2},
                    b::Vector{T}) where {S<:Real,T<:Real}

    # apply the transformation x |--> Ax+b to all coordinates
    for k in keys(X.xy)
        X.xy[k] = A*X.xy[k] + b
    end
    nothing
end

"""
`transform(G,A,b)` applies an affine transformation to all coordinates
in the graph's drawing. Here `A` is 2-by-2 matrix and `b` is a 2-vector.
Each point `p` is mapped to `A*p+b`.
"""
function transform(G::SimpleGraph,
                    A::Array{S,2},
                    b::Vector{T}) where {S<:Real, T<:Real}
  X = get_embedding_direct(G)
  transform!(X,A,b)
end

function rescale!(X::GraphEmbedding, m::T) where T<:Real
    A = zeros(T,2,2)
    A[1,1] = m
    A[2,2] = m
    b = zeros(T,2)
    transform!(X,A,b)
end

function rescale!(X::GraphEmbedding)
    L = mean(edge_length(X))
    if L != 0
        rescale!(X,1/L)
    else
      @warn "Cannot scale by 0. No action taken."
    end
    nothing
end

"""
`scale(G,m)` multiplies all coordinates in the graph's drawing by
`m`. If `m` is omitted, the drawing is rescaled so that the average
length of an edge equals 1.
"""
function scale(G::SimpleGraph, m::T) where T<:Real
  X = get_embedding_direct(G)
  rescale!(X,m)
end

function scale(G::SimpleGraph)
  X = get_embedding_direct(G)
  rescale!(X)
end

function rotate!(X::GraphEmbedding, theta)
    A = zeros(Float64,2,2)
    A[1,1] = cos(theta)
    A[1,2] = -sin(theta)
    A[2,1] = sin(theta)
    A[2,2] = cos(theta)
    b = zeros(Float64,2)
    transform!(X,A,b)
    return
end

"""
`rotate(G,theta)` rotate's the graph's drawing by the angle
`theta`.
"""
rotate(G::SimpleGraph, theta::T) where T<:Real = rotate!(get_embedding_direct(G),theta)

function translate!(X::GraphEmbedding, b::Vector{T}) where T
    #A = eye(T,2)
    A = Matrix{T}(I,2,2)
    transform!(X,A,b)
    return
end

"""
`translate(G,b)` translates the graph's drawing by the vector `b`;
that is, every point `p` in the drawing is replaced by `p+b`.
"""
translate(G::SimpleGraph, b::Vector{T}) where T<:Real = translate!(get_embedding_direct(G),b)

function recenter!(X::GraphEmbedding)
    b = [0;0]

    for v in X.G.V
        b += X.xy[v]
    end
    b /= NV(X.G)
    translate!(X,-b)
end

"""
`recenter(G)` translates the graph's drawing so that the center of mass
of the vertices is at the origin.
"""
recenter(G::SimpleGraph) = recenter!(get_embedding_direct(G))

"""
`draw_labels(G)` will add the vertices names to the drawing
window.

The results are usually ugly. One can try increasing the size of the
vertices (see `set_vertex_size`). The font size of the labels can be
specified using `draw_labels(G,FontSize)`. The default is 10.
"""
function draw_labels(G::SimpleGraph, FontSize::Int=10)
  X = get_embedding_direct(G)
  draw_labels(X,FontSize)
  nothing
end

function draw_labels(X::GraphEmbedding, FontSize::Int)
    for v in X.G.V
        x,y = X.xy[v]
        text(x,y,string(v),fontsize=FontSize)
    end
end

"""
`getxy(G)` returns a copy of the dictionary mapping vertices to their
x,y-coordinates. This is a way of saving an embedding.
"""
getxy(G::SimpleGraph) = deepcopy(get_embedding_direct(G).xy)

function edge_length(X::GraphEmbedding, e)
    p1 = X.xy[e[1]]
    p2 = X.xy[e[2]]
    return norm(p1-p2)
end

function edge_length(X::GraphEmbedding)
    EE = elist(X.G)

    return [ edge_length(X,e) for e in EE ]
end

"""
`edge_length(G)` returns an array containing the
lengths of the edges in the current embedding of `G`.
"""
edge_length(G::SimpleGraph) = edge_length(G.cache[:GraphEmbedding])

"""
`edge_length(G,e)` gives the distance between the
embedded endpoints of the edge `e` in the drawing `G`.
"""
edge_length(G::SimpleGraph,ee) = edge_length(G.cache[:GraphEmbedding],ee)

#include("distxy.jl") # distxy! has been replaced by stress!
include("graffle.jl")
include("spectral.jl")
include("geogebra.jl")
include("tikz.jl")
include("embedded-graphs.jl")


text:
pyplot()
cd(pwd)
dir = "results/"

#Plot parameter recovery
recoveryPlots = plotrecovery(results,(c=0,d=2),(:sampler,:Nd);save=true,dir=dir)

#Plot mean run time as a function of number of data points (Nd) for each sampler
meantimePlot = plotsummary(results,:Nd,:time,(:sampler,);save=true,dir=dir,yscale=:log10)

#Plot mean allocations as a function of number of data points (Nd) for each sampler
meanallocPlot = plotsummary(results,:Nd,:allocations,(:sampler,);save=true,dir=dir,yscale=:log10,
  ylabel="Allocations (log scale)")

#Plot mean ess per second of number of data points (Nd) for each sampler
meanallocPlot = plotsummary(results,:Nd,:ess_ps,(:sampler,);save=true,dir=dir)

#Plot density of effective sample size as function of number of data points (Nd) for each sampler
essPlots = plotdensity(results,:ess,(:sampler,:Nd);save=true,dir=dir)

#Plot density of rhat as function of number of data points (Nd) for each sampler
rhatPlots = plotdensity(results,:r_hat,(:sampler,:Nd);save=true,dir=dir)

#Plot density of time as function of number of data points (Nd) for each sampler
timePlots = plotdensity(results,:time,(:sampler,:Nd);save=true,dir=dir)

#Plot density of gc time percent as function of number of data points (Nd) for each sampler
gcPlots = plotdensity(results,:gcpercent,(:sampler,:Nd);save=true,dir=dir)

#Plot density of memory allocations as function of number of data points (Nd) for each sampler
memPlots = plotdensity(results,:allocations,(:sampler,:Nd);save=true,dir=dir,xscale=:log10,
  xlabel="Allocations (log scale)")

#Plot density of megabytes allocated as function of number of data points (Nd) for each sampler
megPlots = plotdensity(results,:megabytes,(:sampler,:Nd);save=true,dir=dir)

#Scatter plot of epsilon and effective sample size as function of number of data points (Nd) for each sampler
scatterPlots = plotscatter(results,:epsilon,:ess,(:sampler,:Nd);save=true,dir=dir)

# Density plot for tree_depth
depthPlots = plotdensity(results, :tree_depth,(:sampler, :Nd); save=true, dir=dir)
epsPlots = plotdensity(results, :epsilon,(:sampler, :Nd); save=true, dir=dir)

# Plot mean ess per second of number of data points (Nd) for each sampler
ess = plotsummary(results, :Nd, :ess, (:sampler,); save=true, dir=dir)


text:
abstract type AbstractSimData{Y,X} end

"""
    SimData <: AbstractSimData

    SimData(extent::AbstractExtent, ruleset::AbstractRuleset)

Simulation dataset to hold all intermediate arrays, timesteps
and frame numbers for the current frame of the simulation.

A simdata object is accessable in [`applyrule`](@ref) as the first parameter.

Multiple grids can be indexed into using their key if you need to read
from arbitrary locations:

```julia
funciton applyrule(data::SimData, rule::SomeRule{Tuple{A,B}},W}, (a, b), cellindex) where {A,B,W}
    grid_a = data[A]
    grid_b = data[B]
    ...
```

In single grid simulations `SimData` can be indexed directly as if it is a `Matrix`.

## Methods

- `currentframe(data::SimData)`: get the current frame number, an `Int`
- `currenttime(data::SimData)`: the current frame time, which `isa eltype(tspan)`
- `aux(d::SimData, args...)`: get the `aux` data `NamedTuple`, or `Nothing`.
  adding a `Symbol` or `Val{:symbol}` argument will get a field of aux.
- `tspan(d::SimData)`: get the simulation time span, an `AbstractRange`.
- `timestep(d::SimData)`: get the simulaiton time step.
- `radius(data::SimData)` : returns the `Int` radius used on the grid,
  which is also the amount of border padding.
- `boundary(data::SimData)` : returns the [`BoundaryCondition`](@ref) - `Remove` or `Wrap`.
- `padval(data::SimData)` : returns the value to use as grid border padding.

These are available, but you probably shouldn't use them and thier behaviour
is not guaranteed in furture versions. They will mean rule is useful only
in specific contexts.

- `extent(d::SimData)` : get the simulation [`AbstractExtent`](@ref) object.
- `init(data::SimData)` : get the simulation init `AbstractArray`/`NamedTuple`
- `mask(data::SimData)` : get the simulation mask `AbstractArray`
- `ruleset(d::SimData)` : get the simulation [`AbstractRuleset`](@ref).
- `source(data::SimData)` : get the `source` grid that is being read from.
- `dest(data::SimData)` : get the `dest` grid that is being written to.

"""
struct SimData{Y,X,G<:NamedTuple,E,RS,F,A} <: AbstractSimData{Y,X}
    grids::G
    extent::E
    ruleset::RS
    currentframe::F
    auxframe::A
end
# Get the extent, usually from an Output
SimData(x, ruleset::AbstractRuleset) = SimData(extent(x), ruleset)
# Convert grids in extent to NamedTuple
SimData(extent::AbstractExtent, ruleset::AbstractRuleset) =
    SimData(_asnamedtuple(extent), ruleset)
function SimData(extent::AbstractExtent{<:NamedTuple{Keys}}, ruleset::AbstractRuleset) where Keys
    # Calculate the neighborhood radus (and grid padding) for each grid
    y, x = gridsize(extent)
    radii = NamedTuple{Keys}(get(radius(ruleset), key, 0) for key in Keys)
    # Construct the SimData for each grid
    grids = map(init(extent), radii, padval(extent)) do in, r, pv
        ReadableGridData{y,x,r}(
            in, mask(extent), proc(ruleset), opt(ruleset), boundary(ruleset), pv
        )
    end
    SimData(grids, extent, ruleset)
end
function SimData(
    grids::G, extent::AbstractExtent, ruleset::AbstractRuleset
) where {G<:Union{NamedTuple{<:Any,<:Tuple{GridData,Vararg}},GridData}}
    currentframe = 1; auxframe = nothing
    Y, X = gridsize(extent)
    # SimData is isbits-only
    s_extent = StaticExtent(extent)
    s_ruleset = StaticRuleset(ruleset)
    SimData{Y,X,G,typeof(s_extent),typeof(s_ruleset),Int,typeof(auxframe)}(
        grids, s_extent, s_ruleset, currentframe, auxframe
    )
end
# For ConstrutionBase
function SimData{Y,X}(
    grids::G, extent::E, ruleset::RS, currentframe::F, auxframe::A
) where {Y,X,G,E,RS,F,A}
    SimData{Y,X,G,E,RS,F,A}(grids, extent, ruleset, currentframe, auxframe)
end

ConstructionBase.constructorof(::Type{<:SimData{Y,X}}) where {Y,X} = SimData{Y,X}

# Getters
extent(d::SimData) = d.extent
ruleset(d::SimData) = d.ruleset
grids(d::SimData) = d.grids
init(d::SimData) = init(extent(d))
mask(d::SimData) = mask(first(d))
aux(d::SimData, args...) = aux(extent(d), args...)
auxframe(d::SimData, key) = auxframe(d)[_unwrap(key)]
auxframe(d::SimData) = d.auxframe
tspan(d::SimData) = tspan(extent(d))
timestep(d::SimData) = step(tspan(d))
currentframe(d::SimData) = d.currentframe
currenttime(d::SimData) = tspan(d)[currentframe(d)]

# Getters forwarded to data
Base.getindex(d::SimData, key::Symbol) = getindex(grids(d), key)

"""
    Base.get(data::SimData, key::Union{Symbol,Aux,Grid}, I...)

Allows parameters to be taken from a single value, another grid or an aux array.

If aux arrays are a `DimArray` time sequence (with a `Ti` dim) the currect date will be
calculated automatically.

Currently this is cycled by default, but will use Cyclic mode in DiensionalData.jl in future.
"""
@propagate_inbounds Base.get(data::SimData, val, I...) = val
@propagate_inbounds Base.get(data::SimData, key::Grid{K}, I...) where K = data[K][I...]
@propagate_inbounds Base.get(data::SimData, key::Aux, I...) = _auxval(data, key, I...)

@propagate_inbounds Base.setindex!(d::SimData, x, I...) = setindex!(first(grids(d)), x, I...)
@propagate_inbounds Base.getindex(d::SimData, I...) = getindex(first(grids(d)), I...)
Base.keys(d::SimData) = keys(grids(d))
Base.values(d::SimData) = values(grids(d))
Base.first(d::SimData) = first(grids(d))
Base.last(d::SimData) = last(grids(d))

gridsize(d::SimData) = gridsize(first(d))
proc(d::SimData) = proc(ruleset(d))
opt(d::SimData) = opt(ruleset(d))
boundary(d::SimData) = boundary(ruleset(d))
padval(d::SimData) = padval(extent(d))
rules(d::SimData) = rules(ruleset(d))

# Get the actual current timestep, e.g. seconds instead of variable periods like Month
currenttimestep(d::SimData) = currenttime(d) + timestep(d) - currenttime(d)

# Uptate timestamp
function _updatetime(simdata::SimData, f::Integer)
    @set! simdata.currentframe = f
    @set simdata.auxframe = _calc_auxframe(simdata)
end

# When no simdata is passed in, create new SimData
function _initdata!(::Nothing, extent::AbstractExtent, ruleset::AbstractRuleset)
    SimData(extent, ruleset)
end
# Initialise a SimData object with a new `Extent` and `Ruleset`.
function _initdata!(
    simdata::AbstractSimData, extent::AbstractExtent, ruleset::AbstractRuleset
)
    map(copy!, values(simdata), values(init(extent)))
    @set! simdata.extent = StaticExtent(extent)
    @set! simdata.ruleset = StaticRuleset(ruleset)
    simdata
end

# Convert regular index to block index
@inline _indtoblock(x::Int, blocksize::Int) = (x - 1) ÷ blocksize + 1

# Convert block index to regular index
@inline _blocktoind(x, blocksize) = (x - 1) * blocksize + 1


text:
using DataFrames
using StatsPlots
using CSV
# using Plots

function get_number_of_df()
    open("number_of_loops.txt","r") do f
        global number = readlines(f)
    end
    global number = parse(Int,join(number))
    return number
end

function get_exons_symbols(query_name)
    df = DataFrame(symbol = Char[], length = Int[])
    open("$(query_name)_annotated.pir","r") do f
        global lines = readlines(f)
    end
    exons_symbol = split(lines[1])[3]
    return exons_symbol
end

function reconstruct_dataframe(query_name, exons, nb_of_dataframes)
    all_df = DataFrame(nb_df = Int[], exonNumber = Int[],symbol= String[], nbTemp = Int[], maxIdentity = Float64[],color = String[],  legend = String[], first_pos_pir = Int64[], last_pos_pir = Int64[])
    reconstructed_dataframe = DataFrame(exonNumber = Int[],symbol= String[], nbTemp = Int[], maxIdentity = Float64[],color = String[],  legend = String[], nb_df = Int64[], first_pos_pir = Int64[], last_pos_pir = Int64[])
    for i in 1:Int(nb_of_dataframes)
        global df = CSV.read("$(query_name)_$(i)_plot1.csv")
        for j in 1:size(df,1)
            push!(all_df, (i, df.exonNumber[j], string(df.symbol[j]), df.nbTemp[j], df.maxIdentity[j], df.color[j], df.legend[j], df.first_pos_pir[j], df.last_pos_pir[j]))
        end
    end
    subdf = groupby(all_df, :symbol, sort=false)
    for i in 1:length(subdf)
        row_to_take = findmax(subdf[i].maxIdentity)[2]
        # if every maxIdentity == 0, then take the one that has the most number of templates
        if findmax(subdf[i].maxIdentity)[1] ==0
            row_to_take = findmax(subdf[i].nbTemp)[2]
        end
        push!(reconstructed_dataframe, (i, subdf[i].symbol[row_to_take], subdf[i].nbTemp[row_to_take], subdf[i].maxIdentity[row_to_take], subdf[i].color[row_to_take], subdf[i].legend[row_to_take], subdf[i].nb_df[row_to_take], subdf[i].first_pos_pir[row_to_take], subdf[i].last_pos_pir[row_to_take]))
    end
    return reconstructed_dataframe
end

query_name = ARGS[1]
nb_of_dataframes = get_number_of_df()
exons = get_exons_symbols(query_name)
reconstructed_dataframe = reconstruct_dataframe(query_name, exons, nb_of_dataframes)

p1 = @df reconstructed_dataframe bar(:nbTemp,color = :color, group = :legend, bar_width=0.9,
        title = "Number of templates by exons",
        xlabel = "Exons",
        xlim = (0,size(reconstructed_dataframe,1)+1),
        xticks = (1:1:size(reconstructed_dataframe,1)+1, :symbol),
        ylabel="Number of templates",
        ylim = (0,maximum(:nbTemp)+1),
        yticks = 0:5:maximum(:nbTemp),
        legend = :topright)

savefig(p1, "$(query_name)_reconstructed_plots.pdf")
CSV.write("$(query_name)_reconstructed.csv", reconstructed_dataframe)

max_stars_repo_name: elolaine/PhyloSofS

max_issues_repo_name: elolaine/Phylosofs

max_forks_repo_name: PhyloSofS-Team/PhyloSofS


text:
import ..build_algorithm_dict
import Flux: Descent

"""
    Agent for learning online in the four rooms environment.
"""
mutable struct FourRoomsAgent{P<:Resampling.AbstractPolicy} <: JuliaRL.AbstractAgent
    μ::P
    ER::ExperienceReplay
    WER::WeightedExperienceReplay
    gvf::GVF
    training_gap::Int64
    buffer_size::Int64
    batch_size::Int64
    warm_up::Int64
    algo_dict::Dict{String, Resampling.LearningUpdate}
    sample_dict::Dict{String, String}
    value_dict::Dict{String, Array{Resampling.TabularLayer, 1}}
    counter::Int64
    state_t::Array{Int64, 1}
    action_t::Int64
    α_arr::Array{Float64, 1}
    function FourRoomsAgent(μ::P, gvf::GVF, α_arr::Array{Float64, 1},
                            training_gap, buffer_size, batch_size,
                            warm_up, parsed, env_size; max_is_ratio=1.0) where {P<:Resampling.AbstractPolicy}
        algo_dict, sample_dict, value_type_dict = build_algorithm_dict(parsed; max_is=max_is_ratio)
        value_dict = Dict{String, Array{Resampling.TabularLayer, 1}}()
        for key in keys(algo_dict)
            if value_type_dict[key] == "State"
                value_dict[key] = [Tabular(env_size...) for α in α_arr]
            else
                value_dict[key] = [Tabular(4, env_size...) for α in α_arr]
            end
        end
        er_names = [:s_t, :s_tp1, :a_t, :r, :γ_tp1, :terminal, :mu_t, :pi_t, :ρ]
        er_types = [Array{Int64, 1}, Array{Int64, 1}, Int64, Float64, Float64, Bool, Float64, Float64, Float64]
        ER = ExperienceReplay(buffer_size, er_types; column_names=er_names)
        WER = WeightedExperienceReplay(buffer_size, er_types; column_names=er_names)
        new{P}(μ, ER, WER, gvf, training_gap, buffer_size, batch_size, warm_up,
               algo_dict, sample_dict, value_dict, warm_up, [0, 0], 0, α_arr)
    end
end

JuliaRL.get_action(agent::FourRoomsAgent, s; rng=Random.GLOBAL_RNG) =
    StatsBase.sample(rng, agent.μ, s)

function JuliaRL.start!(agent::FourRoomsAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)
    agent.counter = agent.warm_up
    agent.state_t = env_s_tp1
    agent.action_t = JuliaRL.get_action(agent, env_s_tp1; rng=rng)
    return agent.action_t
end

function JuliaRL.step!(agent::FourRoomsAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)

    # After taking action at timestep t
    μ_prob = get(agent.μ, agent.state_t, agent.action_t, nothing, nothing, nothing)
    c, γ, π_prob = get(agent.gvf, agent.state_t, agent.action_t, env_s_tp1, nothing, nothing)

    experience = (agent.state_t, env_s_tp1, agent.action_t, c, γ, terminal, μ_prob, π_prob, π_prob/μ_prob)

    add!(agent.ER, experience)
    add!(agent.WER, experience, π_prob/μ_prob)

    if agent.counter == 0
        agent.counter = agent.training_gap
        train_value_functions(agent; rng=rng)
    end
    agent.counter -= 1

    agent.state_t .= env_s_tp1
    agent.action_t = JuliaRL.get_action(agent, agent.state_t; rng=rng)
    return agent.action_t
end

function train_value_functions(agent::FourRoomsAgent; rng=Random.GLOBAL_RNG)

    α_arr = agent.α_arr

    samp_er = sample(agent.ER, agent.batch_size; rng=rng)
    samp_wer = sample(agent.WER, agent.batch_size; rng=rng)

    action_tp1_er = [StatsBase.sample(rng, agent.μ, s) for s in samp_er[:s_tp1]]
    action_tp1_wer = [StatsBase.sample(rng, agent.μ, s) for s in samp_wer[:s_tp1]]
    # action_tp1 = StatsBase.sample(rng, [1, 2], StatsBase.weights(π), batch_size)

    arg_er = (samp_er[:ρ], samp_er[:s_t], samp_er[:s_tp1],
              samp_er[:r], samp_er[:γ_tp1], samp_er[:terminal],
              samp_er[:a_t], action_tp1_er, agent.gvf.policy)

    # println("ER:", samp_er[:s_t])

    arg_wer = (ones(length(samp_wer[:ρ])), samp_wer[:s_t], samp_wer[:s_tp1],
               samp_wer[:r], samp_wer[:γ_tp1], samp_wer[:terminal],
               samp_er[:a_t], action_tp1_wer, agent.gvf.policy)
    # println("WER:", samp_wer[:s_t])

    for (α_idx, α) in enumerate(α_arr)
        opt = Descent(α)
        for key in keys(agent.algo_dict)
            if agent.sample_dict[key] == "IR"
                corr_term = 1.0
                if key == "BCIR"
                    corr_term = Resampling.total(agent.WER)/size(agent.WER)[1]
                end
                update!(agent.value_dict[key][α_idx], opt, agent.algo_dict[key], arg_wer...; corr_term=corr_term)
            elseif agent.sample_dict[key] == "ER"
                update!(agent.value_dict[key][α_idx], opt, agent.algo_dict[key], arg_er...;)
            elseif agent.sample_dict[key] == "Optimal"
                samp_opt = Resampling.getrow(agent.ER.buffer, 1:size(agent.ER.buffer)[1])
                arg_opt = (samp_opt[:ρ], samp_opt[:s_t], samp_opt[:s_tp1],
                           samp_opt[:r], samp_opt[:γ_tp1], samp_opt[:terminal])
                update!(agent.value_dict[key][α_idx], opt, agent.algo_dict[key], arg_opt...;)
            end
        end
    end
end


text:
__precompile__() # this module is safe to precompile

using Conda: add_channel
using PyCall: PyNULL, PyError, PyObject, PyAny, pyimport_conda, pycall
using VersionParsing: vparse
# Extending methods
import PyCall

export f90nml

const f90nml = PyNULL()

add_channel("conda-forge")

function __init__()
    copy!(f90nml, pyimport_conda("f90nml", "f90nml", "conda-forge"))
    # Code from https://github.com/JuliaPy/PyPlot.jl/blob/caf7f89/src/init.jl#L168-L173
    vers = f90nml.__version__
    global version = try
        vparse(vers)
    catch
        v"0.0.0" # fallback
    end
end

macro pyinterface(T)
    T = esc(T)
    return quote
        # Code from https://github.com/JuliaPy/PyPlot.jl/blob/6b38c75/src/PyPlot.jl#L54-L62
        PyCall.PyObject(f::$T) = getfield(f, :o)
        Base.convert(::Type{$T}, o::PyObject) = $T(o)
        Base.:(==)(f::$T, g::$T) = PyObject(f) == PyObject(g)
        Base.:(==)(f::$T, g::PyObject) = PyObject(f) == g
        Base.:(==)(f::PyObject, g::$T) = f == PyObject(g)
        Base.hash(f::$T) = hash(PyObject(f))
        PyCall.pycall(f::$T, args...; kws...) = pycall(PyObject(f), args...; kws...)
        (f::$T)(args...; kws...) = pycall(PyObject(f), PyAny, args...; kws...)
        Base.Docs.doc(f::$T) = Base.Docs.doc(PyObject(f))
        # Code from https://github.com/JuliaPy/PyPlot.jl/blob/6b38c75/src/PyPlot.jl#L65-L71
        Base.getproperty(f::$T, s::Symbol) = getproperty(PyObject(f), s)
        Base.getproperty(f::$T, s::AbstractString) = getproperty(PyObject(f), s)
        Base.setproperty!(f::$T, s::Symbol, x) = setproperty!(PyObject(f), s, x)
        Base.setproperty!(f::$T, s::AbstractString, x) = setproperty!(PyObject(f), s, x)
        PyCall.hasproperty(f::$T, s::Symbol) = hasproperty(PyObject(f), s)
        Base.propertynames(f::$T) = propertynames(PyObject(f))
        Base.haskey(f::$T, x) = haskey(PyObject(f), x)
        # Common methods
        Base.close(f::$T) = PyObject(f).close()
    end
end # macro pyinterface


text:
@testset "Cellular automata 1D" begin
  using Agents.CA1D
  rules = Dict("111"=>"0", "110"=>"0", "101"=>"0", "100"=>"1", "011"=>"0",
  "010"=>"1", "001"=>"1", "000"=>"0")  # rule 22
  ncols = 101
  model = CA1D.build_model(rules=rules, ncols=ncols)
  model.agents[51].status="1"
  runs = 10
  data = CA1D.ca_run(model, runs);
  @test size(data) == (ncols * (runs + 1), 4)
end

function dummyplot(data; nodesize=2.0, anim=2)
end

@testset "Cellular automata 2D" begin
  using Agents.CA2D
  rules = (2,3,3)
  dims=(10, 10)
  model = CA2D.build_model(rules=rules, dims=dims, Moore=true)
  for i in 1:nv(model)
    if rand() < 0.1
      model.agents[i].status="1"
    end
  end

  runs = 2
  anim = CA2D.ca_run(model, runs, dummyplot);

  @test anim == nothing
end


text:
# common utilities

## convenient type alias
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}
#

typealias RealArray{T<:Real,N} AbstractArray{T,N}
typealias RealVector{T<:Real} AbstractArray{T,1}
typealias RealMatrix{T<:Real} AbstractArray{T,2}

typealias IntegerArray{T<:Integer,N} AbstractArray{T,N}
typealias IntegerVector{T<:Integer} AbstractArray{T,1}
typealias IntegerMatrix{T<:Integer} AbstractArray{T,2}

typealias RealFP Union(Float32, Float64)

## conversion from real to fp types

fptype{T<:Union(Float32,Bool,Int8,Uint8,Int16,Uint16)}(::Type{T}) = Float32
fptype{T<:Union(Float64,Int32,Uint32,Int64,Uint64,Int128,Uint128)}(::Type{T}) = Float64
fptype(::Type{Complex64}) = Complex64
fptype(::Type{Complex128}) = Complex128


text:
using Base: @kwdef

@kwdef struct Rheology
  G::Float64 = 30e9
  ν::Float64 = 0.3
  μ::Float64 = 0.6 # Friction coef
  β::Float64 = 0.1 # Correction factor
  K₁c::Float64 = 1.74e6 # Critical stress intensity factor (Pa.m^(1/2))
  a::Float64 = 1e-3 # Initial flaw size (m)
  ψ::Float64 = atand(0.6)# crack angle to the principal stress (radians)
  D₀::Float64 = 0.1# Initial flaw density
  n::Float64 = 34.0 # Stress corrosion index
  l̇₀::Float64 = 0.24 # Ref. crack growth rate (m/s)
  H::Float64 = 50e3 # Activation enthalpy (J/mol)
  A::Float64 = 5.71 # Preexponential factor (m/s)
end

function Rheology(r::Rheology, kw::NamedTuple)
  values_dict = Dict{Symbol,Float64}()
  for sym in propertynames(r)
    values_dict[sym] = isdefined(kw,sym) ? getproperty(kw,sym) : getproperty(r,sym)
  end
  return Rheology(; G = values_dict[:G],
                    ν = values_dict[:ν],
                    μ = values_dict[:μ],
                    β = values_dict[:β],
                    K₁c = values_dict[:K₁c],
                    a = values_dict[:a],
                    ψ = values_dict[:ψ],
                    D₀ = values_dict[:D₀],
                    n = values_dict[:n],
                    l̇₀ = values_dict[:l̇₀],
                    H = values_dict[:H],
                    A = values_dict[:A])
end

Base.Broadcast.broadcastable(r::Rheology) = Ref(r)

function Base.show(io::IO, ::MIME"text/plain", r::Rheology)
  print(io, "Rheology instance with fields :\n",
  "\t├── G (shear modulus)                             : $(r.G)\n",
  "\t├── ν (poisson ratio)                             : $(r.ν)\n",
  "\t├── μ (flaws friction coefficient)                : $(r.μ)\n",
  "\t├── β (correction factor)                         : $(r.β)\n",
  "\t├── K₁c (fracture toughness)                      : $(r.K₁c)\n",
  "\t├── a (flaws radius)                              : $(r.a)\n",
  "\t├── ψ (flaws angle wrt σ₁ in degree)              : $(r.ψ)\n",
  "\t├── D₀ (initial damage, constrains flaws density) : $(r.D₀)\n",
  "\t├── n (stress corrosion index)                    : $(r.n)\n",
  "\t├── l̇₀ (ref crack growth rate in m/s)             : $(r.l̇₀)\n",
  "\t├── H (Activation enthalpy (J/mol))               : $(r.H)\n",
  "\t└── A (preexponetial factor (m/s))                : $(r.A)\n")
end

Maybe(T)=Union{T,Nothing}

@kwdef struct SolverParams{TI<:Maybe(Integer),TE<:Union{Real,NamedTuple}}
  newton_abstol::Float64 = 1e-12
  newton_maxiter::Int = 100
  time_maxiter::TI = nothing
  e₀::TE = 1e-12
  adaptative_maxrecursions::Int = 50
end

@kwdef struct OutputParams{T1<:Maybe(Real),T2<:Maybe(Integer),T3<:Maybe(Real),T4<:Maybe(Integer)}
  print_frequency::T1=nothing
  print_period::T2=1
  save_frequency::T3=nothing
  save_period::T4=1
  last_tsim_printed::Ref{Float64}=Ref(0.0)
  last_tsim_saved::Ref{Float64}=Ref(0.0)
end

mutable struct Flags
  print::Bool
  save::Bool
  bifurcation::Bool
  acceptable_error::Bool
  nan::Bool
  nan1::Bool
  nan2::Bool
  Flags()=new(false,false,false,false,false,false,false)
end

struct Params{T1,T2,T3,T4,TI,TE}
  solver::SolverParams{TI,TE}
  output::OutputParams{T1,T2,T3,T4}
  flags::Flags
end
Params(s::SolverParams,o::OutputParams) = Params(s,o,Flags())
Params(o::OutputParams,s::SolverParams) = Params(s,o,Flags())
Params() = Params(SolverParams(),OutputParams(),Flags())

@kwdef mutable struct DiffEqParams{R<:Rheology,P<:Params,DU<:AbstractVector{<:Real},A<:AbstractVector{<:Real}}
  r::R = Rheology()
  mp::P = Params()
  du::DU = nothing
  scalars::A = nothing
  allow_Ḋᵒ::Bool = true
  terminate_flag::Bool = false
end

# Model type ? 2points, 1point


text:
run(`git submodule update --init`)

if !isdir("build")
    mkdir("build")
end
cd("build") do
    run(`cmake ..`)
    run(`make`)
end
run(`cp build/deps.jl .`)


text:
using Documenter
using Ju

is_concrete_type_of(t) = x -> x isa Type && x <: t && !isabstracttype(x)

makedocs(
    sitename = "Ju.jl",
    format = Documenter.HTML(),
    assets = ["assets/favicon.ico"],
    analytics = "UA-132847825-1",
    modules = [Ju],
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Interfaces" => "interfaces.md",
        "Components" => "components.md",
        "Utilities" => "utilities.md"
    ]
)

deploydocs(
    repo = "github.com/Ju-jl/Ju.jl.git",
)


text:
import AlgebraicMultigrid
import DPFEHM
import FEHM
import HDF5
import ILUZero
import IterativeSolvers
import LinearAlgebra
import Preconditioners

import Base.eltype
import Base.*

function loaddir(dirname)
	coords, volumes, neighbors, areas, lengths = DPFEHM.load_uge("$dirname/full_mesh_vol_area.uge")
	areasoverlengths = areas ./ lengths
	leftnodes = FEHM.readzone("$dirname/pboundary_left_w.zone")[2][1]
	rightnodes = FEHM.readzone("$dirname/pboundary_right_e.zone")[2][1]
	dirichletnodes = [leftnodes; rightnodes]
	dirichleths = [2e6 * ones(length(leftnodes)); 1e6 * ones(length(rightnodes))]
	Qs = zeros(size(coords, 2))
	Ks = HDF5.h5read("$dirname/dfn_properties.h5", "Permeability")
	Ks2Ks_neighbors(Ks) = sqrt.((Ks[map(p->p[1], neighbors)] .* Ks[map(p->p[2], neighbors)]))#convert from permeabilities at the nodes to permeabilities connecting the nodes
	Ks_neighbors = Ks2Ks_neighbors(Ks)
	dirichleths2 = zeros(length(Qs))
	for (i, j) in enumerate(dirichletnodes)
		dirichleths2[j] = dirichleths[i]
	end
	args = (zeros(length(Qs) - length(dirichletnodes)), Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths2, Qs, ones(length(Qs)), ones(length(Qs)))
	b = -DPFEHM.groundwater_residuals(args...)
	A = DPFEHM.groundwater_h(args...)
	return A, b, (neighbors, areasoverlengths, Ks_neighbors, Qs, dirichletnodes, dirichleths)
end

preconditioners = [("AMGSA", A->AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.smoothed_aggregation(A))),
				   ("ILU", A->ILUZero.ilu0(A)),
				   ("Iden", A->LinearAlgebra.I)]
colors = ["C0", "C1", "C2", "C3", "C4", "C5", "C6", "C7"]
markers = ["", "--", "+", ":", "."]
preconditioned_solvers = [
	("CG", (A, b, pl)->begin
		x, ch = IterativeSolvers.cg(A, b; maxiter=10000, reltol=1e-12, log=true, Pl=pl)
		return x, ch.data[:resnorm]
	end),
	("GMRES", (A, b, pl)->begin
		x, ch = IterativeSolvers.gmres(A, b; maxiter=10000, reltol=1e-12, log=true, Pl=pl)
		return x, ch.data[:resnorm]
	end),
	("BiCGStab", (A, b, pl)->begin
		x, ch = IterativeSolvers.bicgstabl(A, b; max_mv_products=10000, reltol=1e-12, log=true, Pl=pl)
		return x, ch.data[:resnorm]
	end)]
unpreconditioned_solvers = [
	("AMGSA", (A, b)->begin
		ml = AlgebraicMultigrid.smoothed_aggregation(A)
		x, ch = AlgebraicMultigrid.solve(ml, b; maxiter=10000, reltol=1e-12, log=true)
		return x, ch
	end),
	("MINRES", (A, b)->begin
		x, ch = IterativeSolvers.minres(A, b; maxiter=10000, reltol=1e-12, log=true)
		return x, ch.data[:resnorm]
	end),
	("QMR", (A, b)->begin
		x, ch = IterativeSolvers.qmr(A, b; maxiter=10000, reltol=1e-12, log=true)
		return x, ch.data[:resnorm]
	end),
	("Cholesky", (A, b)->begin
		Af = LinearAlgebra.cholesky(A)
		x = Af \ b
		return x, [1.0, sqrt(sum((A * x .- b) .^ 2) / sum(b .^ 2))]
   end)]


text:
module Sudoku

include("puzzle.jl")
include("solver.jl")

end # module

max_stars_repo_name: Joshua-Robison/Sudoku

max_issues_repo_name: Joshua-Robison/Sudoku

max_forks_repo_name: Joshua-Robison/Sudoku


text:
using Genie, Genie.Renderer.Html, Stipple, StipplePlotly, Colors, ColorSchemes

Genie.config.log_requests = false

plotly_palette = ["Greys", "YlGnBu", "Greens", "YlOrRd", "Bluered", "RdBu", "Reds", "Blues", "Picnic", "Rainbow", "Portland", "Jet", "Hot", "Blackbody", "Earth", "Electric", "Viridis", "Cividis"]

function rgb(pix::RGB{Float64})
  R = round(Int, 255 * clamp(pix.r, 0.0, 1.0))
  G = round(Int, 255 * clamp(pix.g, 0.0, 1.0))
  B = round(Int, 255 * clamp(pix.b, 0.0, 1.0))
  return "rgb($R,$G,$B)"
end

# See: https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/
function ColorScale(scheme::Symbol, N = 101)
  x = permutedims(0.0:(1.0/(N - 1)):1.0)
  cs = get(colorschemes[scheme], x, :clamp)
  cs_rgb = rgb.(cs)
  return vcat(x, cs_rgb)
end

function to_plotly_standard(x)
  if isnothing(x)
    m = "null"
  elseif ismissing(x)
    m = "null"
  elseif isa(x, AbstractString)
    m = x
  elseif isinf(x)
    m = "null"
  elseif isnan(x)
    m = "null"
  else
    m = x
  end
  return m
end

pd(name) = PlotData(
  z =  to_plotly_standard.(round.([100.0 0.0  missing;
       30.30235455584154  94.9495352852204 0.7150792978758869;
       63.172349721957175 37 -72.02122607433563;
       51.19876097916769  18.676411380287483 314.15;
       87.03992589163836  88.74217734803698 62.122844059857215]; sigdigits=3)),
  zmin = 0.0,
  zmax = 100.0,
  x = ["L0", "L63", "L128", "L191", "L255"],
  y = ["L255", "L128", "L0"],
  plot = StipplePlotly.Charts.PLOT_TYPE_HEATMAP,
  name = name,
  colorscale = ColorScale(:algae),
  colorbar = ColorBar("Z-data", 18, "right"),
  hoverongaps = false,
  hoverinfo = "x+y+z"
)

pl(title) = PlotLayout(
  plot_bgcolor = "#FFFFFF",
  title = PlotLayoutTitle(text=title, font=Font(24)),
  margin_b = 25,
  margin_t = 80,
  margin_l = 80,
  margin_r = 40,
  xaxis = [PlotLayoutAxis(xy = "x", index = 1,
    title = "From",
    font = Font(18),
    ticks = "outside top",
    side = "top",
    position = 1.0,
    showline = true,
    showgrid = false,
    zeroline = false,
    mirror = "all",
    ticklabelposition = "outside top")],
  yaxis = [PlotLayoutAxis(xy = "y", index = 1,
    showline = true,
    zeroline = false,
    mirror = "all",
    showgrid = false,
    title = "To",
    font = Font(18),
    ticks = "outside",
    scaleanchor = "x",
    scaleratio = 1,
    constrain = "domain",
    constraintoward = "top")],
)

Base.@kwdef mutable struct Model <: ReactiveModel
  data::R{Vector{PlotData}} = [pd("Random 1")], READONLY
  layout::R{PlotLayout} = pl(""), READONLY
  config::R{PlotConfig} = PlotConfig(), READONLY
end

model = Stipple.init(Model())

function ui()
  page(
    vm(model), class="container", [
      plot(:data, layout = :layout, config = :config)
    ]
  ) |> html
end

route("/", ui)

up()


text:
dfg = CloudGraphsDFG{SolverParams}("localhost", 7474, "neo4j", "test",
                            "testUser", "testRobot", "testSession",
                            "Description of test session",
                            solverParams=SolverParams())

# Nuke the user
clearUser!!(dfg)
@test lsSessions(dfg) == []
@test lsRobots(dfg) == []
@test !(Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg)))
# Create sentinel nodes using shortcut
createDfgSessionIfNotExist(dfg)
@test map(s -> s.id, lsSessions(dfg)) == [Symbol(dfg.sessionId)]
@test map(s -> s.id, lsRobots(dfg)) == [Symbol(dfg.robotId)]
@test Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg))
# Test that we can call it again.
createDfgSessionIfNotExist(dfg)
# And nuke it so we can try the longer functions.
clearUser!!(dfg)

# User, robot, and session
# TODO: Make easier ways to initialize these.
user = User(Symbol(dfg.userId), "Bob Zack", "Description", Dict{Symbol, String}())
robot = Robot(Symbol(dfg.robotId), user.id, "Test robot", "Description", Dict{Symbol, String}())
session = Session(Symbol(dfg.sessionId), robot.id, user.id, "Test Session", "Description", Dict{Symbol, String}())

@test createUser(dfg, user) == user
@test createRobot(dfg, robot) == robot
@test createSession(dfg, session) == session
@test map(s -> s.id, lsSessions(dfg)) == [Symbol(dfg.sessionId)]
@test map(s -> s.id, lsRobots(dfg)) == [Symbol(dfg.robotId)]
@test Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg))

# Test errors
dfgError = deepcopy(dfg)
# User/robot/session ID's can't start with numbers and can't have spaces.
dfgError.userId = "1testNope"
user = User(Symbol(dfgError.userId), "Bob Zack", "Description", Dict{Symbol, String}())
@test_throws Exception createUser(dfgError, user)

@test getUserData(dfg) == Dict{Symbol, String}()
@test getRobotData(dfg) == Dict{Symbol, String}()
@test getSessionData(dfg) == Dict{Symbol, String}()

# User/robot/session data
user.data = Dict{Symbol, String}(:a => "Hello", :b => "Goodbye")
robot.data = Dict{Symbol, String}(:c => "Hello", :d => "Goodbye")
session.data = Dict{Symbol, String}(:e => "Hello", :f => "Goodbye")

setUserData!(dfg, user.data)
setRobotData!(dfg, robot.data)
setSessionData!(dfg, session.data)
@test getUserData(dfg) == user.data
@test getRobotData(dfg) == robot.data
@test getSessionData(dfg) == session.data

# Add some nodes.
v1 = addVariable!(dfg, :a, ContinuousScalar, labels = [:POSE])
v2 = addVariable!(dfg, :b, ContinuousScalar, labels = [:POSE])
v3 = addVariable!(dfg, :c, ContinuousScalar, labels = [:LANDMARK])
addFactor!(dfg, [:a], Prior(Normal(0,1)))
f1 = addFactor!(dfg, [:a; :b], LinearConditional(Normal(50.0,2.0)) )
f2 = addFactor!(dfg, [:b; :c], LinearConditional(Normal(50.0,2.0)) )

sessions = lsSessions(dfg)
@test map(s -> s.id, sessions) == [session.id]

dfgLocal = buildSubgraph(LightDFG, dfg, union(ls(dfg), lsf(dfg)))
# Confirm that with sentinels we still have the same graph (doesn't pull in the sentinels)
@test symdiff(ls(dfgLocal), ls(dfg)) == []
@test symdiff(lsf(dfgLocal), lsf(dfg)) == []


text:
function ArcsSet(arg0::BSPTree, arg1::jdouble)
    return ArcsSet((BSPTree, jdouble), arg0, arg1)
end

function ArcsSet(arg0::Collection, arg1::jdouble)
    return ArcsSet((Collection, jdouble), arg0, arg1)
end

function ArcsSet(arg0::jdouble)
    return ArcsSet((jdouble,), arg0)
end

function ArcsSet(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return ArcsSet((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function apply_transform(obj::AbstractRegion, arg0::Transform)
    return jcall(obj, "applyTransform", AbstractRegion, (Transform,), arg0)
end

function as_list(obj::ArcsSet)
    return jcall(obj, "asList", List, ())
end

function build_new(obj::ArcsSet, arg0::BSPTree)
    return jcall(obj, "buildNew", ArcsSet, (BSPTree,), arg0)
end

function check_point(obj::AbstractRegion, arg0::JVector)
    return jcall(obj, "checkPoint", Region_Location, (JVector,), arg0)
end

function check_point(obj::AbstractRegion, arg0::Point)
    return jcall(obj, "checkPoint", Region_Location, (Point,), arg0)
end

function contains(obj::AbstractRegion, arg0::Region)
    return jcall(obj, "contains", jboolean, (Region,), arg0)
end

function copy_self(obj::AbstractRegion)
    return jcall(obj, "copySelf", AbstractRegion, ())
end

function for_each(obj::Iterable, arg0::Consumer)
    return jcall(obj, "forEach", void, (Consumer,), arg0)
end

function get_barycenter(obj::AbstractRegion)
    return jcall(obj, "getBarycenter", Point, ())
end

function get_boundary_size(obj::AbstractRegion)
    return jcall(obj, "getBoundarySize", jdouble, ())
end

function get_size(obj::AbstractRegion)
    return jcall(obj, "getSize", jdouble, ())
end

function get_tolerance(obj::AbstractRegion)
    return jcall(obj, "getTolerance", jdouble, ())
end

function get_tree(obj::AbstractRegion, arg0::jboolean)
    return jcall(obj, "getTree", BSPTree, (jboolean,), arg0)
end

function intersection(obj::AbstractRegion, arg0::SubHyperplane)
    return jcall(obj, "intersection", SubHyperplane, (SubHyperplane,), arg0)
end

function is_empty(obj::AbstractRegion)
    return jcall(obj, "isEmpty", jboolean, ())
end

function is_empty(obj::AbstractRegion, arg0::BSPTree)
    return jcall(obj, "isEmpty", jboolean, (BSPTree,), arg0)
end

function is_full(obj::AbstractRegion)
    return jcall(obj, "isFull", jboolean, ())
end

function is_full(obj::AbstractRegion, arg0::BSPTree)
    return jcall(obj, "isFull", jboolean, (BSPTree,), arg0)
end

function iterator(obj::ArcsSet)
    return jcall(obj, "iterator", JIterator, ())
end

function project_to_boundary(obj::ArcsSet, arg0::Point)
    return jcall(obj, "projectToBoundary", BoundaryProjection, (Point,), arg0)
end

function split(obj::ArcsSet, arg0::Arc)
    return jcall(obj, "split", ArcsSet_Split, (Arc,), arg0)
end

function spliterator(obj::Iterable)
    return jcall(obj, "spliterator", Spliterator, ())
end


text:
# Test the grat circle calculation functions
using Test
using SAC

@testset "Great circles" begin
    let s = SACtr(1, 2)
        s[:evlo], s[:evla], s[:stlo], s[:stla] = 0, 0, 12, 15
        @test s[:gcarc] ≈ 19.047431084635505
        @test s[:az] ≈ 37.992268575139384
        @test s[:baz] ≈ 219.57789440455088
    end
end


text:
# Apply ∇^2 operator
#
# Global vars:
#   Npoints, LF

function op_nabla2( LF, v::Array{Float64,1} )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    #
    nabla2_v = zeros(Float64,Npoints)
    #
    for ip = 1 : Npoints
        i = LF.lin2xyz[1,ip]
        j = LF.lin2xyz[2,ip]
        k = LF.lin2xyz[3,ip]
        #
        for ii = 1 : LF.Nx
            nabla2_v[ip] = nabla2_v[ip] + LF.LFx.D2jl[ii,i] * v[LF.xyz2lin[ii,j,k]]
        end
        #
        for jj = 1 : LF.Ny
            nabla2_v[ip] = nabla2_v[ip] + LF.LFx.D2jl[jj,j] * v[LF.xyz2lin[i,jj,k]]
        end
        #
        for kk = 1 : LF.Nz
            nabla2_v[ip] = nabla2_v[ip] + LF.LFx.D2jl[kk,k] * v[LF.xyz2lin[i,j,kk]]
        end
    end
    #
    return nabla2_v
    #
end  # function


text:
export Device
export nameQueue, getLogicalDevice, getPhysicalDevice, getInstance,
        namedQueue, family, graphicsQueue,
        attachment2D, handleRef, createView

mutable struct Device
    mVkInstance::vk.VkInstance
    mPhysicalDevice::vk.VkPhysicalDevice
    mVkDevice::vk.VkDevice

    mFeatures::Vector{IFeature}
    mQueues::Dict{String, Queue}
    mPools::Dict{UInt32, vk.VkCommandPool}

    mPhyProperties::vk.VkPhysicalDeviceProperties
    mSuballocator::Suballocator

    function Device(vkInstance::vk.VkInstance,
                      features::Vector{IFeature},
          gpuSelectionStrategy::ISelectionStrategy,
                        queues::Vector{QueueRequest})
        this = new()
        this.mVkInstance = vkInstance
        this.mFeatures = features
        this.mQueues = Dict{String, Queue}()
        this.mPools = Dict{UInt32, vk.VkCommandPool}()

        this.pickPhysicalDevice(gpuSelectionStrategy)
        this.createLogicalDevice([this.mPhysicalDevice], queues)

        this.mSuballocator = Suballocator(this, this.mPhyProperties.limits.bufferImageGranularity)
        return this
    end

    function Device(vkInstance::vk.VkInstance,
                      features::Vector{IFeature},
                 groupAssembly::IGroupAssemblyStrategy,
                        queues::Vector{QueueRequest})
        this = new()
        this.mVkInstance = vkInstance
        this.mFeatures = features
        this.mQueues = Dict{String, Queue}()
        this.mPools = Dict{UInt32, vk.VkCommandPool}()

        # pickPhysicalDevice
        groups = VkExt.enumeratePhysicalDeviceGroups(this.mVkInstance)
        group = groupAssembly.assembleFrom(groups)
        this.mPhysicalDevice = Base.first(group)

        for feat in this.mFeatures
            onPhysicalDeviceSelected(feat, this.mPhysicalDevice)
        end
        this.mPhyProperties = VkExt.getProperties(this.mPhysicalDevice)

        this.createLogicalDevice(group, queues);

        this.mSuballocator = Suballocator(this, this.mPhyProperties.limits.bufferImageGranularity)
        return this
    end
end

@class Device

function getLogicalDevice(this::Device)::vk.VkDevice
    return this.mVkDevice
end

function getPhysicalDevice(this::Device)::vk.VkPhysicalDevice
    return this.mPhysicalDevice
end

function getInstance(this::Device)::vk.VkInstance
    return this.mVkInstance
end

function destroy(this::Device)
    for feat in this.mFeatures
        beforeDeviceDestruction(feat)
    end
    for pool in this.mPools
        vk.vkDestroyCommandPool(this.mVkDevice, pool.second, C_NULL)
    end
    empty!(this.mPools)
end

function pickPhysicalDevice(this::Device, gpuSelectionStrategy::ISelectionStrategy)
    devices = VkExt.enumeratePhysicalDevices(this.mVkInstance)
    isGoodDevice = function(device::vk.VkPhysicalDevice)
        for feat in this.mFeatures
            if !supportsDevice(feat, device)
                return false
            end
        end
        return true
    end

    deviceCount::UInt32 = 0
    for device in devices
        if isGoodDevice(device)
            deviceCount += 1
            devices[deviceCount] = device
        end
    end
    resize!(devices, deviceCount)
    this.mPhysicalDevice = gpuSelectionStrategy.selectFrom(devices)

    for feat in this.mFeatures
        onPhysicalDeviceSelected(feat, this.mPhysicalDevice)
    end

    this.mPhyProperties = VkExt.getProperties(this.mPhysicalDevice)
end

mutable struct FamilyInfo
    index::UInt32
    names::Vector{String}
    priorities::Vector{Float32}
    minPriority::Float32

    FamilyInfo() = new(typemax(UInt32), Vector{String}(), Vector{Float32}(), typemax(Float32))
end

function resolveQueueRequest(req::QueueRequest, families::Vector{vk.VkQueueFamilyProperties})
    if (req.index != typemax(UInt32))
        return
    end

    bestBitcount = typemax(UInt32)
    for i = 0 : length(families) - 1 #queueFamilyIndex should start from 0
        family = families[i + 1]
        if (family.queueFlags == req.flags)
            req.index = i
            return
        end

        if (family.queueFlags & req.flags > 0)
            bitcount = count_ones(family.queueFlags)
            if (bitcount < bestBitcount)
                req.index = i
                bestBitcount = bitcount
            end
        end
    end
end

function createLogicalDevice(this::Device, physicalDevices::Vector{vk.VkPhysicalDevice}, queues::Vector{QueueRequest})
    avaliableExts = VkExt.enumerateDeviceExtensionProperties(this.mPhysicalDevice)

    # Throw in the queues requested by Features
    families = VkExt.getQueueFamilyProperties(this.mPhysicalDevice)
    for feat in this.mFeatures
        fams = queueRequests(feat, families)
        append!(queues, fams)
    end

    for q in queues
        resolveQueueRequest(q, families)
    end

    # Group requested queues by family index
    # Combine info that have the same index
    familyInfoDict = Dict{UInt32, FamilyInfo}()
    for q in queues
        info = get!(familyInfoDict, q.index, FamilyInfo())
        info.index = q.index
        push!(info.names, q.name)
        push!(info.priorities, q.priority)
        info.minPriority = min(info.minPriority, q.priority)
    end

    queueCreateInfos = Vector{vk.VkDeviceQueueCreateInfo}()
    for info in values(familyInfoDict)
        queueCreateInfo = vk.VkDeviceQueueCreateInfo(
            vk.VK_STRUCTURE_TYPE_DEVICE_QUEUE_CREATE_INFO, #sType::VkStructureType
            C_NULL, #pNext::Ptr{Cvoid}
            0, #flags::VkDeviceQueueCreateFlags
            info.index, #queueFamilyIndex::UInt32
            # queueCount should be less than or equal to available queue count
            # for this pCreateInfo->pQueueCreateInfos[0].queueFamilyIndex} (=0)
            # obtained previously from vkGetPhysicalDeviceQueueFamilyProperties
            1, #length(info.priorities), #queueCount::UInt32
            pointer([info.minPriority]) #pQueuePriorities::Ptr{Cfloat}
        )
        push!(queueCreateInfos, queueCreateInfo)
    end

    extNames = Vector{String}()
    for feat in this.mFeatures
        add = deviceExtensions(feat)
        append!(extNames, add)
    end

    deviceFatures = VkExt.VkPhysicalDeviceFeatures()
    VkExt.setSamplerAnisotropy(deviceFatures, VkExt.VK_TRUE)
    for feat in this.mFeatures
        addPhysicalDeviceFeatures(feat, deviceFatures)
    end

    # TODO
    # features::IFeature::NextPtr features2 = nullptr
    # for (auto &&feat : mFeatures)
    #     feat->addPhysicalDeviceFeatures2(features2)
    #
    # if (features2) {
    #     // Some extension needs features from the *Features2KHR extension.
    #     // The specific features go in the next-field of the extension struct.

    #     info.setPEnabledFeatures(nullptr)

    #     vk::PhysicalDeviceFeatures2KHR wrapper{features}
    #     wrapper.pNext = features2
    #     info.pNext = &wrapper
    #     mVkDevice = mPhysicalDevice.createDeviceUnique(info)
    # } else {
        # Only level 1 extensions are used, no weird magic here.

        pNext = C_NULL
        phyDeviceCount = length(physicalDevices)
        if (phyDeviceCount > 1)
            groupInfo = vk.VkDeviceGroupDeviceCreateInfo(
                vk.VK_STRUCTURE_TYPE_DEVICE_GROUP_DEVICE_CREATE_INFO, #sType::VkStructureType
                C_NULL, #pNext::Ptr{Cvoid}
                phyDeviceCount, #physicalDeviceCount::UInt32
                pointer(physicalDevices) #pPhysicalDevices::Ptr{VkPhysicalDevice}
                )
            pNext = ref_to_pointer(Ref(groupInfo))
        end
        enabledFeaturesRef = Ref(VkExt.build(deviceFatures))
        createInfo = vk.VkDeviceCreateInfo(
            vk.VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO, #sType::VkStructureType
            pNext, #pNext::Ptr{Cvoid}
            0, #flags::VkDeviceCreateFlags
            length(queueCreateInfos), #queueCreateInfoCount::UInt32
            pointer(queueCreateInfos), #pQueueCreateInfos::Ptr{VkDeviceQueueCreateInfo}
            0, #enabledLayerCount::UInt32, VilidationLayers is disabled
            C_NULL, #ppEnabledLayerNames::Ptr{Cstring}
            length(extNames), #enabledExtensionCount::UInt32
            strings2pp(extNames), #ppEnabledExtensionNames::Ptr{Cstring}
            Base.unsafe_convert(Ptr{vk.VkPhysicalDeviceFeatures}, enabledFeaturesRef)
        )

        GC.@preserve enabledFeaturesRef createInfo begin
            this.mVkDevice = VkExt.createDevice(this.mPhysicalDevice, Ref(createInfo))
        end
    #}

    for pair in familyInfoDict
        info = pair.second
        pool = VkExt.createCommandPool(this.mVkDevice, vk.VkCommandPoolCreateInfo(
            vk.VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO, #sType::VkStructureType
            C_NULL, #pNext::Ptr{Cvoid}
            vk.VK_COMMAND_POOL_CREATE_RESET_COMMAND_BUFFER_BIT |
            vk.VK_COMMAND_POOL_CREATE_TRANSIENT_BIT, #flags::VkCommandPoolCreateFlags
            info.index #queueFamilyIndex::UInt32
        ))
        this.mPools[info.index] = pool
        for i::UInt32 = 1 : length(info.names)
            name = info.names[i]
            family = info.index
            queue = Ref{vk.VkQueue}()
            vk.vkGetDeviceQueue(this.mVkDevice, family, 0, queue)
            this.mQueues[name] = Queue(family, queue[], pool, this.mVkDevice)
        end
    end
end

function createPipelineLayout(this::Device, type::Type, descriptorSets::Vector{DescriptorSetLayout} = Vector{DescriptorSetLayout}())::PipelineLayout
    range = vk.VkPushConstantRange(
        vk.VK_SHADER_STAGE_ALL, #stageFlags::VkShaderStageFlags
        0, #offset::UInt32
        sizeof_obj(type()) #size::UInt32
    )
    return createPipelineLayout(this, [range], descriptorSets);
end

function createPipelineLayout(this::Device, constantRanges::Vector{vk.VkPushConstantRange} = Vector{vk.VkPushConstantRange}(),
                                            descriptorSets::Vector{DescriptorSetLayout} = Vector{DescriptorSetLayout}())::PipelineLayout
    return PipelineLayout(this.mVkDevice, descriptorSets, constantRanges)
end

function createDescriptorSetLayout(this::Device, info::DescriptorSetLayoutCreateInfo, poolSize::UInt32 = UInt32(4))
    return DescriptorSetLayout(this.mVkDevice, info, poolSize)
end

function createDescriptorPool(vkDevice::vk.VkDevice, info::DescriptorPoolCreateInfo)::DescriptorPool
    return DescriptorPool(vkDevice, info)
end

function createRenderPass(this::Device, info::RenderPassCreateInfo)::RenderPass
    return RenderPass(this.mVkDevice, info)
end

function createShader(this::Device, code::String)
    # TODO
end

function createShaderFromFile(this::Device, filePath::String)::ShaderModule
    file = Base.read(filePath)
    size = filesize(filePath)
    stage = identifyShader(filePath)
    return ShaderModule(this.mVkDevice, pointer(file), size, stage)
end

function namedQueue(this::Device, name::String)::Queue
    @assert haskey(this.mQueues, name) "No Queue with this name:" * name * " exists!"
    return this.mQueues[name]
end

function graphicsQueue(this::Device)::Queue
    return namedQueue(this, "graphics")
end

function transferQueue(this::Device)::Queue
    if haskey(this.mQueues, "transfer")
        return this.mQueues["transfer"]
    else
        return graphicsQueue(this)
    end
end

function getSuballocator(this::Device)::Suballocator
    return this.mSuballocator
end


text:
export  E_EXPECTED_EOF, E_UNEXPECTED_EOF, E_UNEXPECTED_CHAR, E_BAD_KEY,
        E_BAD_ESCAPE, E_BAD_CONTROL, E_LEADING_ZERO, E_BAD_NUMBER, E_BAD_HEADER,
        E_BAD_TRAILER, E_NOT_IMPLEMENTED,
        E_INVALID_OBJECT, E_INVALID_PAGE_NUMBER, E_INVALID_PAGE_LABEL,
        E_NOT_TAGGED_PDF

# The following errors may be thrown by the reader
const E_EXPECTED_EOF    = "Expected end of input"
const E_UNEXPECTED_EOF  = "Unexpected end of input"
const E_UNEXPECTED_CHAR = "Unexpected character"
const E_BAD_KEY         = "Invalid object key"
const E_BAD_ESCAPE      = "Invalid escape sequence"
const E_BAD_CONTROL     = "ASCII control character in string"
const E_BAD_NUMBER      = "Invalid number"
const E_BAD_HEADER      = "Invalid file header"
const E_BAD_TRAILER     = "Invalid file trailer"
const E_NOT_IMPLEMENTED = "Not Implemented"

const E_INVALID_OBJECT  = "Invalid Object Found"
const E_INVALID_PAGE_NUMBER = "Page number is invalid"
const E_INVALID_PAGE_LABEL  = "Page label is invalid or page label definitions are not present in the document (Table 28 PDF Specification 1.7)"
const E_NOT_TAGGED_PDF = "PDF file is not tagged"


text:
#supplyimpedance.jl: Template for plotting impedance profiles
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#glyph=set(shape=:circle)
line_Zconst = [#Constituents
	cons(:a, line = set(style=:dash, color=colorant"grey20", width=2)),
	cons(:a, line = set(style=:dashdot, color=colorant"grey20", width=2)),
	cons(:a, line = set(style=:dot, color=colorant"grey20", width=2)),
]
line_Ztot = cons(:a, line = set(style=:solid, color=:blue, width=4))

#==
===============================================================================#
function fnbuild_zsupply(data)
	_title = "Supply Impedance"
	plot = cons(:plot, nstrips=2, title=_title, legend=false,
		ystrip1 = set(scale=:lin, axislabel=LBLAX_MAGZVA, min=0, max=.5),
		ystrip2 = set(scale=:log, axislabel=LBLAX_MAGZVA, min=1e-3, max=1e3),
		xaxis = set(scale=:log, label=LBLAX_FREQ, min=1e3, max=10e9),
	)

	for (i, elem) in enumerate(data[:constituents])
		magZ = abs(elem.Z)
		line = line_Zconst[i]
		idy = value(magZ, x=elem.id𝑓)
		push!(plot,
			cons(:wfrm, magZ, line_Zconst[i], label="", strip=2),
			cons(:atext, elem.id, x=elem.id𝑓, y=idy, align=elem.align, strip=2),
		)
	end

	Z = data[:Z]
	magZ = abs(Z)
	push!(plot,
		cons(:wfrm, magZ, line_Ztot, label="", strip=1),
		cons(:wfrm, magZ, line_Ztot, label="", strip=2),
	)

	pcoll = push!(cons(:plotcoll, title=""), plot)
	return pcoll
end

#==Register plot builder
===============================================================================#
push!(plot_builders, :Zsupply => EasyPlot.EasyPlotBuilder(fnbuild_zsupply))

#Last line


text:
function backprop(net::Network,input::Vector{Float64},label::Vector{Float64})::Vector{Vector{Float64}}
    feedforward!(net,input)

    # δ^L = ∇C_a ⊙ σ'(z^L)
    output_error = net.∇C_a(net.neurons[end],label).*net.act_p(net.z[end])
    errors = [output_error]

    for i=1:(length(net.layers)-1)
        # δ^l=((w^l+1)T δ^l+1) ⊙ σ′(z^l)
        next_error = (net.weights[length(net.layers)-i+1]*errors[1]).*net.act_p(net.z[length(net.layers)-i])
        prepend!(errors,next_error)
    end

    return errors, net.neurons
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
const TESTDIR = dirname(@__FILE__)

test_println(kind, f, pad = "  ") = println(pad, "testing $(kind): `$(f)`...")

println("running TapeTests...")
t = @elapsed include(joinpath(TESTDIR, "TapeTests.jl"))
println("done (took $t seconds).")

println("running TrackedTests...")
t = @elapsed include(joinpath(TESTDIR, "TrackedTests.jl"))
println("done (took $t seconds).")

println("running MacrosTests...")
t = @elapsed include(joinpath(TESTDIR, "MacrosTests.jl"))
println("done (took $t seconds).")

println("running ScalarTests...")
t = @elapsed include(joinpath(TESTDIR, "derivatives/ScalarTests.jl"))
println("done (took $t seconds).")

println("running LinAlgTests...")
t = @elapsed include(joinpath(TESTDIR, "derivatives/LinAlgTests.jl"))
println("done (took $t seconds).")

println("running ElementwiseTests...")
t = @elapsed include(joinpath(TESTDIR, "derivatives/ElementwiseTests.jl"))
println("done (took $t seconds).")

println("running GradientTests...")
t = @elapsed include(joinpath(TESTDIR, "api/GradientTests.jl"))
println("done (took $t seconds).")

println("running JacobianTests...")
t = @elapsed include(joinpath(TESTDIR, "api/JacobianTests.jl"))
println("done (took $t seconds).")

println("running HessianTests...")
t = @elapsed include(joinpath(TESTDIR, "api/HessianTests.jl"))
println("done (took $t seconds).")

println("running ConfigTests...")
t = @elapsed include(joinpath(TESTDIR, "api/ConfigTests.jl"))
println("done (took $t seconds).")


text:
function backward(backend::GPUBackend, state::SoftmaxLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  diff = diffs[1]
  if isa(diff, CuTensorBlob)
    copy!(diff, state.softmax.blobs[1])

    data_type = eltype(diff)
    spatial_dim, channels, num = split_dims(diff, state.logistic.op_dim)
    prob_dim = channels

    x_block = round(Int, ceil(convert(Float64, num)/CUDA.THREADS_PER_BLOCK_X))
    y_block = spatial_dim

    if isa(state.logistic.weights_blob, NullBlob)
      weights = convert(Ptr{data_type}, 0)
    else
      weights = state.logistic.weights_blob.ptr.p
    end

    if data_type == Float32
      kernel = backend.mocha.softmax_loss_backward_float
    elseif data_type == Float64
      kernel = backend.mocha.softmax_loss_backward_double
    else
      error("Unsupported data type $data_type")
    end
    CUDA.launch(kernel, (x_block, y_block), (CUDA.THREADS_PER_BLOCK_X, 1),
        (diff.ptr.p, inputs[2].ptr.p, weights, num, spatial_dim, prob_dim))
    CuBLAS.scal(backend.cublas_ctx, length(diff), convert(data_type, state.layer.weight/(spatial_dim*num)),
        diff.ptr, 1)
  end
end


text:
include("figure_helpers.jl")
include("../src/helpers.jl")

## Edit point embeddings
xy = [
1.3 1.7;
2.5 1.75;
1.2 .7;
2.1 0;
3.6 1.5;
2.6 .8;
3.5 .25
]
xy = [xy[:,2] xy[:,1]]
Edges = [
1 2 3;
2 5 7;
5 6 7;
2 3 4;
4 6 7;
3 4 6;
]
radius = .23
fsize = 24
lw = 2
al = .3
gverts = collect(1:4)
bverts = collect(5:7)
p = plot_hypergraph(xy,Edges,radius,fsize,lw,al,gverts,bverts)
savefig("Figures/ToyHyper.pdf")

## Plot the affinity scores
n = size(xy,1)
Hyp = Vector{Vector{Int64}}()
for i = 1:size(Edges,1)
    push!(Hyp,Edges[i,:])
end
H = elist2incidence(Hyp,n)
classes = zeros(Int64,n)
classes[gverts] .= 1
k = 3
N = get_hyperedge_counts(H,classes)
B1,B2 = n1_n2_baselines(k,sum(classes),n-sum(classes))
Hw, Hm = affinities_from_N(N)
H1 = Hw[k,:]
H2 = Hm[k,:]
R1 = H1./B1
R2 = H2./B2

lw = 2.5
ms = 9
gfs = 15
tfs = 15
titlesize = 14
s1 = 350
s2 = 300
k = 3
color1 = :green
color2 = :blue
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity / Baseline "
p = Plots.plot([1; k],[1; 1],fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, 1:k, 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

plot!(p,1:k,R2[1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2)
plot!(p,1:k,R1[1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1)

savefig("Figures/ToyRatioPlot.pdf")

## Table of types degrees, for LaTeX code
d = zeros(Int64,n,3)
for i = 1:n
    ei = findall(x->x==1,H[:,i])
    di = length(ei)
    for j = 1:di
        J = ei[j]
        edge = findall(x->x==1,H[J,:])
        t = sum(classes[edge])
        if classes[i] == 0
            t = k-t
        end
        # @show i,t, edge
        d[i,t]+=1
    end
end

## Sum differently

print("
{\\footnotesize
\\centering
\\begin{tabular}{lc c c c c c c c}
    \\toprule
    ")
print("& ")
for i = 1:n-1
    if i < 5
        print("{\\color{dgreen} \\bf $i} & ")
    else
        print("{\\color{blue} \\bf $i} & ")
    end
end
print("{\\color{blue}\\bf $n}\\\\ \n")
print("\\midrule\n")
for j = 1:4
    if j == 4
        print("\$ d_{v}\$ & ")
        for i = 1:n-1
            print("$(sum(d[i,:])) & ")
        end
        print("$(sum(d[n,:])) \\\\\n")
    else
        print(" type-$j & ")
        for i = 1:n-1
            print("$(d[i,j]) & ")
        end
        print("$(d[n,j]) \\\\\n")
    end
end
sumd = sum(d,dims =2)

print("\\bottomrule \\\\ \n \\end{tabular} \\\\ }")

## Now print out the array of affinity and baseline scores

G = "{\\color{dgreen} \\text{G}}"
B = "{\\color{blue} \\text{B}}"
hg = ["\\textbf{h}_1($G)","\\textbf{h}_2($G)","\\textbf{h}_3($G)"]
bg = ["\\textbf{b}_1($G)","\\textbf{b}_2($G)","\\textbf{b}_3($G)"]
hb = ["\\textbf{h}_1($B)","\\textbf{h}_2($B)","\\textbf{h}_3($B)"]
bb = ["\\textbf{b}_1($B)","\\textbf{b}_2($B)","\\textbf{b}_3($B)"]

Hw = round.(H1,digits = 2)
Bw = round.(B1,digits = 2)
Hm = round.(H2,digits = 2)
Bm = round.(B2,digits = 2)
print("\\centering \n {\\footnotesize \n \$	\\begin{array}{lll} \n")
print("$(hg[1]) = $(Hw[1]) & $(hg[2]) = $(Hw[2]) & $(hg[3]) = $(Hw[3]) \\\\ \n")
print("$(bg[1]) = $(Bw[1]) & $(bg[2]) = $(Bw[2]) & $(bg[3]) = $(Bw[3]) \\\\ \n")
print("$(hb[1]) = $(Hm[1]) & $(hb[2]) = $(Hm[2]) & $(hb[3]) = $(Hm[3])\\\\ \n")
print("$(bb[1]) = $(Bm[1]) & $(bb[2]) = $(Bm[2]) & $(bb[3]) = $(Bm[3])\n")
print("\\end{array} \$ }")

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily