text:
struct LoadBus <: Bus
    id::Int # for internal use only, don't mess with this.
    coords::LatLon # latitude and longitude
    load::Real #in MW
    voltage::Real #in kV
    population::Integer
    connected_to::Set{Bus}
    connections::Set{TransLine}

    function LoadBus(id, coords, load, volt, pop, connected_to = Set(), connections = Set())
        return new(id, coords, load, volt, pop, connected_to, connections)
    end
end

function LoadBus{T <: Real}(
        id,
        coords::Tuple{T,T},
        load,
        volt,
        pop,
        connected_to = Set(),
        connections = Set()
    )
    return LoadBus(
        id,
        LatLon(coords[1], coords[2]),
        load,
        volt,
        pop,
        connected_to,
        connections
    )
end

function LoadBus(coords::LatLon, load, volt, pop, connected_to = Set(), connections = Set())
    return LoadBus(-1, coords, load, volt, pop, connected_to, connections)
end

function LoadBus{T <: Real}(
        coords::Tuple{T,T},
        load,
        volt,
        pop,
        connected_to = Set(),
        connections = Set()
    )
    return LoadBus(
        -1,
        LatLon(coords[1], coords[2]),
        load,
        volt,
        pop,
        connected_to,
        connections
    )
end

function show(io::IO, bus::LoadBus)
    println(io, "LoadBus(")
    println(io, "\tid=$(bus.id),")
    println(io, "\tcoords=$(bus.coords),")
    println(io, "\tload=$(bus.load)")
    println(io, "\tvoltage=$(bus.voltage),")
    println(io, "\tpopulation=$(bus.population), ")
    if length(bus.connected_to) > 0
        println(io, "\tconnected_to=Set{Bus}(...)")
    else
        println(io, "\tconnected_to=Set{Bus}()")
    end
    if length(bus.connections) > 0
        println(io, "\tconnections=Set{TransLine}(...)")
    else
        println(io, "\tconnections=Set{Transline}()")
    end
    print(io, ")")
end

function hash(a::LoadBus, h::UInt)
    h = hash(a.id, hash(a.coords, hash(a.load, hash(a.voltage, hash(a.population, h)))))
end

struct GenBus <: Bus # this actually represents a power plant
    id::Int # for internal use only, don't mess with this.
    coords::LatLon # latitude and longitude
    generation::Real # MW
    voltage::Vector{Real} # kV
    tech_type::Vector{AbstractString} # for all generators in a plant
    connected_to::Set{Bus}
    connections::Set{TransLine}
    pfactor::Real
    summgen::Real # MW
    wintgen::Real # MW
    gens::Vector{Generator}

    function GenBus(
        id,
        coords,
        gen,
        volt,
        tech,
        connected_to = Set(),
        connections = Set(),
        pfactor = -1,
        summgen = -1,
        wintgen = -1,
        gens = []
    )
        return new(
            id,
            coords,
            gen,
            volt,
            tech,
            connected_to,
            connections,
            pfactor,
            summgen,
            wintgen,
            gens
        )
    end
end

function GenBus{T <: Real}(
    id,
    coords::Tuple{T,T},
    gen,
    volt,
    tech,
    connected_to = Set(),
    connections = Set(),
    pfactor = -1,
    summgen = -1,
    wintgen = -1,
    gens = []
)
    return GenBus(
        id,
        LatLon(coords[1], coords[2]),
        gen,
        volt,
        tech,
        connected_to,
        connections,
        pfactor,
        summgen,
        wintgen,
        gens
    )
end

function GenBus(
    coords::LatLon,
    gen,
    volt,
    tech,
    connected_to = Set(),
    connections = Set(),
    pfactor = -1,
    summgen = -1,
    wintgen = -1,
    gens = []
)
    return GenBus(
        -1,
        coords,
        gen,
        volt,
        tech,
        connected_to = Set(),
        connections = Set(),
        pfactor,
        summgen,
        wintgen,
        gens
    )
end

function GenBus{T <: Real}(
    coords::Tuple{T,T},
    gen,
    volt,
    tech,
    connected_to = Set(),
    connections = Set(),
    pfactor = -1,
    summgen = -1,
    wintgen = -1,
    gens = []
)
    return GenBus{T <: Real}(
        -1,
        coords::Tuple{T,T},
        gen,
        volt,
        tech,
        connected_to,
        connections,
        pfactor,
        summgem,
        wintgen,
        gens = []
    )
end

function GenBus(
    id,
    plant::AbstractArray{Generator},
    on::AbstractArray{Generator},
    pfactorfunc
)
    return GenBus(
        id,
        plant[1].coords,
        sum([g.cap for g in on]),
        plant[1].volt,
        [t for t in Set([g.tech for g in plant])],
        Set(),
        Set(),
        pfactorfunc(plant),
        sum([g.scap for g in on]),
        sum([g.wcap for g in on]),
        plant,
    )
end

function GenBus(id, plant::AbstractArray{Generator}, pfactorfunc)
    on = filter(g -> g.status in ["OP", "SP"], plant)
    return GenBus(id, plant, on, pfactorfunc)
end

function hash(a::GenBus, h::UInt)
    h = hash(a.id, h)
    h = hash(a.gens, h)
    h = hash(a.wintgen, h)
    h = hash(a.summgen, h)
    h = hash(a.pfactor, h)
    h = hash(a.tech_type, h)
    h = hash(a.voltage, h)
    h = hash(a.generation, h)
    return hash(a.coords, h)
end

function show(io::IO, bus::GenBus)
    println(io, "GenBus(")
    println(io, "\tid=$(bus.id)")
    println(io, "\tcoords=$(bus.coords),")
    println(io, "\tgeneration=$(bus.generation)")
    println(io, "\tvoltage=$(bus.voltage),")
    println(io, "\ttech_type=$(bus.tech_type), ")
    if length(bus.connected_to) > 0
        println(io, "\tconnected_to=Set{Bus}(...)")
    else
        println(io, "\tconnected_to=Set{Bus}()")
    end
    if length(bus.connections) > 0
        println(io, "\tconnections=Set{TransLine}(...)")
    else
        println(io, "\tconnections=Set{Transline}()")
    end
    println(io, "\tpfactor=$(bus.pfactor)")
    println(io, "\tsummgen=$(bus.summgen)")
    println(io, "\twintgen=$(bus.wintgen)")
    println(io, "\tgens=$(bus.gens)")
    print(io, ")")
end

function ==(a::Union{LoadBus, GenBus}, b::Union{LoadBus, GenBus})
    typeof(a) != typeof(b) && return false
    compared = Set{Tuple{UInt, UInt}}()

    bus_equal(a::GenBus, b::LoadBus) = false
    bus_equal(a::LoadBus, b::GenBus) = false

    function bus_equal(a::LoadBus, b::LoadBus)
        a.id == b.id &&
        a.coords == b.coords &&
        a.load == b.load &&
        a.voltage == b.voltage &&
        a.population == b.population &&
        hash(a.connections) == hash(b.connections) &&
        length(a.connected_to) == length(b.connected_to) || return false

        pair = (object_id(a), object_id(b))
        if !(pair in compared)
            push!(compared, pair)

            all(bus_equal(a, b) for (a, b) in zip(a.connected_to, b.connected_to)) ||
            return false
        end
        return true
    end

    function bus_equal(a::GenBus, b::GenBus)
        a.id == b.id &&
        a.coords == b.coords &&
        a.generation == b.generation &&
        a.voltage == b.voltage &&
        a.tech_type == b.tech_type &&
        hash(a.connections) == hash(b.connections) &&
        length(a.connected_to) == length(b.connected_to) &&
        a.pfactor == b.pfactor &&
        a.summgen == b.summgen &&
        a.wintgen == b.wintgen &&
        a.gens == b.gens || return false

        pair = (object_id(a), object_id(b))
        if !(pair in compared)
            push!(compared, pair)

            all(bus_equal(a, b) for (a, b) in zip(a.connected_to, b.connected_to)) ||
            return false
        end
        return true
    end

    return bus_equal(a, b)
end

capacity(b::GenBus) = b.generation
technology(b::GenBus) = b.tech_type
power_factor(b::GenBus) = b.pfactor
summer_capacity(b::GenBus) = b.summgen
winter_capacity(b::GenBus) = b.wintgen
generators(b::GenBus) = b.gens
trans_lines(b::Union{LoadBus, GenBus}) = b.connections
coords(b::Bus) = b.coords
load(b::LoadBus) = b.load
voltage(b::Bus) = b.voltage
population(b::LoadBus) = b.population
neighbors(b::Bus) = b.connected_to


text:
"""

    euclideanDistances(dataMatrix)

Calculate Euclidean distances between pairs.

# Arguments
- `dataMatrix::Array{Float64, 1}`: Data matrix with dimensions n x p, where n is the number of observations and p is the number of variables.

# Notes
    This is the helper function for the dataimage() function defined in Marchette & Solka (2003).

# References
Marchette, David J., and Jeffrey L. Solka. "Using data images for outlier detection."
Computational Statistics & Data Analysis 43.4 (2003): 541-552.
"""
function euclideanDistances(dataMatrix::Array{Float64,2})::Array{Float64,2}
    n, p = size(dataMatrix)
    d = zeros(Float64, n, n)
    for i in 1:n
        for j in i:n
            if i != j
                @inbounds d[i, j] = sqrt(sum((dataMatrix[i,:] .- dataMatrix[j,:]).^2.0))
                @inbounds d[j, i] = d[i, j]
            end
        end
    end
    return d
end

