text:
### Beta -- Bernoulli or Binomial

posterior(prior::Beta, ss::BernoulliStats) = Beta(prior.alpha + ss.cnt1, prior.beta + ss.cnt0)

posterior(prior::Beta, ss::BinomialStats) = Beta(prior.alpha + ss.ns, prior.beta + (ss.ne * ss.n - ss.ns))

function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T})
	posterior(prior, suffstats(Binomial, n, x))
end

function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64})
	posterior(prior, suffstats(Binomial, n, x, w))
end

### Dirichlet -- Categorical or Multinomial

function posterior{T<:Real}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
	α = copy(prior.alpha)
	for i in 1:length(x)
		@inbounds xi = x[i]
		α[xi] += 1.0   # cannot put @inbounds here, as no guarantee xi is inbound
	end
	return Dirichlet(α)
end

function posterior{T<:Real}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
	if length(x) != length(w)
		throw(ArgumentError("Inconsistent argument dimensions."))
	end
	α = copy(prior.alpha)
	for i in 1:length(x)
		@inbounds xi = x[i]
		@inbounds wi = w[i]
		α[xi] += wi
	end
	return Dirichlet(α)
end

function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
	return Dirichlet(prior.alpha + x)
end

# each column are counts for one experiment
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
	d::Int = dim(prior)
	if d != size(X, 1)
		throw(ArgumentError("Inconsistent argument dimensions."))
	end
	n = size(X, 2)
	α = copy(prior.alpha)
	o = 0
	for j in 1:n
		for i in 1:d
			@inbounds α[i] += X[o+i]
		end
		o += d
	end
	return Dirichlet(α)
end

function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
	d::Int = dim(prior)
	if d != size(X, 1)
		throw(ArgumentError("Inconsistent argument dimensions."))
	end
	n = size(X, 2)
	if n != length(w)
		throw(ArgumentError("Inconsistent argument dimensions."))
	end
	α = copy(prior.alpha)
	o = 0
	for j in 1:n
		@inbounds wj = w[j]
		for i in 1:d
			@inbounds α[i] += X[o+i] * wj
		end
		o += d
	end
	return Dirichlet(α)
end

### Gamma -- Exponential

function posterior(prior::Gamma, ss::ExponentialStats)
	return Gamma(prior.shape + ss.sw, 1.0 / (rate(prior) + ss.sx))
end

### For Normal distributions

function posterior(prior::Normal, ss::NormalKnownSigmaStats)
	μ0 = prior.μ
	c0 = 1.0 / abs2(prior.σ)
	c1 = 1.0 / abs2(ss.σ)

	τ = c0 + ss.tw * c1
	μ1 = (μ0 * c0 + ss.s * c1) / τ
	σ1 = sqrt(1 / τ)
	return Normal(μ1, σ1)
end

function posterior(prior::InvertedGamma, ss::NormalKnownMuStats)
	α1 = prior.shape + ss.tw / 2
	β1 = rate(prior) + ss.s2 / 2
	return InvertedGamma(α1, 1.0 / β1)
end

function posterior{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T})
	pri_μ::Normal = prior[1]
	σ::Float64 = prior[2]
	posterior(pri_μ, suffstats(NormalKnownSigma(σ), x))
end

function posterior{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T}, w::Array{Float64})
	pri_μ::Normal = prior[1]
	σ::Float64 = prior[2]
	posterior(pri_μ, suffstats(NormalKnownSigma(σ), x, w))
end

function posterior{T<:Real}(prior::(Float64, InvertedGamma), ::Type{Normal}, x::Array{T})
	μ::Float64 = prior[1]
	pri_σ::InvertedGamma = prior[2]
	posterior(pri_σ, suffstats(NormalKnownMu(μ), x))
end

function posterior{T<:Real}(prior::(Float64, InvertedGamma), ::Type{Normal}, x::Array{T}, w::Array{Float64})
	μ::Float64 = prior[1]
	pri_σ::InvertedGamma = prior[2]
	posterior(pri_σ, suffstats(NormalKnownMu(μ), x, w))
end


text:
module ImageReconstruction

export radon, iradon

"""
radon transform

https://en.wikipedia.org/wiki/Radon_transform
"""
radon

"""
inverse radon transform

https://en.wikipedia.org/wiki/Radon_transform
"""
iradon

