text:
println("Welcome to Julia, Julian!")

arithmetic_mean(x) = sum(x) / length(x)

# Arithmetic Mean
avg(x) = (ones(length(x))'x) / length(x)

# Norm of a vector
norm(x) = sqrt(x'x)

# Standard Deviation
std(x) = norm(x .- avg(x)) / sqrt(length(x))

# Quadratic Mean
rms(x) = sqrt(avg(x)^2 + std(x)^2)

# Distance function
distance(x, y) = norm(x - y)

head = first
tail = last

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "LibVPX"
version = v"1.10.0"

# Collection of sources required to build LibVPX
sources = [
    ArchiveSource("https://github.com/webmproject/libvpx/archive/v$(version).tar.gz",
                  "85803ccbdbdd7a3b03d930187cb055f1353596969c1f92ebec2db839fa4f834a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/libvpx-*/
sed -i 's/cp -p/cp/' build/make/Makefile

mkdir vpx_build && cd vpx_build
apk add diffutils yasm

if [[ "${bb_full_target}" == i686-linux-* ]]; then
    export TARGET=x86-linux-gcc
elif [[ "${bb_full_target}" == x86_64-linux-* ]]; then
    export TARGET=x86_64-linux-gcc
elif [[ "${bb_full_target}" == armv7l-linux-* ]]; then
    export TARGET=armv7-linux-gcc
elif [[ "${bb_full_target}" == aarch64-linux-* ]]; then
    export TARGET=arm64-linux-gcc
elif [[ "${bb_full_target}" == powerpc64le-linux-* ]]; then
    export TARGET=ppc64le-linux-gcc
elif [[ "${bb_full_target}" == x86_64-apple-* ]]; then
    export TARGET=x86_64-darwin14-gcc
elif [[ "${bb_full_target}" == aarch64-apple-* ]]; then
    export TARGET=arm64-darwin20-gcc
elif [[ "${bb_full_target}" == i686-w64-mingw32* ]]; then
    export TARGET=x86-win32-gcc
elif [[ "${bb_full_target}" == x86_64-w64-mingw32* ]]; then
    export TARGET=x86_64-win64-gcc
    export CFLAGS="${CFLAGS} -fno-asynchronous-unwind-tables"
elif [[ "${bb_full_target}" == *-freebsd* ]]; then
    export TARGET=generic-gnu
fi

CONFIG_OPTS=()
if [[ "${target}" == aarch64-apple-* ]]; then
    # This feature isn't currently available for this platforms
    # check again in the future.
    CONFIG_OPTS+=(--disable-runtime-cpu-detect)
else
    CONFIG_OPTS+=(--enable-runtime-cpu-detect)
    if [[ "${target}" == *-freebsd* ]]; then
        CONFIG_OPTS+=(--disable-multithread)
    fi
fi

../configure --prefix=$prefix --target=${TARGET} \
    --as=yasm \
    --enable-postproc \
    --enable-pic \
    --enable-vp8 \
    --enable-vp9 \
    --enable-vp9-highbitdepth \
    ${CONFIG_OPTS[@]}
echo "SRC_PATH_BARE=.." >> config.mk
echo "target=libs" >> config.mk
make -j${nproc}
make install

# pkgconfig file on Windows is installed to ${libdir}/pkgconfig,
# we have to move it to ${prefix}/lib/pkgconfig/
if [[ "${target}" == *-mingw* ]] && [[ -d "${libdir}/pkgconfig" ]] ; then
    mkdir -p "${prefix}/lib"
    mv "${libdir}/pkgconfig" "${prefix}/lib/pkgconfig"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    # While we want this, it's not available on windows.  :/
    #LibraryProduct("libvpx", :libvpx),
    FileProduct("\${libdir}/libvpx.a", :libvpx),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8", julia_compat="1.6", lock_microarchitecture=false)

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using ThArrays
using Test

@testset "Simple TorchScript" begin

    @testset "Run simple method" begin
        script = """
def main(a, b):
    return a + b
"""
        cu = ThJIT.compile(script)

        ta = Tensor(rand(3, 2))
        tb = Tensor(rand(3, 2))
        # cu["main"], cu[:main], cu.main
        res = cu.main(ta, tb)

        @test res == ta + tb
    end
end


text:
function motion(v,α,dt,m₀,μ,ξ,Ω,graph)

    n = length(ξ)
    rnoise = diagm(0 => [8e-4,8e-4,1e-6])
    jcb = jacobian_motion(μ[3],v,α,dt)

    ψ = inv(jcb)-I
    ind = [[2*i+2;2*i+3] for i ∈ sort(m₀)]
    ind = vcat(ind...)
    ind = [1:3;ind]
    Ω₁ = Ω[ind,1:3]
    ψ₁ = Ω₁*ψ

    λ = zeros(n,n)
    λ[ind,1:3] = ψ₁
    λ[1:3,ind] .+= ψ₁'
    λ[1:3,1:3] .+= ψ'*Ω[1:3,1:3]*ψ
    Φ = Ω + λ

    Φₓ = Φ[1:3,ind]
    κ = zeros(n,n)
    κ[ind,ind] = Φₓ'*inv(inv(rnoise)+Φ[1:3,1:3])*Φₓ
    Ω .+= λ - κ

    δ = equation_motion(μ[3],v,α,dt)
    ξ[ind] .+= Ω[ind,1:3]*δ + (λ-κ)[ind,:]*μ
    μ[1:3] += δ

    if length(m₀)>0
        edges = [[i+1 j+1] for i ∈ m₀,j ∈ m₀ if i ≠ j]
        edges = vcat(edges...)
        if size(edges,1)>0
            sg = size(graph,1)
            adjacency = sparse(edges[:,1],edges[:,2],true,sg,sg)
            graph .|= adjacency
        end
    end
end


text:
__precompile__()

module MazeSolving

using SimpleImage
using DataStructures

export Maze

mutable struct Maze
	data::Array{UInt32, 2}
	w::Int
	h::Int
	entr::Tuple{Int, Int}
	exit::Tuple{Int, Int}
	function Maze(filepath::String)
		mz = load_maze(filepath)
		h, w = size(mz)
		entr, exit = find_endpoints(mz)
		new(mz, w, h, entr, exit)
	end
end

function load_maze(filepath::String)::Array{UInt32, 2}
	return SimpleImage.load_png(filepath)
end

function save_maze(filepath::String, maze::Array{UInt32, 2})
	SimpleImage.write_png(filepath, maze)
end

function save_maze(filepath::String, maze::Maze)
	save_maze(filepath, maze.data)
end

function find_endpoints(maze::Array{UInt32, 2})::Tuple{Tuple{Int, Int}, Tuple{Int, Int}}
	h, w = size(maze)
	points = []

	for i in 1 : h
		if maze[i, 1] == 0xffffffff
			push!(points, (1, i))
			if length(points) > 1
				return (points[1], points[2])
			end
		end
	end
	for i in 1 : w
		if maze[1, i] == 0xffffffff
			push!(points, (i, 1))
			if length(points) > 1
				return (points[1], points[2])
			end
		end
	end
	for i in 1 : h
		if maze[i, w] == 0xffffffff
			push!(points, (w, i))
			if length(points) > 1
				return (points[1], points[2])
			end
		end
	end
	for i in 1 : w
		if maze[h, i] == 0xffffffff
			push!(points, (i, h))
			if length(points) > 1
				return (points[1], points[2])
			end
		end
	end

	return ((0, 0), (0, 0))
end

function mark_endpoints(maze::Array{UInt32, 2})::Tuple{Tuple{Int, Int}, Tuple{Int, Int}}
	entr, exit = find_endpoints(maze)
	maze[entr[2], entr[1]] = 0xff00ff00
	maze[exit[2], exit[1]] = 0xff0000ff
	return entr, exit
end

function get_dirs(maze::Maze, pos::Tuple{Int, Int})::UInt8
	newdirs::UInt8 = 0

	newdirs |= maze.data[pos[2], pos[1] + 1] == 0xffffffff ? 1 : 0
	newdirs |= maze.data[pos[2] + 1, pos[1]] == 0xffffffff ? 2 : 0
	newdirs |= maze.data[pos[2], pos[1] - 1] == 0xffffffff ? 4 : 0
	newdirs |= maze.data[pos[2] - 1, pos[1]] == 0xffffffff ? 8 : 0

	return newdirs
end

get_mazecolor(maze::Maze, pos::Tuple{Int, Int})::UInt32 = maze.data[pos[2], pos[1]]

function set_mazecolor(maze::Maze, pos::Tuple{Int, Int}, color::UInt32)
	maze.data[pos[2], pos[1]] = color
end

include("depthfirst.jl")

include("bredthfirst.jl")

include("mazemanip.jl")

include("nodify.jl")

end


text:
"""
Each learn step adds to a running sum of gradients, and doesn't update θ
until `k` samples have been collected, at which point it updates with the
average gradient over those `k` samples.
"""
type OnlineGradAvg{PU <: ParamUpdater} <: LearningStrategy
    lr::Float64
    pu::PU
    idx::Int # current step
    k::Int # number of steps to average
    ∇avg::Vector{Float64}
    # OnlineGradAvg(lr::Float64, pu::PU, k::Int) = new(lr, pu, 1, k)
end

function OnlineGradAvg(k::Int; lr::Number = 1e-2, pu::ParamUpdater = RMSProp())
    OnlineGradAvg(lr, pu, 1, k, zeros(0))
end

# initialize
@with ga function pre_hook(ga::OnlineGradAvg, model::Learnable)
    idx = 1
    init(pu, model)
    ∇avg = zeros(length(grad(model)))
end

# initialize
@with ga function pre_hook(ga::OnlineGradAvg, θ::AbstractVector)
    idx = 1
    init(pu, θ)
    ∇avg = zeros(θ)
end

function learn!(model::Learnable, ga::OnlineGradAvg)
    learn!(params(model), ga, grad(model))
end

# one iteration update... we assume θ/∇ are pre-populated
@with ga function learn!(θ, ga::OnlineGradAvg, ∇)
    # add the ∇
    @simd for i in 1:length(∇)
        @inbounds ∇avg[i] += ∇[i]
    end

    # is it time to update θ?
    if idx >= k

        # convert sums to means
        scalar = 1 / k
        for i in eachindex(∇avg)
            @inbounds ∇avg[i] *= scalar
        end

        # update the params using the search direction
        update!(θ, pu, ∇avg, lr)

        # reset
        fill!(∇avg, 0.0)
        idx = 1
    else
        idx += 1
    end
end

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
const Tup = Union{Tuple, NamedTuple}
const EmptyTup = Union{Tuple{}, NamedTuple{(), Tuple{}}}

@generated function staticschema(::Type{T}) where {T}
    name_tuple = Expr(:tuple, [QuoteNode(f) for f in fieldnames(T)]...)
    type_tuple = Expr(:curly, :Tuple, [Expr(:call, :fieldtype, :T, i) for i in 1:fieldcount(T)]...)
    Expr(:curly, :NamedTuple, name_tuple, type_tuple)
end

staticschema(::Type{T}) where {T<:Tup} = T

createinstance(::Type{T}, args...) where {T} = T(args...)
createinstance(::Type{T}, args...) where {T<:Union{Tuple, NamedTuple}} = T(args)


text:
using Distributions
import EvaluationOfRLAlgs

function sample_fourierbasis(dorder_range, iorder_range, state_ranges, rng)
    sdim = size(state_ranges)[1]
    logp = 0.

    maxParams = 10000 # max params for dorder is 10k
    domp = max(floor((Float64(maxParams) + 1.0)^(1.0 / Float64(sdim)) - 1.), 0.)  # largest dorder that produces <= maxParams features
    maxDorder = min(dorder_range[2], domp)

    dorder = rand(rng, dorder_range[1]:dorder_range[2])
    logp += log(1. / (maxDorder - dorder_range[1] + 1))
    iorder = rand(rng, iorder_range[1]:iorder_range[2])
    logp += log(1. / (iorder_range[2] - iorder_range[1] + 1))

    # full = rand(rng, [true, false])
    # logp += log(0.5)
    full = false
    ϕ = EvaluationOfRLAlgs.FourierBasis(state_ranges, dorder, iorder, full)
    return ϕ, [dorder, iorder,full], logp
end

function logRand(low, high, rng)
    X = exp(rand(rng, Uniform(log(low), log(high))))
    logp = -log(log(high) - log(high) - log(X))
    return X, logp
end

function log2RandDisc(low, high, rng)
    X = round(Int,2^(rand(rng, Uniform(log2(low), log2(high)))))
    lxlow, lxhigh = log2(X), log2(X+1)
    Fxhigh = (lxhigh - log2(low)) / (log2(high) - log2(low))
    Fxlow = (lxlow - log2(low)) / (log2(high) - log2(low))
    logp = log(Fxhigh-Fxlow)
    return X, logp
end

function CreateSarsaLambda(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    if discrete_states
        qarange = [1e-3, 0.1]
    else
        qarange = [1e-6, 0.001]
    end
    αq, qalogp = logRand(qarange[1], qarange[2], rng)
    logp += qalogp
    push!(hyps, αq)
    qopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(qf, αq, γ, λ)

    agent = EvaluationOfRLAlgs.Sarsa(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateSarsaLambdaScaled(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    if discrete_states
        qarange = [1e-3, 0.1]
    else
        phi_dim = EvaluationOfRLAlgs.get_num_outputs(ϕ)
        qarange = [1e-3, 1.0] ./ phi_dim
    end
    αq, qalogp = logRand(qarange[1], qarange[2], rng)
    logp += qalogp
    push!(hyps, αq)
    qopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(qf, αq, γ, λ)

    agent = EvaluationOfRLAlgs.Sarsa(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateSarsaParl2(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    qopt = EvaluationOfRLAlgs.Parl2AccumulatingTraceOptimizer(qf, γ, λ)

    agent = EvaluationOfRLAlgs.Sarsa(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateQLambda(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    if discrete_states
        qarange = [1e-3, 0.1]
    else
        qarange = [1e-6, 0.001]
    end
    αq, qalogp = logRand(qarange[1], qarange[2], rng)
    logp += qalogp
    push!(hyps, αq)
    qopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(qf, αq, γ, λ)

    agent = EvaluationOfRLAlgs.QLearning(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateQLambdaScaled(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    if discrete_states
        qarange = [1e-3, 0.1]
    else
        phi_dim = EvaluationOfRLAlgs.get_num_outputs(ϕ)
        qarange = [1e-3, 1.0] ./ phi_dim
    end
    αq, qalogp = logRand(qarange[1], qarange[2], rng)
    logp += qalogp
    push!(hyps, αq)
    qopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(qf, αq, γ, λ)

    agent = EvaluationOfRLAlgs.QLearning(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateQParl2(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
    else
        qf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    ϵrange = [0.0, 0.1]
    ϵ = rand(rng, Uniform(ϵrange[1], ϵrange[2]))
    logp += log(1.0 / (ϵrange[2] - ϵrange[1]))
    push!(hyps, ϵ)

    qopt = EvaluationOfRLAlgs.Parl2AccumulatingTraceOptimizer(qf, γ, λ)

    agent = EvaluationOfRLAlgs.QLearning(qf, qopt, γ, ϵ)
    return agent, hyps, logp
end

function CreateActorCritic(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    else
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    if discrete_states
        varange = [1e-3, 0.1]
        parange = [1e-3, 0.1]
    else
        varange = [1e-6, 1e-3]
        parange = [1e-6, 1e-3]
    end
    αp, palogp = logRand(parange[1], parange[2], rng)
    αv, valogp = logRand(varange[1], varange[2], rng)
    logp += palogp
    logp += valogp
    push!(hyps, αp)
    push!(hyps, αv)
    popt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(policy, αp, γ, λ)
    vopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(vf, αv, γ, λ)

    agent = EvaluationOfRLAlgs.ActorCritic(policy, vf, popt, vopt, γ)
    return agent, hyps, logp
end

function CreateActorCriticScaled(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    else
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    if discrete_states
        varange = [1e-3, 0.1]
        parange = [1e-3, 0.1]
    else
        phi_dim = EvaluationOfRLAlgs.get_num_outputs(ϕ)
        varange = [1e-3, 1.0] ./ phi_dim
        parange = [1e-3, 1.0] ./ (phi_dim*num_actions)
    end
    αp, palogp = logRand(parange[1], parange[2], rng)
    αv, valogp = logRand(varange[1], varange[2], rng)
    logp += palogp
    logp += valogp
    push!(hyps, αp)
    push!(hyps, αv)
    popt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(policy, αp, γ, λ)
    vopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(vf, αv, γ, λ)

    agent = EvaluationOfRLAlgs.ActorCritic(policy, vf, popt, vopt, γ)
    return agent, hyps, logp
end

function CreateActorCriticParl2(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    else
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    if discrete_states
        parange = [1e-3, 0.1]
    else
        phi_dim = EvaluationOfRLAlgs.get_num_outputs(ϕ)
        parange = [1e-3, 1.0] ./ (phi_dim*num_actions)
    end
    αp, palogp = logRand(parange[1], parange[2], rng)
    logp += palogp
    push!(hyps, αp)
    popt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(policy, αp, γ, λ)
    vopt =  EvaluationOfRLAlgs.Parl2AccumulatingTraceOptimizer(vf, γ, λ)

    agent = EvaluationOfRLAlgs.ActorCritic(policy, vf, popt, vopt, γ)
    return agent, hyps, logp
end

function CreateNACTD(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    else
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    end
    wf = EvaluationOfRLAlgs.LinearFunction(Float64, EvaluationOfRLAlgs.get_num_params(policy), 1)

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    if discrete_states
        varange = [1e-3, 0.1]
        warange = [1e-3, 0.1]
        parange = [1e-3, 0.1]
    else
        varange = [1e-6, 1e-3]
        warange = [1e-6, 1e-3]
        parange = [1e-6, 1e-3]
    end
    αp, palogp = logRand(parange[1], parange[2], rng)
    αv, valogp = logRand(varange[1], varange[2], rng)
    αw, walogp = logRand(warange[1], warange[2], rng)
    logp += palogp
    logp += valogp
    logp += walogp
    push!(hyps, αp)
    push!(hyps, αw)
    push!(hyps, αv)
    popt = EvaluationOfRLAlgs.SGA(policy, αp)
    wopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(wf, αw, γ, λ)
    vopt = EvaluationOfRLAlgs.AccumulatingTraceOptimizer(vf, αv, γ, λ)

    # normalizew = rand(rng, [true, false])
    # logp += log(0.5)
    normalizew = true
    push!(hyps, normalizew)
    agent = EvaluationOfRLAlgs.NACTD(policy, vf, wf, popt, vopt, wopt, γ, Bool(normalizew))
    return agent, hyps, logp
end

function CreatePPO(env, rng)
    sdesc = EvaluationOfRLAlgs.get_state_desc(env)
    num_actions = EvaluationOfRLAlgs.get_action_desc(env)
    hyps = []
    logp = 0.0
    if typeof(sdesc) <: AbstractArray
        discrete_states = false
    else
        discrete_states = true
    end
    if !discrete_states
        dorder_range = [0,9]
        iorder_range = [1,9]
        ϕ, bhyps, blogp = sample_fourierbasis(dorder_range, iorder_range, sdesc, rng)
        append!(hyps, bhyps)
        logp += blogp
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    else
        vf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, 1)
        pf = EvaluationOfRLAlgs.LinearFunction(Float64, sdesc, num_actions)
        policy = EvaluationOfRLAlgs.LinearSoftmaxPolicy(pf)
    end

    # sample gamma
    Γ = EvaluationOfRLAlgs.get_gamma(env)
    γrange = (1e-4, 0.05)
    gamma, glogp = logRand(γrange[1], γrange[2], rng)
    logp += glogp
    γ = Γ - gamma
    push!(hyps, gamma)

    λrange = [0.0, 1.0]
    λ = rand(rng, Uniform(λrange[1], λrange[2]))
    logp += log(1.0 / (λrange[2] - λrange[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, λ)

    clip_range = [0.1, 0.3]
    clip_param = rand(rng, Uniform(clip_range[1], clip_range[2]))
    logp += log(1.0 / (clip_range[2] - clip_range[1]))  # logdensity for uniform distribution on [a,b]
    push!(hyps, clip_param)

    entropy_range = [1e-8, 1e-2]
    entropy_coef, elogp = logRand(entropy_range[1], entropy_range[2], rng)
    logp += elogp
    push!(hyps, entropy_coef)

    steps_range = [64, 256]
    steps_per_batch, sblogp = log2RandDisc(steps_range[1], steps_range[2], rng)
    logp += sblogp
    push!(hyps, steps_per_batch)

    epochs_range = [1,10]
    num_epochs = rand(rng, epochs_range[1]:epochs_range[2])
    logp += log(1. / (epochs_range[2] - epochs_range[1] + 1))
    push!(hyps, num_epochs)

    bsize_range = [16, 64]
    batch_size, bslogp = log2RandDisc(bsize_range[1], min(bsize_range[2], steps_per_batch), rng)
    logp += bslogp
    push!(hyps, batch_size)

    ϵ = 1e-5
    push!(hyps, ϵ)

    if discrete_states
        parange = [1e-3, 0.1]
    else
        parange = [1e-6, 1e-3]
    end
    αp, palogp = logRand(parange[1], parange[2], rng)
    logp += palogp
    push!(hyps, αp)

	popt = EvaluationOfRLAlgs.Adam(policy, αp, ϵ=ϵ)
	vopt = EvaluationOfRLAlgs.Adam(vf, αp, ϵ=ϵ)

    agent = EvaluationOfRLAlgs.PPO(policy, vf, popt, vopt, γ, λ, clip_param, entropy_coef, steps_per_batch, num_epochs, batch_size, typeof(env.state), Int)

    return agent, hyps, logp
end

function sample_alg(name, env, rng, num_episodes)

    if name == "actorcritic"
        agent, hyps, logp = CreateActorCritic(env, rng)
    elseif name == "actorcritic-scaled"
        agent, hyps, logp = CreateActorCriticScaled(env, rng)
    elseif name == "actorcritic-parl2"
        agent, hyps, logp = CreateActorCriticParl2(env, rng)
    elseif name == "nactd"
        agent, hyps, logp = CreateNACTD(env, rng)
    elseif name == "sarsalambda"
        agent, hyps, logp = CreateSarsaLambda(env, rng)
    elseif name == "sarsalambda-scaled"
        agent, hyps, logp = CreateSarsaLambdaScaled(env, rng)
    elseif name == "sarsa-parl2"
        agent, hyps, logp = CreateSarsaParl2(env, rng)
    elseif name == "qlambda"
        agent, hyps, logp = CreateQLambda(env, rng)
    elseif name == "qlambda-scaled"
        agent, hyps, logp = CreateQLambdaScaled(env, rng)
    elseif name == "q-parl2"
        agent, hyps, logp = CreateQParl2(env, rng)
    elseif name == "ppo"
        agent, hyps, logp = CreatePPO(env, rng)
    else
        println("algorithm ", name, " not yet supported")
    end

    return agent, hyps, logp
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
export run!, collect_agent_data!, collect_model_data!,
       init_agent_dataframe, init_model_dataframe, aggname,
       should_we_collect

###################################################
# Definition of the data collection API
###################################################
get_data(a, s::Symbol, obtainer::Function) = obtainer(getproperty(a, s))
get_data(a, f::Function, obtainer::Function) = obtainer(f(a))

should_we_collect(s, model, when::AbstractVector) = s ∈ when
should_we_collect(s, model, when::Bool) = when
should_we_collect(s, model, when) = when(model, s)

"""
    run!(model, agent_step! [, model_step!], n::Integer; kwargs...) → agent_df, model_df
    run!(model, agent_step!, model_step!, n::Function; kwargs...) → agent_df, model_df

Run the model (step it with the input arguments propagated into [`step!`](@ref)) and collect
data specified by the keywords, explained one by one below. Return the data as
two `DataFrame`s, one for agent-level data and one for model-level data.

## Data-deciding keywords
* `adata::Vector` means "agent data to collect". If an entry is a `Symbol`, e.g. `:weight`,
  then the data for this entry is agent's field `weight`. If an entry is a `Function`, e.g.
  `f`, then the data for this entry is just `f(a)` for each agent `a`.
  The resulting dataframe columns are named with the input symbol (here `:weight, :f`).

* `adata::Vector{<:Tuple}`: if `adata` is a vector of tuples instead,
  data aggregation is done over the agent properties.

  For each 2-tuple, the first entry is the "key" (any entry like the ones mentioned above,
  e.g. `:weight, f`). The second entry is an aggregating function that aggregates the key,
  e.g. `mean, maximum`. So, continuing from the above example, we would have
  `adata = [(:weight, mean), (f, maximum)]`.

  It's also possible to provide a 3-tuple, with the third entry being a conditional
  function (returning a `Bool`), which assesses if each agent should be included in the
  aggregate. For example: `x_pos(a) = a.pos[1]>5` with `(:weight, mean, x_pos)` will result
  in the average weight of agents conditional on their x-position being greater than 5.

  The resulting data name columns use the function [`aggname`](@ref), and create something
  like `:mean_weight` or `:maximum_f_x_pos`.
  This name doesn't play well with anonymous functions, but you can simply use
  `DataFrames.rename!` to change the returned dataframe's column names.

  **Notice:** Aggregating only works if there are agents to be aggregated over.
  If you remove agents during model run, you should modify the aggregating functions.
  *E.g.* instead of passing `mean`, pass `mymean(a) = isempty(a) ? 0.0 : mean(a)`.

* `mdata::Vector` means "model data to collect" and works exactly like `adata`.
  For the model, no aggregation is possible (nothing to aggregate over).

By default both keywords are `nothing`, i.e. nothing is collected/aggregated.

### Other keywords
* `when=true` : at which steps `s` to perform the data collection and processing.
  A lot of flexibility is offered based on the type of `when`. If `when::Vector`,
  then data are collect if `s ∈ when`. Otherwise data are collected if `when(model, s)`
  returns `true`. By default data are collected in every step.
* `when_model = when` : same as `when` but for model data.
* `obtainer = identity` : method to transfer collected data to the `DataFrame`.
  Typically only change this to [`copy`](https://docs.julialang.org/en/v1/base/base/#Base.copy)
  if some data are mutable containers (e.g. `Vector`) which change during evolution,
  or [`deepcopy`](https://docs.julialang.org/en/v1/base/base/#Base.deepcopy) if some data are
  nested mutable containers.
  Both of these options have performance penalties.
* `replicates=0` : Run `replicates` replicates of the simulation.
* `parallel=false` : Only when `replicates>0`. Run replicate simulations in parallel.
* `agents_first=true` : Whether to update agents first and then the model, or vice versa.
"""
function run! end

run!(model::ABM, agent_step!, n::Int = 1; kwargs...) =
run!(model::ABM, agent_step!, dummystep, n; kwargs...)

function run!(model::ABM, agent_step!, model_step!, n;
  replicates::Int=0, parallel::Bool=false, kwargs...)

  r = replicates
  if r > 0
    if parallel
      return parallel_replicates(model, agent_step!, model_step!, n, r; kwargs...)
    else
      return series_replicates(model, agent_step!, model_step!, n, r; kwargs...)
    end
  else
    return _run!(model, agent_step!, model_step!, n; kwargs...)
  end
end

###################################################
# Core data collection loop
###################################################
"""
  _run!(model, agent_step!, model_step!, n; kwargs...)
Core function that loops over stepping a model and collecting data at each step.
"""
function _run!(model, agent_step!, model_step!, n;
               when = true, when_model = when,
               mdata=nothing, adata=nothing, obtainer = identity,
               agents_first=true)

    df_agent = init_agent_dataframe(model, adata)
    df_model = init_model_dataframe(model, mdata)
    if n isa Integer
        if when == true; for c in eachcol(df_agent); sizehint!(c, n); end; end
        if when_model == true; for c in eachcol(df_model); sizehint!(c, n); end; end
    end

    s = 0
    while until(s, n, model)
        if should_we_collect(s, model, when)
            collect_agent_data!(df_agent, model, adata, s; obtainer = obtainer)
        end
        if should_we_collect(s, model, when_model)
            collect_model_data!(df_model, model, mdata, s; obtainer = obtainer)
        end
        step!(model, agent_step!, model_step!, 1, agents_first)
        s += 1
    end
    if should_we_collect(s, model, when)
        collect_agent_data!(df_agent, model, adata, s; obtainer = obtainer)
    end
    if should_we_collect(s, model, when_model)
        collect_model_data!(df_model, model, mdata, s; obtainer = obtainer)
    end
    return df_agent, df_model
end

###################################################
# core data collection functions per step
###################################################
"""
    init_agent_dataframe(model, adata) → agent_df
Initialize a dataframe to add data later with [`collect_agent_data!`](@ref).
"""
init_agent_dataframe(model, properties::Nothing) = DataFrame()

"""
    collect_agent_data!(df, model, properties, step = 0; obtainer = identity)
Collect and add agent data into `df` (see [`run!`](@ref) for the dispatch rules
of `properties` and `obtainer`). `step` is given because the step number information
is not known.
"""
collect_agent_data!(df, model, properties::Nothing, step::Int=0; kwargs...) = df

function init_agent_dataframe(model::ABM, properties::AbstractArray)
    nagents(model) < 1 && throw(ArgumentError("Model must have at least one agent to "*
    "initialize data collection"))

    headers = Vector{Symbol}(undef, 2+length(properties))
    headers[1] = :step
    headers[2] = :id
    for i in 1:length(properties); headers[i+2] = Symbol(properties[i]); end

    types = Vector{Vector}(undef, 2+length(properties))
    types[1] = Int[]
    types[2] = Int[]
    a = random_agent(model)
    for (i, k) in enumerate(properties)
        current_type = typeof(get_data(a, k, identity))
        isconcretetype(current_type) || warn("Type is not concrete when using $(k)"*
        "on agents. Consider narrowning the type signature of $(k).")
        types[i+2] = current_type[]
    end
    DataFrame(types, headers)
end

function collect_agent_data!(df, model, properties::Vector, step::Int=0; obtainer = identity)
    alla = sort!(collect(values(model.agents)), by=a->a.id)
    dd = DataFrame()
    dd[!, :step] = fill(step, length(alla))
    dd[!, :id] = map(a->a.id, alla)
    for fn in properties
        dd[!, Symbol(fn)] = collect(get_data(a, fn, obtainer) for a in alla)
    end
    append!(df, dd)
    return df
end

# Aggregating version
function init_agent_dataframe(model::ABM, properties::Vector{<:Tuple})
    nagents(model) < 1 && throw(ArgumentError("Model must have at least one agent to "*
    "initialize data collection"))
    headers = Vector{String}(undef, 1+length(properties))
    types = Vector{Vector}(undef, 1+length(properties))
    alla = allagents(model)

    headers[1] = "step"
    types[1] = Int[]

    for (i, property) in enumerate(properties)
        k, agg = property
        headers[i+1] = aggname(property)
        # This line assumes that `agg` can work with iterators directly
        current_type = typeof( agg( get_data(a, k, identity) for a in Iterators.take(alla,1) ) )
        isconcretetype(current_type) || warn("Type is not concrete when using function $(agg) "*
            "on key $(k). Consider using type annotation, e.g. $(agg)(a)::Float64 = ...")
        types[i+1] = current_type[]
    end
    DataFrame(types, headers)
end

"""
    aggname(k) → name
    aggname(k, agg) → name
    aggname(k, agg, condition) → name

Return the name of the column of the `i`-th collected data where `k = adata[i]`
(or `mdata[i]`).
`aggname` also accepts tuples with aggregate and conditional values.
"""
aggname(k, agg) = string(agg)*"_"*string(k)
aggname(k, agg, condition) = string(agg)*"_"*string(k)*"_"*string(condition)
aggname(x::Tuple{K,A}) where {K,A} = aggname(x[1], x[2])
aggname(x::Tuple{K,A,C}) where {K,A,C} = aggname(x[1], x[2], x[3])
aggname(x::Union{Function, Symbol, String}) = string(x)

function collect_agent_data!(df, model, properties::Vector{<:Tuple}, step::Int=0; kwargs...)
    alla = allagents(model)
    push!(df[!, 1], step)
    for (i, prop) in enumerate(properties)
        _add_col_data!(df[!, i+1], prop, alla; kwargs...)
    end
    return df
end

# Normal aggregates
function _add_col_data!(col::AbstractVector{T}, property::Tuple{K,A}, agent_iter; obtainer = identity) where {T,K,A}
    k, agg = property
    res::T = agg(get_data(a, k, obtainer) for a in agent_iter)
    push!(col, res)
end

# Conditional aggregates
function _add_col_data!(col::AbstractVector{T}, property::Tuple{K,A,C}, agent_iter; obtainer = identity) where {T,K,A,C}
    k, agg, condition = property
    res::T = agg(get_data(a, k, obtainer) for a in Iterators.filter(condition, agent_iter))
    push!(col, res)
end

# Model data
"""
    init_model_dataframe(model, mdata) → model_df
Initialize a dataframe to add data later with [`collect_model_data!`](@ref).
"""
function init_model_dataframe(model::ABM, properties::Vector)
    headers = Vector{Symbol}(undef, 1+length(properties))
    headers[1] = :step
    for i in 1:length(properties); headers[i+1] = Symbol(properties[i]); end

    types = Vector{Vector}(undef, 1+length(properties))
    types[1] = Int[]
    for (i,k) in enumerate(properties)
        types[i+1] =
            if typeof(k) <: Symbol
                typeof(model.properties[k])[]
            else
                current_type = typeof(k(model))
                isconcretetype(current_type) || warn("Type is not concrete when using $(k)"*
                "on the model. Considering narrowing the type signature of $(k).")
                current_type[]
            end
    end
    DataFrame(types, headers)
end

init_model_dataframe(model::ABM, properties::Nothing) = DataFrame()

"""
    collect_model_data!(df, model, properties, step = 0, obtainer = identity)
Same as [`collect_agent_data!`](@ref) but for model data instead.
"""
function collect_model_data!(df, model, properties::Vector, step::Int=0; obtainer = identity)
  push!(df[!, :step], step)
  for fn in properties
    push!(df[!, Symbol(fn)], get_data(model, fn, obtainer))
  end
  return df
end

collect_model_data!(df, model, properties::Nothing, step::Int=0; kwargs...) = df

###################################################
# Parallel / replicates
###################################################
function replicate_col!(df, rep)
  df[!, :replicate] = [rep for i in 1:size(df, 1)]
end

"Run replicates of the same simulation"
function series_replicates(model, agent_step!, model_step!, n, replicates; kwargs...)

  df_agent, df_model = _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...)
  replicate_col!(df_agent, 1)
  replicate_col!(df_model, 1)

  for rep in 2:replicates
    df_agentTemp, df_modelTemp = _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...)
    replicate_col!(df_agentTemp, rep)
    replicate_col!(df_modelTemp, rep)

    append!(df_agent, df_agentTemp)
    append!(df_model, df_modelTemp)
  end
  return df_agent, df_model
end

"Run replicates of the same simulation in parallel"
function parallel_replicates(model::ABM, agent_step!, model_step!, n, replicates; kwargs...)

  all_data = pmap(j -> _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...),
                  1:replicates)

  df_agent = DataFrame()
  df_model = DataFrame()
  for (rep, d) in enumerate(all_data)
    replicate_col!(d[1], rep)
    replicate_col!(d[2], rep)
    append!(df_agent, d[1])
    append!(df_model, d[2])
  end

  return df_agent, df_model
end


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

# aws ce get-cost-and-usage --time-period Start=2019-12-01,End=2019-12-28     --granularity DAILY    --metrics "BlendedCost" "UnblendedCost" "UsageQuantity"     --group-by Type=DIMENSION,Key=SERVICE Type=TAG,Key=Environment --region us-east-1
import JSON
function getTheData(
    from::DateTime,
    to::DateTime
    )

    local s = open("hidden.json") do file
        read(file, String)
    end
    return JSON.parse(s)

    # TODO: implement something here!
    # TODO: Prevent executing this multiple times!

end

function fetchData!(
    store::DataStore,
    attribute::Aws,
    ctx::DataAttributeContext,
    from::DateTime,
    to::DateTime,
)::Tuple{Array{DateTime,1},Array{Float64,1},Tuple{String,String}}

    local dataLines = Dict{String, Array{Tuple{DateTime,Float64}, 1}}()
    for event in getTheData(from, to)["ResultsByTime"]
        local start = Dates.Date(event["TimePeriod"]["Start"])
        for group in event["Groups"]
            for (metric, value) in group["Metrics"]
                local key = metric * " " * group["Keys"][1]
                if !(key in keys(dataLines))
                    dataLines[key] = []
                end
                push!(dataLines[key], (start, parse(Float64, value["Amount"])))
            end
        end
    end

    # fill it
    for event in getTheData(from, to)["ResultsByTime"]
        local start = Dates.Date(event["TimePeriod"]["Start"])
        for (k,v) in dataLines
            local found = false
            for (d,x) in v
                if d == start
                    found = true
                    break
                end
            end
            if found
                continue
            end

            push!(dataLines[k], (start, .0))
        end
    end
    for (k,v) in dataLines
        dataLines[k] = sort(v, by=x->x[1])
    end

    return pseudoColumn(dataLines[attribute.property],1), pseudoColumn(dataLines[attribute.property],2), ("Aws", attribute.property)
end

function fetchContext!(
    store::DataStore,
    attribute::Aws,
    from::DateTime,
    to::DateTime,
)::Dict{String, DataAttributeContext}
    local ctx = DataAttributeContext(Set())
    return Dict("$ctx" => ctx)
end

# @show fetchData!(DataStore(), Aws("BlendedCost Amazon DynamoDB"), DataAttributeContext(Set()), Dates.now(), Dates.now())

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
module LorenzMapDefinition

using ..PwDynamicDefinition: PwMap
using ValidatedNumerics

export LorenzMap

struct LorenzMap
    D::PwMap
end

LorenzMap(θ, α) = PwMap([x->θ*abs(x-0.5)^α, x->1-θ*abs(x-0.5)^α],
                    [@interval(0), @interval(0.5), @interval(1)])

import ..DynamicDefinition: derivative

end


text:
######################################################################
# c_interface.jl
######################################################################

# internal pointer for ccall stuff
mutable struct rtlsdr_dev
end

mutable struct RTLSDRError <: Exception
    message::String
end

# returns a pointer to
function rtlsdr_open()
    index = 0
    rd = Array{Ptr{rtlsdr_dev},1}(undef,1)
    ret = ccall( (:rtlsdr_open, "librtlsdr"),
                 Cint,
                 (Ptr{Ptr{rtlsdr_dev}}, UInt32),
                 rd,
                 index
               )

    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error opening device.")); end

    rf = rd[1]

    # resetting buffers is critical to reading bytes
    rtlsdr_reset_buffer(rf)

    return rf
end

function rtlsdr_reset_buffer(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_reset_buffer, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},),
                 rf
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error resetting buffer (error code $ret).")); end
end

function rtlsdr_close(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_close, "librtlsdr"), Cint, (Ref{rtlsdr_dev},), rf)
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error closing device (error code $ret).")); end
end

# center frequency
function rtlsdr_set_center_freq(rf::Ref{rtlsdr_dev}, freq)
    ret = ccall( (:rtlsdr_set_center_freq, "librtlsdr"), Cint, (Ref{rtlsdr_dev}, UInt32), rf, Cint(freq))
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting center frequency (error code $ret).")); end
end

function rtlsdr_get_center_freq(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_get_center_freq, "librtlsdr"), UInt32, (Ref{rtlsdr_dev},), rf)
    if ret < 0
        throw(RTLSDRError("RTLSDR.jl reports: Error getting center frequency (error code $ret)."));
    end
    return ret
end

# sample rate
function rtlsdr_get_sample_rate(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_get_sample_rate, "librtlsdr"),
                 UInt32,
                 (Ref{rtlsdr_dev},),
                 rf
               )
    if ret < 0
        throw(RTLSDRError("RTLSDR.jl reports: Error getting sample rate (error code $ret)."));
    end
    return ret
end

function rtlsdr_set_sample_rate(rf::Ref{rtlsdr_dev}, sample_rate)
    ret = ccall( (:rtlsdr_set_sample_rate, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},UInt32),
                 rf,
                 sample_rate
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error satting sample rate (error code $ret).")); end
end

# gain
function rtlsdr_set_agc_mode(rf::Ref{rtlsdr_dev}, on)
    ret = ccall( (:rtlsdr_set_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 on
                )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting AGC mode (error code $ret).")); end
end

function rtlsdr_set_tuner_gain_mode(rf::Ref{rtlsdr_dev}, manual)
    ret = ccall( (:rtlsdr_set_tuner_gain_mode, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 manual
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting tuner gain mode (error code $ret).")); end
end

function rtlsdr_set_tuner_gain(rf::Ref{rtlsdr_dev}, gain)
    ret = ccall( (:rtlsdr_set_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 gain
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting tuner gain (error code $ret).")); end
end

function rtlsdr_get_tuner_gain(rf::Ref{rtlsdr_dev})
    gain = Array{Cint,1}(undef,1)
    ret = ccall( (:rtlsdr_get_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Ptr{Cint}),
                 rf,
                 gain
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error getting tuner gain (error code $ret).")); end
    return gain[1]
end

function read_bytes(rf::Ref{rtlsdr_dev}, num_bytes)
    buf = Vector{Cuchar}(undef,num_bytes)
    bytes_read = Ref{Cint}(0)
    ret = ccall( (:rtlsdr_read_sync, "librtlsdr"),
                 Int32,
                 (Ref{rtlsdr_dev}, Ref{Cuchar}, Cint, Ref{Cint}),
                 rf,
                 buf,
                 num_bytes,
                 bytes_read
                )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error while reading from device (error code $ret).")); end

    return buf
end


text:
# This file reimplements gmsh/tutorial/t5.geo in Julia.

import gmsh

model = gmsh.model
factory = model.geo

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)

model.add("t5")

lcar1 = .1
lcar2 = .0005
lcar3 = .055

factory.addPoint(0.5,0.5,0.5, lcar2, 1)
factory.addPoint(0.5,0.5,0, lcar1, 2)
factory.addPoint(0,0.5,0.5, lcar1, 3)
factory.addPoint(0,0,0.5, lcar1, 4)
factory.addPoint(0.5,0,0.5, lcar1, 5)
factory.addPoint(0.5,0,0, lcar1, 6)
factory.addPoint(0,0.5,0, lcar1, 7)
factory.addPoint(0,1,0, lcar1, 8)
factory.addPoint(1,1,0, lcar1, 9)
factory.addPoint(0,0,1, lcar1, 10)
factory.addPoint(0,1,1, lcar1, 11)
factory.addPoint(1,1,1, lcar1, 12)
factory.addPoint(1,0,1, lcar1, 13)
factory.addPoint(1,0,0, lcar1, 14)

factory.addLine(8,9, 1);   factory.addLine(9,12, 2)
factory.addLine(12,11, 3); factory.addLine(11,8, 4)
factory.addLine(9,14, 5);  factory.addLine(14,13, 6)
factory.addLine(13,12, 7); factory.addLine(11,10, 8)
factory.addLine(10,13, 9); factory.addLine(10,4, 10)
factory.addLine(4,5, 11);  factory.addLine(5,6, 12)
factory.addLine(6,2, 13);  factory.addLine(2,1, 14)
factory.addLine(1,3, 15);  factory.addLine(3,7, 16)
factory.addLine(7,2, 17);  factory.addLine(3,4, 18)
factory.addLine(5,1, 19);  factory.addLine(7,8, 20)
factory.addLine(6,14, 21);

factory.addCurveLoop([-11,-19,-15,-18], 22)
factory.addPlaneSurface([22], 23)
factory.addCurveLoop([16,17,14,15], 24)
factory.addPlaneSurface([24], 25)
factory.addCurveLoop([-17,20,1,5,-21,13], 26)
factory.addPlaneSurface([26], 27)
factory.addCurveLoop([-4,-1,-2,-3], 28)
factory.addPlaneSurface([28], 29)
factory.addCurveLoop([-7,2,-5,-6], 30)
factory.addPlaneSurface([30], 31)
factory.addCurveLoop([6,-9,10,11,12,21], 32)
factory.addPlaneSurface([32], 33)
factory.addCurveLoop([7,3,8,9], 34)
factory.addPlaneSurface([34], 35)
factory.addCurveLoop([-10,18,-16,-20,4,-8], 36)
factory.addPlaneSurface([36], 37)
factory.addCurveLoop([-14,-13,-12,19], 38)
factory.addPlaneSurface([38], 39)

shells = []

# When the tag is not specified, a new one is automatically provided
sl = factory.addSurfaceLoop([35,31,29,37,33,23,39,25,27])
push!(shells, sl)

function cheeseHole(x, y, z, r, lc, shells)
    p1 = factory.addPoint(x,  y,  z,   lc)
    p2 = factory.addPoint(x+r,y,  z,   lc)
    p3 = factory.addPoint(x,  y+r,z,   lc)
    p4 = factory.addPoint(x,  y,  z+r, lc)
    p5 = factory.addPoint(x-r,y,  z,   lc)
    p6 = factory.addPoint(x,  y-r,z,   lc)
    p7 = factory.addPoint(x,  y,  z-r, lc)

    c1 = factory.addCircleArc(p2,p1,p7)
    c2 = factory.addCircleArc(p7,p1,p5)
    c3 = factory.addCircleArc(p5,p1,p4)
    c4 = factory.addCircleArc(p4,p1,p2)
    c5 = factory.addCircleArc(p2,p1,p3)
    c6 = factory.addCircleArc(p3,p1,p5)
    c7 = factory.addCircleArc(p5,p1,p6)
    c8 = factory.addCircleArc(p6,p1,p2)
    c9 = factory.addCircleArc(p7,p1,p3)
    c10 = factory.addCircleArc(p3,p1,p4)
    c11 = factory.addCircleArc(p4,p1,p6)
    c12 = factory.addCircleArc(p6,p1,p7)

    l1 = factory.addCurveLoop([c5,c10,c4])
    l2 = factory.addCurveLoop([c9,-c5,c1])
    l3 = factory.addCurveLoop([c12,-c8,-c1])
    l4 = factory.addCurveLoop([c8,-c4,c11])
    l5 = factory.addCurveLoop([-c10,c6,c3])
    l6 = factory.addCurveLoop([-c11,-c3,c7])
    l7 = factory.addCurveLoop([-c2,-c7,-c12])
    l8 = factory.addCurveLoop([-c6,-c9,c2])

    s1 = factory.addSurfaceFilling([l1])
    s2 = factory.addSurfaceFilling([l2])
    s3 = factory.addSurfaceFilling([l3])
    s4 = factory.addSurfaceFilling([l4])
    s5 = factory.addSurfaceFilling([l5])
    s6 = factory.addSurfaceFilling([l6])
    s7 = factory.addSurfaceFilling([l7])
    s8 = factory.addSurfaceFilling([l8])

    sl = factory.addSurfaceLoop([s1, s2, s3, s4, s5, s6, s7, s8])
    v = factory.addVolume([sl])
    push!(shells, sl)
    return v
end

x = 0
y = 0.75; z = 0; r = 0.09
for t in 1:5
    global x, z
    x += 0.166
    z += 0.166
    v = cheeseHole(x, y, z, r, lcar3, shells)
    model.addPhysicalGroup(3, [v], t)
end

factory.addVolume(shells, 186);

model.addPhysicalGroup(3, [186], 10);
factory.synchronize()
gmsh.fltk.run()

model.mesh.generate(3)
gmsh.write("t5.msh")

gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
#Context and current_context should be overloaded by users of the library! They are standard Symbols
abstract type AbstractContext end

struct Context <: AbstractContext
    id::Symbol
end
#=
We need to track the current OpenGL context.
Since we can't do this via pointer identity  (OpenGL may reuse the same pointers)
We go for this slightly ugly version.
In the future, this should probably be part of GLWindow.
=#
const GLOBAL_CONTEXT = Base.RefValue{AbstractContext}(Context(:none))

function current_context()
    GLOBAL_CONTEXT[]
end
function is_current_context(x)
    x == GLOBAL_CONTEXT[]
end
function clear_context!()
    GLOBAL_CONTEXT[] = Context(:none)
end
#this should remain here, maybe, it uses a glframebuffer
# mutable struct Context <: AbstractContext
#     context
# end
function set_context!(x)
    GLOBAL_CONTEXT[] = x
end

function exists_context()
    if current_context().id == :none
        @error "Couldn't find valid OpenGL Context. OpenGL Context active?"
    end
end

#These have to get overloaded for the pipeline to work!
swapbuffers(c::AbstractContext) = return
clear!(c::AbstractContext) = return


text:
abstract SpinModel{T}

#---
# building SpinModel from SpinInfo
##### include quantum spins
function buildSpinSystem{T}(spininfo::SpinInfo{T}, trotterinfo::TrotterInfo{T} = TrotterInfo(0.01,10))
	if isClassical(spininfo)
		if spininfo.modelname == "ising"
			return IsingModel(spininfo)
		elseif spininfo.modelname == "potts"
			return PottsModel(spininfo)
		elseif spininfo.modelname == "clock"
			return ClockModel(spininfo)
		else
			println("warning: no matching modelname")
			return SpinModel
		end
	else
		if spininfo.modelname == "ising"
			println("set trotterparameter 0.01, trotteriteration 10")
			return QuantumIsingModel(spininfo, trotterinfo)
		else
			error("model name error")
		end
	end
end

#---
# functions for SpinModel
function getMeasureOperator{T}(spinmodel::SpinModel{T})
	println("SpinModel was used")
	return [-one(T), one(T)]
end

function getFactorW(spinmodel::SpinModel)
	return spinmodel.factorW
end

function isPottsModel(spinModel::SpinModel)
	return false
end

#---
# functions for SpinModel
# inherit from SpinInfo
single_args_non_preserving = [
							  :isZeroTemperature,
							  :isClassical,
							  :getModelname,
							  :getStates,
							  :isSymmetricFactorization,
							  :getTemperature,
							  :getExternalfield,
							  :getEnvParameters
							 ]

for f in single_args_non_preserving
	eval(quote

		function $(f)(spinmodel::SpinModel)
			return $(f)(spinmodel.spininfo)
		end

	end)
end

set_single_parameter = [:setTemperature!,
						:setExternalfield!]
for f in set_single_parameter
	eval(quote

	  function $(f){T}(spinmodel::SpinModel{T}, parameter::T)
		 $(f)(spinmodel.spininfo, parameter)
		 initialize!(spinmodel)
	  end

	end)
end

function setEnvParameters!{T}(spinmodel::SpinModel{T}, temperature::T, externalfield::T)
	setEnvParameters!(spinmodel.spininfo, temperature, externalfield)
	initialize!(spinmodel)
end

#---
# functions to make factorW from Hamiltonian
function getHamiltonian{T}(spinmodel::SpinModel{T})
	println("return random hamiltonian")
	return rand(getStates(spinmodel), getStates(spinmodel))
end

function initializeFactorWasym{T}(spinmodel::SpinModel{T})
	hamiltonian = getHamiltonian(spinmodel)
	factorW = getFactorWfromHamiltonian(getTemperature(spinmodel), hamiltonian)
	return factorW
end

function getFactorWfromHamiltonian{T}(temperature::T, hamiltonian::Array{T,2})
	numberOfState = size(hamiltonian,1)
	boltzmannWeight = getBoltzmannWeight(temperature, hamiltonian)
	D, P = eig(boltzmannWeight)
	factorW = zeros(T, numberOfState, numberOfState)
	for sigma = 1:numberOfState, s = 1:numberOfState
		factorW[sigma,s] = P[sigma,s] * sqrt(D[s])
	end
	return factorW
end

function getBoltzmannWeight{T}(temperature::T, hamiltonian::Array{T,2})
	BoltzmannWeight = exp( -hamiltonian ./temperature )
	return BoltzmannWeight
end

#---
# auxilary functions

function getExpGamma{T}(revTemp::T,externalfield::T)
	expterm = exp(revTemp*externalfield/4.0)
end

#---
# testFunction
function testSpinModel()
	return testIsingModel()
end


text:
import Base.show, Base.rand
import POMDPs:
        isterminal,
        pdf,
        observation,
        generate_s,
        generate_o,
        initialstate_distribution,
        discount,
        reward,
        rand
using Random
import Random: rand
import LPDM.isterminal
using StatsFuns
using LinearAlgebra

const Vec2 = SVector{2,Float64}
Vec2() = Vec2(0.0,0.0)
Vec2Iter(S::Array{Array{Float64,1},1}) = [Vec2(s) for s in S]

# Typealias appropriately
const LD2State  = Vec2
const LD2Action = Vec2
const LD2Obs    = Vec2
const LD2Belief = LPDMBelief

# Base.show(io::IO, x::Vec2) = print(io,"[$(@sprintf("%.2f", x[1])),$(@sprintf("%.2f",x[2]))]")
Base.show(io::IO, x::Vec2) = print(io,"[$(@sprintf("%.2f", x[1])),$(@sprintf("%.2f",x[2]))]")

abstract type AbstractLD2 <: POMDP{LD2State, LD2Action, LD2Obs} end

"""
2-Dimensional light-dark problem

    min_noise::Float64
        minimum standard deviation of observation noise

    min_noise_loc::Float64
        x location of minimum observation noise

    Q::Matrix{Float64}
    R::Matrix{Float64}
        cost function matrices: reward = -s'Qs -a'Ra

    term_radius::Float64
        problem will terminate if the state is within this radius of zero

    init_dist::Any
        initial state distribution

    discount::Float64
        discount factor
"""
@with_kw mutable struct LightDark2D <: AbstractLD2
    min_noise::Float64      = 0.0
    min_noise_loc::Float64  = 5.0
    Q::Matrix{Float64}      = diagm(0=>[0.5, 0.5])
    R::Matrix{Float64}      = diagm(0=>[0.5, 0.5])
    term_radius::Float64    = 0.05
    init_dist::Any          = SymmetricNormal2([2.0, 2.0], 0.5)
    discount::Float64       = 1.0
    count::Int              = 0
end

# POMDPs.jl API functions:
POMDPs.generate_s(p::AbstractLD2, s::Vec2, a::Vec2, rng::AbstractRNG) = generate_s(p, s, a, rng, false)
POMDPs.generate_o(p::AbstractLD2, s::Vec2, a::Vec2, sp::Vec2, rng::AbstractRNG) = generate_o(p, sp, rng)
POMDPs.observation(p::AbstractLD2, a::Vec2, sp::Vec2) = observation(p, sp)
# POMDPs.isterminal(p::AbstractLD2, s::Vec2) = LinearAlgebra.norm(s) <= p.term_radius
function POMDPs.isterminal(p::AbstractLD2, s::Vec2) #DEBUG
    #NOTE: yes, this is not strictly accurate, but makes things easier
    if abs(s[1]) < p.term_radius && abs(s[2]) < p.term_radius
        return true
    else
        return false
    end
end

# Replaces the default call
function LPDM.isterminal(pomdp::AbstractLD2, particles::Vector{LPDMParticle{LD2State}})
    expected_state = Vec2(0.0,0.0)
    weight_sum = 0.0
    # LPDM.normalize!(particles)
    for p in particles
        expected_state += p.state*p.weight
        weight_sum += p.weight
    end
    expected_state /= weight_sum
    # println("expected state: $expected_state")
    return isterminal(pomdp,expected_state)
end

POMDPs.initialstate_distribution(p::AbstractLD2) = p.init_dist

function POMDPs.reward(p::AbstractLD2, s::Vec2, a::Vec2)
    if isterminal(p,s)
        r = 0.0
    else
        r = (p.reward_mode == :quadratic) ?  -(dot(s, p.Q*s) + dot(a, p.R*a)) : -1
    end
    return r
end

POMDPs.reward(p::AbstractLD2, s::Vec2, a::Vec2, sp::Vec2) = POMDPs.reward(p, s, a)
POMDPs.discount(p::AbstractLD2) = p.discount

LPDM.default_action(p::AbstractLD2) = Vec2(0.0,0.0)
LPDM.default_action(p::AbstractLD2, ::Vector{LPDMParticle{Vec2}}) = LPDM.default_action(p)

struct SymmetricNormal2D
    mean::Vec2
    std::Float64
end

# Compose action spaces via permutations of available moves
function permute(moves::Vector{Float64})::Vector{LD2Action}
    actions = Vector{Vec2}(undef,0)
    all_moves = vcat(-moves, [0.0], moves)
    for i in 1:length(all_moves)
        for j in 1:length(all_moves)
            if all_moves[i] != 0.0 || all_moves[j] != 0.0 # don't include 'idle'
                push!(actions,Vec2(all_moves[i], all_moves[j]))
            end
        end
    end
    return actions
end

# Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{SymmetricNormal2}) = d[].mean + d[].std*Vec2(rand(rng)-0.5,rand(rng)-0.5)
# Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{SymmetricNormal2}) = Vec2(rand(rng)-0.5,rand(rng)-0.5)
# POMDPs.pdf(d::SymmetricNormal2, s::Vec2) = exp(-0.5*sum((s-d.mean).^2)/d.std^2)/(2*pi*d.std^2)
POMDPs.pdf(d::SymmetricNormal2D, o::LD2Obs) = Distributions.pdf(Distributions.MvNormal([d.mean[1],d.mean[2]],[d.std,d.std]), o)
mean(d::SymmetricNormal2D) = d.mean
mode(d::SymmetricNormal2D) = d.mean
Base.eltype(::Type{SymmetricNormal2D}) = Vec2
# POMDPs.rand(p::AbstractLD2, s::LD2State, rng::LPDM.RNGVector) = rand(rng, SymmetricNormal2D(s,p.resample_std)) # for resampling
POMDPs.rand(p::AbstractLD2, s::LD2State, rng::LPDM.RNGVector) = LD2State(rand(rng, Distributions.MvNormal([s[1],s[2]],[p.resample_std,p.resample_std]))) # for resampling
# POMDPs.rand(p::AbstractLD2, a::LD2State, rng::LPDM.RNGVector) = LD2State(rand(rng, Distributions.MvNormal([a[1],a[2]],[p.action_std,p.action_std]))) # for resampling

# chose this on 2/6/17 because I like the bowtie particle patterns it produces
# unclear which one was actually used in the paper
# Masters Thesis by Pas assumes the sqrt version
obs_std(p::AbstractLD2, x::Float64) = sqrt(0.5*(p.min_noise_loc-x)^2 + p.min_noise)
# obs_std(p::AbstractLD2, x::Float64) = 0.5*(p.min_noise_loc-x)^2 + p.min_noise

function generate_s(p::AbstractLD2, s::Vec2, a::Vec2, rng::AbstractRNG, det::Bool = false)
    p.count += 1
    if p.action_std == 0.0 || det
        return s + a
    else
        return s + LPDM.rand(rng, Distributions.MvNormal([a[1],a[2]], [p.action_std,p.action_std]))
    end
end
POMDPs.observation(p::AbstractLD2, sp::Vec2) = SymmetricNormal2D(sp, obs_std(p, sp[1]))
POMDPs.observation(p::AbstractLD2, s::Vec2, a::Vec2, sp::Vec2) = observation(p,sp)

# generate_o(p::AbstractLD2, sp::Vec2, rng::AbstractRNG) = rand(rng, observation(p, sp))
function generate_o(p::AbstractLD2, sp::Vec2, rng::AbstractRNG)
    d = observation(p, sp)
    # println("sp: $sp")
    # println("d: $d")
    o = rand(rng, Distributions.MvNormal([d.mean[1],d.mean[2]],[d.std,d.std]))
    if p.obs_mode == :discrete
        o_disc = Vec2(p.bin_centers[encode(p.lindisc,o[1])],
                      p.bin_centers[encode(p.lindisc,o[2])])
        # println("generate_o: discrete")
        return o_disc
    elseif p.obs_mode == :continuous
        # println("generate_o: continuous")
        o[1]=clamp(o[1], -p.max_xy, p.max_xy)
        o[2]=clamp(o[2], -p.max_xy, p.max_xy)
        return Vec2(o)
    else
        error("Invalid obs_mode = $(p.obs_mode)")
    end
    # return obs_index(p,o_disc) # return a single combined obs index
end

# Initial distribution corresponds to how the observations would look
#TODO: possibly move to LightDarkPOMDPs
function state_distribution(pomdp::AbstractLD2, s0::LD2State, config::LPDMConfig, rng::RNGVector)
    states = Vector{LPDMParticle{LD2State}}();
    weight = 1/(config.n_particles^2) # weight of each individual particle
    particle = LPDMParticle{LD2State}(LD2State(0.0,0.0), 1, weight)

    d = observation(pomdp,s0)
    for i = 1:config.n_particles^2 #TODO: Inefficient, possibly improve. Maybe too many particles
        s = Vec2(rand(rng, Distributions.MvNormal([d.mean[1],d.mean[2]],[d.std,d.std])))
        particle = LPDMParticle{LD2State}(s,
                                          i,
                                          weight)
        push!(states, particle)
    end
    return states
end

function POMDPs.generate_sor(p::AbstractLD2, s::Vec2, a::Vec2, rng::AbstractRNG, det::Bool = false)

    s = generate_s(p, s, a, rng, det)
    o = generate_o(p, s, rng)
    r = reward(p, s, a)

    return s, o, r
end

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
@testset "1137.n-th-tribonacci-number.jl" begin
    @test tribonacci_1137(0) == 0
    @test tribonacci_1137(3) == 2
    @test tribonacci_1137(4) == 4
    @test tribonacci_1137(37) == 2082876103
end


text:
using JuMP
using CPLEX
include("parser.jl")

"""
Fonction permettant de tester si un callback est appelé car une solution entière a été obtenue
"""
function isIntegerPoint(cb_data::CPLEX.CallbackContext, context_id::Clong)

    # context_id == CPX_CALLBACKCONTEXT_CANDIDATE si le callback est appelé dans un des deux cas suivants :
    # cas 1 - une solution entière a été obtenue; ou
    # cas 2 - une relaxation non bornée a été obtenue
    if context_id != CPX_CALLBACKCONTEXT_CANDIDATE
        return false
    end

    # Pour déterminer si on est dans le cas 1 ou 2, on essaie de récupérer la
    # solution entière courante
    ispoint_p = Ref{Cint}()
    ret = CPXcallbackcandidateispoint(cb_data, ispoint_p)

    # S'il n'y a pas de solution entière
    if ret != 0 || ispoint_p[] == 0
        return false
    else
        return true
    end
end

function branch_and_cut(MyFileName::String)
  n, s, t, S, d1, d2, p, ph, d, D = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  # Il est imposé d'utiliser 1 seul thread en Julia avec CPLEX pour utiliser les callbacks
  MOI.set(m, MOI.NumberOfThreads(), 1)

  # probleme maitre
  ## Variables
  @variable(m, x[1:n,1:n], Bin) # certains arcs n'existent pas, ils ne sont pas utilises dans les contraintes
  @variable(m, y[1:n], Bin)
  @variable(m, z >= 0)

  ## Constraints
  @constraint(m,sum(x[s,j] for j in 1:n if d[s,j] != 0) == 1)
  @constraint(m,sum(x[i,t] for i in 1:n if d[i,t] != 0) == 1)
  @constraint(m,[v in 1:n; v != s && v != t], sum(x[i,v] for i in 1:n if d[i,v] != 0) == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[v in 1:n; v != t], y[v] == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[i in 1:n; i != t], y[t] == sum(x[i,t] for i in 1:n if d[i,t] != 0))
  @constraint(m, sum(p[i]*y[i] for i in 1:n) <= S)
  @constraint(m, z >= sum(d[i,j]*x[i,j] for i in 1:n, j in 1:n if d[i,j] != 0)) # objectif robuste reformule

  ## Objective
  @objective(m, Min, z)

  # Fonction exécutée à chaque callback
  #
  # Input
  # - context_id : permet de déterminer pour quelle raison le callback à
  #   été appelé (solution entière, nouvelle relaxation, ...)
  # - cb_data  :   permet  d'obtenir   diverses  autres   informations
  #   (valeur des bornes inférieures et supérieures , meilleure solution
  #   connue, ...)
  function my_callback_function(cb_data::CPLEX.CallbackContext, context_id::Clong)

      # On teste d'abord si le callback est appelé car une solution entière a été obtenue
      # (cette fonction est définie ci-dessous mais son contenu n'est pas très important)
      if isIntegerPoint(cb_data, context_id)

          # Cette ligne doit être appelée avant de pouvoir récupérer la solution entière
          CPLEX.load_callback_variable_primal(cb_data, context_id)

          # On récupère la valeur de x & y
          x_val = JuMP.callback_value.(cb_data, x)
          y_val = JuMP.callback_value.(cb_data, y)
          z_val = JuMP.callback_value.(cb_data, z)

          #### sous-probleme 1 ####
          # Create the model
          m1 = Model(CPLEX.Optimizer)

          ## Variables
          @variable(m1, delta1[1:n,1:n] >= 0)

          ## Constraints
          @constraint(m1,[i in 1:n, j in 1:n], delta1[i,j] <= D[i,j])
          @constraint(m1, sum(delta1[i,j] for i in 1:n, j in 1:n) <= d1)

          ## Objective
          @objective(m1, Max, sum(d[i,j]*(1 + delta1[i,j])*x_val[i,j] for i in 1:n, j in 1:n if d[i,j] != 0))

          #resolution
          optimize!(m1)
          z1 = JuMP.objective_value.(m1)
          delta_1 = JuMP.getvalue.( m1[:delta1] )

          #### sous-probleme 2 ####
          # Create the model
          m2 = Model(CPLEX.Optimizer)

          ## Variables
          @variable(m2, delta2[1:n] >= 0)

          ## Constraints
          @constraint(m2,[v in 1:n], delta2[v] <= 2)
          @constraint(m2, sum(delta2[v] for v in 1:n) <= d2)

          ## Objective
          @objective(m2, Max, sum((p[i] + delta2[i]*ph[i])*y_val[i] for i in 1:n))

          #resolution
          optimize!(m2)
          z2 = JuMP.objective_value.(m2)
          delta_2 = JuMP.getvalue.( m2[:delta2] )

          if z2 >= S + 1e-3 || z1 > z_val + 1e-3 || z1 < z_val - 1e-3
            con1 = @build_constraint(z >= sum(x[i,j]*d[i,j]*(1+ delta_1[i,j]) for i in 1:n, j in 1:n if d[i,j] != 0))
            MOI.submit(m, MOI.LazyConstraint(cb_data), con1)
            con2 = @build_constraint(sum(y[v]*(p[v] + ph[v]*delta_2[v]) for v in 1:n) <= S)
            MOI.submit(m, MOI.LazyConstraint(cb_data), con2)
            println("Add constraints")
          end
      end
  end

  MOI.set(m, CPLEX.CallbackFunction(), my_callback_function)
  optimize!(m)

end

max_stars_repo_name: boumpteryx/ECMA

max_issues_repo_name: boumpteryx/ECMA

max_forks_repo_name: boumpteryx/ECMA


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 095703b0-294a-4dcf-be5f-79e4c64f956e
using PlutoUI

# ╔═╡ 5a94afc0-0fc6-11ec-3c01-030a879b4c06
md"""# What do you know about data, information, and knowledge?

Before we proceed, I want to let you know that this is a safe-space for the use of the word _data_. If you hae been in academics for very long you will likely have discovered that the use of _data_ exclusively as a plural noun and the word _datum_ as the singular version is the _sine qua non_ of academic speech. However, I disagree. I quote Bill Walsh

>To many, the use of [_data_ and _media_] in their original plural form is a timus test for ltieracy, so it takes some bravery to deviate from the norm. But I think the words are now singular in most cases.
>
>_Data_ was a plural word when its singular form was still alive. When was the last time you saw _datum_ anywhere but in a discussion of the issue? The purists are trying to keep the singular form on lie support, but I say it's time to pull the plug and acknowledge that _data_ is a collective noun, like _information_.
>
>_Media_ is a plural when it refers to the plural of _medium_, but...usually it doesn't. When. an artist works in _mixed media_, that's a plural usage. Oil painting is a medium, collage art is a medium, sculpture is a medium.
>
>_The news media_ began as a similar plural term---it referred to print, radio and television. The plural usage still exists, with the online worlds as another medium, but usually when people refer to _the media_ (meaning the enws media), they're using it as a collective singular. That singular often overlaps with the print-radio-television-and-Internet plural that I described, but listen more closely and it's clear that those media are not what people have in mind. (Bill Walsh, _The Elephants of Style_)

I no longer hear gasps in the crowd when a well respected person says _data_ in that collective noun sense.

![data/datum ngram frequency in Google Books](https://github.com/chapmanbe/raw-stuff/blob/main/GoogleNgrams/data_datum.png?raw=true)

"""

# ╔═╡ 9c96a6c0-42aa-4985-84ad-0e7ccd4974a7
md"""

### Data in the news

I recently (7 September 2021) went to the [abc.com.au search site](https://search-beta.abc.net.au/index.html?siteTitle=news#/) and searched for `data`. Here are four health-related news stories that come up on the first page of the results:

- [NSW got 45 per cent of the Pfizer vaccines allocated to GPs last month, data shows](https://www.abc.net.au/news/2021-09-06/nsw-gps-got-45-per-cent-of-pfizer-covid-vaccines-last-month/100380908)
- [VIDEO: Queensland Premier accused of misrepresenting data on reopening](https://www.abc.net.au/news/2021-09-02/queensland-premier-accused-of-misrepresenting-data/13525132)
- [VIDEO: New data finds thousands of NSW healthcare staff unvaccinated](https://www.abc.net.au/news/2021-09-02/new-data-finds-thousands-of-healthcare-staff-unvaccinated/13523916)
- [Vaccination rates 'set to rise' in Mildura, as data highlights rocky start to rollout](https://www.abc.net.au/news/2021-08-31/vaccination-rates-expected-to-rise-in-mildura/100422418)

### Information in the news

- [VIDEO: Epidemiologist says authorities need to provide more information on how the virus is spreading](https://www.abc.net.au/news/2021-08-16/authorities-need-to-provide-more-info-on-how-virus-spreading/13497556)
- [Australia's official COVID-19 health information translations left to become eight weeks out of date](https://www.abc.net.au/news/2021-08-12/covid-19-information-weeks-out-of-date/100369794)
- [As Tasmania's homeless population receives coronavirus vaccines , there are calls for more accessible public health information](https://www.abc.net.au/news/2021-08-02/mobile-clinics-vaccinating-tasmanias-vulnerable-covid19/100340842)

### Knowledge in the news

Repeating the steps with knowledge, I had to go to the third page of search results to find a knowledge story related to health.

- [Local knowledge vital to contact tracing success in fight against coronavirus](https://www.abc.net.au/news/2020-08-15/local-contact-tracing-teams-needed-to-combat-coronavirus/12557504)
- [COVID-19 vaccine researchers draw on knowledge from history's 'huge disasters' and success stories](https://www.abc.net.au/news/2020-07-15/coronavirus-covid-19-vaccine-researchers-drawing-on-history/12448736)
- [How neuro-knowledge can help us find happiness and improve mental health](https://www.abc.net.au/news/2020-07-05/how-neuro-knowledge-can-help-us-find-happiness-coronavirus/12402152)

## Activity

Based on thise headlines (no need to read the full articles), how would you describe data, information, and knowledge?

$(@bind dik_news TextField((80,7)))

"""

# ╔═╡ d1ada150-cc7e-42fc-aaaa-fc9d7426e87e
let dik_news = md"""

Here are some of the things that I see in the usage.

First, _data_ is much more commonly used than either _information_ or _knowledge_. This is somewhat reflected in the Google Books ngram frequencies for the three words.

![Google N-gram frequencies for data, informationk, and knowledge](https://github.com/chapmanbe/raw-stuff/blob/main/GoogleNgrams/dik_ngrams.png?raw=true)

Second, look at how frequently _data_ seems to be an actor in the universe. _Data_ "shows" and "finds" and "highlights" seemingly without human help. By contrast _information_ and _knwoledge_ are passive.
"""

HTML("""My Answers$(html(dik_news))""")

end

# ╔═╡ d37b4cfd-0745-4e41-8813-3752ab9db1f3
md"""
## Your definitions

Before I tell you how I'm going to define data, information, and knowledge, I'd like to gather your definitions. I don't want you opening a dictionary or going to wikipedia. Imamgine you are in a car going for a long drive when suddenly your four-year-old self appears in the back seat and asks, "Mom/Dad what's 'data'?" and the same for information and knowledge. Don't think too long; you might get in an accident, or worse you'll ask yourself again.

$(@bind dik_self TextField((80,7)))

#### Type your answers here.

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─5a94afc0-0fc6-11ec-3c01-030a879b4c06
# ╟─9c96a6c0-42aa-4985-84ad-0e7ccd4974a7
# ╟─d1ada150-cc7e-42fc-aaaa-fc9d7426e87e
# ╟─d37b4cfd-0745-4e41-8813-3752ab9db1f3
# ╟─095703b0-294a-4dcf-be5f-79e4c64f956e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
include("support.jl")

function remove_dups!(ll :: LinkedList)
  record :: Set = Set()
  push!(record, ll.head.item)

  prev :: Node = ll.head      # pointer to the first occurance of some item
  curr = ll.head.next
  while curr != nothing
    if curr.item in record
      prev.next = curr.next
    else
      push!(record, curr.next.item)
      prev = curr             # update when new item appears
    end
    curr = curr.next
  end
end

# Time  : O(n)
# Space : O(n)

using Base: Test

ll = LinkedList()
append_to_tail!(ll, 1)
append_to_tail!(ll, 2)
append_to_tail!(ll, 2)
append_to_tail!(ll, 3)
append_to_tail!(ll, 3)
append_to_tail!(ll, 3)
remove_dups!(ll)
delete_node!(ll, 1)
delete_node!(ll, 2)
delete_node!(ll, 3)
@test ll.head == nothing

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
"""
Describes the reservoir bathymetry, i.e. the relationship between the reservoir level and the total storage.
The function is only valid within the bounds of the reservoir (i.e. between its minimum and maximum volumes).
The parameter is the volume and its output is the corresponding reservoir level. The function must be one-way.

`isempty(bathymetry)` allows checking whether the bathymetry is actually available, with `NullReservoirBathymetry` being
a nonexistent bathymetry (i.e. no information available).

`computeLevel(bathymetry, v)` gets the reservoir level from the given volume `v` (10^6 m^3).
`computeVolume(bathymetry, h)` get the volume from the given reservoir level `h` (m).
"""
abstract ReservoirBathymetry

immutable NullReservoirBathymetry <: ReservoirBathymetry end
Base.isempty(g::ReservoirBathymetry) = false
Base.isempty(g::NullReservoirBathymetry) = true

computeVolume(g::ReservoirBathymetry, h::Float64) =
  error("Operation not implemented: volume = f(reservoir level), for bathymetry " * string(typeof(g)))
volume(g::ReservoirBathymetry, h::Float64) = computeVolume(g, h)
computeLevel(g::ReservoirBathymetry, v::Float64) =
  error("Operation not implemented: volume = f(reservoir level), for bathymetry " * string(typeof(g)))
level(g::ReservoirBathymetry, h::Float64) = computeLevel(g, h)

computeVolume(g::ReservoirBathymetry, h::Height) = computeLevel(g, _from_unitful(h)) * 1.m^3
computeLevel(g::ReservoirBathymetry, v::Volume) = computeLevel(g, _from_unitful(v)) * 1.m

"""
This special kind of relationship is completely transparent, and can directly be used by an optimisation solver
(if the mathematical structure is allowable).
"""
abstract WhiteBoxReservoirBathymetry <: ReservoirBathymetry

iswhitebox(g::ReservoirBathymetry) = false
iswhitebox(g::WhiteBoxReservoirBathymetry) = true
isblackbox(g::ReservoirBathymetry) = ! iswhitebox(g)

"""
Among the white-box geometries are polynomial-based ones, which can be represented as follows:
    level(volume) = \sum_{i=0}^N a_{i+1} volume^i

These geometries are thus completely defined by a series of coefficients, which *should* be stored
in the `coefficients` field with increasing order (publicly accessible with the `coefficient` function).
In this case, the functions `coefficients` and `level` are automatically provided (not `volume`, as it requires
solving a polynomial equation).
"""
abstract PolynomialReservoirBathymetry <: WhiteBoxReservoirBathymetry

ispolynomial(g::ReservoirBathymetry) = false
ispolynomial(g::PolynomialReservoirBathymetry) = true

isPolynomial(g::ReservoirBathymetry) = ispolynomial(g)

"""
Retrieves the coefficient of the polynomial in `g` at the given order, i.e. power of the volume.
For example, the constant part of the polynomial has order zero, and the linear one has order one.
"""
coefficient(g::PolynomialReservoirBathymetry, order::Int) = g.coefficients[order + 1]

"""
Evaluates the given polynomial reservoir bathymetry at the given volume(s).
"""
function computeLevel(g::PolynomialReservoirBathymetry, v::Float64)
  if length(g.coefficients) == 0
    return zero(volume)
  end

  # Inefficient implementation of polynomial evaluation. Works for single numbers and vectors.
  level = g.coefficients[1]
  for i = 2:length(g.coefficients)
      level += g.coefficients[i] * v .^ (i - 1)
  end
  return level
end

immutable LinearReservoirBathymetry <: PolynomialReservoirBathymetry
  coefficients::Tuple{Float64, Float64}
  LinearReservoirBathymetry(c::Tuple{Float64, Float64}) = new(c)
  LinearReservoirBathymetry(a::Float64, b::Float64) = new((a, b))
end

getIntercept(g::LinearReservoirBathymetry) = g.coefficients[1]
getSlope(g::LinearReservoirBathymetry) = g.coefficients[2]
islinear(g::ReservoirBathymetry) = false
islinear(g::LinearReservoirBathymetry) = true

computeVolume(g::LinearReservoirBathymetry, h::Float64) = (h - getIntercept(g)) / getSlope(g)

isLinear(g::ReservoirBathymetry) = islinear(g)
intercept(g::LinearReservoirBathymetry) = getIntercept(g)
slope(g::LinearReservoirBathymetry) = getSlope(g)

"""
Quadratic reservoir bathymetry, whose form is:
    level = a + b volume + c volume^2
Both convex and concave polynomials make sense, depending on the position of the inflexion point with respect to the
part of the domain that is used:
        level                     level
          ^                         ^
          |    /                    |  ___
          |   /                     | /
          |_ /                      |/
          |                         |
          +------> vol              +------> vol
          Convex polynomial         Concave polynomial
It is the user's responsibility to ensure that this polynomial makes sense, i.e. is strictly increasing between the
reservoir capacities
"""
immutable QuadraticReservoirBathymetry <: PolynomialReservoirBathymetry
  coefficients::Tuple{Float64, Float64, Float64}
  QuadraticReservoirBathymetry(c::Tuple{Float64, Float64, Float64}) = new(c)
  QuadraticReservoirBathymetry(a::Float64, b::Float64, c::Float64) = new((a, b, c))
end

isquadratic(g::ReservoirBathymetry) = false
isquadratic(g::QuadraticReservoirBathymetry) = true

function computeVolume(g::QuadraticReservoirBathymetry, h::Float64)
  # Solve a quadratic equation:
  #    a v^2 + b v + c == h
  # But the well-known formulae only work for == 0, hence the c.
  a = coefficient(g, 2)
  b = coefficient(g, 1)
  c = coefficient(g, 0) - h
  Δ = b^2 - 4 * a * c

  if Δ < 0.
    error("The bathymetry has no real root for h = " * string(h) * ".")
  end

  v1 = (- b + sqrt(Δ)) / (2 * a)
  v2 = (- b - sqrt(Δ)) / (2 * a)

  if v1 < 0. && v2 < 0.
    error("The bathymetry has no acceptable root for h = " * string(h) * ".")
  elseif v1 != v2 && v1 > 0. && v2 > 0.
    error("The bathymetry has two acceptable roots for h = " * string(h) * ".")
  elseif (v1 == 0. && v2 <= 0.) || (v2 == 0. && v1 <= 0.)
    error("The bathymetry has one root for h = " * string(h) * ", but it is zero.")
  elseif v1 == v2
    return v1
  elseif v1 > 0. && v2 <= 0.
    return v1
  elseif v2 > 0. && v1 <= 0.
    return v2
  end
end

isQuadratic(g::ReservoirBathymetry) = isquadratic(g)

immutable PiecewiseLinearReservoirBathymetry <: WhiteBoxReservoirBathymetry
  control_points::Array{Float64, 1} # Strictly increasing.
  values::Array{Float64, 1} # Strictly increasing (reservoir level increases with volume of water).

  function PiecewiseLinearReservoirBathymetry(control_points::Array{Float64, 1}, values::Array{Float64, 1})
    if length(control_points) < 2
      error("Must provide at least a segment, i.e. two control points.")
    end

    if length(control_points) - 1 > length(values)
      error("Inconsistent set of points: too many control points for the values (exactly one more than the values is needed).")
    elseif length(control_points) - 1 < length(values)
      error("Inconsistent set of points: too few control points for the values (exactly one more than the values is needed).")
    end

    for i in 2:length(control_points)
      if control_points[i] <= control_points[i - 1]
        error("Control points must be strictly increasing: problem found near index " * string(i) * ".")
      end
    end

    for i in 2:length(values)
      if values[i] <= values[i - 1]
        error("Values must be strictly increasing: problem found near index " * string(i) * ".")
      end
    end

    return new(control_points, values)
  end
end

ispiecewiselinear(g::ReservoirBathymetry) = false
ispiecewiselinear(g::PiecewiseLinearReservoirBathymetry) = true

isPiecewiseLinear(g::ReservoirBathymetry) = ispiecewiselinear(g)
isPWL(g::ReservoirBathymetry) = ispiecewiselinear(g)

"""
Among the white-box geometries is the square root, which can be represented as follows:
    level(volume) = a_0 + \sqrt{ a_1 \times volume }
As a consequence, the volume is given by:
    volume(level) = \frac{(level - a_0)^2}{a_1}
"""
immutable SquareRootReservoirBathymetry <: WhiteBoxReservoirBathymetry
  a0::Float64
  a1::Float64
  SquareRootReservoirBathymetry(c::Tuple{Float64, Float64}) = new(c[1], c[2])
  SquareRootReservoirBathymetry(a::Float64, b::Float64) = new(a, b)
end

issquareroot(g::ReservoirBathymetry) = false
issquareroot(g::SquareRootReservoirBathymetry) = true

isSquareRoot(g::ReservoirBathymetry) = issquareroot(g)

computeLevel(g::SquareRootReservoirBathymetry, v::Float64) = g.a0 + sqrt(g.a1 * v)
computeVolume(g::SquareRootReservoirBathymetry, h::Float64) = (h - g.a0)^2 / g.a1


text:
using Ion
Ion.command_main(String[])
Ion.command_main(["-h"])
Ion.command_main(["release", "-h"])
Ion.command_main(["create", "-h"])
Ion.command_main(["clone", "-h"])
Ion.command_main(["package", "-h"])
Ion.command_main(["format", "-h"])

cd(tempdir()) do
    Ion.command_main(["clone", "Bloqade", "-f", "-y"])
    Ion.command_main(["format", "Bloqade"])
    Ion.command_main(["create", "FakePkg", "-f"])
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, classification, Datasets

Documenter.makedocs(
    sitename = "DDL Documentation",
    repo = "https://github.com/SubhadityaMukherjee/DDL.jl"
)


text:
using Snappy
using Base.Test
using Compat

const SnappyOK = Snappy.SnappyOK
const SnappyInvalidInput = Snappy.SnappyInvalidInput
const SnappyBufferTooSmall = Snappy.SnappyBufferTooSmall

# Low-level Interfaces

function originals()
    ["", "\x00", "\x00\x00", "foo", "foobarbaz", "x"^50]
end

let
    # Compress & uncompress small data using pre-allocated buffer
    buffer_size = 100
    for original in originals()
        # compress
        compressed = Array(UInt8, buffer_size)
        olen, st = Snappy.snappy_compress(original.data, compressed)
        @test st == SnappyOK
        resize!(compressed, olen)

        # uncompress
        uncompressed = Array(UInt8, buffer_size)
        olen, st = Snappy.snappy_uncompress(compressed, uncompressed)
        @test st == SnappyOK
        resize!(uncompressed, olen)

        restored = ASCIIString(uncompressed)
        @test original == restored
    end
end

let
    # Prepare too small buffer and try to compress and uncompress data
    original = "orig"

    # prepare too small buffer to compress
    compressed = Array(UInt8, 1)
    _, st = Snappy.snappy_compress(original.data, compressed)
    @test st == SnappyBufferTooSmall

    # now enough buffer size
    compressed = Array(UInt8, 100)
    olen, st = Snappy.snappy_compress(original.data, compressed)
    @test st == SnappyOK
    resize!(compressed, olen)

    # again the prepared buffer is too small to uncompress
    uncompressed = Array(UInt8, 1)
    olen, st = Snappy.snappy_uncompress(compressed, uncompressed)
    @test st == SnappyBufferTooSmall
end

let
    # Break compressed data and detect it
    original = "orig"
    compressed = Array(UInt8, 100)
    olen, st = Snappy.snappy_compress(original.data, compressed)
    @test st == SnappyOK
    resize!(compressed, olen)
    @test Snappy.snappy_validate_compressed_buffer(compressed) == SnappyOK
    # perturb the compressed buffer
    compressed[1] = 0x00
    @test Snappy.snappy_validate_compressed_buffer(compressed) == SnappyInvalidInput
end

let
    # Estimate compressed size
    for original in originals()
        maxlen = Snappy.snappy_max_compressed_length(@compat(UInt(length(original))))
        compressed = Array(UInt8, 100)
        olen, st = Snappy.snappy_compress(original.data, compressed)
        @test st == SnappyOK
        #@show Int(olen), Int(maxlen)
        @test olen <= maxlen
    end
end

let
    # Estimate uncompressed size
    for original in originals()
        compressed = Array(UInt8, 100)
        olen, st = Snappy.snappy_compress(original.data, compressed)
        @test st == SnappyOK
        resize!(compressed, olen)
        olen, st = Snappy.snappy_uncompressed_length(compressed)
        @test st == SnappyOK
        @test olen == length(original)
    end
end

# High-level Interfaces

let
    # QuickCheck-like property satisfaction tests (compress ○ uncompress = id)
    srand(2014)

    # byte arrays
    randbytes(n) = rand(UInt8, n)
    for original in map(randbytes, 0:100:10000)
        @test uncompress(compress(original)) == original
    end

    # strings
    for original in map(randstring, 0:100:1000)
        @test ASCIIString(uncompress(compress(original.data))) == original
    end
end

let
    # Large data
    # 128MiB
    srand(2014)
    original = randstring(128 * 1024^2)
    @test uncompress(compress(original.data)) == original.data
end


text:
using DelimitedFiles

# Count number of qs in a group
function count_occur(group, qs)
    n = 0

    # Loop over qs
    for q in qs

        # n++ if q in group
        if q in join(group)
            n +=1
        end
    end
    return n
end

# Count qs answered by each person
function count_each_person(group, qs)
    n = 0 # total number of questions answered by each member

    # Loop over questions answered by first member
    for q in group[1]
        dne = true
        for p in range(2,stop = size(group)[1])
            entry = group[p]
            if !(q in entry)
                # Setting to false signals this question did not exist in this response
                dne = false
            end
        end

        # If dne is still true, then increment i
        if dne == true
            n += 1
        end
    end
    return n
end

# Read file. ('new=true' for part b)
function read_file(file, qs; new=false)
    # using old or new instructions?
    if !new
        count = count_occur
    else
        count = count_each_person
    end

    open(file) do f
        ntot = 0
        group = []
        while !eof(f)
            line = readline(f)

            # Empty line denotes end of group
            if line == ""
                ntot += count(group, qs)
                group = []
            else
                group = vcat(group, line)
            end
        end

        # Count last group before EOF
        ntot += count(group, qs)
        return ntot
    end
end

##############
# START HERE #
##############

file = "./input.txt"

qs = [
    'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l',
    'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'
]
println("Old instructions:\t", read_file(file, qs))
println("New instructions:\t", read_file(file, qs, new=true))