"""

    mahalanobisBetweenPairs(dataMatrix)

Calculate Mahalanobis distances between pairs.

# Arguments
- `dataMatrix::Array{Float64, 1}`: Data matrix with dimensions n x p, where n is the number of observations and p is the number of variables.

# Notes
    Differently from Mahalabonis distances, this function calculates Mahalanobis distances between
    pairs, rather than the distances to center of the data. This is the helper function for the
    dataimage() function defined in Marchette & Solka (2003).

# References
Marchette, David J., and Jeffrey L. Solka. "Using data images for outlier detection."
Computational Statistics & Data Analysis 43.4 (2003): 541-552.
"""
function mahalanobisBetweenPairs(dataMatrix::Array{Float64,2})::Array{Float64,2}
    n, p = size(dataMatrix)
    d = zeros(Float64, n, n)
    covmat = cov(dataMatrix)
    if det(covmat) == 0.0
        @warn "Covariance matrix is singular, mahalanobis distances can not be calculated."
    end
    covinv = inv(covmat)
    for i in 1:n
        for j in i:n
            if i != j
                @inbounds d[i, j] = sqrt((dataMatrix[i,:] .- dataMatrix[j,:]) * covinv * (dataMatrix[i,:] .- dataMatrix[j,:])')
                @inbounds d[j, i] = d[i, j]
            end
        end
    end
    return d
end

"""

    dataimage(dataMatrix; distance = :mahalanobis)

Generate the Marchette & Solka (2003) data image for a given data matrix.

# Arguments
- `dataMatrix::Array{Float64, 1}`: Data matrix with dimensions n x p, where n is the number of observations and p is the number of variables.
- `distance::Symbol`: Optional argument for the distance function.

# Notes
    distance is :mahalanobis by default, for the Mahalanobis distances.
    use

        dataimage(mat, distance = :euclidean)

    to use Euclidean distances.

# Examples
```julia-repl
julia> x1 = hbk[:,"x1"];
julia> x2 = hbk[:,"x2"];
julia> x3 = hbk[:,"x3"];
julia> mat = hcat(x1, x2, x3);
julia> dataimage(mat)
```

# References
Marchette, David J., and Jeffrey L. Solka. "Using data images for outlier detection."
Computational Statistics & Data Analysis 43.4 (2003): 541-552.
"""
function dataimage(dataMatrix::Array{Float64,2}; distance=:mahalanobis)
    d = nothing
    if distance == :mahalanobis
        d = mahalanobisSquaredBetweenPairs(dataMatrix)
    elseif distance == :euclidean
        d = euclideanDistances(dataMatrix)
    else
        @error "Distance function unknown: " distance
        @error "Using mahalanobis instead"
        d = mahalanobisSquaredBetweenPairs(dataMatrix)
    end
    colours = 1.0 .- d / maximum(d)
    n, _ = size(d)
    colormatrix = Array{RGB{Float64}}(undef, n, n)
    for i in 1:n
        for j in 1:n
            @inbounds colormatrix[i, j] = RGB(colours[i, j])
        end
    end
    plot(colormatrix)
end

ext: jl

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
"""
    child_size(::Type{T}, args...) where {T}

The number of bytes needed to allocate children of `T`, not including `self_size(T)`.

Defaults to 0.
"""
child_size(::Type{T}) where T = 0

"""
    init(blob::Blob{T}, args...) where T

Initialize `blob`.

Assumes that `blob` it at least `self_size(T) + child_size(T, args...)` bytes long.
"""
function init(blob::Blob{T}, args...) where T
    init(blob, Blob{Nothing}(blob + self_size(T)), args...)
end

"""
    init(blob::Blob{T}, free::Blob{Nothing}, args...)::Blob{Nothing} where T

Initialize `blob`, where `free` is the beginning of the remaining free space. Must return `free + child_size(T, args...)`.

The default implementation where `child_size(T) == 0` does nothing. Override this method to add custom initializers for your types.
"""
function init(blob::Blob{T}, free::Blob{Nothing}) where T
    @assert child_size(T) == 0 "Default init cannot be used for types for which child_size(T) != 0"
    # TODO should we zero memory?
    free
end

child_size(::Type{Blob{T}}, args...) where T = self_size(T) + child_size(T, args...)

function init(blob::Blob{Blob{T}}, free::Blob{Nothing}, args...) where T
    nested_blob = Blob{T}(free)
    blob[] = nested_blob
    init(nested_blob, free + self_size(T), args...)
end

child_size(::Type{BlobVector{T}}, length::Int64) where {T} = self_size(T) * length

function init(blob::Blob{BlobVector{T}}, free::Blob{Nothing}, length::Int64) where T
    blob.data[] = Blob{T}(free)
    blob.length[] = length
    free + child_size(BlobVector{T}, length)
end

child_size(::Type{BlobBitVector}, length::Int64) = self_size(UInt64) * Int64(ceil(length / 64))

function init(blob::Blob{BlobBitVector}, free::Blob{Nothing}, length::Int64)
    blob.data[] = Blob{UInt64}(free)
    blob.length[] = length
    free + child_size(BlobBitVector, length)
end

child_size(::Type{BlobString}, length::Int64) = length

function init(blob::Blob{BlobString}, free::Blob{Nothing}, length::Int64)
    blob.data[] = Blob{UInt8}(free)
    blob.len[] = length
    free + child_size(BlobString, length)
end

child_size(::Type{BlobString}, string::Union{String, BlobString}) = sizeof(string)

function init(blob::Blob{BlobString}, free::Blob{Nothing}, string::Union{String, BlobString})
    free = init(blob, free, sizeof(string))
    unsafe_copyto!(blob[], string)
    free
end

"""
    malloc(::Type{T}, args...)::Blob{T} where T

Allocate an uninitialized `Blob{T}`.
"""
function malloc(::Type{T}, args...)::Blob{T} where T
    size = self_size(T) + child_size(T, args...)
    Blob{T}(Libc.malloc(size), 0, size)
end

"""
    malloc_and_init(::Type{T}, args...)::Blob{T} where T

Allocate and initialize a new `Blob{T}`.
"""
function malloc_and_init(::Type{T}, args...)::Blob{T} where T
    size = self_size(T) + child_size(T, args...)
    blob = Blob{T}(Libc.malloc(size), 0, size)
    used = init(blob, args...)
    @assert used - blob == size
    blob
end

"""
    free(blob::Blob)

Free the underlying allocation for `blob`.
"""
function free(blob::Blob)
    Libc.free(getfield(blob, :base))
end


text:
using AES
using MD5
using JSON

function padMessage(message)
    cur_nr_of_blocks = length(message)÷16
    nr_of_fill_zeros = (cur_nr_of_blocks+1)*16-length(message)
    fill_zeros = zeros(UInt8, nr_of_fill_zeros)
    return vcat(message, fill_zeros)
end

function encryptMessage(p::Player, mode::Int, message::String, cryptomaterial::Array{Array{UInt8,1},1})
    message = Array{UInt8, 1}(message)
    #println("Prepad message: ", message)
    #println("Prepad message length: ", length(message))
    message = padMessage(message)
    #println("Postpad message: ", message)
    #println("Postpad message length: ", length(message))
    if mode == 1
        AESCBC(message, cryptomaterial[1], cryptomaterial[2], true)
    elseif mode == 2
        AESCFB(message, cryptomaterial[1], cryptomaterial[2], true)
    elseif mode == 3
		AESOFB(message, cryptomaterial[1], cryptomaterial[2])
    elseif mode == 4
        AESECB(message, cryptomaterial[1], true)
    else
		AESCTR(message, cryptomaterial[1], cryptomaterial[2])
    end

end

function generate_cryptomaterial(p::Player)
    cryptomaterial = Array{Array{UInt8,1}, 1}(undef, 2)

	global aesSeed
    key = md5(aesSeed)
    cryptomaterial[2] = rand(UInt8, 16)
	cryptomaterial[2][9] = UInt8(0x00)
    cryptomaterial[1] = vcat(key, md5(string(p.dimension)))
    return cryptomaterial
end

function sendSecret(p::Player, mode::Int, customMessage::String, tokenize::Bool=false)

	if !tokenize
		print_dict(p, "cryptomat_sender_1")
	end

	if customMessage == ""
		if mode == 1
			codetype = "cheesecode"
		elseif mode == 3
			codetype = "bombcode"
		end

		path = string("data/.", codetype, "_", p.dimension)
		if !(isfile(path))
			customMessage = string(rand(Int))
		else
			f = open_file_try(path, "r")
			customMessage = read(f, String)
			close(f)
		end
	end

	if tokenize
		messages = [customMessage]
	else
	    messages = ["ATOM:-:BOMB:-:CODE:-:START::SUPER:-:SAFE:-:CRYPTOMAT:-:PROTOCOL",
	                customMessage,
	                "ATOM:-:BOMB:-:CODE:-:END::SUPER:-:SAFE:-:CRYPTOMAT:-:PROTOCOL"]
	end
    cryptomaterial = generate_cryptomaterial(p::Player)

    for cur_message in messages
        #println("\n", cur_message)
        enc_Msg = encryptMessage(p, mode, cur_message, cryptomaterial)

		crypto_json = JSON.json(cryptomaterial[1])
		out = Pipe()
		try
			run(pipeline(`./cryptomat/rsa.py $crypto_json`, stdout=out))
		catch
			close(out.in)
			write(p.socket, "Error\n")
			return
		end
		close(out.in)
		enc_key = JSON.parse(String(read(out)), inttype=UInt8)

        full_Msg = [enc_Msg, cryptomaterial[2], enc_key]
		if tokenize
			return JSON.json(full_Msg)
		end
		write(p.socket, "Message:\n$(JSON.json(full_Msg))\n")
    end

    write(p.socket, "Sending finished\n")

end


text:
function computeStepByCrankNicholson(A,B,u,dt,f1,f2)
    LHS = B + 0.5*dt*A
    RHS = (B - 0.5*dt*A)*u + 0.5*dt*(f1+f2)

    return LHS\RHS
end

function computeStepByBackwardEuler(A,B,u,dt,f)
    LHS = B + dt*A
    RHS = B*u + dt*f

    return LHS\RHS
end

function computeStepByForwardEuler(A,B,u,dt,f)
    LHS = B
    RHS = (B-dt*A)*u+ dt*f

    return LHS\RHS
end

function initializeChemical(sMesh,CHEMICAL::Symbol)
    ψinit = zeros(length(sMesh.xy))

    if CHEMICAL == :ψA
        ENTRANCE = :lower
    elseif CHEMICAL == :ψB
        ENTRANCE = :upper
    elseif CHEMICAL == :ψC
        ENTRANCE = :none
    end

    for i=1:length(sMesh.xy)
        ψinit[i] = scalarBCf_ydom(sMesh.xy[i].x,sMesh.xy[i].y,ENTRANCE)
    end

    return ψinit
end

function updateChemical(M,Stiff,u,p;solver=:backwarddeuler)
    A = Stiff.Op
    B = M
    uold = u
    dt = p.Δt
    f1 = Stiff.rhs
    f2 = copy(f1)   # change this for general case w/ reaction

    if solver==:forwardeuler
        return computeStepByForwardEuler(A,B,uold,dt,f1)
    elseif solver==:backwardeuler
        return computeStepByBackwardEuler(A,B,uold,dt,f1)
    elseif solver==:cranknicholson
        # NOTE: forcing will change with time in general, so this will need to be modified
        return computeStepByCrankNicholson(A,B,uold,dt,f1,f2)
    else
        error("solver $solver not recognized")
    end
end

function computeReaction!(p::precipParam,soln::precipSolution)
    N = length(soln.ψA.u)

    for i=1:N
        ψA = soln.ψA.u[i]
        ψB = soln.ψB.u[i]
        θs = soln.θs.u[i]
        a,b,c,d = QSSRstep(ψA,ψB,θs,p.Ma,p.Mb,p.ρf,p.ρs,p.a,p.b)
        soln.ψA.u[i] = a
        soln.ψB.u[i] = b
        soln.θs.u[i] = c
        soln.θf.u[i] = d
    end

    nothing
end

function QSSRstep(ψA::T,ψB::T,θs::T,Ma::T,Mb::T,ρf::T,ρs::T,a,b) where T<:Real
    θf = one(T) - θs
    J = (ρs-ρf)*θs + (Ma*ψA + Mb*ψB)*θf
    if isAlimitingReagent(ψA,ψB,a,b)
        ψAnew = zero(T)
        ϕB = (ψB - b*ψA/a)*θf
        θsnew = (J - Mb*ϕB)/(ρs-ρf)
        ψBnew = ϕB/(1-θsnew)
    else
        ψBnew = zero(T)
        ϕA = (ψA - a*ψB/b)*θf
        θsnew = (J - Ma*ϕA)/(ρs-ρf)
        ψAnew = ϕA/(1-θsnew)
    end
    θfnew = one(T)-θsnew
    return ψAnew,ψBnew,θsnew,θfnew
end

isAlimitingReagent(ψA,ψB,a,b) = (ψA/a <= ψB/b ? true : false)
isBlimitingReagent(ψA,ψB,a,b) = !(isAlimitingReagent(ψA,ψB,a,b))

function QSSRstep_test()
    @testset "QSSR reaction" begin
        a,b,Ma,Mb,Mc,ρs,ρf = defineChemistry()
        θfT = 0.738; θsT = 1.0 - θfT
        ψAT = 0.5*1e-2
        ψBT = 0.5*1e-2
        ψCT = 0.0
        ψA = 0.0
        ψB = 0.0
        θf = 0.0

        @testset "A limiting reagent" begin
            a=2; b=1
            initialMass = computeDensity(ψAT,ψBT,ψCT,θfT,ρf,ρs,Ma,Mb,Mc)
            ψA,ψB,θs,θf = QSSRstep(ψAT,ψBT,θsT,Ma,Mb,ρf,ρs,a,b)
            finalMass = computeDensity(ψA,ψB,ψCT,θf,ρf,ρs,Ma,Mb,Mc)
            @test isAlimitingReagent(ψAT,ψBT,a,b)
            @test abs(finalMass - initialMass)<3*eps()
            @test !(θfT < 0.0)
        end

        a,b,Ma,Mb,Mc,ρs,ρf = defineChemistry()
        θf = 1.0; θs = 1.0-θf
        ψA = 0.5*1e-2
        ψB = 0.5*1e-2
        ψC = 0.0

        @testset "B limiting reagent" begin
            a=1; b=2

            initialMass = computeDensity(ψA,ψB,ψC,θf,ρf,ρs,Ma,Mb,Mc)
            ψAnew,ψBnew,θsnew,θfnew = QSSRstep(ψA,ψB,θs,Ma,Mb,ρf,ρs,a,b)
            finalMass = computeDensity(ψA,ψB,ψC,θf,ρf,ρs,Ma,Mb,Mc)

            @test isBlimitingReagent(ψA,ψB,a,b)
            @test abs(finalMass - initialMass)<3*eps()
            @test !(θfnew < 0.0)
        end
    end

    nothing
end

# assumes that ψC==0
function computeTotalMass(p::precipParam,mesh::eFEM.ScalarMesh,soln::precipSolution)
    N = length(soln.ψA.u)
    density = zeros(N)
    # compute mass density function
    for i=1:N
        ψA = soln.ψA.u[i]
        ψB = soln.ψB.u[i]
        ψC = soln.ψC.u[i]
        θf = soln.θf.u[i]
        θs = soln.θs.u[i]
        density[i] = computeDensity(ψA,ψB,ψC,θf,p.ρf,p.ρs,p.Ma,p.Mb,p.Mc)
    end

    return eFEM.DomainQuad(mesh,density)
end

function computeDensity(ψA::T,ψB::T,ψC::T,θf::T,ρf::T,ρs::T,Ma::T,Mb::T,Mc::T) where T
    return ρf*θf + ρs*(1-θf) + Ma*ψA*θf + Mb*ψB*θf + Mc*ψC*θf
end

function defineChemistry()
    ρs = 4.1        # density (g/cm^3) of Ni(OH)2
    ρf = 1.0        # density (g/cm^3) of water
    Ma = 58.693     # molar mass (g/mol) of Ni2+
    Mb = 17.008     # molar mass (g/mol) of OH-
    Mc = 92.724     # molar mass (g/mol) of (Ni)OH_2
    a = 1
    b = 2
    c = 1

    return a,b,c,Ma,Mb,Mc,ρs,ρf
end

function checkValidChemistry(Ma,Mb,Mc,a,b,c)
    lhs = a*Ma + b*Mb
    rhs = c*Mc
    if abs(lhs - rhs) > 2*eps()
        return false
    else
        return true
    end
end

function computePureAqeuousReaction(ψAn::T,ψBn::T,ψCn::T,Ma::T,Mb::T,Mc::T,a::Int,b::Int,c::Int) where T<:Real
    ψnp1 = zero(typeof(ψAn))
    if isAlimitingReagent(ψAn,ψBn,a,b)
        ψBnp1 = ψBn - (b/a)*ψAn
        ψnp1 = (Ma*ψAn + Mb*ψBn + Mc*ψCn - Mb*ψBnp1)/Mc
    elseif isBlimitingReagent(ψAn,ψBn,a,b)
        ψAnp1 = ψAn - (a/b)*ψBn
        ψnp1 = (Ma*ψAn + Mb*ψBn + Mc*ψCn - Ma*ψAnp1)/Mc
    end
    return ψnp1
end

function checkψCExceedsThreshold(ψC,ψCthreshold)
    return (ψC > ψCthreshold) ? true : false
end


text:
include("nls_solver.jl")
include("ls_solver.jl")


text:
using HTTP

r = HTTP.request("GET", "https://ipecho.net/plain")
println(String(r.body))


text:
# под вопросом. пока не подключаем

"strftime(\"%F/%H\", now())"
strftime{T<:DateTime}(fmt::AbstractString, t::T) =
    Libc.strftime(fmt,Dates.datetime2unix(t))

"strftime(\"%F\", Dates.today())"
strftime{T<:Date}(fmt::AbstractString, d::T) =
    strftime(fmt,DateTime(d))

"""
[Dates.Date(2016,4,1), Dates.Date(2016,4,2)] |>
    strftime("../%F/total.gz") |> first
"""
strftime(fmt::AbstractString) =
    i->imap(t->strftime(fmt, t), i)

max_stars_repo_name: closescreen/Many

max_issues_repo_name: closescreen/Many

max_forks_repo_name: closescreen/Many


text:
# Taken from: https://github.com/JuliaML/OpenAIGym.jl/blob/master/src/OpenAIGym.jl

# --------------------------------------------------------------

"A simple wrapper around the OpenAI gym environments to add to the Reinforce framework"
mutable struct GymEnv{T}
    name::Symbol
    ver::Symbol
    pyenv::PyObject   # the python "env" object
    pystep::PyObject  # the python env.step function
    pyreset::PyObject # the python env.reset function
    pystate::PyObject # the state array object referenced by the PyArray state.o
    pystepres::PyObject # used to make stepping the env slightly more efficient
    info::PyObject    # store it as a PyObject for speed, since often unused
    state::T
    reward::Float64
    total_reward::Float64
    actions
    done::Bool
    function GymEnv{T}(name, ver, pyenv, pystate, state) where T
        env = new{T}(name, ver, pyenv, pyenv."step", pyenv."reset",
                                 pystate, PyNULL(), PyNULL(), state)
        reset!(env)
        env
    end
end

use_pyarray_state(envname::Symbol) = !(envname ∈ (:Blackjack,))

function GymEnv(name::Symbol, ver::Symbol = :v0;
                stateT = ifelse(use_pyarray_state(name), PyArray, PyAny))
    if PyCall.ispynull(pysoccer) && name ∈ (:Soccer, :SoccerEmptyGoal)
        copy!(pysoccer, pyimport("gym_soccer"))
    end

    GymEnv(name, ver, pygym.make("$name-$ver"), stateT)
end

GymEnv(name::AbstractString; kwargs...) =
    GymEnv(Symbol.(split(name, '-', limit = 2))...; kwargs...)

function GymEnv(name::Symbol, ver::Symbol, pyenv, stateT)
    pystate = pycall(pyenv."reset", PyObject)
    state = convert(stateT, pystate)
    T = typeof(state)
    GymEnv{T}(name, ver, pyenv, pystate, state)
end

function Base.show(io::IO, env::GymEnv)
  println(io, "GymEnv $(env.name)-$(env.ver)")
  if hasproperty(env.pyenv, :class_name)
    println(io, "  $(env.pyenv.class_name())")
  end
  println(io, "  r  = $(env.reward)")
  print(  io, "  ∑r = $(env.total_reward)")
end

# --------------------------------------------------------------

"""
    close(env::GymEnv)
"""
Base.close(env::GymEnv) =
	!ispynull(env.pyenv) && env.pyenv.close()

# --------------------------------------------------------------

"""
    render(env::GymEnv; mode = :human)
# Arguments
- `mode`: `:human`, `:rgb_array`, `:ansi`
"""
render(env::GymEnv, args...; kwargs...) =
    pycall(env.pyenv.render, PyAny; kwargs...)

# --------------------------------------------------------------

function actionset(A::PyObject)
    if hasproperty(A, :n)
        # choose from n actions
		collect(1:A.n)
        # DiscreteSet(0:A.n-1)
    elseif hasproperty(A, :spaces)
        # a tuple of action sets
        sets = [actionset(a) for a in A.spaces]
        TupleSet(sets...)
    elseif hasproperty(A, :high)
        # continuous interval
        [A.low, A.high]
        # if A[:shape] == (1,)  # for now we only support 1-length vectors
        #     IntervalSet{Float64}(A[:low][1], A[:high][1])
        # else
        #     # @show A[:shape]
        #     lo,hi = A[:low], A[:high]
        #     # error("Unsupported shape for IntervalSet: $(A[:shape])")
        #     [IntervalSet{Float64}(lo[i], hi[i]) for i=1:length(lo)]
        # end
    elseif hasproperty(A, :actions)
        # Hardcoded
        collect(1:length(A))
    else
        @show A
        @show keys(A)
        error("Unknown actionset type: $A")
    end
end

pyaction(a::Vector) = Any[pyaction(ai) for ai=a]
pyaction(a) = a

"""
`reset!(env::GymEnv)` reset the environment
"""
function reset!(env::GymEnv)
    pycall!(env.pystate, env.pyreset, PyObject)
    convert_state!(env)
    env.reward = 0.0
    env.total_reward = 0.0
    env.actions = actionset(env.pyenv.action_space)
    env.done = false
    return env.state
end

function PyObject2Dict(obj)
	ks = collect(obj.keys())
	Dict(k=>obj.get(k) for k in ks)
end

"""
    step!(env::GymEnv, a)
take a step in the enviroment
"""
function step!(env::GymEnv, a)
    pyact = pyaction(a)
    pycall!(env.pystepres, env.pystep, PyObject, pyact)

    env.pystate, r, env.done, env.info =
        convert(Tuple{PyObject, Float64, Bool, PyObject}, env.pystepres)

    convert_state!(env)

    env.total_reward += r
    return (r, env.state, PyObject2Dict(env.info))
end

@inline step!(env::GymEnv, s, a) = step!(env, a)

convert_state!(env::GymEnv{T}) where T =
    env.state = convert(T, env.pystate)

convert_state!(env::GymEnv{<:PyArray}) =
    env.state = PyArray(env.pystate)

finished(env::GymEnv)     = env.done
finished(env::GymEnv, s′) = env.done

# --------------------------------------------------------------

const pygym    = PyNULL()
const pysoccer = PyNULL()

function __init__()
    # the copy! puts the gym module into `pygym`, handling python ref-counting
    copy!(pygym, pyimport("gym"))
end

max_stars_repo_name: ancorso/POMDPGym

max_issues_repo_name: ancorso/POMDPGym

max_forks_repo_name: ancorso/POMDPGym


text:
"""
How to add a type:
- Add an abstract type with the same name but starting with an "A"
- Add the type with no arguments
- Add a funcion with the same name arguments
- export the function
- Add `_show` function
- Add `to_latex` function
"""

include("Baseline.jl")
include("SetTheory.jl")
include("Logic.jl")


text:
mutable struct Unknowns
    n_unknowns::Int
    unknown_indexes::Dict{String, Vector{CartesianIndex}}

    function Unknowns(ssm::StateSpaceModel)
        dict = Dict{String, Vector{CartesianIndex}}()
        dict["Z"] = find_nans_Z(ssm.Z)
        dict["T"] = find_nans(ssm.T)
        dict["R"] = find_nans(ssm.R)
        dict["d"] = find_nans(ssm.d)
        dict["c"] = find_nans(ssm.c)
        dict["sqrtH"] = find_nans_cov(ssm.H)
        dict["sqrtQ"] = find_nans_cov(ssm.Q)

        n_unknowns = num_unknowns(dict)

        if !isempty(dict["d"]) || !isempty(dict["c"])
            error("StateSpaceModels does not currently support estimating parameters in neither c and d.")
        end

        return new(n_unknowns, dict)
    end
end

function find_nans(mat::Matrix{T}) where T
    return findall(isnan, mat)
end

function find_nans_Z(Z::Array{T, 3}) where T
    t_dim = size(Z, 3)
    cart_idxs = Vector{Vector{CartesianIndex}}(undef, t_dim)
    for i in 1:t_dim
        cart_idxs[i] = findall(isnan, Z[:, :, i])
    end

    # Guarantee that the unknown is the same across all t
    @assert all_indexes_are_equal(cart_idxs)

    return cart_idxs[1]
end

function find_nans_cov(cov::Matrix{T}) where T
    return findall(isnan, tril(cov))
end

function all_indexes_are_equal(cart_vec::Vector{Vector{CartesianIndex}})
    for i in eachindex(cart_vec)
        if cart_vec[i] != cart_vec[1]
            return false
        end
    end
    return true
end

function num_unknowns(dict::Dict{String, Vector{CartesianIndex}})
    n = 0
    for (k, v) in dict
        n += length(v)
    end
    return n
end

function fill_model_with_parameters!(model::StateSpaceModel{Typ}, psitilde::Vector{Typ}, unknowns::Unknowns) where Typ
    # The order to fill the model is Z => T => R => sqrtH => sqrtQ
    offset = 0

    # Fill Z
    for un_Z in unknowns.unknown_indexes["Z"]
        offset += 1
        for t in axes(model.Z, 3)
            model.Z[un_Z, t] = psitilde[offset]
        end
    end

    # Fill T
    for un_T in unknowns.unknown_indexes["T"]
        offset += 1
        model.T[un_T] = psitilde[offset]
    end

    # Fill R
    for un_R in unknowns.unknown_indexes["R"]
        offset += 1
        model.R[un_R] = psitilde[offset]
    end

    # Fill d
    for un_d in unknowns.unknown_indexes["d"]
        offset += 1
        model.d[un_d] = psitilde[offset]
    end

    # Fill c
    for un_c in unknowns.unknown_indexes["c"]
        offset += 1
        model.c[un_c] = psitilde[offset]
    end

    # Fill sqrtH in the H matrix and do a gram (XX') after
    for un_H in unknowns.unknown_indexes["sqrtH"]
        offset += 1
        model.H[un_H] = psitilde[offset]
    end
    !isempty(unknowns.unknown_indexes["sqrtH"]) && fill_H!(model, gram(tril(model.H)))

    # Fill sqrtQ in the Q matrix and do a gram (XX') after
    for un_Q in unknowns.unknown_indexes["sqrtQ"]
        offset += 1
        model.Q[un_Q] = psitilde[offset]
    end
    !isempty(unknowns.unknown_indexes["sqrtQ"]) && fill_Q!(model, gram(tril(model.Q)))

    return
end

function fill_H!(model::StateSpaceModel{T}, H::Matrix{T}) where T
    for i in axes(H, 1), j in axes(H, 2)
        model.H[i, j] = H[i, j]
    end
end

function fill_Q!(model::StateSpaceModel{T}, Q::Matrix{T}) where T
    for i in axes(Q, 1), j in axes(Q, 2)
        model.Q[i, j] = Q[i, j]
    end
end


text:
if nprocs() < 2
    addprocs(1)
    using MessageUtils
end

function foo()
    while true
        m = take!()
#        println("recd : ", m)
        put!((:response, m))
    end
end

ct = ctask(foo)

(1,2,3) |> ct
(1,3) |> ct

remotecall_fetch(2, ct -> take!(ct), ct)
take!(ct)

ctn = ctask(foo; name="testct")

("Hello", ) |> "testct"
("World", 1, 2) |> "testct"

remotecall_fetch(2, ()->take!("testct"))
take!("testct")


text:
module RoundedNumbersForHumans

export RoundedNumber

abstract type SuffixType end
struct Scientific <: SuffixType end
struct Common <: SuffixType end

const suffixlistcommon = Dict(
    -1=>"m", # one-thousandth
    -2=>"μ", # one-millionth
    -3=>"n", # one-billionth
    -4=>"p", # one-trillionth
    -5=>"f", # one-quadrillionth
    -6=>"a", # one-quintillionth
    1=>"k", # kilo
    2=>"M", # millions
    3=>"B", # billions
    4=>"T", # trillions
    5=>"Qd", # quadrillions
    6=>"Qn", # quintillions
    )

const suffixlistscientific = Dict(
    -1=>"m", # milli
    -2=>"μ", # micro
    -3=>"n", # nano
    -4=>"p", # pico
    -5=>"f", # femto
    -6=>"a", # atto
    -7=>"z", # zepto
    -8=>"y", # yocto
    1=>"k", # kilo
    2=>"M", # mega
    3=>"G", # giga
    4=>"T", # tera
    5=>"P", # peta
    6=>"E", # exa
    7=>"Z", # zetta
    8=>"Y", # yotta
    )

suffixlist(::Type{Common}) = suffixlistcommon
suffixlist(::Type{Scientific}) = suffixlistscientific
suffixlist(s::SuffixType) = suffixlist(typeof(s))

struct RoundedNumber{S<:SuffixType,T} <: Number
    n :: T
    exp1000 :: Int
end

RoundedNumber{S}(n, exp1000) where {S<:SuffixType} = RoundedNumber{S, typeof(n)}(n, exp1000)

exponent1000(n) = iszero(n) ? 0 : floor(Int, log(1000,abs(n)))
minexp1000(s::SuffixType) = minimum(keys(suffixlist(s)))
maxexp1000(s::SuffixType) = maximum(keys(suffixlist(s)))

function RoundedNumber(n, suffixtype::SuffixType = Common())
    exp1000 = exponent1000(n)
    exp1000 = min(exp1000, maxexp1000(suffixtype))
    exp1000 = max(exp1000, minexp1000(suffixtype))
    RoundedNumber{typeof(suffixtype)}(n, exp1000)
end

prefixnumdigits(r::RoundedNumber) = prefixnumdigits(r.n,r.exp1000)
function prefixnumdigits(n, exp1000)
    ceil(Int,log10(n)) - 3exp1000
end

function numberprefix(r::RoundedNumber; kwargs...)
    numberprefix(r.n, r.exp1000; kwargs...)
end
function numberprefix(n, exp1000; sigdigits=3)
    y = sign(n) * round(10^(log10(abs(n)) - 3exp1000), sigdigits=sigdigits)
    maybeInt(y)
end

maybeInt(y) = isinteger(y) ? Integer(y) : y

function Base.string(r::RoundedNumber{S}; kwargs...) where {S}
    if r.exp1000 != 0
        suffix = suffixlist(S)[r.exp1000]
    else
        suffix = ""
    end
    prefixnum = numberprefix(r.n, r.exp1000; kwargs...)
    prefix = string(prefixnum)
    prefix*suffix
end

(::Type{T})(r::RoundedNumber) where {T<:Number} = T(r.n)

Base.show(io::IO, r::RoundedNumber) = print(io, string(r))

end # module


text:
# Copyright (c) 2012: John Myles White and other contributors.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

@testset "Test simple function" begin

    function f_lm!(x, r)
        r[1] = x[1]
        r[2] = 2.0 - x[2]
    end

    function g_lm!(x, j)
        j[1, 1] = 1.0
        j[1, 2] = 0.0
        j[2, 1] = 0.0
        j[2, 2] = -1.0
    end

    initial_x = [100.0, 100.0]

    m = 2

    results = ParallelTrainingNN.levenberg_marquardt(f_lm!, g_lm!, initial_x, m)
    @test norm(Optim.minimizer(results) - [0.0, 2.0]) < 0.01
end

@testset "Test Rosenbrock" begin

    function rosenbrock_res!(x, r)
        r[1] = 10.0 * (x[2] - x[1]^2 )
        r[2] =  1.0 - x[1]
    end

    function rosenbrock_jac!(x, j)
        j[1, 1] = -20.0 * x[1]
        j[1, 2] =  10.0
        j[2, 1] =  -1.0
        j[2, 2] =   0.0
    end

    m = 2

    initial_xrb = [-1.2, 1.0]

    results = ParallelTrainingNN.levenberg_marquardt(rosenbrock_res!,
                                        rosenbrock_jac!,
                                        initial_xrb, m)

    @test norm(Optim.minimizer(results) - [1.0, 1.0]) < 0.01

    # check estimate is within the bound PR #278
    result = ParallelTrainingNN.levenberg_marquardt(rosenbrock_res!,
                                       rosenbrock_jac!,
                                       [150.0, 150.0], m;
                                       lower = [10.0, 10.0],
                                       upper = [200.0, 200.0])
    @test Optim.minimizer(result)[1] >= 10.0
    @test Optim.minimizer(result)[2] >= 10.0
end


text:
# @testset "Experimental AD"  begin

##

using ACE, ACEbase, Zygote, ChainRules
using Printf, Test, LinearAlgebra, ACE.Testing, Random
using ACE: evaluate, evaluate_d, SymmetricBasis, NaiveTotalDegree, PIBasis
using ACEbase.Testing: fdtest
using ACE: EVAL, EVAL_D
import ChainRulesCore: rrule
##

@info("Basic test of LinearACEModel construction and evaluation")

# construct the 1p-basis
D = NaiveTotalDegree()
maxdeg = 12
ord = 4

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, D = D)