"""
Radon transform of an image using a pixel-driven algorithm.
"""
function radon(image::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
    P = zeros(eltype(image), length(t), length(θ))
    nr, nc = size(image)

    for i = 1:nr, j = 1:nc
        x = j - nc / 2 + 0.5
        y = i - nr / 2 + 0.5
        @inbounds for (k, θₖ) in enumerate(θ)
            t′ = x * cos(θₖ) + y * sin(θₖ)

            a = convert(Int, round((t′ - t.start) / step(t) + 1))

            (a < 1 || a > length(t)) && continue
            α = abs(t′ - t[a])
            P[a, k] += (1 - α) * image[i, j]

            (a > length(t) + 1) && continue
            P[a+1, k] += α * image[i, j]
        end
    end

    P
end

end # module

ext: jl


text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type UserInformationProto
    effectiveUser::AbstractString
    realUser::AbstractString
    UserInformationProto() = (o=new(); fillunset(o); o)
end #type UserInformationProto
hash(v::UserInformationProto) = ProtoBuf.protohash(v)
isequal(v1::UserInformationProto, v2::UserInformationProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::UserInformationProto, v2::UserInformationProto) = ProtoBuf.protoeq(v1, v2)

type IpcConnectionContextProto
    userInfo::UserInformationProto
    protocol::AbstractString
    IpcConnectionContextProto() = (o=new(); fillunset(o); o)
end #type IpcConnectionContextProto
const __fnum_IpcConnectionContextProto = Int[2,3]
meta(t::Type{IpcConnectionContextProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_IpcConnectionContextProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::IpcConnectionContextProto) = ProtoBuf.protohash(v)
isequal(v1::IpcConnectionContextProto, v2::IpcConnectionContextProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::IpcConnectionContextProto, v2::IpcConnectionContextProto) = ProtoBuf.protoeq(v1, v2)

export UserInformationProto, IpcConnectionContextProto


text:
using Random, Test, HDF5

@testset "reference" begin
    data = rand(100)

    fn = tempname()

    f = h5open(fn, "w")
    f["data"] = data
    f["group_data"] = data
    f["attr_data"] = data
    # reference attached to file
    f["file_ref"] = HDF5.Reference(f, "data")
    # reference attached to group
    g = create_group(f, "sub")
    g["group_ref"] = HDF5.Reference(f, "group_data")
    # reference attached to dataset
    f["data"]["attr_ref"] = HDF5.Reference(f, "attr_data")

    close(f)

    f = h5open(fn, "r")
    # read back file-attached reference
    ref = read(f["file_ref"])
    @test ref isa HDF5.Reference
    @test data == read(f[ref])
    # read back group-attached reference
    gref = read(f["sub"]["group_ref"])
    @test gref isa HDF5.Reference
    @test data == read(f["sub"][gref])
    # read back dataset-attached reference
    aref = read(f["data"]["attr_ref"])
    @test aref isa HDF5.Reference
    @test data == read(f["data"][aref])

    close(f)
    rm(fn)
end


text:
import Base.transpose

#=
 transpose : return the transpose of Given Assoc
=#

function transpose(A::Assoc)
    return Assoc(A.col,A.row,A.val,A.A')
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using ArgParse
using Literate

function parse_commandline()
  s = ArgParseSettings()
  @add_arg_table! s begin
    "--input-file"
    help = "path to the source Julia to script to be converted into a notebook"
    arg_type = String
    required = true
    metavar = "<path>"
    "--output-dir"
    help = "path to the destination folder for the created notebook"
    arg_type = String
    required = true
    metavar = "<path>"
    "--markdown"
    help = "make a markdown file instead of a notebook"
    action = :store_true
    "--no-execution"
    help = "make a markdown or notebook file but without executing it"
    action = :store_false
  end
  return parse_args(s)
end

function main()
  parsed_args = parse_commandline()
  if parsed_args["markdown"]
    Literate.markdown(parsed_args["input-file"], outputdir=parsed_args["output-dir"]; execute=parsed_args["no-execution"])
  else
    Literate.notebook(parsed_args["input-file"], outputdir=parsed_args["output-dir"]; execute=parsed_args["no-execution"])
  end
end

main()


text:
using Statistics
using Random
using Distributions
using JSON
using POMDPs
using POMDPPolicies
using MCTS
using SCoBA
using SCoBA.SCoBADomains
using SCoBA.SCoBASolver

# Vary the number of objects as START:SKIP:STOP
const STOP = 200
const SKIP = 40
const START = 40

# For random box generation reproducibility
rng = MersenneTwister(1234)

function main(outfile::String, baseline::String="scoba", n_trials::Int64=100)

    results = Dict()

    for n_boxes = START:SKIP:STOP

        @show n_boxes

        comp_times = Float64[]
        num_nodes = Int64[]

        # Do one more trial as we will ignore the first (Julia)
        for t = 1:n_trials+1

            server, cbelt = setup_conveyor_scenario()

            if baseline == "scoba"

                # Add a batch of boxes
                for _ = 1:n_boxes
                    # choose a point in the first arm's workspace
                    start_point = rand(rng, Uniform(0.05, 0.35))
                    add_box_to_belt!(cbelt, start_point)
                end # 1:n_boxes

                # Default steps to run SCoBA
                update_time_windows!(cbelt, server)

                pickup_util(ie) = 1.0
                success_prob_fn_first(rt, sf) = cumul_grasp_prob(server.agent_set["arm1"], rt, sf)

                # Measure the time to run SCoBA tree generation
                time = @elapsed generate_search_tree!(server, "arm1",
                                                    server.agent_prop_set["arm1"].workspace_objs,
                                                    success_prob_fn_first, pickup_util,
                                                    server.agent_set["arm1"].downtime)

                # Update the statistics
                push!(comp_times, time)
                push!(num_nodes, length(server.agent_prop_set["arm1"].tree.nodes))

            elseif baseline == "hungarian"

                for _ = 1:n_boxes
                    # choose a point in the first arm's workspace
                    start_point = rand(rng, Uniform(0.05, 0.95))
                    add_box_to_belt!(cbelt, start_point)
                end # 1:n_boxes

                # Default steps to run SCoBA
                update_time_windows!(cbelt, server)

                time = @elapsed expected_hungarian!(server, cbelt)
                push!(comp_times, time)
            end

        end # for t in n_trials

        meantime = mean(comp_times[2:end])
        meannodes = Inf
        if ~isempty(num_nodes)
            meannodes = mean(num_nodes[2:end])
        end
        results[n_boxes] = Dict("meantime"=>meantime, "semtime"=>std(comp_times[2:end])/sqrt(n_trials), "medtime"=>median(comp_times[2:end]),
                               "meannodes"=>meannodes)
    end

    # Write to JSON outfile
    open(outfile, "w") do f
        JSON.print(f, results, 2)
    end
end # end main

# Example call to main
# main("conveyor_belt_scaling_scoba.json", "scoba", 100)


text:
using CSV
using Dates
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using DataFrames

data = CSV.read("input/history.csv", DataFrame)
days = Dates.value.(data[:,"time_start"] .- Date("1900-01-01"))
#dates = Date.(Dates.UTD.(days .+ Dates.value(Date("1900-01-01"))))

composition = [ r.double == "t" ? [[r.w1_id,r.w1_id*r.ground,r.w2_id,r.w2_id*r.ground],[r.l1_id,r.l1_id*r.ground,r.l2_id,r.l2_id*r.ground]] : [[r.w1_id,r.w1_id*r.ground],[r.l1_id,r.l1_id*r.ground]] for r in eachrow(data) ]

players = Set(vcat(([ r.double == "t" ? [[r.w1_id,r.w2_id],[r.l1_id,r.l2_id]] : [[r.w1_id],[r.l1_id]] for r in eachrow(data) ]...)...))

priors = Dict([(p, ttt.Player(ttt.Gaussian(0., 1.6), 1.0, 0.036) ) for p in players])

function fit(s,g)
    h = ttt.History(composition=composition, times = days, sigma = s, gamma = g, beta = 0.0, priors = priors)
    ttt.convergence(h,epsilon=0.01, iterations=10)
    return h
end
@time h = fit(1.0,0.01)
ttt.log_evidence(h)
lc = ttt.learning_curves(h)

federer = "f324"
nadal = "n409"
djokovic = "d643"

dict = Dict{String,Vector}()
dict["federer_mu"] = [tp[2].mu for tp in lc["f324"]]
dict["federer_sigma"]= [tp[2].sigma for tp in lc["f324"]]
dict["federer_time"] = [tp[1] for tp in lc["f324"]]
dict["federer_mu_hard"] = [tp[2].mu for tp in lc["f324Hard"]]
dict["federer_sigma_hard"]= [tp[2].sigma for tp in lc["f324Hard"]]
dict["federer_time_hard"]= [tp[1] for tp in lc["f324Hard"]]
dict["federer_mu_clay"] = [tp[2].mu for tp in lc["f324Clay"]]
dict["federer_sigma_clay"]= [tp[2].sigma for tp in lc["f324Clay"]]
dict["federer_time_clay"]= [tp[1] for tp in lc["f324Clay"]]
dict["federer_mu_carpet"] = [tp[2].mu for tp in lc["f324Carpet"]]
dict["federer_sigma_carpet"]= [tp[2].sigma for tp in lc["f324Carpet"]]
dict["federer_time_carpet"]= [tp[1] for tp in lc["f324Carpet"]]
dict["federer_mu_grass"] = [tp[2].mu for tp in lc["f324Grass"]]
dict["federer_sigma_grass"]= [tp[2].sigma for tp in lc["f324Grass"]]
dict["federer_time_grass"]= [tp[1] for tp in lc["f324Grass"]]

dict["nadal_mu"] = [tp[2].mu for tp in lc["n409"]]
dict["nadal_sigma"] = [tp[2].sigma for tp in lc["n409"]]
dict["nadal_time"] = [tp[1] for tp in lc["n409"]]
dict["nadal_mu_hard"] = [tp[2].mu for tp in lc["n409Hard"]]
dict["nadal_sigma_hard"] = [tp[2].sigma for tp in lc["n409Hard"]]
dict["nadal_time_hard"] = [tp[1] for tp in lc["n409Hard"]]
dict["nadal_mu_clay"] = [tp[2].mu for tp in lc["n409Clay"]]
dict["nadal_sigma_clay"] = [tp[2].sigma for tp in lc["n409Clay"]]
dict["nadal_time_clay"] = [tp[1] for tp in lc["n409Clay"]]
dict["nadal_mu_carpet"] = [tp[2].mu for tp in lc["n409Carpet"]]
dict["nadal_sigma_carpet"] = [tp[2].sigma for tp in lc["n409Carpet"]]
dict["nadal_time_carpet"] = [tp[1] for tp in lc["n409Carpet"]]
dict["nadal_mu_grass"] = [tp[2].mu for tp in lc["n409Grass"]]
dict["nadal_sigma_grass"] = [tp[2].sigma for tp in lc["n409Grass"]]
dict["nadal_time_grass"] = [tp[1] for tp in lc["n409Grass"]]

dict["djokovic_mu"] = [tp[2].mu for tp in lc["d643"]]
dict["djokovic_sigma"] = [tp[2].sigma for tp in lc["d643"]]
dict["djokovic_time"] = [tp[1] for tp in lc["d643"]]
dict["djokovic_mu_hard"] = [tp[2].mu for tp in lc["d643Hard"]]
dict["djokovic_sigma_hard"] = [tp[2].sigma for tp in lc["d643Hard"]]
dict["djokovic_time_hard"] = [tp[1] for tp in lc["d643Hard"]]
dict["djokovic_mu_clay"] = [tp[2].mu for tp in lc["d643Clay"]]
dict["djokovic_sigma_clay"] = [tp[2].sigma for tp in lc["d643Clay"]]
dict["djokovic_time_clay"] = [tp[1] for tp in lc["d643Clay"]]
dict["djokovic_mu_carpet"] = [tp[2].mu for tp in lc["d643Carpet"]]
dict["djokovic_sigma_carpet"] = [tp[2].sigma for tp in lc["d643Carpet"]]
dict["djokovic_time_carpet"] = [tp[1] for tp in lc["d643Carpet"]]
dict["djokovic_mu_grass"] = [tp[2].mu for tp in lc["d643Grass"]]
dict["djokovic_sigma_grass"] = [tp[2].sigma for tp in lc["d643Grass"]]
dict["djokovic_time_grass"] = [tp[1] for tp in lc["d643Grass"]]

maxlen = maximum([length(value) for (key, value) in dict])

using DataFrames
df = DataFrame(Dict(key => [value;repeat([missing],maxlen-length(value))] for (key, value) in dict))

CSV.write("output/atp_ground_learning_curves.csv", df; header=true)

@time h_nadal_optimum = fit(0.35,0.0)
lc = ttt.learning_curves(h_nadal_optimum )
lc[nadal*"Clay"][end][2].mu - lc[nadal*"Hard"][end][2].mu
lc[nadal*"Clay"][end][2].mu - lc[nadal*"Grass"][end][2].mu
lc[djokovic*"Clay"][end][2].mu - lc[djokovic*"Grass"][end][2].mu
lc[djokovic*"Clay"][end][2].mu - lc[djokovic*"Hard"][end][2].mu


text:
# Julia wrapper for header: /usr/local/include/libavcodec/vdpau.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_alloc_vdpaucontext,
    av_vdpau_hwaccel_get_render2,
    av_vdpau_hwaccel_set_render2,
    av_vdpau_bind_context,
    av_vdpau_get_surface_parameters,
    av_vdpau_alloc_context,
    av_vdpau_get_profile

function av_alloc_vdpaucontext()
    ccall((:av_alloc_vdpaucontext,libavcodec),Ptr{AVVDPAUContext},())
end

function av_vdpau_hwaccel_get_render2(arg1)
    ccall((:av_vdpau_hwaccel_get_render2,libavcodec),AVVDPAU_Render2,(Ptr{AVVDPAUContext},),arg1)
end

function av_vdpau_hwaccel_set_render2(arg1,arg2::AVVDPAU_Render2)
    ccall((:av_vdpau_hwaccel_set_render2,libavcodec),Void,(Ptr{AVVDPAUContext},AVVDPAU_Render2),arg1,arg2)
end

function av_vdpau_bind_context(avctx,device::Integer,get_proc_address,flags::Integer)
    ccall((:av_vdpau_bind_context,libavcodec),Cint,(Ptr{AVCodecContext},Cint,Ptr{Cint},UInt32),avctx,device,get_proc_address,flags)
end

function av_vdpau_get_surface_parameters(avctx,_type,width,height)
    ccall((:av_vdpau_get_surface_parameters,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Cint},Ptr{UInt32},Ptr{UInt32}),avctx,_type,width,height)
end

function av_vdpau_alloc_context()
    ccall((:av_vdpau_alloc_context,libavcodec),Ptr{AVVDPAUContext},())
end

function av_vdpau_get_profile(avctx,profile)
    ccall((:av_vdpau_get_profile,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{Cint}),avctx,profile)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
function optical_path_bending(target::AbstractTarget)
end


text:
module NineTTT
include("./game.jl")

end


text:
# This file is a part of DataDeps.jl. License is MIT.

"""
    splitpath(path)

The opposite of `joinpath`,
splits a path unto each of its directories names / filename (for the last).
"""
function splitpath(path::AbstractString)
   ret=String[]
   prev_path = path
   while(true)
       path, lastpart = splitdir(path)
       length(lastpart)>0 && pushfirst!(ret, lastpart)
       length(path)==0 && break
       if prev_path==path
            # catch the casewhere path begins with a root
            pushfirst!(ret, path)
            break
       end
       prev_path = path
   end
   return ret
end

########################################
# Enviroment variable stuff

"""
    env_bool(key)

Checks for an enviroment variable and fuzzy converts it to a bool
"""
env_bool(key, default=false) = haskey(ENV, key) ? lowercase(ENV[key]) ∉ ["0","","false", "no"] : default

"""
    env_list(key)

Checks for an enviroment variable and converts it to a list of strings, sperated with a colon
"""
env_list(key, default=String[]) = haskey(ENV, key) ? split(ENV[key], ":") : default

#########################################
# User input stuff

"""
    bool_input

Prompted the user for a yes or no.
"""
function input_bool(prompt="")::Bool
    input_choice(prompt, 'y','n')=='y'
end

"""
    input_choice

Prompted the user for one of a list of options
"""
function input_choice(prompt, options::Vararg{Char})::Char
    while(true)
        println(prompt)
        println("["*join(options, '/')*"]")
        response = readline()
        length(response)==0 && continue
        reply = lowercase(first(response))
        for opt in lowercase.(options)
            reply==opt && return opt
        end
    end
end

"""
    input_choice

Prompts the user for one of a list of options.
Takes a vararg of tuples of Letter, Prompt, Action (0 argument function)

Example:
```
input_choice(
    ('A', "Abort -- errors out", ()->error("aborted")),
    ('X', "eXit -- exits normally", ()->exit()),
    ('C', "Continue -- continues running", ()->nothing)),
)

```
"""
function input_choice(options::Vararg{Tuple{Char, <:AbstractString, Any}})
    acts = Dict{Char, Any}()
    prompt = ""
    chars = Char[]
    for (cc, prmt, act) in options
        prompt*="\n [$cc] $prmt"
        push!(chars, cc)
        acts[lowercase(cc)] = act
    end
    prompt*="\n"

    acts[input_choice(prompt, chars...)]()
end


text:
## Makes a random walk Metropolis kernel
function makeRWMKernel(logTargetDensity::F,
  Σ::SMatrix{d, d, Float64}, rng::RNG=Random.GLOBAL_RNG) where
  {F<:Function, d, RNG<:AbstractRNG}

  A::SMatrix{d, d, Float64} = cholesky(Σ).L
  scratchv::MVector{d, Float64} = MVector{d, Float64}(undef)
  prevx::MVector{d, Float64} = MVector{d, Float64}(undef)
  ldprevx = Ref(-Inf)
  accepts = Ref(0)
  calls = Ref(0)

  function P(x::SVector{d, Float64})
    calls.x += 1
    randn!(rng, scratchv)
    scratchv .= A * scratchv
    z::SVector{d, Float64} = x + scratchv
    if x == prevx
      lpi_x = ldprevx.x
    else
      lpi_x = logTargetDensity(x)
      prevx .= x
    end
    lpi_z = logTargetDensity(z)
    if -randexp(rng) < lpi_z - lpi_x
      prevx .= z
      ldprevx.x = lpi_z
      accepts.x += 1
      rval = z
    else
      rval = x
    end
    return rval
  end

  function P(s::Symbol)
    s == :acceptanceRate && return accepts.x / calls.x
  end

  return P
end


text:
# This dictionary maps easy to remember names to Youtube video IDs
# after adding an ID here, you can use the {{youtube <shortname>}}
# syntax in your markdown files to embed the video into the page!
videos = Dict(
              "installation" => "OOjKEgbt8AI",
              "jupyter-notebooks" => "MFgMO8Xwx-k",
              "session1" => "x5o0nTpDbT8",
              "session2" => "sHCqosMRfJ8",
              "push" => "iQJQR0bsnDc",
              "session3" => "bDiBuRLxdSY",
              "solverpaths" => "PUlOQdgWlFE",
              "constrained-banana" => "kqjr_NYUomU",
              "pr-video" => "kKjmrAExsNQ"
)


text:
function turing_inference(prob::DiffEqBase.DEProblem,alg,t,data,priors;
  likelihood_dist_priors = [InverseGamma(2, 3)],
  likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
  num_samples=1000, sampler = Turing.NUTS(0.65),
  syms = [Symbol("theta$i") for i in 1:length(priors)],
  kwargs...)

  N = length(priors)
  mf(vi,sampler,model) = begin
    # If x is provided, use the provided values.
    # Otherwise, treat x as an empty vector with
    # two entries.
    x = model.defaults.x

    N = length(priors)
    _theta = Vector(undef,length(priors))

    for i in 1:length(priors)
      _theta[i], __lp = Turing.assume(sampler,
                        priors[i],
                        Turing.VarName([:mf, syms[i]], ""),
                        vi)
      vi.logp += __lp
    end

    theta = convert(Array{typeof(first(_theta))},_theta)
    _σ = Vector(undef,length(likelihood_dist_priors))

    for i in 1:length(likelihood_dist_priors)
      _σ[i], __lp = Turing.assume(sampler,
                      likelihood_dist_priors[i],
                      Turing.VarName([:mf, Symbol("σ$i")], ""),
                      vi)
      vi.logp += __lp
    end
    σ = convert(Array{typeof(first(_σ))},_σ)

    p_tmp = remake(prob, u0 = convert.(eltype(theta), (prob.u0)), p = theta)
    _saveat = t === nothing ? Float64[] : t
    sol_tmp = solve(p_tmp, alg; saveat = _saveat, kwargs...)

    if sol_tmp isa DiffEqBase.AbstractEnsembleSolution
      failure = any((s.retcode != :Success for s in sol_tmp)) && any((s.retcode != :Terminated for s in sol_tmp))
    else
      failure = sol_tmp.retcode != :Success && sol_tmp != :Terminated
    end

    if failure
      vi.logp -= Inf
    else
      if sol_tmp isa DiffEqBase.AbstractNoTimeSolution
        res = sol_tmp.u
        __lp = Turing.observe(
          sampler,
          likelihood(res,theta,Inf,σ),   # Distribution
          x,    # Data point
          vi
        )
        vi.logp += __lp
      else
        for i = 1:length(t)
          res = sol_tmp.u[i]
          _t   = sol_tmp.t[i]
          __lp = Turing.observe(
            sampler,
            likelihood(res,theta,_t,σ),   # Distribution
            x[:,i],    # Data point
            vi
          )
          vi.logp += __lp
        end
      end
    end
    vi
  end

  bigtup = Tuple{([[syms[i] for i in 1:length(priors)]; [Symbol("σ$i") for i in 1:length(likelihood_dist_priors)]])...}

  # Define the default value for x when missing
  defaults = (x = data,)

  # Instantiate a Model object.
  model = Turing.Model{bigtup, Tuple{:x}}(mf, data, defaults)

  chn = sample(model, sampler, num_samples)
end


text:
#!/usr/bin/env julia

#Dependencies
using JuliaDB
using HTTP
using CodecZlib
using EzXML
using Gumbo
using Cascadia
using Dates

#Functions
function PareDownIO(MyIOSTREAM) #Reduce Size of XML Files (keep RAM usage as small as possible)
	ParedIO=IOBuffer(append=true);

		for LocLine in eachline(MyIOSTREAM)
			if occursin("credit",LocLine) || occursin("id",LocLine) || occursin("model",LocLine) || occursin("host",LocLine) || occursin("coproc",LocLine) || occursin("xml",LocLine)
				println(ParedIO, string(LocLine));
			end
		end

	return ParedIO
	close(ParedIO)
end

function MyStreamXMLparse(XMLstream,OutFile) #Read XML files line by line and extract desired values (avoids memory leak from libXML2)
reader = EzXML.StreamReader(XMLstream)
	#println("Loaded Reader")
	#Arrays to accumulate needed data
	LocHostID = [];
	LocTotCred = [];
	LocRAC = [];
	LocPModel = [];
	LocGModel = [];

	#Flags to indicate if all values where found for each host
	foundID=true
	foundTotCred=true
	foundRAC=true
	foundPModel=true
	foundGModel=true

	for line in reader
		if (reader.type==1) && (reader.name=="host") # If we find new host reset flags for new node
			if foundID==false
				push!(LocHostID,"0")
			end
			if foundTotCred==false
				push!(LocTotCred,"0")
			end
			if foundRAC==false
				push!(LocRAC,"0")
			end
			if foundPModel==false
				push!(LocPModel,"NONE")
			end
			if foundGModel==false
				push!(LocGModel,"NONE")
			end
			foundID=false
			foundTotCred=false
			foundRAC=false
			foundPModel=false
			foundGModel=false
		end
		#Locate which data type we found and save value
		if (reader.type==1) && (reader.name=="id") # Type=1 is an element/node
			push!(LocHostID,reader.content)
			foundID=true
		end
		if (reader.type==1) && (reader.name=="total_credit") # Type=1 is an element/node
			push!(LocTotCred,reader.content)
			foundTotCred=true
		end
		if (reader.type==1) && (reader.name=="expavg_credit") # Type=1 is an element/node
			push!(LocRAC,reader.content)
			foundRAC=true
		end
		if (reader.type==1) && (reader.name=="p_model") # Type=1 is an element/node
			#println(string(reader.content))
			push!(LocPModel,uppercase(string(reader.content)))
			foundPModel=true
		end
		if (reader.type==1) && (reader.name=="coprocs") # Type=1 is an element/node
			push!(LocGModel,uppercase(string(reader.content)))
			foundGModel=true
		end
	end
	#Final cleanup in case last host had missing data
	if foundID==false
		push!(LocHostID,"0")
	end
	if foundTotCred==false
		push!(LocTotCred,"0")
	end
	if foundRAC==false
		push!(LocRAC,"0")
	end
	if foundPModel==false
		push!(LocPModel,"NONE")
	end
	if foundGModel==false
		push!(LocGModel,"NONE")
	end

	LocHostID=parse.([Int64],LocHostID)
	LocTotCred=parse.([Float64],LocTotCred)
	LocRAC=parse.([Float64],LocRAC)

	LocalTable=table(LocHostID,LocPModel,LocGModel,LocTotCred,LocRAC; names = [:ID, :CPUmodel, :GPUmodel, :TotCred, :RAC]);
	LocalTable=filter(host -> host.RAC > 1.0 , LocalTable)		#Remove any inactive hosts from database
	save(LocalTable,OutFile)

end

###
### Start Main run
###

WhiteListFile=joinpath(".","WhiteList.csv");# Import Gridcoin WhiteList from CSV file
println("Reading $WhiteListFile")
WhiteListTable=loadtable(WhiteListFile);
WLlength=size(select(WhiteListTable,1),1);

#Remove old host data
FullHostFilePath=joinpath(pwd(),"HostFiles");
if isdir(FullHostFilePath)
	if Sys.iswindows()
		run(`cmd /C rmdir /Q /S $FullHostFilePath`)	#Windows File Permisions issue (workaround)
	else
		rm("HostFiles"; force=true, recursive=true);
	end
end
mkdir("HostFiles")					#Make new folder to store host data

#Check with block explorer to verify greylist/TeamRAC
statsURL="https://www.gridcoinstats.eu/project";
statsHTML=joinpath(tempdir(),"stats.html")

if Sys.iswindows()
	run(`cmd /C curl $statsURL -s -o $statsHTML`)
else
	run(`wget $statsURL -q -O $statsHTML`);
end

HTMLdat=Gumbo.parsehtml(read(statsHTML,String));
HTMLtab=eachmatch(Selector("tr"), HTMLdat.root);
TableLines=size(HTMLtab,1);
CurrentWLsize=TableLines-1;
WLTab_RACvect=[ Inf for ind=1:WLlength];
for lineNum = 2:TableLines
	line=HTMLtab[lineNum];
	ProjName=string(line[1][1][1]);
	TeamRAC=parse(Int64,replace(string(line[7][1][1]),' ' => ""));
	locIndex=findall(x-> x==ProjName, select(WhiteListTable,:FullName));

	if ~isempty(locIndex)
		WLTab_RACvect[locIndex[1]]=TeamRAC;
	else
		#Print notice if QuickMag whitelist disagrees with block explorer
		#Einstein does not publish host data, so it is not listed in WhiteList.csv
		println("    Not building host database for: $ProjName")
	end

end

rm(statsHTML);

WhiteListTable=JuliaDB.pushcol(WhiteListTable, :TeamRAC, WLTab_RACvect)
WhiteListTable=JuliaDB.pushcol(WhiteListTable, :NumWL_Proj, [CurrentWLsize for ind=1:WLlength]);
WhiteListTable=JuliaDB.pushcol(WhiteListTable, :TimeStamp, [Dates.now() for ind=1:WLlength]);
GreyList=findall(x-> x==Inf,select(WhiteListTable,:TeamRAC))	#Print notice if projects are on greylist
for line in GreyList
	println("    Project on greylist: $(WhiteListTable[line].FullName)")
end
println("")

println("Downloading Host Data")

FailedDownloads=[];		#Vector to keep track of any failed downloads

Threads.@threads for ind=1:WLlength	#Process projects in WhiteList.csv (Runs in parallel if julia started with multiple threads)
		row=WhiteListTable[ind];
		if (row.TeamRAC!=Inf)	#Skip if project if there was no current team data available
			println("    Starting to download project: $(row.Project), $(row.Type), $(row.URL)")

			LocFilePath=joinpath(".","HostFiles","$(row.Type)"*"_"*"$(row.Project).jldb") #Path to saved data

			try
				#Alternate between downloading & processing host data
				CompressedFileStream = Base.BufferStream();
				@async while !eof(CompressedFileStream)
					MyStreamXMLparse(PareDownIO(GzipDecompressorStream(CompressedFileStream)),LocFilePath) #Remove most unnecessary elements from XML to save RAM & Convert XML to binary JuliaDB file
				end

				#Download host data and send to CompressedFileStream
				Request = HTTP.get(row.URL, response_stream=CompressedFileStream,verbose=0, connect_timeout=30, retries=3)

				#Wait until there is no data left in CompressedFileStream
				close(CompressedFileStream)
				println("    Finished downloading project: $(row.Project)")

			catch e					#catch errors that occur if a project website is down
				println("Error: Unable to download data for $(row.Project)")
				push!(FailedDownloads,ind)
			end

		end
end

# Finalize WhiteListTable.jldb noting missing data (Host data & Team data from block explorer)
TeamRacVect=select(WhiteListTable,:TeamRAC);
for jnd in FailedDownloads
	TeamRacVect[jnd] = Inf ;
end
WhiteListTable=JuliaDB.popcol(WhiteListTable, :TeamRAC)
WhiteListTable=JuliaDB.pushcol(WhiteListTable, :TeamRAC, TeamRacVect)
save(WhiteListTable,joinpath(".","HostFiles","WhiteList.jldb")) #Save checked and parsed WhiteListTable for quicker access


text:
#! /usr/bin/julia

# Rosetta Code, Problem of Apollonius

using Polynomials

immutable Point{T<:Real}
    x::T
    y::T
end

immutable Circle{T<:Real}
    c::Point{T}
    r::T
end
Circle{T<:Real}(x::T, y::T, r::T) = Circle(Point(x,y), r)

function apollonius{T<:Real}(ap::Array{Circle{T},1},
                             enc::Array{Int,1}=Int[])
    length(ap) == 3 || error("This Apollonius problem needs 3 circles.")
    x = map(u->u.c.x, ap)
    y = map(u->u.c.y, ap)
    r = map(u->u in enc ? -1 : 1, 1:3) .* map(u->u.r, ap)
    a = 2x[1] - 2x[2:3]
    b = 2y[1] - 2y[2:3]
    c = 2r[1] - 2r[2:3]
    d = (x[1]^2 + y[1]^2 - r[1]^2)
    d -= (map(u->u^2, x) + map(u->u^2, y) - map(u->u^2, r))[2:3]
    u = Poly([-det(hcat(b,d)), det(hcat(b,c))]/det(hcat(a,b)))
    v = Poly([det(hcat(a,d)), -det(hcat(a,c))]/det(hcat(a,b)))
    w = Poly([r[1], 1.0])^2
    s = (u - x[1])^2 + (v - y[1])^2 - w
    r = filter(x->imag(x)==0 && 0<x, roots(s))
    length(r) < 2 || error("The solution is not unique.")
    length(r) == 1 || error("There is no solution.")
    r = r[1]
    Circle(polyval(u,r), polyval(v,r), r)
end

function detail(c::Circle)
    x = @sprintf "%.4f" c.c.x
    y = @sprintf "%.4f" c.c.y
    r = @sprintf "%.4f" c.r
    @sprintf "centered at (%s, %s) with radius %s" x y r
end

test = [Circle(0.,0.,1.),
        Circle(4.,0.,1.),
        Circle(2.,4.,2.)]

println("The defining circles are:")
for c in test
    println("    ", detail(c))
end
println("The internal circle is:")
println("    ", detail(apollonius(test)))
println("The external circle is:")
println("    ", detail(apollonius(test, [1:3])))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using SupervisedLearning

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

symbols = [:error, :uninitialized, :initialized, :training, :trystop, :trained, :converged]
@test SupervisedLearning.STATES == symbols

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

msg("State: constructors")

@test_throws MethodError State(1, :error)
@test_throws KeyError State(:foo)
@test_throws BoundsError State(0)
@test_throws BoundsError State(8)

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

msg(string("State: execute show: ", State(2)))

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

msg("State: construction from symbols")

for sym in symbols
  @test State(sym) == convert(State, sym)
  @test State(sym) == State(sym)
  @test State(sym) <= State(sym)
  @test State(sym) >= State(sym)
  @test State(sym) == sym
  @test State(sym) <= sym
  @test State(sym) >= sym
  @test sym == State(sym)
  @test sym <= State(sym)
  @test sym >= State(sym)
end

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

msg("State: construction from numbers")

for i in 1:length(symbols)
  @test State(i) == convert(State, i)
  @test State(i) == State(i)
  @test State(i) <= State(i)
  @test State(i) >= State(i)
  @test State(i) == i
  @test State(i) <= i
  @test State(i) >= i
  @test i == State(i)
  @test i <= State(i)
  @test i >= State(i)
end

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

msg("State: comparing with numbers and symbols")

for (i, sym) in enumerate(symbols)
  @test !(State(i) != State(sym))
  @test !(State(i) < State(sym))
  @test !(State(i) > State(sym))
  @test State(i) == State(sym)
  @test State(i) <= State(sym)
  @test State(i) >= State(sym)
  @test !(State(sym) != State(i))
  @test !(State(sym) < State(i))
  @test !(State(sym) > State(i))
  @test State(sym) == State(i)
  @test State(sym) <= State(i)
  @test State(sym) >= State(i)
  @test !(State(i) != sym)
  @test !(State(i) < sym)
  @test !(State(i) > sym)
  @test State(i) == sym
  @test State(i) <= sym
  @test State(i) >= sym
  @test !(sym != State(i))
  @test !(sym < State(i))
  @test !(sym > State(i))
  @test sym == State(i)
  @test sym <= State(i)
  @test sym >= State(i)
  @test State(sym) == i
  @test State(sym) <= i
  @test State(sym) >= i
  @test !(State(sym) != i)
  @test !(State(sym) < i)
  @test !(State(sym) > i)
  @test sym == State(i)
  @test sym <= State(i)
  @test sym >= State(i)
  @test !(sym != State(i))
  @test !(sym < State(i))
  @test !(sym > State(i))

  for j = 1:(i-1)
    @test State(j) != sym
    @test State(j) < sym
    @test State(j) <= sym
    @test !(State(j) > sym)
    @test !(State(j) == sym)
    @test !(State(j) >= sym)
    @test sym != State(j)
    @test sym > State(j)
    @test sym >= State(j)
    @test !(sym < State(j))
    @test !(sym == State(j))
    @test !(sym <= State(j))
    @test State(sym) != j
    @test State(sym) > j
    @test State(sym) >= j
    @test !(State(sym) < j)
    @test !(State(sym) == j)
    @test !(State(sym) <= j)
    @test sym != State(j)
    @test sym > State(j)
    @test sym >= State(j)
    @test !(sym == State(j))
    @test !(sym < State(j))
    @test !(sym <= State(j))
  end

  for j = (i+1):length(symbols)
    @test State(j) != sym
    @test State(j) > sym
    @test State(j) >= sym
    @test !(State(j) < sym)
    @test !(State(j) == sym)
    @test !(State(j) <= sym)
    @test sym != State(j)
    @test sym < State(j)
    @test sym <= State(j)
    @test !(sym > State(j))
    @test !(sym == State(j))
    @test !(sym >= State(j))
    @test State(sym) != j
    @test State(sym) < j
    @test State(sym) <= j
    @test !(State(sym) > j)
    @test !(State(sym) == j)
    @test !(State(sym) >= j)
    @test sym != State(j)
    @test sym < State(j)
    @test sym <= State(j)
    @test !(sym == State(j))
    @test !(sym > State(j))
    @test !(sym >= State(j))
  end
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLUtils.jl" begin

   # function lar4mesh(verts,cells) # cells are triangles
   @testset "lar4mesh" begin
      @test
      @test
      @test
      @test
   end

   # function two2three(points)
   @testset "two2three" begin
      @test
      @test
      @test
      @test
   end

   # function glGenBuffer()
   @testset "glGenBuffer" begin
      @test
      @test
      @test
      @test
   end

   # function glGenVertexArray()
   @testset "glGenVertexArray" begin
      @test
      @test
      @test
      @test
   end

   # function glLineWidth()
   @testset "glLineWidth" begin
      @test
      @test
      @test
      @test
   end

   # function glPointSize()
   @testset "glPointSize" begin
      @test
      @test
      @test
      @test
   end

   # function glCheckError(actionName="")
		  @testset "glCheckError" begin
      @test
      @test
      @test
      @test
   end

   # function glErrorMessage()
   @testset "glErrorMessage" begin
      @test
      @test
      @test
      @test
   end

   # function glDeleteLater(fun::Function)
   @testset "glDeleteLater" begin
      @test
      @test
      @test
      @test
   end

   # function glDeleteNow()
   @testset "glDeleteNow" begin
      @test
      @test
      @test
      @test
   end

   # function normalize2(V::Lar.Points; flag=true, fold=-1)
   @testset "normalize2" begin
      @test
      @test
      @test
      @test
   end

   # function normalize3(V::Lar.Points, flag=true)
   @testset "normalize3" begin
      @test
      @test
      @test
      @test
   end

end


text:
using BlackBoxOptim: epsilon_box_dominates, epsilon_box_progress

describe("epsilon_box dominance and progress") do
  describe("with scalar epsilon") do
    test("vectors of length 2") do
      @check  epsilon_box_dominates([0.5, 0.5],  [1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([0.5, 1.5],  [1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([1.5, 0.5],  [1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([1.1, 1.3],  [1.5, 1.5], 1.0)
      @check !epsilon_box_dominates([1.6, 1.5],  [1.5, 1.5], 1.0)
      @check !epsilon_box_dominates([1.5, 1.9],  [1.5, 1.5], 1.0)
      @check !epsilon_box_dominates([1.5, 1.5],  [1.5, 1.5], 1.0)

      @check  epsilon_box_progress([0.5, 0.5],  [1.5, 1.5], 1.0)
      @check  epsilon_box_progress([0.5, 1.5],  [1.5, 1.5], 1.0)
      @check  epsilon_box_progress([1.5, 0.5],  [1.5, 1.5], 1.0)
      @check !epsilon_box_progress([1.1, 1.3],  [1.5, 1.5], 1.0)
      @check !epsilon_box_progress([1.6, 1.5],  [1.5, 1.5], 1.0)
      @check !epsilon_box_progress([1.5, 1.9],  [1.5, 1.5], 1.0)
      @check !epsilon_box_progress([1.5, 1.5],  [1.5, 1.5], 1.0)

      @check epsilon_box_dominates([0.5, 0.5],  [1.0, 0.99], 1.0)
      @check epsilon_box_dominates([0.5, 0.75], [0.5, 1.0], 1.0)
      @check epsilon_box_dominates([0.5, 1.75], [2.5, 1.0], 1.0)

      @check !epsilon_box_dominates([0.5,  0.5],  [0.5, 0.5], 1.0)
      @check  epsilon_box_dominates([0.49, 0.5],  [0.5, 0.5], 1.0)
      @check  epsilon_box_dominates([0.5,  0.49], [0.5, 0.5], 1.0)

      @check epsilon_box_dominates([0.5, 0.5], [0.99, 0.99], 1.0)
    end

    test("vectors of length 3") do
      @check  epsilon_box_dominates([0.5, 0.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([1.5, 0.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([0.5, 1.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([0.5, 0.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([1.5, 1.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([1.5, 0.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_dominates([0.5, 1.5, 1.5], [1.5, 1.5, 1.5], 1.0)

      @check !epsilon_box_dominates([1.5, 1.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check !epsilon_box_dominates([1.6, 1.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check !epsilon_box_dominates([0.4, 2.7, 1.5], [1.5, 1.5, 1.5], 1.0)

      @check  epsilon_box_progress([0.5, 0.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([1.5, 0.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([0.5, 1.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([0.5, 0.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([1.5, 1.5, 0.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([1.5, 0.5, 1.5], [1.5, 1.5, 1.5], 1.0)
      @check  epsilon_box_progress([0.5, 1.5, 1.5], [1.5, 1.5, 1.5], 1.0)
    end
  end

  describe("with vector epsilon") do
    test("symmetric epsilon vector") do
      @check  epsilon_box_dominates([0.5, 0.5],  [1.5, 1.5], [1.0, 1.0])
    end

    test("asymmetric epsilon vector") do
      @check  epsilon_box_dominates([0.5, 0.5],  [1.5, 1.5], [1.0, 2.0])
      @check  epsilon_box_dominates([1.5, 0.5],  [1.5, 1.5], [1.0, 2.0])

      @check !epsilon_box_dominates([1.5, 1.5],  [1.5, 1.5], [1.0, 2.0])
      @check !epsilon_box_dominates([1.6, 1.5],  [1.5, 1.5], [1.0, 2.0])

      @check  epsilon_box_progress([0.5, 0.5],  [1.5, 1.5], [1.0, 2.0])
      @check !epsilon_box_progress([1.5, 0.5],  [1.5, 1.5], [1.0, 2.0])
      @check !epsilon_box_progress([0.5, 0.5],  [1.5, 1.5], [2.0, 2.0])
      @check !epsilon_box_progress([1.6, 1.5],  [1.5, 1.5], [1.0, 2.0])
    end
  end
end


text:
set_stepping_mode!(mode) = ()->GLOBAL_STEPPING_MODE[]=mode
const actions = OrderedDict([
   :CC => (desc="Continue",  act=set_stepping_mode!(StepContinue)),
   :SI => (desc="Step In",   act=set_stepping_mode!(StepIn)),
   :SN => (desc="Step Next", act=set_stepping_mode!(StepNext)),
   :SO => (desc="Step Out",  act=set_stepping_mode!(StepOut)),
   :XX => (desc="Abort",     act=()->throw(UserAbortedException())),
])

function print_commands()
   printstyled("Commands: "; color=:green)
   format_command((txt, cmd)) = "$txt ($(cmd.desc))"
   println(join(format_command.(collect(actions)), ", "))
end
################################

function breadcrumbs(meth, statement_ind)
   printstyled("\nBreakpoint Hit: "; color=:blue)
   printstyled(string(meth, "\n"); color=:light_blue)
   line_num = statement_ind2src_linenum(meth, statement_ind)
   breadcrumbs(string(meth.file), line_num)
   println()
end

function breadcrumbs(file::AbstractString, line_num; nbefore=2, nafter=2)
   return breadcrumbs(stdout, file, line_num; nbefore=nbefore, nafter=nafter)
end

function breadcrumbs(io, file::AbstractString, line_num; nbefore=2, nafter=2)
   @assert(nbefore >= 0)
   @assert(nafter >= 0)

   all_lines = loc_for_file(file)
   first_line_num = max(1, line_num - nbefore)
   last_line_num = min(length(all_lines), line_num + nafter)

   for ln in first_line_num:last_line_num
      line = all_lines[ln]
      if ln == line_num
         line = "➧" * line
         color = :cyan
      else
         line = " " * line
         color = :light_green
         if ln ∈ (first_line_num, last_line_num)
            color = :light_black
         end
      end
      printstyled(io, line, "\n"; color=color)
   end
end

# this function exists only for mocking so we can test it.
breakpoint_hit(meth, statement_ind, variables) = nothing

function iron_repl(meth, statement_ind, variables)
    breadcrumbs(meth, statement_ind)

    printstyled("Vars: "; color=:light_yellow)
    println(join(keys(variables), ", "))
    print_commands()

    run_repl(variables, Main)
end

"""
    should_break
Determines if we should actualy break at a potential breakpoint
"""
function should_break(meth, statement_ind)
    return GLOBAL_STEPPING_MODE[] === StepNext ||
        should_breakon(GLOBAL_BREAKPOINT_RULES, meth, statement_ind)
end

"""
    break_action

What to do when a breakpoint is hit
"""
function break_action(meth, statement_ind, slotnames, slotvals)
   variables = LittleDict(slotnames, slotvals)
   breakpoint_hit(meth, statement_ind, variables)

   code_word = iron_repl(meth, statement_ind, variables)
   actions[code_word].act()
end


text:
using Dash

app = dash()
app.layout = html_div(children = "Hello world!", id = "hello-div")
run_server(app)


text:
struct YieldCurve
    rate::Float64
end

function discount(self::YieldCurve, t::Float64)
    return exp(-self.rate*t)
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
module SpecialFuncs

using PyCall
using SymPy
using SpecialFunctions
import SpecialFunctions: besseli, besselj, besselk, bessely
import SpecialFunctions: airyai, airyaiprime, airybi, airybiprime, airyaix,airyaiprimex,
airybix, airybiprimex, besselh, besselhx, besselix, besselj0, besselj1, besseljx, besselkx,
bessely0, bessely1, besselyx, beta, dawson, erf, erfc, erfcinv, erfcx, erfi, erfinv,
eta, digamma, gamma, invdigamma, lgamma, polygamma, trigamma,
hankelh1, hankelh1x, hankelh2, hankelh2x, zeta

## https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/SpecialFunctions.jl
## julia -> sympy mapping.
julia_sympy_map = (
                   :airyai       => :airyai,
                   :airyaiprime  => :airyaiprime,
                   :airybi       => :airybi,
                   :airybiprime  => :airybiprime,
                   :airyaix      => :nothing,
                   :airyaiprimex => :nothing,
                   :airybix      => :nothing,
                   :airybiprimex => :nothing,
                   :besselh      => :nothing,
                   :besselhx     => :nothing,
#                  :besseli      => :besseli,
                   :besselix     => :nothing,
#                  :besselj      => :besselj,
                   :besselj0     => :nothing,
:besselj1   => :nothing,
:besseljx   => :nothing,
## :besselk    => :besselk,
:besselkx   => :nothing,
## :bessely    => :bessely,
:bessely0   => :nothing,
:bessely1   => :nothing,
:besselyx   => :nothing,
#:beta       => :beta,
:dawson     => :nothing,
:erf        => :erf,
:erfc       => :erfc,
:erfcinv    => :erfcinv,
:erfcx      => :erfcx,
:erfi       => :erfi,
:erfinv     => :erfinv,
:eta        => :nothing,
:digamma    => :digamma,
:gamma      => :gamma,
:invdigamma => :nothing,
:lgamma     => :lgamma,
:polygamma  => :polygamma,
:trigamma   => :trigamma,
:hankelh1   => :nothing,
:hankelh1x  => :nothing,
:hankelh2   => :nothing,
:hankelh2x  => :nothing,
:zeta       => :zeta,
#
#  SymPy only
#
:nothing    => :uppergamma,
:nothing    => :lowergamma,
:nothing    => :erf2,
:nothing    => :erf2inv,
:nothing    => :fresnels,
:nothing    => :fresnelc,
:nothing    => :Ei,
:nothing    => :expint,
:nothing    => :li,
:nothing    => :Li,
:nothing    => :Si,
:nothing    => :Ci,
:nothing    => :Shi,
:nothing    => :Chi,
:nothing    => :hankel1,
:nothing    => :hankel2,
:nothing    => :jn,
:nothing    => :yn,
:nothing    => :hn1,
:nothing    => :hn2,
:nothing    => :jn_zeros,
:nothing    => :dirichlet_eta,
:nothing    => :polylog,
:nothing    => :lerchphi,
:nothing    => :stieltjes,
:nothing    => :elliptic_k,
:nothing    => :elliptic_f,
:nothing    => :elliptic_e,
:nothing    => :elliptic_pi,
:nothing    => :mathieus,
:nothing    => :mathieuc,
:nothing    => :mathieusprime,
:nothing    => :mathieucprime,
:nothing    => :jacobi,
:nothing    => :jacobi_normalized,
:nothing    => :gegenbauer,
:nothing    => :chebyshevt,
:nothing    => :chebyshevu,
:nothing    => :chebyshevt_root,
:nothing    => :chebyshevu_root,
:nothing    => :legendre,
:nothing    => :assoc_legendre,
:nothing    => :hermite,
:nothing    => :laguerre,
:nothing    => :assoc_laguerre,
:nothing    => :Ynm,
:nothing    => :Ynm_c

)

for (jmeth, smeth) in [(j,s) for (j,s) in julia_sympy_map if s !== :nothing && j!== :nothing]
    meth_name = string(smeth)
    @eval import SpecialFunctions: $jmeth
    @eval begin
        # @doc """
        # `$($meth_name)`: a SymPy function.
        #     The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
        #     """ ->
        ($jmeth)(x::Sym, xs...;kwargs...) = SymPy._sympy_meth($meth_name, x, xs...; kwargs...)
    end
    eval(Expr(:export,smeth))
end

## only in SymPy
for meth in [v for (k,v) in julia_sympy_map if k == :nothing]
    meth_name = string(meth)
    @eval begin
        # @doc """
        #     `$($meth_name)`: a SymPy function.
        #         The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
        #         """ ->
        ($meth)(xs...; kwargs...) = SymPy._sympy_meth($meth_name, xs...; kwargs...)
    end
    eval(Expr(:export, meth))
end

#gamma(x::Sym, args...; kwargs...) = sympy_meth(:gamma, x, args...; kwargs...)
#lgamma(x::Sym, args...; kwargs...) = sympy_meth(:loggamma, x, args...; kwargs...)

beta(x::Sym, y::Number) = sympy_meth(:beta, x, y)
export beta
lgamma(x::Sym) = log(gamma(x))

## these have (parameter, x) signature. Use symbolic x to call sympy version, othewise
## should dispatch to julia version.
for fn in (:besselj, :bessely, :besseli, :besselk)
    meth = string(fn)
    @eval import SpecialFunctions: $fn
    @eval ($fn)(nu::Number, x::Sym; kwargs...) = SymPy._sympy_meth($meth, nu, x; kwargs...)
    @eval ($fn)(nu::Number, a::AbstractArray{Sym}) = map(x ->$fn(nu, x), a)
end

## Hyper and friends don't really have symbolic use...
"""

Evaluates the generalized hypergeometric function:
[hyper](http://docs.sympy.org/dev/modules/mpmath/functions/hypergeometric.html#hyper)

"""
hyper(as::Vector{T}, bs::Vector{S}, z::Number) where {T<:Number, S<:Number} = sympy_meth(:hyper, as, bs, z)
export hyper

"""

[docs](http://docs.sympy.org/dev/modules/mpmath/functions/hypergeometric.html#meijer-g-function)

"""
function meijerg(a1s::Vector{T}, a2s::Vector{T}, b1s::Vector{S}, b2s::Vector{S}, z::Number, r=1;kwargs...) where {T<:Number, S<:Number}
    sympy_meth(:meijerg, a1s, a2s, b1s, b2s, z;  kwargs...)
end
export meijerg

end


text:
"""
    ScalarFunction <: AbstractFunction

Represents a scalar function of the form:

```math
c = f(x,u)
```
where ``c \\in \\mathbb{R}``.

# Evaluation
Since the function return a scalar, both `evaluate` and `evaluate!` call the
same function methods. To avoid confusion, `evaluate` should always be preferred
when working with a `ScalarFunction`. To use, simply implement one of the following
methods:

    evaluate(fun, x, u, p)
    evaluate(fun, x, u)

where `p` is tuple of parameters.

# Differentiation
First and second-order derivatives of scalar functions are commonly referred to as
gradients and Hessians. We use the convention that a gradient is a 1-dimensional array
(i.e. an `AbstractVector` with size `(n,)`) while the Jacobian of a scalar function is a
row vector (i.e. an `AbstractMatrix` with size `(1,n)`). Theses methods can be
called using:

    gradient!(::DiffMethod, fun, grad, z)
    hessian!(::DiffMethod, fun, hess, z)

Which allows the user to dispatch on the [`DiffMethod`](@ref). These methods can
also be called by calling the more generic `jacobian` and `∇jacobian!` methods:

    jacobian!(sig, diff, fun, J, y, z)
    ∇jacobian!(sig, diff, fun, H, b, y, z)

where the length of `y`, and `b` is 1, and `b[1] == one(eltype(b))`.

To implement `UserDefined` methods, implement any one of the following gradient methods:

    gradient!(::UserDefined, fun, grad, z)
    gradient!(fun, grad, z)
    gradient!(fun, grad, x, u, p)
    gradient!(fun, grad, x, u)

and any one of the following Hessian methods:

    hessian!(::UserDefined, fun, hess, z)
    hessian!(fun, hess, z)
    hessian!(fun, hess, x, u, p)
    hessian!(fun, hess, x, u)

"""
abstract type ScalarFunction <: AbstractFunction end
output_dim(fun::ScalarFunction) = 1

# Inplace reverts to scalar return
evaluate!(fun::ScalarFunction, y, x, u, p) = evaluate(fun, x, u, p)

# Gradients
function jacobian!(sig::FunctionSignature, diff::DiffMethod, fun::ScalarFunction, J, y, z)
    @assert length(y) == 1
    gradient!(sig, diff, fun, J', z)
end

gradient!(::FunctionSignature, diff::DiffMethod, fun::ScalarFunction, grad, z) =
    gradient!(diff, fun, grad, z)
gradient!(diff::UserDefined, fun::ScalarFunction, grad, z) =
    gradient!(fun, grad, z)
gradient!(fun::AbstractFunction, grad, z::AbstractKnotPoint) =
    gradient!(fun, grad, state(z), control(z), getparams(z))
gradient!(fun::AbstractFunction, grad, x, u, p) =
    gradient!(fun, grad, x, u)
gradient!(fun::AbstractFunction, grad, x, u) =
    throw(NotImplementedError("Gradient not implemented for scalar function $(typeof(fun))"))

# Hessian
function ∇jacobian!(sig::FunctionSignature, diff::DiffMethod, fun::ScalarFunction, H, b, y, z)
    @assert b[1] ≈ one(eltype(b))
    @assert length(y) == 1
    hessian!(sig, diff, fun, H, z)
end

hessian!(::FunctionSignature, diff::DiffMethod, fun::ScalarFunction, hess, z) =
    hessian!(diff, fun, hess, z)
hessian!(diff::UserDefined, fun::ScalarFunction, hess, z) = hessian!(fun, hess, z)
hessian!(fun::ScalarFunction, hess, z::AbstractKnotPoint) =
    hessian!(fun, hess, state(z), control(z), getparams(z))
hessian!(fun::ScalarFunction, hess, x, u, p) =
    hessian!(fun, hess, x, u)
hessian!(fun::ScalarFunction, hess, x, u) =
    throw(NotImplementedError("Hessian not implemented for scalar function $(typeof(fun))"))


text:
using Winston
large = [i^4 + 1e12 for i in 1:10^3]
small = [1:10^3]
p = plot(small,large)


text:
using AbstractLogic
using Test
using Suppressor

@test LogicalCombo() |> nfeasible == 0

logicset = @suppress logicalparse("a,b,c,d,e in 1:2")
@test logicset[1,1] == 1
@test logicset[4^3,3] == 2
@test logicset[:,:a] == logicset[:,1]

logicset = @suppress logicalparse("a.1,a.2,a.3,b.1,b.2,b.3 in 1:3 || {{j}}.1 != {{j}}.2, {{j}}.3 &&& {{j}}.2 != {{j}}.3")
@test nfeasible(logicset) == 36
@test @suppress nfeasible(logicalparse("a,b,c in 1:3 || {{i}}!={{>i}}")) == nfeasible(logicalparse("a,b,c in unique"))

logicset = @suppress logicalparse("a,b,c in 1:3 || {{i}}!={{!i}}")
@test AbstractLogic.operatoreval("a^=b",logicset) |> nfeasible == 0
@test AbstractLogic.operatoreval("a|=b",logicset) |> nfeasible == 0
@test AbstractLogic.operatoreval("a!|=b",logicset) |> nfeasible == 6
@test AbstractLogic.operatoreval("a<=b",logicset) |> nfeasible == 3
@test AbstractLogic.operatoreval("a<<b",logicset) |> nfeasible == 1
@test AbstractLogic.operatoreval("a>=b",logicset) |> nfeasible == 3
@test AbstractLogic.operatoreval("a>>b",logicset) |> nfeasible == 1
@test AbstractLogic.operatoreval("a % 3 ==0",logicset) |> nfeasible == 2

@test @suppress logicalparse("{{<3|>3}}=1", logicset) |> nfeasible  == 0

@test @suppress search("{{i+4}}=1", logicset)[1] === missing

@test @suppress logicalparse("a,b ∈ 0:1; a==a === b==b; a == b !=> b != a") |> nfeasible == 4

x = @suppress logicalparse("a,b in 1:5; true &&& true ||| false ==> true <== false ||| true ^^^ false")
@test nfeasible(x) == 25

x = @suppress logicalparse("a,b in 1:5; true &&&& true |||| false ===> true <=== false |||| true ^^^^ false")
@test nfeasible(x) == 25

logicset = @suppress logicalparse("a, b in 1; a != b")
@test @suppress( search("=1", logicset)) === missing
@test  @suppress(checkfeasible("a=1", logicalparse("a, b in 1; a != b")) === missing)

@test AbstractLogic.definelogicalset("a ∈ 1") |> nfeasible == 1
@test AbstractLogic.definelogicalset("a,b,c ∈ 1,2,3 || unique") |> nfeasible == 6

@test expand(LogicalCombo(), [:a => 1:2, :b => 1:3]) |> nfeasible == 6

logicset = @suppress logicalparse("a, b, c in 1:3; a > b")
@test @suppress !dependenton("a ⊥ b", logicset)
@test @suppress dependenton("a !⊥ b", logicset)

@test @suppress !dependenton("a ⊥ b", logicset)
@test @suppress dependenton("a !⊥ b", logicset)
@test @suppress dependenton("a ⊥ c", logicset)

@test @suppress !setcompare("a ⊂ b", logicset)
@test @suppress setcompare("a !⊂ b", logicset)
@test @suppress setcompare("c ⊃ b", logicset)
@test @suppress setcompare("b !⊃ c", logicset)
@test @suppress setcompare("a,b ⋂ b,c", logicset)
@test @suppress !setcompare("a !⋂ c", logicset)

@test logicset |> showfeasible == logicset[:,:,true]

@test @suppress(logicalparse("a,b,c,d,e in 1:2; {{i}}!={{i+1}}") |> nfeasible) == 2
@test @suppress(logicalparse("a,b,c,d,e in 1:2; {{i}}!={{i-1}}") |> nfeasible) == 2
@test @suppress(logicalparse("a,b in 1:2; {{i}}!={{>i}}") |> nfeasible) == 2

@test @capture_out(help()) |> length > 500
@test @capture_out(help("")) |> length > 500
@test @capture_out(help("in")) |> length > 250
@test @capture_out(help("back")) |> length > 50
@test @capture_out(help("!=")) |> length > 50
@test @capture_out(help("!==")) |> length > 430
@test @capture_out(help("!===")) |> length > 440
@test @capture_out(help("{{i}}")) |> length > 370

@test @suppress(checkfeasible("{{i}}=1", logicalparse("a,b in 1:2")))[1] == .25
@test @capture_out(logicalparse("a,b in 1:2; error"))[(end-16):end] == "not interpreted!\n"

x = @suppress logicalparse("a,b in 1:2; a==b")
@test @suppress(AbstractLogic.back(x) |> nfeasible) == 4

@test @suppress nfeasible(logicalparse(["a,b in 1:2","a==b"])) == nfeasible(x)

@test @suppress(logicalparse("a, b ∈ 0:1; a | b ||| a & b")) |> nfeasible == 3


text:
### LINEAR MAP CONSTRUCTION

#=
If a function is already in-place, then `linear_map` will preserve this, so, e.g,

f_lm = linear_map(f,input,output)
mul!(output_vec,f_lm,input_vec)

will produce the same effect as f(output,input), but on vectors
=#

# for a given function of function-like object A, which acts upon data of type `input`
# and returns data of type `output`
# return a LinearMap that acts upon a vector form of u
linear_map(A,input,output;eltype=Float64) = _linear_map(A,input,output,eltype,Val(length(input)),Val(length(output)))

linear_map(A,input;eltype=Float64) = _linear_map(A,input,eltype,Val(length(input)))

linear_map(A::AbstractMatrix{T},input::AbstractVector{T};eltype=Float64) where {T} = LinearMap{eltype}(A)

linear_map(A::SaddleSystem{T,Ns,Nc},::Any;eltype=Float64) where {T,Ns,Nc} = LinearMap{T}(x->A*x,Ns+Nc)

linear_map(A::UniformScaling,input;eltype=Float64) = LinearMap{eltype}(x->A*x,length(input))

function linear_inverse_map(A,input;eltype=Float64)
    hasmethod(\,Tuple{typeof(A),typeof(input)}) || error("No such backslash operator exists")
    return LinearMap{eltype}(_create_vec_backslash(A,input),length(input))
end
function linear_inverse_map!(A,input;eltype=Float64)
    hasmethod(ldiv!,Tuple{typeof(input),typeof(A),typeof(input)}) || error("No such ldiv! operator exists")
    return LinearMap{eltype}(_create_vec_backslash!(A,input),length(input))
end

linear_inverse_map(A::SaddleSystem{T,Ns,Nc},::Any;eltype=Float64) where {T,Ns,Nc} = LinearMap{T}(x->A\x,Ns+Nc)

# Square operators. input of zero length
_linear_map(A,input,eltype,::Val{0}) =
      LinearMap{eltype}(x -> (),0,0)

# Square operators. input of non-zero length
_linear_map(A,input,eltype,::Val{M}) where {M} =
      LinearMap{eltype}(_create_fcn(A,input),length(input))

# input and output have zero lengths
_linear_map(A,input,output,eltype,::Val{0},::Val{0}) =
      LinearMap{eltype}(x -> (),0,0)

# input is 0 length, output is not
_linear_map(A,input,output,eltype,::Val{0},::Val{M}) where {M} =
      LinearMap{eltype}(x -> _unwrap_vec(zero(output)),length(output),0)

# output is 0 length, input is not
_linear_map(A,input,output,eltype,::Val{N},::Val{0}) where {N} =
      LinearMap{eltype}(x -> (),0,length(input))

# non-zero lengths of input and output
_linear_map(A,input,output,eltype,::Val{N},::Val{M}) where {N,M} =
      LinearMap{eltype}(_create_fcn(A,output,input),length(output),length(input))

ismultiplicative(A,input) = hasmethod(*,Tuple{typeof(A),typeof(input)})

# Create a function for operator A that can act upon an input of type AbstractVector
# and return an output of type AbstractVector. It should wrap the input vector
# in the input data type associated with A and it should then unwrap its
# output back into vector form
#=function _create_fcn(A,output,input)
    # if A has an associated * operation, then use this
    if _ismultiplicative(A,input)
        fcn = _create_vec_multiplication(A,input)
    # or if A is a function or function-like object, then use this
    elseif hasmethod(A,Tuple{typeof(input)})
        fcn = _create_vec_function(A,input)
    # or just quit
    else
        error("No function exists for this operator to act upon this type of data")
    end
    return fcn
end
=#
_create_fcn(A,input) = _create_fcn(A,nothing,input)
_create_fcn(A,output,input) = _create_fcn(A,output,input,Val(ismultiplicative(A,input)))

_create_fcn(A,output,input,::Val{true}) = _create_vec_multiplication(A,input)
_create_fcn(A,output,input,::Val{false}) = _create_fcn_function(A,output,input)
_create_fcn_function(A,output,input) = _create_fcn_function(A,output,input,Val(isinplace(A,2)))

# out of place
function _create_fcn_function(A,output,input,::Val{false})
    hasmethod(A,Tuple{typeof(input)}) || error("No function exists for this operator to act upon this type of data")
    _create_vec_function(A,input)
end

function _create_fcn_function(A,output,input,::Val{true})
    hasmethod(A,Tuple{typeof(output),typeof(input)}) || error("No function exists for this operator to act upon this type of data")
    _create_vec_function!(A,output,input)
end

# In each of these, u, outp, inp only provide the templates/sizes for the wrapping.
@inline _create_vec_multiplication(A,u::TU) where {TU} = (x -> _unwrap_vec(A*_wrap_vec(x,u)))
@inline _create_vec_function(A,u::TU) where {TU} = (x -> _unwrap_vec(A(_wrap_vec(x,u))))
@inline _create_vec_function!(A,outp::TO,inp::TI) where {TO,TI} = ((y,x) -> (_outpwrap = _wrap_vec(y,outp); A(_outpwrap,_wrap_vec(x,inp))))
@inline _create_vec_backslash(A,u::TU) where {TU} = (x -> _unwrap_vec(A\_wrap_vec(x,u)))
@inline _create_vec_backslash!(A,u::TU) where {TU} = (y,x) -> (_ywrap = _wrap_vec(y,u); ldiv!(_ywrap,A,_wrap_vec(x,u)); y)
@inline _create_vec_backslash!(A::UniformScaling,u::TU) where {TU} = (y,x) -> (_wrap_vec(y,u) .= A\_wrap_vec(x,u))
@inline _create_vec_backslash!(A::AbstractMatrix,u::TU) where {TU} = (Afact = factorize(A); (y,x) -> (_ywrap = _wrap_vec(y,u); _ywrap .= Afact\_wrap_vec(x,u)))

#### WRAPPERS ####
# wrap the vector x in type u, unless u is already a subtype of AbstractVector,
# in which case it just keeps it as is.
_wrap_vec(x::AbstractVector{T},u::TU) where {T,TU} = TU(reshape(x,size(u)...))
#_wrap_vec(x::AbstractVector{T},u::AbstractVector{U}) where {T,U} = x
_wrap_vec(x::Vector{T},u::Vector{T}) where {T} = x

# if the vector x is simply a reshaped form of type u, then just get the
# parent of x
_wrap_vec(x::Base.ReshapedArray,u::TU) where {TU} = parent(x)

#### UNWRAPPERS ####
# Usually vec suffices
_unwrap_vec(x) = vec(x)
_unwrap_vec(x::Tuple) = x


text:
using Formatting
using Base.Test

# with positional arguments

@test format("{1}", 10) == "10"
@test format("abc {1}", 10) == "abc 10"
@test format("{1} efg", 10) == "10 efg"
@test format("abc {1} efg", 10) == "abc 10 efg"
@test format("{1} + {2}", 10, "xyz") == "10 + xyz"
@test format("abc {1} + {2}", 10, "xyz") == "abc 10 + xyz"
@test format("{1} + {2} efg", 10, "xyz") == "10 + xyz efg"
@test format("abc {1} + {2} efg", 10, "xyz") == "abc 10 + xyz efg"
@test format("abc {1}{2} efg", 10, "xyz") == "abc 10xyz efg"

@test format("{1:d} + {2:s}", 10, "xyz") == "10 + xyz"
@test format("{1:04d} + {2:*>5}", 10, "xyz") == "0010 + **xyz"
@test format("let {2:<5} := {1:.4f};", 12.3, "var") == "let var   := 12.3000;"

@test format("{}", 10) == "10"
@test format("{} + {}", 10, 20) == "10 + 20"
@test format("{} + {:04d}", 10, 20) == "10 + 0020"
@test format("{:03d} + {}", 10, 20) == "010 + 20"
@test format("{:03d} + {:04d}", 10, 20) == "010 + 0020"
@test_throws(ErrorException, format("{1} + {}", 10, 20) )
@test_throws(ErrorException, format("{} + {1}", 10, 20) )

# escape {{ and }}

@test format("{{}}") == "{}"
@test format("{{{1}}}", 10) == "{10}"
@test format("v: {{{2}}} = {1:.4f}", 1.2, "ab") == "v: {ab} = 1.2000"

# with filter
@test format("{1|>abs2} + {2|>abs2:.2f}", 2, 3) == "4 + 9.00"


text:
using VisualCircuits, Compose

begin
    width, height = 12.5, 5.3

    nand3 = Gate(4.5, 1.1, :NAND)
    nand4 = Gate(4.5, 3.1, :NAND)
    nand1 = Gate(2.8, nand3[:in1][2], :NAND, [:in1, :in2, :in3])
    nand2 = Gate(2.8, nand4[:in2][2], :NAND, [:in1, :in2, :in3])

    nand7 = Gate(9, 1.1, :NAND)
    nand8 = Gate(9, 3.1, :NAND)
    nand5 = Gate(7.3, nand7[:in1][2], :NAND)
    nand6 = Gate(7.3, nand8[:in2][2], :NAND)

    not = Gate(5, 4.3, :NOT)

    box = compose(
        context(),
        rectangle(0, 0, width, height),
        stroke("cyan"), fill(nothing)
    )

    feedforward = compose(
        context(),
        connect((1, nand1[:in2][2]), nand1[:in2]),
        connect((1, nand2[:in2][2]), nand2[:in2]),
        connect(nand1[:out], nand3[:in1]),
        connect(nand3[:out], nand5[:in2]),
        connect(nand5[:out], nand7[:in1]),
        connect(nand2[:out], nand4[:in2]),
        connect(nand4[:out], nand6[:in1]),
        connect(nand6[:out], nand8[:in2]),
        connect(nand7[:out], (width-1, nand7[:out][2])),
        connect(nand8[:out], (width-1, nand8[:out][2]))
    )
    # gate low/high pos
    loop(p1, h1, h2, p2) = connect(
        p1, p1 - (0.3, 0), (p1[1] - 0.3, h1), (p2[1], h2), p2
    )
    loops = compose(
        context(),
        loop(nand3[:in2], 1.7, 2.5, nand4[:out] + (0.3, 0)),
        loop(nand4[:in1], 2.5, 1.7, nand3[:out] + (0.3, 0)),
        loop(nand7[:in2], 1.7, 2.5, nand8[:out] + (0.3, 0)),
        loop(nand8[:in1], 2.5, 1.7, nand7[:out] + (0.3, 0)),
        dots(
            nand3[:out] + (0.3, 0), nand4[:out] + (0.3, 0),
            nand7[:out] + (0.3, 0), nand8[:out] + (0.3, 0)
        )
    )

    clock = compose(
        context(),
        connect((1, 2.1), (nand1[:in1][1]-0.3, 2.1)),
        connect(
            nand1[:in3], nand1[:in3] - (0.3, 0),
            nand2[:in1] - (0.3, 0), nand2[:in1]
        ),
        connect(
            nand5[:in1], nand5[:in1] - (0.3, 0),
            nand6[:in2] - (0.3, 0), nand6[:in2]
        ),
        connect(nand2[:in1]-(0.3,0), (nand2[:in1][1]-0.3, not[:in][2]), not[:in]),
        connect(nand6[:in2]-(0.3,0), (nand6[:in2][1]-0.3, not[:out][2]), not[:out]),
        dots((nand1[:in1][1]-0.3, 2.1), nand2[:in1]-(0.3,0), nand6[:in2]-(0.3,0))
    )

    feedback = compose(
        context(),
        connect(
            nand8[:out] + (0.6, 0),
            (nand8[:out][1] + 0.6, 0.2),
            (nand1[:in1][1] - 0.3, 0.2),
            nand1[:in1] - (0.3, 0),
            nand1[:in1]
        ),
        connect(
            nand7[:out] + (0.9, 0),
            (nand7[:out][1] + 0.9, height - 0.2),
            (nand2[:in3][1] - 0.6, height - 0.2),
            nand2[:in3] - (0.6, 0),
            nand2[:in3]
        ),
        dots(nand7[:out] + (0.9, 0), nand8[:out] + (0.6, 0))
    )

    labels = compose(
        context(),
        text(0.9, nand1[:in2][2], "J",   hright, vcenter),
        text(0.9, 2.1,            "CLK", hright, vcenter),
        text(0.9, nand2[:in2][2], "K",   hright, vcenter),
        text(width-0.9, nand7[:out][2], "Q",   hleft, vcenter),
        text(width-0.9, nand8[:out][2], "~Q",   hleft, vcenter),
        fontsize(30pt), font("Helvetica-Bold"), linewidth(0mm)
    )

    circuit = compose(
        context(units = UnitBox(0, 0, width, height)),
        map(render, [eval(Symbol("nand$n")) for n in 1:8])...,
        render(not),

        # box,
        feedforward, loops, clock, feedback, labels,

        stroke("black"), linewidth(2mm)
    )

    circuit |> SVG("MS_JK_FlipFlop.svg", autosize(width, height)...)
    # circuit |> PDF("MS_JK_FlipFlop.pdf", autosize(width, height)...)
end


text:
import PerronFrobenius.RectangularInvariantMeasure

@recipe function f(riv::RectangularInvariantMeasure;
                    boxfillfactor::Int = 3, linesegments = true,
                    lw = 0.8, lc = :black, lα = 0.5, ls = :do,
                    ms = 0.2, mc = :black, mα = 0.3)

    pts = riv.pts
    ϵ = riv.ϵ
    axisminima, stepsizes = minima_and_stepsizes(pts, ϵ)

    # All bins are assigned a bin origin.
    v = unique(riv.visited_bins_coordinates, dims = 2)

    # Bins may be visited by several times, so we'll get rid
    # of the repetitions.
    n_visited_boxes = size(v, 2)
    D = size(v, 1)

    # Plot the partition grid over the points of the reconstructed
    # orbit.
    legend --> false

    for i = 1:n_visited_boxes
        origin = riv.visited_bins_coordinates[:, i]
        # Get the measure of the box, and fill the box with a
        # number of points scaled to the measure.
        μ = riv.measure.dist[i]
        fillpts = fill_hyperrectangle_3D(origin,
                    stepsizes,
                    ceil(Int, μ*100)*boxfillfactor)
        @series begin
            seriestype := :scatter3d
            markercolor --> mc
            markeralpha --> μ*5
            markersize --> μ*10
            linecolor --> mc
            fillcolor --> mc
            fillalpha --> μ*2
            splitaxes(fillpts)
        end

        if linesegments
            linesegments = connectvertices(origin, stepsizes)

            for segment in linesegments
                @series begin
                    seriestype := :path
                    linealpha --> lα
                    linewidth --> lw
                    linecolor --> lc
                    splitaxes(segment)
                end
            end
        end
    end
end


text:
using StaticArrays, ThreadingUtilities
struct MulStaticArray{P} end
function (::MulStaticArray{P})(p::Ptr{UInt}) where {P}
  _, (ptry,ptrx) = ThreadingUtilities.load(p, P, 2*sizeof(UInt))
  unsafe_store!(ptry, unsafe_load(ptrx) * 2.7)
  nothing
end
@generated function mul_staticarray_ptr(::A, ::B) where {A,B}
    c = MulStaticArray{Tuple{A,B}}()
    :(@cfunction($c, Cvoid, (Ptr{UInt},)))
end

@inline function setup_mul_svector!(p, y::Base.RefValue{T}, x::Base.RefValue{T}) where {T}
    py = Base.unsafe_convert(Ptr{T}, y)
    px = Base.unsafe_convert(Ptr{T}, x)
    fptr = mul_staticarray_ptr(py, px)
    offset = ThreadingUtilities.store!(p, fptr, sizeof(UInt))
    ThreadingUtilities.store!(p, (py,px), offset)
end

@inline function launch_thread_mul_svector(tid, y, x)
    ThreadingUtilities.launch(setup_mul_svector!, tid, y, x)
end

function waste_time(a, b)
  s = a * b'
  for i ∈ 1:00
    s += a * b'
  end
  s
end

function mul_svector_threads(a::T, b::T, c::T) where {T}
    ra = Ref(a)
    rb = Ref(b)
    rc = Ref(c)
    rx = Ref{T}()
    ry = Ref{T}()
    rz = Ref{T}()
    GC.@preserve ra rb rc rx ry rz begin
        launch_thread_mul_svector(1, rx, ra)
        launch_thread_mul_svector(2, ry, rb)
        launch_thread_mul_svector(3, rz, rc)
        w = waste_time(a, b)
        ThreadingUtilities.wait(1)
        ThreadingUtilities.wait(2)
        ThreadingUtilities.wait(3)
    end
    rx[], ry[], rz[], w
end

@testset "SVector Test" begin
  a = @SVector rand(16);
  b = @SVector rand(16);
  c = @SVector rand(16);
  w,x,y,z = mul_svector_threads(a, b, c)
  @test @allocated(mul_svector_threads(a, b, c)) == 0
  @test w == a*2.7
  @test x == b*2.7
  @test y == c*2.7
  @test z ≈ waste_time(a, b)
  A = @SMatrix rand(7,9);
  B = @SMatrix rand(7,9);
  C = @SMatrix rand(7,9);
  Wans = A*2.7; Xans = B*2.7; Yans = C*2.7; Zans = waste_time(A, B)
  for i ∈ 1:100 # repeat rapdily
    W,X,Y,Z = mul_svector_threads(A, B, C)
    @test W == Wans
    @test X == Xans
    @test Y == Yans
    @test Z ≈ Zans
  end
end


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames
using Revise
import KenyaCoV
using LinearAlgebra:eigen,normalize

"""
Load age structured data
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
N = sum(u0[:,:,1],dims = 1)

@load "data/agemixingmatrix_china.jld2" M_China
@load "data/agemixingmatrix_Kenya_norestrictions.jld2" M_Kenya
@load "data/agemixingmatrix_Kenya_homeonly.jld2" M_Kenya_ho

"""
Can adjust β to match a desired R₀ value by evaluating the leading eigenvalue
The idea is to match to the chinese epidemic R₀ -- it will be different in Kenya
"""

P.ϵ = 0.1
P.χ = ones(KenyaCoV.n_a)
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale
P.β = 2*P.γ

KenyaCoV.calculate_R₀(P)
KenyaCoV.calculate_R₀_homeonly(P)

P.dt = 0.25

"""
Can vary the spatial contact assumptions as well
"""
# P.ρ = zeros(20)
# KenyaCoV.transportstructure_params!(P,P.ρ,P_dest)

"""
Run model

"""

u0[KenyaCoV.ind_nairobi_as,5,4] = 5#10 diseased
P.ϵ_D = 1
function ramp_down(t)
    if t < 60.
        return (1-t/60) + 0.5*t/60
    else
        return 0.5
    end
end
P.c_t = ramp_down
P.c_t = t -> 1.
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,365.),P)
@time sol = solve(prob,FunctionMap(),dt = P.dt)

cum_inc = [sum(sol(t)[:,:,8]) for t = 0.:1.:365]
plot(diff(cum_inc))
plot!(diff(cum_inc).+1,yscale = :log10)

sum(sol[end][:,:,8])/sum(u0)
sum(sol[end][:,:,8])

prob_ode = KenyaCoV.create_KenyaCoV_ode_prob(u0,(0.,365.),P)
@time sol_ode = solve(prob_ode,Tsit5())

times = 0:1:365
I_area = zeros(Int64,20,length(sol.t))
for i = 1:20,(ind,t) in enumerate(sol.t)
    I_area[i,ind] = sum(sol(t)[i,:,3:4])
end
# I = [sum(sol(t)[:,:,3:4]) for t in sol.t]
plt = plot(sol.t,I_area[1,:], lab = 1,xlims = (0.,30),ylims = (0.,100));
for i = 2:20
    plot!(plt,sol.t,I_area[i,:],lab = i);
end
display(plt)

"""
Example of adding an event
"""
function isolation_limit(u,t,integrator) # Event when event_f(u,t) == 0
  integrator.p.isolating_detecteds && sum(u[:,:,8]) < 10000
end
function affect_isolation_limit!(integrator)
  integrator.p.τ = 0
  integrator.p.isolating_detecteds = false
end
cb = DiscreteCallback(isolation_limit,affect_isolation_limit!)

P.isolating_detecteds = true
P.τ = 1/3
@time sol_cb = solve(prob,FunctionMap(),dt = P.dt,callback = cb)
sum(sol_cb[end][:,:,7:8])

"""
Example of adding a constant rate jump to the method
"""

#1. Define the jump rate and affect on integrator
jumprate(u,p,t) = 1.
function jumpaffect!(integrator)
    integrator.u[4,5,3] += 1
end
#2. Declare the ConstantRateJump
example_jump = ConstantRateJump(jumprate,jumpaffect!)
#3. Make a jump problem that BOTH inherits the DiscreteProblem (prob) AND is given an "aggregator" method
#In this case the aggregator method "Direct()" is the Gillespie alg

Jmp_prob = JumpProblem(prob,Direct(),example_jump)

#Solve the whole problem using FunctionMap() this time steps forward discretely BUT if a ConstantRateJump
# event occurs in the interval it includes that as well
sol = solve(Jmp_prob,FunctionMap(),dt = P.dt)
sol.t


text:
using DynamicAxisWarping
#using DynamicAxisWarping.WarpPlots
using DynamicAxisWarping.Datasets

# UCI data repository must be downloaded first, run DynamicAxisWarping.Datasets.download_ucr()
data, labels = ucr_traindata("MedicalImages");

# c = [ cl==1 ? :blue : :red for cl in class ]'
# plot(y, linecolor=c, legend=false)
nclust = 3
init_centers = DynamicAxisWarping.dbaclust_initial_centers(data, nclust, FastDTW(10))
centers, clustids, result =
  dbaclust(data, nclust, FastDTW(10); init_centers = deepcopy(init_centers), iterations = 1)

plot(data)
plot(centers)
plot(init_centers)


text:
using Colors, GLVisualize, GLAbstraction, GLWindow

if !isdefined(:runtests)
    window = glscreen()
end
description = """
Example of how to record a video from GLVisualize
"""
kitty = visualize(loadasset("cat.obj"))
_view(kitty, window)

# save video to report dir, or in some tmp dir we'll delete later
path = if haskey(ENV, "CI_REPORT_DIR")
    ENV["CI_REPORT_DIR"] * "videorecord"
else
    homedir()
end

name = ""
isdir(path) || mkdir(path)
while true # for some reason, folder retured by mktempdir isn't usable -.-
    name = path * "/$(randstring()).mkv"
    isfile(name) || break
end

# create a stream to which we can add frames
io = GLVisualize.create_video_stream(name, window)
for i=1:10 # record 10 frames
    # do something
    GLAbstraction.set_arg!(kitty, :color, RGBA{Float32}(1, 0, 1-(i/10), i/10))
    #render current frame
    # if you call @async renderloop(window) you can replace this part with yield
    GLWindow.render_frame(window)
    GLWindow.swapbuffers(window)
    GLWindow.reactive_run_till_now()

    # add the frame from the current window
    GLVisualize.add_frame!(io)
end
# closing the stream will trigger writing the video!
close(io.io)

# clean up, only when we're not recording this!
if !haskey(ENV, "CI_REPORT_DIR")
    rm(name)
end


text:
######################################################################
# Extra general purpose functions & stubs
# -----
# Copyright (c) Kiruse 2021. Licensed under MIT License
module ExtraFun
using Reexport
export use

include("./Stubs.jl")
include("./Errors.jl")
include("./Macros.jl")

include("./Types.jl")
include("./Functions.jl")
include("./Functionals.jl")
include("./XCopy.jl")
include("./CancellableTasks.jl")
include("./Events.jl")

@reexport using .CancellableTasks

end # module

max_stars_repo_name: Kirusifix/ExtraFun

max_issues_repo_name: Kirusifix/ExtraFun

max_forks_repo_name: Kirusifix/ExtraFun