# generate a configuration
nX = 10
Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

φ = ACE.Invariant()
pibasis = PIBasis(B1p, ord, maxdeg; property = φ)
basis = SymmetricBasis(pibasis, φ)

BB = evaluate(basis, cfg)
c = rand(length(BB)) .- 0.5
model = ACE.LinearACEModel(basis, c, evaluator = :standard)
@show length(basis)

##

# set_params!(model, parms); evaluate(model, cfg)
# EVAL(model, cfg)(parms)
# AD w.r.t. parms
#  vs. evaluate(model, cfg, parms)  ????
e = EVAL(model, cfg)
e(c)
Zygote.gradient(e, c)[1]
# ACEbase.Testing.fdtest(e, c -> Zygote.gradient(e, c)[1], c)

##

norm2(x::AbstractVector) = sum(abs2, x)

#  vs. evaluate_d(model, cfg, parms)  ????
de = EVAL_D(model, cfg)
de(c)
fd = c -> sum(norm2, de(c))
# fd = c -> norm2(de(c)[1])
fd(c)
Zygote.gradient(fd, c)[1]
# ACEbase.Testing.fdtest(fd, c -> Zygote.gradient(fd, c)[1], c)

##

cfgs = [ ACEConfig(rand(EuclideanVectorState, B1p.bases[1], nX))
         for _ = 1:10 ]
loss = params -> sum( ( EVAL(model, cfg)(params)
                         + sum(norm2, EVAL_D(model, cfg)(params)) )
                      for cfg in cfgs ) / length(cfgs) / nX
loss(c)
Zygote.gradient(loss, c)[1]
# ACEbase.Testing.fdtest(loss, c -> Zygote.gradient(loss, c)[1], c)

##

naive = ACE.LinearACEModel(basis, c, evaluator = :naive)
loss_naive = params -> sum( ( EVAL(naive, cfg)(params)
                           + sum(norm2, EVAL_D(naive, cfg)(params)) )
                           for cfg in cfgs ) / length(cfgs) / nX
loss_naive(c)
Zygote.gradient(loss_naive, c)[1]

##

@time loss(c)
@time loss_naive(c)
@time Zygote.gradient(loss, c)[1]
@time Zygote.gradient(loss_naive, c)[1]

##

# end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.NoRenderer

No renderer (so no animations are shown).

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module NoRenderer

import Modia3D
import Modia3D.Composition
import Modia3D.Solids

const NOTHING = Nothing

include("renderer.jl")
include("handler.jl")

end


text:
Base.@kwdef struct RenderState
    blending_mode = nothing
    enable_depth_testing::Bool = false
    enable_depth_bias::Bool = false
    enable_depth_clamp::Bool = false
    color_write_mask::Vk.ColorComponentFlag = Vk.COLOR_COMPONENT_R_BIT | Vk.COLOR_COMPONENT_G_BIT | Vk.COLOR_COMPONENT_B_BIT
end

Base.@kwdef struct ProgramInvocationState
    face_culling::Vk.CullModeFlag = Vk.CULL_MODE_BACK_BIT
    triangle_orientation::Vk.FrontFace = Vk.FRONT_FACE_CLOCKWISE
    enable_primitive_restart::Bool = false
    primitive_topology::Vk.PrimitiveTopology = Vk.PRIMITIVE_TOPOLOGY_TRIANGLE_LIST
    vertex_input_rate::Vk.VertexInputRate = Vk.VERTEX_INPUT_RATE_VERTEX
    polygon_mode::Vk.PolygonMode = Vk.POLYGON_MODE_FILL
end

"""
Set of buffer handles for loading per-material and per-vertex data, along with global camera data.
"""
struct PushConstantData
    camera_data::UInt64
    vertex_data::UInt64
    material_data::UInt64
end

PushConstantData() = PushConstantData(0, 0, 0)

struct DrawState
    render_state::RenderState
    program_state::ProgramInvocationState
    push_data::PushConstantData
end

DrawState() = DrawState(RenderState(), ProgramInvocationState(), PushConstantData())


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DiscriminativeCircuits

include("../helper/little_circuits.jl")

@testset "probabilistic circuit nodes" begin

    c1 = little_3var()
    classes = 2
    @test all(isleaf, intersect(linearize(LogisticCircuit(c1, classes)), linearize(LogisticCircuit(c1, classes))))
    p1 = LogisticCircuit(c1, classes)
    lit3 = children(children(p1)[1])[1]

    # traits
    @test p1 isa LogisticCircuit
    @test p1 isa LogisticInnerNode
    @test p1 isa Logistic⋁Node
    @test children(p1)[1] isa Logistic⋀Node
    @test lit3 isa LogisticLiteralNode
    @test is⋁gate(p1)
    @test is⋀gate(children(p1)[1])
    @test isliteralgate(lit3)
    @test length(or_nodes(p1)) == 5

    # extension methods
    @test literal(lit3) === literal(children(children(c1)[1])[1])
    @test variable(left_most_descendent(p1)) == Var(3)
    @test ispositive(left_most_descendent(p1))
    @test !isnegative(left_most_descendent(p1))
    @test num_nodes(p1) == 15
    @test num_edges(p1) == 18
    @test num_parameters(p1) == 10*2
    @test num_parameters_per_class(p1) == 10

    r1 = fully_factorized_circuit(LogisticCircuit,10; classes=classes)
    @test num_parameters(r1) == 2 * (2*10+1)
    @test length(or_nodes(r1)) == 11

end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
using LightDarkPOMDPs
using Base.Test
using POMDPs
using POMDPToolbox
using StaticArrays
using Plots
using ParticleFilters

#### LightDark2D Problem ####

p = LightDark2D()

# simulate with a policy that moves -0.01 times the observation
hr = HistoryRecorder(max_steps=10)
h = sim(p, simulator=hr) do o
    @show o
    return -0.01*o
end
@show p.count

# simulate with a kalman filter
bpol = FunctionPolicy(b -> -0.5*mean(b))
up = LightDark2DKalman(p)
h = simulate(hr, p, bpol, up)

plotly()
plot(p)
plot!(h)
# gui()

#### LightDark2DTarget Problem ####

p = LightDark2DTarget()

plot(p)
up = SIRParticleFilter(p, 1000)
b = initialize_belief(up, initial_state_distribution(p))
plot!(b)
# gui()

hist = simulate(hr, p, bpol, up)
println("state trajectory")
for s in state_hist(hist)
    @show s
end


text:
using CUDA: CuArray
using Oceananigans.Architectures: AbstractCPUArchitecture, AbstractGPUArchitecture

import Base: zeros

zeros(FT, ::AbstractCPUArchitecture, N...) = zeros(FT, N...)

function zeros(FT, ::AbstractGPUArchitecture, N...)
    a = CuArray{FT}(undef, N...)
    a .= 0
    return a
end

zeros(arch::AbstractArchitecture, grid, N...) = zeros(eltype(grid), arch, N...)


text:
module StaleC

using StaleA

StaleA.stale(x::String) = length(x)
call_buildstale(x) = StaleA.build_stale(x)

call_buildstale("hey")

end # module


text:
"""
compute_scc(m::Model=get_model(); year::Union{Int, Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta::Float64=1.5, pulse_size=1e10)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiDICE2010 model.
If no model is provided, the default model from MimiDICE2010.get_model() is used.
Ramsey discounting is used with a pure rate of time preference of `prtp` and inequality aversion of `eta`.
`pulse_size` controls the size of the marginal emission pulse.
"""
function compute_scc(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta::Float64=1.5, pulse_size=1e10)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2015)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = get_marginal_model(m; year=year, pulse_size=pulse_size)

    return _compute_scc(mm, year=year, last_year=last_year, prtp=prtp, eta=eta)
end

"""
compute_scc_mm(m::Model=get_model(); year::Union{Int, Nothing} = nothing, last_year::Int = model_years[end], prtp::Float64 = 0.015, eta::Float64=1.5, pulse_size=1e10)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiDICE2010 model.
If no model is provided, the default model from MimiDICE2010.get_model() is used.
Ramsey discounting is used with a pure rate of time preference of `prtp` and inequality aversion of `eta`.
`pulse_size` controls the size of the marginal emission pulse.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta::Float64=1.5, pulse_size=1e10)
    year === nothing ? error("Must specify an emission year. Try `compute_scc_mm(m, year=2015)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = get_marginal_model(m; year=year, pulse_size=pulse_size)
    scc = _compute_scc(mm; year=year, last_year=last_year, prtp=prtp, eta=eta)

    return (scc = scc, mm = mm)
end

# helper function for computing SCC from a MarginalModel, not to be exported or advertised to users
function _compute_scc(mm::MarginalModel; year::Int, last_year::Int, prtp::Float64, eta::Float64)
    ntimesteps = findfirst(isequal(last_year), model_years)     # Will run through the timestep of the specified last_year
    run(mm, ntimesteps=ntimesteps)

    marginal_damages = -1 * mm[:neteconomy, :C][1:ntimesteps] * 1e12 * 12 / 44 # convert from trillion $/ton C to $/ton CO2; multiply by -1 to get positive value for damages

    cpc = mm.base[:neteconomy, :CPC]

    year_index = findfirst(isequal(year), model_years)

    df = [zeros(year_index - 1)..., ((cpc[year_index] / cpc[i])^eta * 1 / (1 + prtp)^(t - year) for (i, t) in enumerate(model_years) if year <= t <= last_year)...]
    scc = sum(df .* marginal_damages * 10)  # currently implemented as a 10year step function; so each timestep of discounted marginal damages is multiplied by 10
    return scc
end

"""
get_marginal_model(m::Model=get_model(); year::Union{Int, Nothing} = nothing, pulse_size::Float64=1e10)

Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from MimiDICE2010.get_model() is used as the base model.
`pulse_size` controls the size of the marginal emission pulse.
"""
function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing, pulse_size::Float64=1e10)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2015)`.") : nothing
    !(year in model_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = create_marginal_model(m, pulse_size) # Pulse has a value of 1GtC per year for ten years
    add_marginal_emissions!(mm.modified, year, pulse_size)

    return mm
end

"""
Adds a marginal emission component to `m` which adds `pulse_size` of additional C emissions over ten years starting in the specified `year`.
"""
function add_marginal_emissions!(m::Model, year::Int, pulse_size::Float64)
    add_comp!(m, Mimi.adder, :marginalemission, before=:co2cycle)

    time = Mimi.dimension(m, :time)
    addem = zeros(length(time))
    addem[time[year]] = pulse_size / 1e10     # Unit of pulse_size is tons, but units of emissions in DICE are GtC, so we convert to GtC by dividing by 1e9, and then divide by 10 again because that pulse is emitted for ten years.

    set_param!(m, :marginalemission, :add, addem)
    connect_param!(m, :marginalemission, :input, :emissions, :E)
    connect_param!(m, :co2cycle, :E, :marginalemission, :output)
end

# Old available function
function getmarginal_dice_models(;emissionyear=2015)

    DICE = get_model()
    run(DICE)

    mm = MarginalModel(DICE)
    m1 = mm.base
    m2 = mm.modified

    add_comp!(m2, Mimi.adder, :marginalemission, before=:co2cycle)

    time = Mimi.dimension(m1, :time)
    addem = zeros(length(time))
    addem[time[emissionyear]] = 1.0

    set_param!(m2, :marginalemission, :add, addem)
    connect_param!(m2, :marginalemission, :input, :emissions, :E)
    connect_param!(m2, :co2cycle, :E, :marginalemission, :output)

    run(m1)
    run(m2)

    return m1, m2
end


text:
# Pointer/type aliases
const QueueInterceptPacketWriter = Ptr{Cvoid}
const QueueInterceptHandler = Ptr{Cvoid}
const RuntimeQueueNotifier = Ptr{Cvoid}
const SignalHandler = Ptr{Cvoid}
const DeallocationCallback = Ptr{Cvoid}
const AQLProfileDataCallback = Ptr{Cvoid}
const SystemEventCallback = Ptr{Cvoid}
const ToolsEvent = UInt64
const SignalValue = Int64
#const SignalValue = Union{Int64,UInt64} #TODO: check me on use of signal_store_relaxed
const MachineVersion = UInt16
const File = Cint
const ExtModule = Ptr{Cvoid}

# Composite types

struct AMDGPUSamplerDescriptor
    size::UInt16
    kind::MetadataKind
    coord::SamplerCoord
    filter::SamplerFilter
    addressing::SamplerAddressing
    reserved1::UInt8
end

struct AMDGPUHSAImageDescriptor
    size::UInt16
    kind::MetadataKind
    geometry::ImageGeometry
    channel_order::ImageChannelOrder
    channel_type::ImageChannelType
    reserved1::UInt8
    width::UInt64
    height::UInt64
    depth::UInt64
    array::UInt64
end

struct NoteCodeObjectVersion
    major_version::UInt32
    minor_version::UInt32
end

struct NoteHSAIL
    hsail_major_version::UInt32
    hsail_minor_version::UInt32
    profile::UInt8
    machine_model::UInt8
    default_float_round::UInt8
end

struct NoteISA
    vendor_name_size::UInt16
    architecture_name_size::UInt16
    major::UInt32
    minor::UInt32
    stepping::UInt32
    vendor_and_architecture_name::NTuple{1, UInt8}
end

struct NoteProducer
    producer_name_size::UInt16
    reserved::UInt16
    producer_major_version::UInt32
    producer_minor_version::UInt32
    producer_name::NTuple{1, UInt8}
end

struct NoteProducerOptions
    producer_options_size::UInt16
    producer_options::NTuple{1, UInt8}
end

struct Dim3
    x::UInt32
    y::UInt32
    z::UInt32
end

struct Agent
    handle::UInt64
end

struct Cache
    handle::UInt64
end

struct Signal
    handle::UInt64
end

struct SignalGroup
    handle::UInt64
end

struct Region
    handle::UInt64
end

struct Queue
    type::QueueType
    features::UInt32
    base_address::Ptr{Cvoid}
    doorbell_signal::Signal
    size::UInt32
    reserved1::UInt32
    id::UInt64
end

struct KernelDispatchPacket
    header::UInt16
    setup::UInt16
    workgroup_size_x::UInt16
    workgroup_size_y::UInt16
    workgroup_size_z::UInt16
    reserved0::UInt16
    grid_size_x::UInt32
    grid_size_y::UInt32
    grid_size_z::UInt32
    private_segment_size::UInt32
    group_segment_size::UInt32
    kernel_object::UInt64
    kernarg_address::Ptr{Cvoid}
    reserved2::UInt64
    completion_signal::Signal
end

struct AgentDispatchPacket
    header::UInt16
    type::UInt16
    reserved0::UInt32
    return_address::Ptr{Cvoid}
    arg::NTuple{4, UInt64}
    reserved2::UInt64
    completion_signal::Signal
end

struct BarrierAndPacket
    header::UInt16
    reserved0::UInt16
    reserved1::UInt32
    dep_signal::NTuple{5, Signal}
    reserved2::UInt64
    completion_signal::Signal
end

struct BarrierOrPacket
    header::UInt16
    reserved0::UInt16
    reserved1::UInt32
    dep_signal::NTuple{5, Signal}
    reserved2::UInt64
    completion_signal::Signal
end

struct ISA
    handle::UInt64
end

struct Wavefront
    handle::UInt64
end

struct CodeObjectReader
    handle::UInt64
end

struct Executable
    handle::UInt64
end

struct LoadedCodeObject
    handle::UInt64
end

struct ExecutableSymbol
    handle::UInt64
end

struct CodeObject
    handle::UInt64
end

struct CallbackData
    handle::UInt64
end

struct CodeSymbol
    handle::UInt64
end

struct AMDControlDirectives
    enabled_control_directives::EnabledControlDirective
    enable_break_exceptions::UInt16
    enable_detect_exceptions::UInt16
    max_dynamic_group_size::UInt32
    max_flat_grid_size::UInt64
    max_flat_workgroup_size::UInt32
    required_dim::UInt8
    reserved1::NTuple{3, UInt8}
    required_grid_size::NTuple{3, UInt64}
    required_workgroup_size::NTuple{3, UInt32}
    reserved2::NTuple{60, UInt8}
end

struct AMDKernelCode
    AMD_kernel_code_version_major::KernelCodeVersion
    AMD_kernel_code_version_minor::KernelCodeVersion
    AMD_machine_kind::MachineKind
    AMD_machine_version_major::MachineVersion
    AMD_machine_version_minor::MachineVersion
    AMD_machine_version_stepping::MachineVersion
    kernel_code_entry_byte_offset::Int64
    kernel_code_prefetch_byte_offset::Int64
    kernel_code_prefetch_byte_size::UInt64
    max_scratch_backing_memory_byte_size::UInt64
    compute_pgm_rsrc1::ComputePgmRsrcOne
    compute_pgm_rsrc2::ComputePgmRsrcTwo
    kernel_code_properties::KernelCodeProperties
    workitem_private_segment_byte_size::UInt32
    workgroup_group_segment_byte_size::UInt32
    gds_segment_byte_size::UInt32
    kernarg_segment_byte_size::UInt64
    workgroup_fbarrier_count::UInt32
    wavefront_sgpr_count::UInt16
    workitem_vgpr_count::UInt16
    reserved_vgpr_first::UInt16
    reserved_vgpr_count::UInt16
    reserved_sgpr_first::UInt16
    reserved_sgpr_count::UInt16
    debug_wavefront_private_segment_offset_sgpr::UInt16
    debug_private_segment_buffer_sgpr::UInt16
    kernarg_segment_alignment::PowerTwo
    group_segment_alignment::PowerTwo
    private_segment_alignment::PowerTwo
    wavefront_size::PowerTwo
    call_convention::Int32
    reserved1::NTuple{12, UInt8}
    runtime_loader_kernel_symbol::UInt64
    control_directives::AMDControlDirectives
end

struct RuntimeLoaderDebugInfo
    elf_raw::Ptr{Cvoid}
    elf_size::UInt64
    kernel_name::Cstring
    owning_segment::Ptr{Cvoid}
end

struct AMDQueue
    queue::Queue
    reserved1::NTuple{4, UInt32}
    write_dispatch_id::UInt64
    group_segment_aperture_base_hi::UInt32
    private_segment_aperture_base_hi::UInt32
    max_cu_id::UInt32
    max_wave_id::UInt32
    max_legacy_doorbell_dispatch_id_plus_1::UInt64
    legacy_doorbell_lock::UInt32
    reserved2::NTuple{9, UInt32}
    read_dispatch_id::UInt64
    read_dispatch_id_field_base_byte_offset::UInt32
    compute_tmpring_size::UInt32
    scratch_resource_descriptor::NTuple{4, UInt32}
    scratch_backing_memory_location::UInt64
    scratch_backing_memory_byte_size::UInt64
    scratch_wave64_lane_byte_size::UInt32
    queue_properties::QueueProperties
    reserved3::NTuple{2, UInt32}
    queue_inactive_signal::Signal
    reserved4::NTuple{14, UInt32}
end

struct AMDSignal
    kind::SignalKind
    event_mailbox_ptr::UInt64
    event_id::UInt32
    reserved1::UInt32
    start_ts::UInt64
    end_ts::UInt64
    reserved3::NTuple{2, UInt32}
end

struct APITableVersion
    major_id::UInt32
    minor_id::UInt32
    step_id::UInt32
    reserved::UInt32
end

#=
struct HSAAPITable
    version::APITableVersion
    core::Ptr{CoreAPITable}
    AMD_ext::Ptr{AMDExtTable}
    finalizer_ext::Ptr{FinalizerExtTable}
    image_ext::Ptr{ImageExtTable}
end
=#
struct Image
    handle::UInt64
end

struct ImageFormat
    channel_type::ImageChannelType
    channel_order::ImageChannelOrder
end

struct ExtImageDescriptor
    geometry::ImageGeometry
    width::UInt64
    height::UInt64
    depth::UInt64
    array_size::UInt64
    format::ImageFormat
end

struct ImageDataInfo
    size::UInt64
    alignment::UInt64
end

struct ImageRegion
    offset::Dim3
    range::Dim3
end

struct Sampler
    handle::UInt64
end

struct ExtSamplerDescriptor
    coordinate_mode::SamplerCoordinateMode
    filter_mode::SamplerFilterMode
    address_mode::SamplerAddressingMode
end

struct Images_1_00Pfn
    image_get_capability::Ptr{Cvoid}
    image_data_get_info::Ptr{Cvoid}
    image_create::Ptr{Cvoid}
    image_destroy::Ptr{Cvoid}
    image_copy::Ptr{Cvoid}
    image_import::Ptr{Cvoid}
    image_export::Ptr{Cvoid}
    image_clear::Ptr{Cvoid}
    sampler_create::Ptr{Cvoid}
    sampler_destroy::Ptr{Cvoid}
end

struct Images_1Pfn
    image_get_capability::Ptr{Cvoid}
    image_data_get_info::Ptr{Cvoid}
    image_create::Ptr{Cvoid}
    image_destroy::Ptr{Cvoid}
    image_copy::Ptr{Cvoid}
    image_import::Ptr{Cvoid}
    image_export::Ptr{Cvoid}
    image_clear::Ptr{Cvoid}
    sampler_create::Ptr{Cvoid}
    sampler_destroy::Ptr{Cvoid}
    image_get_capability_with_layout::Ptr{Cvoid}
    image_data_get_info_with_layout::Ptr{Cvoid}
    image_create_with_layout::Ptr{Cvoid}
end

struct AMDHDPFlush
    HDP_MEM_FLUSH_CNTL::Ptr{UInt32}
    HDP_REG_FLUSH_CNTL::Ptr{UInt32}
end

struct ProfilingDispatchTime
    start::UInt64
    _end::UInt64
end

struct ProfilingAsyncCopyTime
    start::UInt64
    _end::UInt64
end

struct MemoryPool
    handle::UInt64
end

struct PitchedPtr
    base::Ptr{Cvoid}
    pitch::UInt64
    slice::UInt64
end

struct MemoryPoolLinkInfo
    min_latency::UInt32
    max_latency::UInt32
    min_bandwidth::UInt32
    max_bandwidth::UInt32
    atomic_support_32bit::Bool
    atomic_support_64bit::Bool
    coherent_support::Bool
    link_type::LinkInfoType
    numa_distance::UInt32
end

struct AMDImageDescriptor
    version::UInt32
    deviceID::UInt32
    data::NTuple{1, UInt32}
end

struct PointerInfo
    size::UInt32
    type::PointerType
    agentBaseAddress::Ptr{Cvoid}
    hostBaseAddress::Ptr{Cvoid}
    sizeInBytes::UInt64
    userData::Ptr{Cvoid}
    agentOwner::Agent
end

struct IPCMemory
    handle::NTuple{8, UInt32}
end

const IPCSignal = IPCMemory

struct GpuMemoryFaultInfo
    agent::Agent
    virtual_address::UInt64
    fault_reason_mask::UInt32
end

struct Event
    event_type::AMDEventType
end

struct ExtToolsExceptionPolicy
    exception_mask::UInt32
    wave_action::WaveAction
    host_action::HostAction
    wave_mode::WaveMode
end

struct DispatchDebugInfo
    scratch_memory_descriptor::NTuple{4, UInt32}
    global_memory_descriptor::NTuple{4, UInt32}
end

struct WaveMsgAmdGen2
    ui32::UInt32
    reserved0::UInt32
end

struct WaveControlMessage
    memory_va::Ptr{Cvoid}
end

struct KernelExecutionMode
    ui32::UInt32
end

struct ExtProgram
    handle::UInt64
end

struct ExtControlDirectives
    control_directives_mask::UInt64
    break_exceptions_mask::UInt16
    detect_exceptions_mask::UInt16
    max_dynamic_group_size::UInt32
    max_flat_grid_size::UInt64
    max_flat_workgroup_size::UInt32
    reserved1::UInt32
    required_grid_size::NTuple{3, UInt64}
    required_workgroup_size::Dim3
    required_dim::UInt8
    reserved2::NTuple{75, UInt8}
end

struct ExtFinalizer_1_00Pfn
    EXT_program_create::Ptr{Cvoid}
    EXT_program_destroy::Ptr{Cvoid}
    EXT_program_add_module::Ptr{Cvoid}
    EXT_program_iterate_modules::Ptr{Cvoid}
    EXT_program_get_info::Ptr{Cvoid}
    EXT_program_finalize::Ptr{Cvoid}
end

struct AQLProfileEvent
    block_name::AQLProfileBlockName
    block_index::UInt32
    counter_id::UInt32
end

struct AQLProfileParameter
    parameter_name::AQLProfileParameterName
    value::UInt32
end

struct AQLProfileDescriptor
    ptr::Ptr{Cvoid}
    size::UInt32
end

struct AQLProfile
    agent::Agent
    type::AQLProfileEventType
    events::Ptr{AQLProfileEvent}
    event_count::UInt32
    parameters::Ptr{AQLProfileParameter}
    parameter_count::UInt32
    output_buffer::AQLProfileDescriptor
    command_buffer::AQLProfileDescriptor
end

struct AQLPm4Packet
    header::UInt16
    pm4_command::NTuple{27, UInt16}
    completion_signal::Signal
end

struct AQLProfileInfoData
    sample_id::UInt32
end

struct AQLProfileIDQuery
    name::Cstring
    id::UInt32
    instance_count::UInt32
end

struct AQLProfile_1_00Pfn
    version_major::Ptr{Cvoid}
    version_minor::Ptr{Cvoid}
    error_string::Ptr{Cvoid}
    validate_event::Ptr{Cvoid}
    start::Ptr{Cvoid}
    stop::Ptr{Cvoid}
    read::Ptr{Cvoid}
    legacy_get_pm4::Ptr{Cvoid}
    get_info::Ptr{Cvoid}
    iterate_data::Ptr{Cvoid}
end

struct SegmentDescriptor
    agent::Agent
    executable::Executable
    code_object_storage_type::CodeObjectStorageType
    code_object_storage_base::Ptr{Cvoid}
    code_object_storage_size::UInt64
    code_object_storage_offset::UInt64
    segment_base::Ptr{Cvoid}
    segment_size::UInt64
end

struct AMDLoader_1_00Pfn
    host_address::Ptr{Cvoid}
    segment_descriptors::Ptr{Cvoid}
    executable::Ptr{Cvoid}
end

struct AMDLoader_1_01Pfn
    query_host_address::Ptr{Cvoid}
    query_segment_descriptors::Ptr{Cvoid}
    query_executable::Ptr{Cvoid}
    executable_iterate_loaded_code_objects::Ptr{Cvoid}
    loaded_code_object_get_info::Ptr{Cvoid}
end


text:
using MeshCat
using GeometryBasics
using CoordinateTransformations
using Rotations
import Plots as pl

function geom_init()
    vis = Visualizer()
    render(vis)
    delete!(vis)
    return vis
end

function geom_vis(vis, qhist, dt, speed=1)

    green_material = MeshPhongMaterial(color=pl.RGBA(0, 1, 0, 0.8))
    red_material = MeshPhongMaterial(color=pl.RGBA(1, 0, 0, 0.8))

    cylinderb = Cylinder(Point(0, -0.02, l_cb[2]), Point(0, 0.02, l_cb[2]), 0.04)  #float(l_cb)
    cylinder0 = Cylinder(Point(-l_c0[1], -l_c0[2], -l_c0[3]), Point((l0-l_c0)[1], (l0-l_c0)[2], (l0-l_c0)[3]), 0.008)
    cylinder1 = Cylinder(Point(-l_c1[1], -l_c1[2], -l_c1[3]), Point((l1-l_c1)[1],(l1-l_c1)[2],(l1-l_c1)[2]), 0.008)
    cylinder2 = Cylinder(Point(-l_c2[1], -l_c2[2], -l_c2[3]), Point((l2-l_c2)[1], (l2-l_c2)[2], (l2-l_c2)[3]), 0.008)
    cylinder3 = Cylinder(Point(-l_c3[1], -l_c3[2], -l_c3[3]), Point((lee-l_c3)[1], (lee-l_c3)[2],  (lee-l_c3)[3]), 0.008)

    setobject!(vis["cylinderb"],cylinderb,red_material)
    setobject!(vis["cylinder0"],cylinder0,green_material)
    setobject!(vis["cylinder1"],cylinder1,green_material)
    setobject!(vis["cylinder2"],cylinder2,green_material)
    setobject!(vis["cylinder3"],cylinder3,green_material)

    for k = 1:N

        q = qhist[:, k]

        # set position and attitude
        positionb = Translation(q[1:3]...)
        attitudeb = LinearMap(UnitQuaternion(q[4:7]))
        position0 = Translation(q[8:10]...)
        attitude0 = LinearMap(UnitQuaternion(q[11:14]))
        position1 = Translation(q[15:17]...)
        attitude1 = LinearMap(UnitQuaternion(q[18:21]))
        position2 = Translation(q[22:24]...)
        attitude2 = LinearMap(UnitQuaternion(q[25:28]))
        position3 = Translation(q[29:31]...)
        attitude3 = LinearMap(UnitQuaternion(q[32:35]))

        settransform!(vis["cylinderb"], compose(positionb,attitudeb))
        settransform!(vis["cylinder0"], compose(position0,attitude0))
        settransform!(vis["cylinder1"], compose(position1,attitude1))
        settransform!(vis["cylinder2"], compose(position2,attitude2))
        settransform!(vis["cylinder3"], compose(position3,attitude3))
        sleep(dt/speed)
    end

    return nothing
end

function urdf_init()
    curdir = pwd()
    urdfpath = joinpath(curdir, "res/flyhopper_robot/urdf/flyhopper_robot_base_call.urdf")
    leg = parse_urdf(urdfpath, floating=true)
    # state = MechanismState(leg)
    # mvis = MechanismVisualizer(doublependulum, Skeleton(randomize_colors=true, inertias=false));
    mvis = MechanismVisualizer(leg, URDFVisuals(urdfpath));

    render(mvis)
    return mvis
end

function urdf_vis(mvis, qhist, dt, speed=1)

    for k = 1:N
        q = copy(qhist[:, k])
        pb = q[1:3]
        Qb = q[4:7]
        Q0 = q[11:14]
        Q1 = q[18:21]
        Q2 = q[25:28]
        Q3 = q[32:35]

        # convert quaternions to joint angles
        a0 = -angle_y(Qb, Q0) +30*(pi/180)
        a1 = -angle_y(Q0, Q1) +120*(pi/180)
        a2 = -angle_y(Qb, Q2) +150*(pi/180)
        a3 = -angle_y(Q2, Q3) -120*(pi/180)

        #-0.5235987755982988, -2.6179938779914944, -2.0943951023931953, 2.0943951023931953
        q_array = vcat(Qb, pb, [a0, a2, a1, a3])
        # @show [a0, a2, a1, a3]
        set_configuration!(mvis, q_array)
        sleep(dt/speed)
    end
end

#=
function a_urdf(a)
    # adjust angles for urdf animation
    a0 = a[1] +30*(pi/180)
    a1 = a[2] +120*(pi/180)
    a2 = a[3] +150*(pi/180)
    a3 = a[4] -120*(pi/180)

    return [a0, a2, a1, a3]
end
=#

max_stars_repo_name: bbokser/hopper-sim

max_issues_repo_name: bbokser/hopper-sim

max_forks_repo_name: bbokser/hopper-sim


text:
# # Example 5.4 of [AJPR14]

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/AJPR14e54.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/AJPR14e54.ipynb)
# **Contributed by**: Benoît Legat

# In the Example 5.4 of [AJPR14], the authors study the following system whose JSR is 3.917384715148.
# In this example, we reproduce their numerical results.
# The system is characterized by the matrices $\mathcal{A} = \{A_1, A_2\}$ where the matrices are given below:
#
# [AJPR14] A. Ahmadi, R. Jungers, P. Parrilo and M. Roozbehani,
# *Joint spectral radius and path-complete graph Lyapunov functions.*
# SIAM J. CONTROL OPTIM 52(1), 687-717, **2014**.

using Test #src
using SwitchOnSafety
A1 = [-1 -1
      -4  0]
A2 = [ 3  3
      -2  1]
s = discreteswitchedsystem([A1, A2])

# Pick an SDP solver from [this list](http://jump.dev/JuMP.jl/v0.21.3/installation/#Getting-Solvers-1).

import CSDP
optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true);

# We first try with a Common Quadratic Lyapunov Function (CQLF) and obtain $\rho_{\mathcal{V}^2}(\mathcal{A}) \approx 3.980613$ and the corresponding lower bound $\rho_{\mathcal{V}^2}(\mathcal{A}) / \sqrt{2} ≈ 2.814547$.

lb, ub = soslyapb(s, 1, optimizer_constructor=optimizer_constructor, tol=1e-6, verbose=1)
@test ub ≈ 3.9805 rtol=1e-4 #src

# From the infeasibility certificate of the last infeasible SemiDefinite Program (SDP) solved by in order to obtain the upper bound, we find the lower bound $\rho(A_2A_1)^{1/2} \approx 3.917384715148$.

seq = sosbuildsequence(s, 1, p_0=:Primal)
psw = findsmp(seq)
@test psw.growthrate ≈ 3.917384715148 #src

# Now with a quartic lyapunov function, we get the tighter upper bound
# $\rho_{\mathcal{V}^{\text{SOS},4}}(\mathcal{A}) \approx 3.9241$.

lb, ub = soslyapb(s, 2, optimizer_constructor=optimizer_constructor, tol=1e-6, verbose=1)
@test ub ≈ 3.9241 rtol=1e-4 #src

# We find the same lower bound from the infeasibility certificate.

seq = sosbuildsequence(s, 2, p_0=:Primal)
psw = findsmp(seq)
@test psw.growthrate ≈ 3.917384715148 #src

# This upper bound is further improved by computing quadratic lyapunov functions
# for the the system obtained with the graph $G_1$ of [Fig. 3.1, AJPR14].

G1 = mdependentlift(s, 1, false)
lb, ub = soslyapb(G1, 1, optimizer_constructor=optimizer_constructor, tol=1e-6, verbose=1)
@test ub ≈ 3.9224 rtol=1e-4 #src

# We again find the same lower bound from the infeasibility certificate.

seq = sosbuildsequence(s, 2, p_0=:Primal)
psw = findsmp(seq)
@test psw.growthrate ≈ 3.917384715148 #src


text:
# mriview.jl
# Simple MRI viewer

# This script displays a NIfTI volume using ImageView.

using NIfTI, Tk, Images, ImageView

function mriview(ni::NIVolume, xy=["x", "z"], zoom::Int=5)
    img = Image(ni.raw, ["spatialorder" => ["x", "y", "z"], "timedim" => 4, "colorspace" => "Gray",
                         "pixelspacing"=>voxel_size(ni.header)])
    img = scale(scaleminmax(img), img)
    imgc, imgslice = display(img, xy=xy)
    tksize = get_size(toplevel(imgc))
    controlheight = tksize[2] - size(img, xy[2])
    set_size(toplevel(imgc), max(tksize[1], size(img, xy[1])*zoom),
                             size(img, xy[2])*zoom+controlheight)
    (imgc, imgslice)
end
mriview(path::String) = mriview(niread(path, mmap=true))


text:
using MLDatasets
using ImageFiltering
using Images
using DelimitedFiles

cifar_train,_  = CIFAR10.traindata()

function load_96_unlabeled_stl10()
    train_images = Array{UInt8}(undef, 96, 96, 3, 100000)
    read!(joinpath("data","STL10", "unlabeled_X.bin"), train_images)
    return (train_images)
end

function load_stl10()
    train_images = Array{UInt8}(undef, 96, 96, 3, 8000)
    read!(joinpath("data","STL10", "test_X.bin"), train_images)
    return (train_images)
end

stl10_train = zeros(10,10,10,10) #load_stl10()

function generate_kappa!(n, m; type=0, smooth=false, threshold=50, kernel=3)
    if type == 0 # Uniform Model
        return ones(r_type, n-1, m-1)
    end
    if type > 4 # Geographical Model
        if type == 5
            full_kappa = readdlm("SEGmodel2Dsalt.dat")
        elseif type == 6
            full_kappa = readdlm("overthrust_XZ_slice.dat")
        elseif type == 7
            full_kappa = readdlm("MarmousiModelWithoutPad.dat")
        end
        sample = full_kappa[n+ceil(Int64,(size(full_kappa,1)-(n-1))/2-1):-1:ceil(Int64,(size(full_kappa,1)-(n-1))/2)+1,
                         ceil(Int64,(size(full_kappa,2)-(m-1))/2)+1:m+ceil(Int64,(size(full_kappa,2)-(m-1))/2)-1]
        sample = 1.0 ./ sample
    else
        if type == 2 # STL10 Model
            index = rand(1:size(stl10_train,4))
            sample = Gray.(colorview(RGB, stl10_train[:,:,1,index] / 255,stl10_train[:,:,2,index] / 255,stl10_train[:,:,3,index] / 255))
        else # CIFAR10 Model
            index = 7706 # Const Model - 4
            if type == 1
                index = rand(1:size(cifar_train,4))
            end
            sample = Gray.(colorview(RGB, cifar_train[:,:,1,index],cifar_train[:,:,2,index],cifar_train[:,:,3,index]))
        end

        # Resize
        sample = imresize(sample, (n-1,m-1))

        # Smooth
        if smooth == true
            sample = imfilter(sample, Kernel.gaussian(kernel))
        end
    end

    # Sample ∈ [random threshold, 1]
    threshold = 0.01 * threshold #  rand(threshold:90)
    sample_normal = threshold .+ (((1.0 - threshold) .* (sample .- minimum(sample))) ./ (maximum(sample) - minimum(sample)))
    sample_normal1 = r_type.(sample_normal)

    return sample_normal1
end


text:
using Distributions, LinearAlgebra
using Random: seed!
using AugmentedGaussianProcesses
using CSV, DataFrames
seed!(42)
N = 3000
D = 20
K = 4
k = RBFKernel(1.0)
X = rand(N, D)
y = rand(MvNormal(kernelmatrix(X, k) + 1e-3I))
df = DataFrame(X)
df.y_reg = y
df.y_class = sign.(y)
width = maximum(y) - minimum(y);
normy = (y .- minimum(y)) / width * K;
df.y_multi = floor.(Int64, normy)
CSV.write(joinpath(@__DIR__, "benchmarkdata.csv"), df)


text:
"""

In the 20×20 grid below, four numbers along a diagonal line have been marked with >

08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 >26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 >63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 >78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 >14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48

The product of these numbers is 26 × 63 × 78 × 14 = 1788696.

What is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally) in the 20×20 grid?
"""

grid =[08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48]

diag(A) = [A[i,i] for i in 1:size(A)[1]]

function maxgridf(grid::AbstractArray{Int},ksz::Int,f=prod)
    maxval = 0
    # horizontally right
    for i in 1:size(grid)[1]
        for j in 1:size(grid)[2]-ksz
            t = f(grid[i,j:j+ksz-1])
            if t > maxval
                maxval=t
            end
        end
    end

    # vertical down
    for i in 1:size(grid)[1]-ksz
        for j in 1:size(grid)[2]
            t = f(grid[i:i+ksz-1,j])
            if t > maxval
                maxval=t
            end
        end
    end

    # # right diagonal

    for i in 1:size(grid)[1]-ksz
        for j in 1:size(grid)[2]-ksz
            t = f(diag(grid[i:i+ksz-1,j:j+ksz-1]))
            if t > maxval
                maxval=t
            end
        end
    end

    # # left diagonal

    for i in 1:size(grid)[1]-ksz
        for j in 1:size(grid)[2]-ksz
            t = f(diag(grid[i:i+ksz-1,j:j+ksz-1][end:-1:1,:]))
            if t > maxval
                maxval=t
            end
        end
    end

    return maxval

end

maxgridf(grid,4)


text:
using StanModels

df = filter(row -> row[:age] >= 18,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "Howell1.csv"), delim=';'))

max_height = maximum(df[!, :height])
min_height = minimum(df[!, :height])

m4_2s = "
data {
  int N;
  real h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=$(min_height),upper=$(max_height)> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 0.1);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

sm = SampleModel("m4.2s", m4_2s);

m4_2_data = Dict("N" => size(df, 1), "h" => df[!, :height]);

(sample_file, log_file) = stan_sample(sm, data=m4_2_data);

if !(sample_file == nothing)
  chn = read_samples(sm)
  chn = set_names(chn, Dict("mu" => "μ", "sigma" => "σ"))
  describe(chn)
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Objective #
# ========================== #
struct FunctionizeDecisionObjectiveBridge{T} <: MOIB.Objective.AbstractBridge end

function MOIB.Objective.bridge_objective(::Type{FunctionizeDecisionObjectiveBridge{T}}, model::MOI.ModelLike,
                                         f::SingleDecision) where T
    F = AffineDecisionFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), convert(F, f))
    return FunctionizeDecisionObjectiveBridge{T}()
end

function MOIB.Objective.supports_objective_function(
    ::Type{<:FunctionizeDecisionObjectiveBridge}, ::Type{SingleDecision})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:FunctionizeDecisionObjectiveBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:FunctionizeDecisionObjectiveBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.set_objective_function_type(::Type{FunctionizeDecisionObjectiveBridge{T}}) where T
    return AffineDecisionFunction{T}
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FunctionizeDecisionObjectiveBridge, ::MOI.NumberOfVariables)
    return 0
end
function MOI.get(bridge::FunctionizeDecisionObjectiveBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# No variables or constraints are created in this bridge so there is nothing to
# delete.
function MOI.delete(model::MOI.ModelLike, bridge::FunctionizeDecisionObjectiveBridge) end

function MOI.set(::MOI.ModelLike, ::MOI.ObjectiveSense,
                 ::FunctionizeDecisionObjectiveBridge, ::MOI.OptimizationSense)
    # `FunctionizeDecisionObjectiveBridge` is sense agnostic, therefore, we don't need to change
    # anything.
    return nothing
end
function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{SingleDecision},
                 bridge::FunctionizeDecisionObjectiveBridge{T}) where T
    F = AffineDecisionFunction{T}
    return MOI.get(model, MOIB.ObjectiveFunctionValue{F}(attr.result_index))
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ObjectiveFunction{SingleDecision},
                 bridge::FunctionizeDecisionObjectiveBridge{T}) where T
    F = AffineDecisionFunction{T}
    func = MOI.get(model, MOI.ObjectiveFunction{F}())
    return convert(SingleDecision, func)
end


text:
export MtlDevice, devices

## construction

const MTLDevice = Ptr{MtDevice}

struct MtlDevice
    handle::MTLDevice
end

Base.convert(::Type{MTLDevice}, dev::MtlDevice) = dev.handle
Base.unsafe_convert(::Type{MTLDevice}, d::MtlDevice) = convert(MTLDevice, d.handle)

Base.:(==)(a::MtlDevice, b::MtlDevice) = a.handle == b.handle
Base.hash(dev::MtlDevice, h::UInt) = hash(dev.handle, h)

"""
    devices()

Get an iterator for the compute devices.
"""
function devices()
    handles = mtCopyAllDevices()
    devices = Vector{MtlDevice}()
    while true
        handle = unsafe_load(handles, length(devices)+1)
        handle == C_NULL && break
        push!(devices, MtlDevice(handle))
    end
    Libc.free(handles)

    return devices
end

"""
    MtlDevice(i::Integer)

Get a handle to a compute device.
"""
MtlDevice(i::Integer) = devices()[i]

function Base.show(io::IO, d::MtlDevice)
    print(io, "MtlDevice($(name(d)))")
end

function Base.show(io::IO, ::MIME"text/plain", d::MtlDevice)
    println(io, "MtlDevice:")
    println(io, " name :             ", name(d))
    println(io, " lowpower :         ", is_lowpower(d))
    println(io, " headless :         ", is_headless(d))
    println(io, " removable :        ", is_removable(d))
    println(io, " unified memory :   ", has_unified_memory(d))
    println(io, " id :               ", registry_id(d))
    print(io,   " transfer rate :    ", max_transfer_rate(d))
end

## properties

name(d::MtlDevice) = unsafe_string(mtDeviceName(d))

is_lowpower(d::MtlDevice) = mtDeviceLowPower(d)
is_headless(d::MtlDevice) = mtDeviceHeadless(d)
is_removable(d::MtlDevice) = mtDeviceRemovable(d)
has_unified_memory(d::MtlDevice) = mtDeviceHasUnifiedMemory(d)
registry_id(d::MtlDevice) = mtDeviceRegistryID(d)
max_transfer_rate(d::MtlDevice) = mtDeviceMaxTransferRate(d)

max_workingsetsize(d::MtlDevice) = mtDeviceRecommendedMaxWorkingSetSize(d)
max_threadgroupmemorylength(d::MtlDevice) = mtDeviceMaxThreadgroupMemoryLength(d)
max_threadspergroup(d::MtlDevice) = mtMaxThreadsPerThreadgroup(d)
max_bufferlength(d::MtlDevice) = mtDeviceMaxBufferLength(d)

allocatedsize(d::MtlDevice) = mtDeviceCurrentAllocatedSize(d)


text:
module ConjugateComputationVI

using AbstractGPs
using ChainRulesCore
using FastGaussQuadrature
using LinearAlgebra
using Optim
using Zygote

using AbstractGPs: AbstractGP, LatentFiniteGP, Normal

include("quadrature.jl")
include("parametrisation_bijections.jl")
include("approximate_inference.jl")

end


text:
using Mux, AssetRegistry

orchid_img = AssetRegistry.register("./assets/orchid.jpg")

function index()
    res = "

        <title>Welcome to Angkaew Data Lake service</title>

        Angkaew Data Lake

"

	# now looping thru items in our database
    for i=1:5
        res *= "

Description: asdfasfasdfadfdsfdfdf
 Contributor:
 Date:

"
    end

    res *=  "
Data science research center, Faculty of Science, Chiang Mai University"
    return res
end

@app test = (
  Mux.defaults,
  page(respond(index())),
  page("/about",
       probabilty(0.1, respond("
Boo!
")),
       respond("
About Me
")),
  page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
  Mux.notfound())

#fetch(serve(test,parse(Int,ARGS[1])))
serve(test)

max_stars_repo_name: jakramate/jlab

max_issues_repo_name: jakramate/jlab

max_forks_repo_name: jakramate/jlab


text:
"""
indegree(tree::Tree,
         node::Int64)

Determine the in degree of a `Node`
"""
function indegree(tree::Tree,
                  node::Int64)
  if !haskey(tree.nodes, node)
    error("Node does not exist")
  end
  return length(tree.nodes[node].in)
end

"""
outdegree(tree::Tree,
          node::Int64)

Determine the out degree of a `Node`
"""
function outdegree(tree::Tree,
                   node::Int64)
  if !haskey(tree.nodes, node)
    error("Node does not exist")
  end
  return length(tree.nodes[node].out)
end

"""
isroot(tree::Tree,
       node::Int64)

Determine if a `Node` is a root `Node`
"""
function isroot(tree::Tree,
                node::Int64)
  if outdegree(tree, node) > 0 && indegree(tree, node) == 0
    return true
  else
    return false
  end
end

"""
isleaf(tree::Tree,
       node::Int64)

Determine if a `Node` is a leaf `Node`
"""
function isleaf(tree::Tree,
                node::Int64)
  if outdegree(tree, node) == 0 && indegree(tree, node) == 1
    return true
  else
    return false
  end
end

"""
isinternal(tree::Tree,
           node::Int64)

Determine if a `Node` is an internal `Node`
"""
function isinternal(tree::Tree,
                    node::Int64)
  if outdegree(tree, node) > 0 && indegree(tree, node) == 1
    return true
  else
    return false
  end
end

"""
nodetype(tree::Tree,
         node::Int64)

Determine if a `Node` is an internal, root, or leaf `Node`
"""
function nodetype(tree::Tree,
                  node::Int64)
  ins = indegree(tree, node)
  outs = outdegree(tree, node)
  if ins == 0
    if outs == 0
      return "Unattached"
    else
      return "Root"
    end
  elseif ins == 1
    if outs == 0
      return "Leaf"
    elseif outs > 0
      return "Internal"
    else
      error("Unknown node type")
    end
  else
    error("Unknown node type")
  end
end

"""
findroots(tree::Tree)

Find the root `Node`s of a `Tree`
"""
function findroots(tree::Tree)
  roots = Int64[]
  for i in keys(tree.nodes)
    if isroot(tree, i)
      push!(roots, i)
    end
  end
  if length(roots) == 0
    warn("No roots detected")
  end
  return roots
end

"""
findleaves(tree::Tree)

Find the leaf `Node`s of a `Tree`
"""
function findleaves(tree::Tree)
  leaves = Int64[]
  for i in keys(tree.nodes)
    if isleaf(tree, i)
      push!(leaves, i)
    end
  end
  if length(leaves) == 0
    warn("No leaves detected")
  end
  return leaves
end

"""
findinternals(tree::Tree)

Find the internal `Node`s of a `Tree`
"""
function findinternals(tree::Tree)
  nodes = Int64[]
  for i in keys(tree.nodes)
    if isinternal(tree, i)
      push!(nodes, i)
    end
  end
  if length(nodes) == 0
    warn("No internal nodes detected")
  end
  return nodes
end

"""
findnonroots(tree::Tree)

Find the non-root `Node`s of a `Tree`
"""
function findnonroots(tree::Tree)
  nonroots = Int64[]
  for i in keys(tree.nodes)
    if !isroot(tree, i)
      push!(nonroots, i)
    end
  end
  if length(nonroots) == 0
    warn("No non-roots detected")
  end
  return nonroots
end

"""
findnonleaves(tree::Tree)

Find the non-leaf `Node`s of a `Tree`
"""
function findnonleaves(tree::Tree)
  nonleaves = Int64[]
  for i in 1:keys(tree.nodes)
    if !isleaf(tree, i)
      push!(nonleaves, i)
    end
  end
  if length(nonleaves) == 0
    warn("No non-leaves detected")
  end
  return nonleaves
end

"""
findexternals(tree::Tree)

Find the external `Node`s of a `Tree`
"""
function findexternals(tree::Tree)
  externals = Int64[]
  for i in keys(tree.nodes)
    if !isinternal(tree, i)
      push!(externals, i)
    end
  end
  if length(externals) == 0
    warn("No external nodes detected")
  end
  return externals
end

"""
parentnode(tree::Tree,
           node::Int64)

Find parent `Node`
"""
function parentnode(tree::Tree,
                    node::Int64)
  if indegree(tree, node) == 1
    return tree.branches[tree.nodes[node].in[1]].source
  else
    error("In degree of specified node != 1")
  end
end

"""
childnodes(tree::Tree,
           node::Int64)

Find child `Node`s
"""
function childnodes(tree::Tree,
                    node::Int64)
  if !haskey(tree.nodes, node)
    error("Node does not exist")
  end
  nodes = Int64[]
  for i in tree.nodes[node].out
    push!(nodes, tree.branches[i].target)
  end
  return nodes
end

"""
descendantnodes(tree::Tree,
                node::Int64)

Find descendant `Node`s
"""
function descendantnodes(tree::Tree,
                         node::Int64)
  nodecount = [0]
  nodelist = [node]
  while nodecount[end] < length(nodelist)
    push!(nodecount, length(nodelist))
    for i in nodelist[(nodecount[end-1]+1):nodecount[end]]
      append!(nodelist, childnodes(tree, i))
    end
  end
  return nodelist[2:end]
end

"""
descendantcount(tree::Tree,
                node::Int64)

Find the number of descendant `Nodes`
"""
function descendantcount(tree::Tree,
                         node::Int64)
  return length(descendantnodes(tree, node))
end

"""
descendantcount(tree::Tree,
                nodes::Array{Int64})

Find the number of descendant `Nodes`
"""
function descendantcount(tree::Tree,
                         nodes::Array{Int64})
  count = fill(0, size(nodes))
  for i in eachindex(nodes)
    count[i] += descendantcount(tree, nodes[i])
  end
  return count
end

"""
nodepath(tree::Tree,
         node::Int64)

`Node` pathway through which a specified `Node` connects to a root
"""
function nodepath(tree::Tree,
                  node::Int64)
  path = [node]
  while isleaf(tree, path[end]) || isinternal(tree, path[end])
    push!(path, parentnode(tree, path[end]))
  end
  return path
end

"""
ancestornodes(tree::Tree,
              node::Int64)

Find ancestral `Node`s
"""
function ancestornodes(tree::Tree,
                       node::Int64)
  return nodepath(tree, node)[2:end]
end

"""
ancestorcount(tree::Tree,
              node::Int64)

Number of ancestral `Node`s
"""
function ancestorcount(tree::Tree,
                       node::Int64)
  return length(ancestornodes(tree, node))
end

"""
ancestorcount(tree::Tree,
              node::Array{Int64})

Number of ancestral `Node`s
"""
function ancestorcount(tree::Tree, nodes::Array{Int64})
  count = fill(0, size(nodes))
  for i in eachindex(nodes)
    count[i] += ancestorcount(tree, nodes[i])
  end
  return count
end

"""
noderoot(tree::Tree,
         node::Int64)

The root associated with a specified `Node`
"""
function noderoot(tree::Tree,
                  node::Int64)
  return nodepath(tree, node)[end]
end

"""
areconnected(tree::Tree,
             node1::Int64,
             node2::Int64)

Check for connectedness of two `Node`s
"""
function areconnected(tree::Tree,
                      node1::Int64,
                      node2::Int64)
  return noderoot(tree, node1) == noderoot(tree, node2)
end

"""
nodepath(tree::Tree,
         node1::Int64,
         node2::Int64)

`Node` pathway through which two specified `Node`s connect
"""
function nodepath(tree::Tree,
                  node1::Int64,
                  node2::Int64)
  if !areconnected(tree, node1, node2)
    error("Nodes are not connected")
  end
  path1 = reverse(nodepath(tree, node1))
  path2 = reverse(nodepath(tree, node2))
  minlength = minimum([length(path1), length(path2)])
  mrcnode_index = findlast(path1[1:minlength] .== path2[1:minlength])
  return [reverse(path1[(mrcnode_index+1):end]); path2[mrcnode_index:end]]
end

"""
branchpath(tree::Tree,
           node::Int64)

Branch pathway through which a specified node connects to a root
"""
function branchpath(tree::Tree,
                    node::Int64)
  path = Int64[]
  while isleaf(tree, node) || isinternal(tree, node)
    push!(path, tree.nodes[node].in[1])
    node = tree.branches[path[end]].source
  end
  return path
end

"""
branchpath(tree::Tree,
           node1::Int64,
           node2::Int64)

Branch pathway through which two specified nodes connect
"""
function branchpath(tree::Tree,
                    node1::Int64,
                    node2::Int64)
  if !areconnected(tree, node1, node2)
    error("Nodes are not connected")
  end
  path1 = reverse(branchpath(tree, node1))
  path2 = reverse(branchpath(tree, node2))
  minlength = minimum([length(path1), length(path2)])
  if minlength == 0
    mrcbranch_index = 0
  else
      branch_match = path1[1:minlength] .== path2[1:minlength]
      if any(branch_match)
          mrcbranch_index = findlast(branch_match)
      else
          mrcbranch_index = 0
      end
  end
  return [reverse(path1[(mrcbranch_index+1):end]); path2[mrcbranch_index+1:end]]
end