text:
using IterTools
using Distributed

N = 5
op = [:+, :-, :*, :/, :^]
fun = [:sqrt, :factorial]

function perform(all_values)
    push!(all_values, N)
    push!(all_values, 10*N + N)
    ll = length(all_values)
    size_op = length(op) * ll^2
    size_fun = ll * length(fun)

    all_op = Iterators.product(op, all_values, all_values)
    all_fun = Iterators.product(fun, all_values)

    all_exprs = Vector{Any}(undef, size_op + size_fun)

    #for (i, t) in enumerate(all_values)
    #    all_exprs[i] = all_values[i]
    #end

    for (i, t) in enumerate(all_op)
        all_exprs[i] = Expr(:call, t...)
    end

    for (i, t) in enumerate(all_fun)
        all_exprs[size_op + i] = Expr(:call, t...)
    end

    return all_exprs
end

function expands(expr::Expr)
    args = []
    for a in expr.args
        if typeof(a) == Expr
            args = [args..., expands(a)...]
        else
            push!(args, a)
        end
    end
    return args
end

count_val(val, expr::Expr) = count(i -> i == val, expands(expr))
count_N(expr::Expr) = count_val(N, expr) + 2*count_val(N*10 + N, expr)

lessThanN(expr::Expr) = count_N(expr) <= N

equal_N(expr::Expr) = count_N(expr) == N

function factorial_div_sqrt(expr::Expr)
    if expr.args[1] != :factorial
        return true
    else
        expr_vals = expands(expr)
        return !(any(f -> f == :/, expr_vals) || any(f -> f == :sqrt, expr_vals))
    end
end

function pow_negative_big(expr::Expr)
    if expr.args[1] != :^
        return true
    else
        return (eval(expr.args[2]) >= 0) && (0 <= eval(expr.args[3]) <= 16)
    end
end

function sqrt_fact_negative(expr::Expr)
    if !in(expr.args[1], [:sqrt, :factorial])
        return true
    else
        return eval(expr.args[2]) >= 0
    end
end

function big_factorials(expr::Expr)
    if expr.args[1] != :factorial
        return true
    else
        return eval(expr.args[2]) < 20
    end
end

function reduce_equivalents(expr_vals)
    dd = Dict()
    for (value, expr) in expr_vals
        if haskey(dd, value)
            push!(dd[value], expr)
        else
            dd[value] = [expr]
        end
    end

    new_expr_vals = []
    for (val, li_expr) in dd
        nb_n_done = []
        for expr in li_expr
            nb_n = count_N(expr)
            if !in(nb_n, nb_n_done)
                push!(new_expr_vals, (val, expr))
                push!(nb_n_done, nb_n)
            end
        end
    end

    return new_expr_vals
end

no_pow_mul(expr::Expr) = !in(expr.args[2], [:*, :^])
no_sqrt(expr::Expr) = !in(expr.args[1], [:sqrt])
two_add(expr::Expr) = count(i -> i == :+, expands(expr)) == 2
no_55(expr::Expr) = count(i -> i == 55, expands(expr)) == 0

function reduce_expr(all_values)
    for expr in all_values
        for (idx, sub) in enumerate(expr.args)
            if typeof(sub) == Expr
                expr.args[idx] = eval(sub)
                print( expr.args[idx])
            end
        end
    end
    return  all_values
end
# (factorial(5) - 5)
function main(number)
    all_values = []
    for i = 1:4
        all_values = perform(all_values)
        #println(all_values)

        #all_values = reduce_expr(all_values)

        all_values = Iterators.filter(lessThanN, all_values)
        all_values = Iterators.filter(factorial_div_sqrt, all_values)
        all_values = Iterators.filter(big_factorials, all_values)
        all_values = Iterators.filter(sqrt_fact_negative, all_values)
        all_values = Iterators.filter(pow_negative_big, all_values)

        ### to remove
        all_values = Iterators.filter(no_pow_mul, all_values)
        all_values = Iterators.filter(no_55, all_values)
        all_values = Iterators.filter(no_sqrt, all_values)

        if i == 4
            all_values = Iterators.filter(two_add, all_values)
        end
        ##

        expr_vals = map(expr -> (eval(expr), expr), all_values)
        expr_vals = filter(item -> -1000 < item[1] < 1000, expr_vals)
        expr_vals = reduce_equivalents(expr_vals)

        all_exprs = []
        for (value, expr) in expr_vals
            if eval(expr) == number && equal_N(expr)
                return expr
            end
            push!(all_exprs, expr)
        end
        all_values = all_exprs # map(item -> item[2], expr_vals)
        if i == 4
            println(all_values)
        else
            println(length(all_values))
        end
    end
end

#for i = 34:60
#    expr = main(i)
#    println(eval(expr),": ", expr)
#end

max_stars_repo_name: atthom/CheckFourFours

max_issues_repo_name: atthom/CheckFourFours

max_forks_repo_name: atthom/CheckFourFours


text:
workspace()

#generally helpful things
include("posit_cgen_helpers.jl")

#conversion headers
include("generate_posit.h.jl")
include("generate_posit_conv.h.jl")
include("generate_posit_ops.h.jl")

#c library source files
include("generate_posit_conv.cpp.jl")
include("generate_posit_basics.cpp.jl")
include("generate_posit_full.cpp.jl")
include("generate_posit_extended.cpp.jl")

#c++ code
include("generate_posit_class.h.jl")
include("generate_posit_class.cpp.jl")

#makefile stuff
include("generate_makefile.jl")

#sample file
include("generate_testfile.jl")

doc"""
  generate_fastsigmoid_c_libary(directory)

  generates the fastsigmoid library in a particular directory.
"""
function generate_fastsigmoid_c_library(posit_defs, d::AbstractString = normpath(Pkg.dir("FastSigmoid"),"c-src"))

  ##############################################################################
  # headers

  #create the include directory
  include_dir = normpath(d, "include")
  isdir(include_dir) || mkdir(include_dir)

  object_dir = normpath(d, "obj")
  isdir(object_dir) || mkdir(object_dir)

  #posit definitions
  posits_h = normpath(include_dir, "posit.h")
  open(posits_h, "w") do io
    generate_posit_h(io, posit_defs)
  end

  #posit conversions
  posit_conv_h = normpath(include_dir, "posit_conv.h")
  open(posit_conv_h, "w") do io
    generate_posit_conv_h(io, posit_defs)
  end

  #posit operations
  posit_ops_status_h = normpath(include_dir, "posit_ops.h")
  open(posit_ops_status_h, "w") do io
    generate_posit_ops_status_h(io, posit_defs)
  end

  posit_ops_errno_h = normpath(include_dir, "posit_ops_errno.h")
  open(posit_ops_errno_h, "w") do io
    generate_posit_ops_errno_h(io, posit_defs)
  end

  ##############################################################################
  # "c" code files.

  filelist = Dict("posit_conv"     => generate_posit_conv_cpp,
                  "posit_basics"   => generate_posit_basics_cpp,
                  "posit_full"     => generate_posit_full_cpp,
                  "posit_extended" => generate_posit_extended_cpp)

  for file in keys(filelist)
    posit_file_path = normpath(d, string(file,".cpp"))
    posit_object_path = normpath(object_dir, string(file,".o"))
    open(posit_file_path, "w") do io
      filelist[file](io, posit_defs)
    end

    #next, compile the files.
    cc = run(`gcc -c -Wall -Werror -fpic $posit_file_path -o $posit_object_path`)
  end

  ##############################################################################
  # c++ code files.

  for n in sort(collect(keys(posit_defs)))
    for es in posit_defs[n]
      posit_class_h = normpath(include_dir, "P$(n)e$(es).h")
      open(posit_class_h, "w") do io
        generate_posit_class_h(io, n, es)
      end

      posit_class_cpp = normpath(d, "P$(n)e$(es).cpp")
      open(posit_class_cpp, "w") do io
        generate_posit_class_cpp(io, n, es)
      end

      posit_class_o = normpath(object_dir, "P$(n)e$(es).o")

      #next, compile the files.
      cc = run(`gcc -c -Wall -Werror -fpic $posit_class_cpp -o $posit_class_o`)
    end
  end

  ##############################################################################
  # make makefile.

  posit_makefile = normpath(d, "Makefile")
  open(posit_makefile, "w") do io
    generate_makefile(io, posit_defs)
  end

  posit_testfile = normpath(d, "test.cpp")
  open(posit_testfile, "w") do io
    generate_testfile(io)
  end

  #link all the files and combine them into libfastposit.so
  object_files = ["$object_dir/$f" for f in readdir(object_dir) if f[end-1:end] == ".o"]
  library_path = normpath(d, "libfastposit.so")
  ln = run(`gcc -shared -o $library_path $object_files`)
end

generate_fastsigmoid_c_library(Dict(8=>[0,1,2],16=>[0,1,2],32=>[0,1,2,3]))

ext: jl


text:
@testset "model2fit/" begin

    include("model_sum.jl")
    include("const_parameters.jl")

    include("transformed_parameters.jl")

    include("peak/Peak.jl")

    include("solver_wrappers/Solver_Wrappers.jl")
end


text:
"""
NTFV fileformat
See EUROCONTROL NEST Manual for Ntfv fileformat description

Example Ntfv fileformat
#;TRAFFIC_VOLUME;2;383;20140206;20140305;7799;EAR_P
T;ABSON;ABSON FLOWS;_;UKLVCTA;AS;G;5
F;BUK>AB;IN
F;BUK>AG;IN
F;DIB>AB;IN
F;INR>AB;IN
F;INR>AG;IN
T;AERODS1;EDNY ARRS VIA EDGGADS;_;EDGGADS;AS;G;1
F;>EDNY;IN
T;AERODS2;EDNY DEPS VIA EDGGADS;_;EDGGADS;AS;G;1
F;EDNY>;IN
;EX
"""
module Ntfv

export read

struct Airblock
    name::String
    upperFL::String
end

struct Trafficvolume
    name::String
    category::String
    reflocname::String
    refloctype::String
    reflocrole::String
    airblocks::Vector{Airblock}
end

function read(file)
    trafficvolumes = Dict{String, Trafficvolume}()
    airblocks = Vector{Airblock}()
    tv_id = ""
    tv_name = ""
    category = ""
    reflocname = ""
    refloctype = ""
    reflocrole = ""

    for line in eachline(file)
        line_elements = split(line, ';')
        if line_elements[1] == "T"
            if tv_id != ""
                trafficvolumes[tv_id] = Trafficvolume(tv_name, category,
                reflocname, refloctype, reflocrole, airblocks)
                airblocks = Vector{Airblock}()
            end
                tv_id = line_elements[2]
                tv_name = line_elements[3]
                category = line_elements[4]
                reflocname = line_elements[5]
                refloctype = line_elements[6]
                reflocrole = line_elements[7]
        elseif line_elements[1] == "F"
            airblocks = vcat(airblocks, Airblock(line_elements[2],
            line_elements[3]))
        end
    end
    trafficvolumes[tv_id] = Trafficvolume(tv_name, category,
    reflocname, refloctype, reflocrole, airblocks)
    return trafficvolumes
end
#
# function splitintodict(iterator)
#     dict = Dict{String, Trafficvolume}()
#     airblocks = Vector{Airblock}(undef, 0)
#     j = 1
#     n_flowelements = j+2
#     trafficvolumeid = ""
#     trafficvolumename = ""
#     trafficvolumecategory = ""
#     reflocname = ""
#     refloctype = ""
#     reflocrole = ""
#     for (i, row) in enumerate(iterator)
#         if row.TorF == "F"
#             airblocks = vcat(airblocks, Airblock(row.ID, row.NAME))
#             if n_flowelements == j
#                 dict[trafficvolumeid] = Trafficvolume(trafficvolumeid,
#                 trafficvolumename, trafficvolumecategory, reflocname,
#                 refloctype, reflocrole, n_flowelements, airblocks)
#             end
#             j += 1
#         elseif row.TorF == "T"
#             trafficvolumeid = row.ID === missing ? "" : row.ID
#             trafficvolumename = row.NAME === missing ? "" : row.NAME
#             trafficvolumecategory = row.CATEGORY === missing ? "" : row.CATEGORY
#             reflocname = row.REFLOCNAME === missing ? "" : row.REFLOCNAME
#             refloctype = row.REFLOCTYPE === missing ? "" : row.REFLOCTYPE
#             reflocrole = row.REFLOCROLE === missing ? "" : row.REFLOCROLE
#             n_flowelements = parse(Int64, row.N_FLOWELEMENTS)
#             airblocks = Vector{Airblock}(undef, 0)
#             j = 1
#         else
#         end
#     end
#     return dict
# end

end # module"


text:
using StingerGraphs
using BenchmarkTools
using JLD

function setupgraph(
    scale::Int64,
    edgefactor::Int64;
    a::Float64=0.57,
    b::Float64=0.19,
    c::Float64 = 0.19
    )
    graph = kronecker(scale, edgefactor, a=a, b=b, c=c)
    s = Stinger()
    for i in 1:size(graph, 2)
        insert_edge!(s, 0, graph[1, i], graph[2, i], 1, 1)
    end
    return s
end

function bfsbenchutil(s::Stinger, nv::Int64)
    for i in 0:1000
        bfs(s, i, nv)
    end
end

function bench(
    scale::Int64,
    edgefactor::Int64;
    a::Float64=0.57,
    b::Float64=0.19,
    c::Float64 = 0.19,
    filename::String="bfs_bench.jld"
    )
    nv = 2^scale
    bfs_bench = @benchmarkable bfsbenchutil(s, $nv) seconds=10 setup=(s=setupgraph($scale, $edgefactor))
    info("Running BFS benchmark")
    bfs_trial = run(bfs_bench)
    jldopen(filename, "w") do f
        write(f, "bfs_trial_$(scale)_$(edgefactor)", bfs_trial)
    end
    bfs_trial
end

function benchgroup(
        scales::Range{Int64},
        edgefactor::Int64;
        a::Float64=0.57,
        b::Float64=0.19,
        c::Float64 = 0.19,
        filename::String="bfs_bench_group.jld"
    )
    jldopen(filename, "w") do f
        for scale in scales
            bfs_trial = bench(scale, edgefactor, a=a, b=b, c=c)
            write(f, "bfs_trial_$(scale)_$(edgefactor)", bfs_trial)
        end
    end
    filename
end


text:
module ClimateCalculator

using Parameters, Optim

export runscen, solveclimate, getscenario, calibrateforcing!, calibratefertilization!, printresults, iyear,
		ClimateState, ClimateParams, annualresults, makecalibrationcache, startserver, importGlobalCarbonProject

const YEARS				= 1765:2500
const GAS				= [:CO2, :CH4, :N2O, :H2O, :O3]
const GAS3 				= [:CO2, :CH4, :N2O]

const y0 = YEARS[1]

iyear(y::Int) = y - y0 + 1
iyear(t::Float64, f::Function) = f(Int,t) - y0 + 1
iyear(t::Float64) = iyear(t,round)

interpolate(t, annualvector) =  (1-(t-floor(t))) * annualvector[iyear(t,floor)] +
									(t-floor(t)) * annualvector[iyear(t,ceil)]

mutable struct ClimateParams
	timestep				::Float64
	lambda					::Float64
	aerosolforcingfactor	::Float64
	fertilization			::Float64
	oceantempfeedback		::Float64
	bioQ10factor			::Float64
	equilibriumCO2			::Float64
end

mutable struct ClimateState
	Concentration			::Dict{Symbol, Float64}
	RadiativeForcing		::Dict{Symbol, Float64}
	TotalRadiativeForcing	::Float64
	DeltaCO2ocean			::Float64
	NetFluxOcean			::Float64
	DeltaDIC				::Float64
	DeltaDICbox				::Vector{Float64}
	NetFluxBiosphere		::Float64
	BioReservoir			::Vector{Float64}
	Temp_global				::Float64
	Temp_land				::Float64
	Temp					::Vector{Float64}
end

include("oceancarbon.jl")
include("biocarbon.jl")
include("concentrations.jl")
include("radiativeforcing.jl")
include("temperatures.jl")

ClimateState(v::Vector{Float64}) = ClimateState(Dict(GAS3[i] => v[i] for i=1:3), Dict(GAS[i] => v[i+3] for i=1:5),
		v[9], v[10], v[11], v[12], v[13:18], v[19], v[20:23], v[24], v[25], v[26:75])

ClimateState() = ClimateState(zeros(75))::ClimateState

state2vector(s::ClimateState)::Array{Float64} = [
		[s.Concentration[g] for g in GAS3]; [s.RadiativeForcing[g] for g in GAS]; s.TotalRadiativeForcing;
		s.DeltaCO2ocean; s.NetFluxOcean; s.DeltaDIC; s.DeltaDICbox; s.NetFluxBiosphere; s.BioReservoir;
		s.Temp_global; s.Temp_land; s.Temp]

Base.copy(s::ClimateState) = ClimateState(copy(s.Concentration), copy(s.RadiativeForcing), s.TotalRadiativeForcing,
											s.DeltaCO2ocean, s.NetFluxOcean, s.DeltaDIC,
											copy(s.DeltaDICbox), s.NetFluxBiosphere, copy(s.BioReservoir),
											s.Temp_global, s.Temp_land, copy(s.Temp))

include("scenario.jl")
include("climate.jl")
include("historicdata.jl")
include("calibrate.jl")
include("cachehistory.jl")
include("webserver.jl")
include("importGlobalCarbonProject.jl")

end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text: remap(x, from, to) = (x - from[1]) * (to[2] - to[1]) / (from[2] - from[1]) + to[1]


text:
module Musical
export TIME_TYPE, AMP_TYPE, totimetype, toamptype,
    MCAudio, addto!, appendzeros!,
    nchannels, nsamples,
    ismono, isstereo,
    left, right,
    centerrange, leftrange, rightrange,
    seconds2saferange, unsafe_seconds2range,
    mono, stereo,
    db, dB,
    importsound, exportsound,
    plotall,
    Envelope, ADSREnvelope, ENVELOPE_BOX, ENVELOPE_STANDARD, apply!,
    sinewave, squarewave, sawtooth, trianglewave, whitenoise, binarynoise, chirpf, ichirpf,
    WaveForm, waveform,
    Note, Notes, scale!, render,
    Tuning, TUNING_STANDARD, notefreq, bar

using WAV, DSP, PyPlot

typealias TIME_TYPE Float64
typealias AMP_TYPE Float32
totimetype(x) = convert(TIME_TYPE, x)
totimetype(x::Array) = convert(Array{TIME_TYPE}, x)
toamptype(x) = convert(AMP_TYPE, x)
toamptype(x::Array) = convert(Array{AMP_TYPE}, x)

# multi channel sound
type MCAudio
    v::Array{AMP_TYPE}
    fs::Int
    function MCAudio(v,fs)
        new(v,fs)
    end
end
function Base.getindex(s::MCAudio, op::Symbol)
    if op == :left
        @assert isstereo(s)
        return s.v[leftrange(s)]
    elseif op == :right
        @assert isstereo(s)
        return s.v[rightrange(s)]
    elseif op == :center
        @assert ismono(s)
        return s.v[centerrange(s)]
    else
        throw_unknown_symbol(op)
    end
end
function Base.getindex(s::MCAudio, op::Symbol, a::Real, b::Real)
    if op == :sec
        i1 = int(a*s.fs) + 1
        i1 = max(1, i1)
        i2 = int(b*s.fs)
        i2 = min(nsamples(s), i2)
        return s.v[i1:i2,:]
    else
        throw_unknown_symbol(op)
    end
end

function Base.setindex!(s::MCAudio, val::Real, op::Symbol)
    if op == :left
        @assert isstereo(s)
        range = leftrange(s)
    elseif op == :right
        @assert isstereo(s)
        range = rightrange(s)
    elseif op == :center
        @assert ismono(s)
        range = centerrange(s)
    else
        throw_unknown_symbol(op)
    end
    @inbounds for i in range
        s.v[i] = val
    end
    return s
end
function Base.setindex!(s::MCAudio, x::AbstractArray, op::Symbol)
    if nchannels(x) != 1 || nsamples(s) != nsamples(x)
        Base.throw_setindex_mismatch(x, (op,))
    end
    if op == :left
        @assert isstereo(s)
        si = leftrange(s)[1]
    elseif op == :right
        @assert isstereo(s)
        si = rightrange(s)[1]
    elseif op == :center
        @assert ismono(s)
        si = centerrange(s)[1]
    else
        throw_unknown_symbol(op)
    end
    copy!(s.v, si, x, 1, nsamples(x))
    return s
end

function Base.setindex!(s::MCAudio, val::Real, op::Symbol, a::Real, b::Real)
    if op == :sec
        @inbounds for j in 1:nchannels(s)
            for i in seconds2saferange(s, a, b)
                s.v[i,j] = val
            end
        end
        return s
    else
        throw_unknown_symbol(op)
    end
end
function Base.setindex!(s::MCAudio, x::AbstractArray, op::Symbol, a::Real, b::Real)
    if nchannels(s) != nchannels(x)
        Base.throw_setindex_mismatch(x, (op,))
    end
    if op == :sec
        r = seconds2saferange(s, a, b)
        if length(r) != nsamples(x)
            Base.throw_setindex_mismatch(x, (op,))
        end
        for j in 1:nchannels(s)
            copy!(s.v, r[1] + (j-1)*nsamples(s), x, 1 + (j-1)*nsamples(x), nsamples(x))
        end
        return s
    else
        throw_unknown_symbol(op)
    end
end

throw_unknown_symbol(x::Symbol) = error("uknown symbol :", string(x))

nchannels(s::MCAudio) = nchannels(s.v)
nchannels(v::AbstractArray) = size(v, 2)
nsamples(s::MCAudio) = nsamples(s.v)
nsamples(v::AbstractArray) = size(v, 1)
timelength(s::MCAudio) = totimetype(nsamples(s)/s.fs)
# mono
centerrange(s::MCAudio) = 1:nsamples(s)
ismono(s::MCAudio) = (nchannels(s) == 1)
# stereo
leftrange(s::MCAudio) = 1:nsamples(s)
rightrange(s::MCAudio) = nsamples(s)+1:2*nsamples(s)
isstereo(s::MCAudio) = (nchannels(s) == 2)

seconds2saferange(s::MCAudio, a::Real, b::Real) = seconds2saferange(s.fs,nsamples(s),a,b)
function seconds2saferange(fs::Integer, nsamples::Integer, a::Real, b::Real)
    i1 = int(totimetype(a)*fs) + 1
    i1 = max(1, i1)
    i2 = int(totimetype(b)*fs)
    i2 = min(nsamples, i2)
    return i1:i2
end
unsafe_seconds2range(s::MCAudio, a::Real, b::Real) = unsafe_seconds2range(s.fs,nsamples(s),a,b)
function unsafe_seconds2range(fs::Integer, nsamples::Integer, a::Real, b::Real)
    i1 = int(totimetype(a)*fs) + 1
    i2 = int(totimetype(b)*fs)
    return i1:i2
end

#Base.convert(::Type(MCAudio), x::Array) = MCAudio(x,fs)
function left(s::MCAudio)
    @assert isstereo(s)
    return MCAudio(s[:left], s.fs)
end
function right(s::MCAudio)
    @assert isstereo(s)
    return MCAudio(s[:right], s.fs)
end
function mono(s::MCAudio)
    @assert nchannels(s) <= 2
    if ismono(s)
        return deepcopy(s)
    else
        return MCAudio((s[:left] + s[:right])/2, s.fs)
    end
end
function stereo(s::MCAudio)
    @assert nchannels(s) <= 2
    if ismono(s)
        return MCAudio(repmat(s.v,1,2), s.fs)
    else
        return deepcopy(s)
    end
end
function stereo(s1::MCAudio, s2::MCAudio)
    @assert ismono(s1) && ismono(s2)
    @assert s1.fs == s2.fs
    return MCAudio([s1.v s2.v], s1.fs)
end

# + - * / ^
+(s::MCAudio, a::Union(Real, AbstractArray)) = MCAudio(s.v .+ a, s.fs)
+(a::Union(Real, AbstractArray), s::MCAudio) = MCAudio(a .+ s.v, s.fs)
-(s::MCAudio, a::Union(Real, AbstractArray)) = MCAudio(s.v .- a, s.fs)
-(a::Union(Real, AbstractArray), s::MCAudio) = MCAudio(a .- s.v, s.fs)
*(s::MCAudio, a::Union(Real, AbstractArray)) = MCAudio(s.v .* a, s.fs)
*(a::Union(Real, AbstractArray), s::MCAudio) = MCAudio(a .* s.v, s.fs)
/(s::MCAudio, a::Union(Real, AbstractArray)) = MCAudio(s.v ./ a, s.fs)
/(a::Union(Real, AbstractArray), s::MCAudio) = MCAudio(a ./ s.v, s.fs)
#(^)(s::MCAudio, a::Real) = MCAudio(s.v .^ a, s.fs) #defined in Base
(^)(a::Union(Real, AbstractArray), s::MCAudio) = MCAudio(a .^ s.v, s.fs)
+(s::MCAudio, d::MCAudio) = MCAudio(s.v .+ d.v, s.fs)
-(s::MCAudio, d::MCAudio) = MCAudio(s.v .- d.v, s.fs)
*(s::MCAudio, d::MCAudio) = MCAudio(s.v .* d.v, s.fs)
/(s::MCAudio, d::MCAudio) = MCAudio(s.v ./ d.v, s.fs)
(^)(s::MCAudio, d::MCAudio) = MCAudio(s.v .^ d.v, s.fs)

function addto!(s::MCAudio, time::Real, a::MCAudio)
    @assert nchannels(s) == nchannels(a)
    @assert time >= 0

    starti = unsafe_seconds2range(s, time, time).start
    ir = starti:(starti + nsamples(a) - 1)
    if ir[end] > nsamples(s)  #append s with zeros
        appendzeros!(s, ir[end])
    end
    addarrays!(s.v, ir, a.v)
    return s
end
function appendzeros!(s::MCAudio, n::Integer)
    @assert n > nsamples(s)
    sv = s.v
    s.v = zeros(AMP_TYPE, n, nchannels(s))
    for j in 1:nchannels(s)
        copy!(s.v, 1 + (j-1)*nsamples(s.v), sv, 1 + (j-1)*nsamples(sv), nsamples(sv))
    end
    return s
end
function addarrays!(a::Array, ir::Range, b::Array)
    @assert nchannels(a) == nchannels(b)
    @assert length(ir) == nsamples(b)
    @assert ir[end] <= nsamples(a)
    @inbounds for c in 1:nchannels(a)
        for i = 1:length(ir)
            a[ir[i],c] += b[i,c]
        end
    end
    return a
end

# utils

type dBconvert
end
const dB = dBconvert()
*(a::Real, d::dBconvert) = db(a)
db(a::Real) = toamptype(10^(a*0.05))

function importsound(file::String; subrange=Any)
    x = wavread(file; subrange=subrange, format="double")
    return MCAudio(toamptype(x[1]), convert(Int,x[2]))
end
function exportsound(s::MCAudio, file::String)
    wavwrite(s.v, file; Fs=s.fs, nbits=32, compression=WAVE_FORMAT_IEEE_FLOAT)
end

include("plot.jl")

include("synth.jl")

end #module


text:
"""
    outneighbors_looped(graph, vertex)

Return lists of outneighbors of `vertex` in `graph`. If outdegree of `vertex` is
zero, then [v] is returned.
"""
function outneighbors_looped(graph::AbstractGraph, vertex::Int)
   outns = outneighbors(graph, vertex)
   if outns == []
      return [vertex]
   end
   outns
end

"""
    isstochastic(stochastic)

Checks whether `stochastic` is column-stochastic matrix. Returns nothing.
"""
function stochasticcheck(stochastic::SparseMatrixCSC{<:Number})
   @assert size(stochastic, 1) == size(stochastic, 2) "Matrix is not square"
   @assert all(sum(stochastic[:,i]) ≈ 1 for i=1:size(stochastic, 1)) "Columns sums do not equal one"
   @assert all(findnz(stochastic)[3] .> 0 ) "Matrix elements are not nonnegative"
   nothing
end

"""
    stochastic_preserves_graph_check(graph, stochastic)

Checks whether `stochastic` matrix presreves graph structure. Returns nothing.
"""
function graphpreservationcheck(graph::AbstractGraph,
                                stochastic::SparseMatrixCSC{<:Number})
   @assert size(stochastic, 1) == nv(graph) "Orders of matrix and graph do not match"
   @assert all([ i ∈ outneighbors_looped(graph, j) for (i, j, _)=zip(findnz(stochastic)...)]) "Nonzero elements of stochastic do not coincide with graph edges"
   nothing
end

"""
    graphstochasticcheck(graph, stochastic)

Check, whether `stochastic` is a column-stochatsic matrix preserving `graph`
structure. Three properties are verified:
* `stochastic` is column-stochastic matrix,
* nonzero element of `stochastic` are located on edges of `graph`.

Function returns `nothing` if `stochastic` satisfies properties above, otherwise
raise error.
"""
function graphstochasticcheck(graph::AbstractGraph,
                              stochastic::SparseMatrixCSC{<:Number})
   stochasticcheck(stochastic)
   graphpreservationcheck(graph, stochastic)
   nothing
end

"""
    default_stochastic(graph)

Generates default column-stochastic matrix, which represents random walk with
uniform spreading. If outdegree of vertex is zero, additional 1 is added on the
diagonal.
"""
function default_stochastic(graph::AbstractGraph)
   result = adjacency_matrix(graph, Float64, dir=:in)
   outdegrees = outdegree(graph)
   for i = 1:size(result, 1)
      if outdegrees[i] == 0
         result[i,i] = 1.
      else
         result[:,i] /= outdegrees[i]
      end
   end
   result
end


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

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

# ╔═╡ 9a0cec14-08db-11eb-3cfa-4d1c327c63f1
begin
    import Pkg

    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="StatsBase", version="0.33"),
    ])

    using Plots, PlutoUI, StatsBase, Statistics
end

# ╔═╡ fb6cdc08-8b44-11eb-09f5-43c167aa53fd
PlutoUI.TableOfContents(aside=true)

# ╔═╡ f3aad4f0-8cc2-11eb-1a25-535297327c65
md"""
# Modeling with stochastic simulation
"""

# ╔═╡ b6b055b6-8cae-11eb-29e5-b507c1a2b9bf
md"""
## Julia features
"""

# ╔═╡ bcfaedfa-8cae-11eb-10a1-cb7be7dc2e6b
md"""
- Extending a function from another package (not `Base`)

- Why create a new type? Abstractions and concentrating information around objects (organisation!)

- Plotting shapes

- String interpolation

- Conventions on capitalization of function, variable and type names
"""

# ╔═╡ d2c19564-8b44-11eb-1077-ddf6d1395b59
md"""
## Individual-based ("microscopic") models
"""

# ╔═╡ ac4f6944-08d9-11eb-0b3c-f5e0a8b8c17e
md"""

In **individual-based** models, we literally specify the behaviour, or actions, of each individual in a set of rules. But often we are interested in a global picture of how the whole system, consisting of many individuals, evolves in time: how many infectious individuals there are in total at a given time, for example, or the behaviour of the whole stock market.

In this notebook we will see how we can start from an individual-based probabilistic  (probabilistic) model, and how sometimes we can find **deterministic** equations for a **macroscopic** (system-level) description.

Those macroscopic equations can either be in discrete time (**recurrence relations** or **difference equations**), or we can take a **continuous limit** and convert them into **ordinary differential equations**.

"""

# ╔═╡ 4ca399f4-8b45-11eb-2d2b-8189e04fc804
md"""
## Modelling time to success (or time to failure)
"""

# ╔═╡ 57080632-8b45-11eb-1003-05afb2331b25
md"""
Let's start with a very simple model of **time to success**. Suppose we are playing a game in which we have a probability $p$ of success on each turn. How many turns do we need until we succeed? For example, how many rolls of a die do we need until we roll a 6? How many rolls of 2 dice until we get a double 6?
"""

# ╔═╡ 139ecfec-8b46-11eb-2649-2f77833d749a
md"""
This can be used as a model in many different situations, for example: Time to failure of a light bulb or a piece of machinery; time for a radioactive nucleus to decay; time to recover from an infection, etc.
"""

# ╔═╡ f9a75ac4-08d9-11eb-3167-011eb698a32c
md"""

The basic question is:

> Suppose we have $N$ light bulbs that are working correctly on day $0$.
>
> - If each bulb has probability $p$ to fail on each day, how many are still working at day number $t$?
> - How long, on average, will a given bulb last?
> - And, do light bulbs really fail exactly at midnight each night? Can you imagine a more realistic model?

As usual we will take a computational thinking point of view: let's code up the simulation and plot the results. Then we can step back and
"""

# ╔═╡ 17812c7c-8cac-11eb-1d0a-6512415f6938
md"""
## Visualizing component failure
"""

# ╔═╡ 178631ec-8cac-11eb-1117-5d872ba7f66e
function simulate(N, p)
	v = fill(0, N, N)
	t = 0

	while any( v .== 0 ) && t < 100
		t += 1

		for i= 1:N, j=1:N
			if rand() < p && v[i,j]==0
				v[i,j] = t
			end
		end

	end

	return v
end

# ╔═╡ 179a4db2-8cac-11eb-374f-0f24dc81ebeb
md"""
M= $(@bind M Slider(2:20, show_value=true, default=8))
"""

# ╔═╡ 17cf895a-8cac-11eb-017e-c79ffcab60b1
md"""
p= $(@bind prob Slider(0.01:.01:1, show_value=true, default=.1))

t = $(@bind tt Slider(1:100, show_value=true, default=1))
"""

# ╔═╡ 17bbf532-8cac-11eb-1e3f-c54072021208
simulation = simulate(M, prob)

# ╔═╡ a38fe2b2-8cae-11eb-19e8-d563e82855d3
gr()

# ╔═╡ 18da7920-8cac-11eb-07f4-e109298fd5f1
begin
	rectangle(w, h, x, y) = (x .+ [0, w, w, 0], y .+ [0, 0, h, h])

	circle(r, x, y) = (θ = range(0, 2π, length=30); (x .+ r .* cos.(θ), y .+ r .* sin.(θ)))
end

# ╔═╡ 17e0d142-8cac-11eb-2d6a-fdf175f5d419
begin
	w = .9
	h = .9
	c = [RGB(0,1,0), RGB(1,0,0), :purple][1 .+ (simulation .< tt) .+ (simulation .<  (tt.-1))]

	plot(ratio=1, legend=false, axis=false, ticks=false)

	for i=1:M, j=1:M
		plot!(rectangle(w, h, i, j), c=:black, fill=true, alpha=0.5)
		plot!(circle(0.3, i+0.45, j+0.45), c = c[i, j], fill=true, alpha=0.5)
	end

	for i=1:M, j=1:M
		if simulation[i,j] < tt
	       annotate!(i+0.45, j+0.5, text("$(simulation[i, j])", 7, :white))
		end
	end

	plot!(lims=(0.5, M+1.1), title="time = $(tt-1);  failed count: $(sum(simulation.<tt))")

end

# ╔═╡ 17fe87a0-8cac-11eb-2938-2d9cd19ecc0f
begin

	plot(size=(500, 300))
	cdf= [ count(simulation .≤ i) for i=0:100]
	bar!(cdf, c=:purple, legend=false, xlim=(0,tt),alpha=0.8)
end

# ╔═╡ 1829091c-8cac-11eb-1b77-c5ed7dd1261b
begin
	newcdf = [ count(simulation .> i) for i=0:100]
	bar!( newcdf, c=RGB(0,1,0), legend=false, xlim=(0,tt),alpha=0.8)
end

# ╔═╡ 1851dd6a-8cac-11eb-18e4-87dbe1714be0
bar(countmap(simulation[:]), c=:red, legend=false, xlim=(0, tt+.5), size=(500, 300))

# ╔═╡ a9447530-8cb6-11eb-38f7-ff69a640e3c4
md"""
## String interpolation
"""

# ╔═╡ c1fde6ba-8cb6-11eb-2170-af6bc84c01a7
md"""
As an aside, how could we display a picture of Daniel Bernoulli and resize it in Pluto? To do so we use a piece of HTML, which we represent as a string. We need to substitute the *value* of a Julia variable into the string, which we do with string interpolation, with the syntax `$(variable)` inside the string.

Then we convert the string to HTML with the `HTML(...)` constructor:
"""

# ╔═╡ 18755e3e-8cac-11eb-37bf-1dfa5fbe730a
@bind bernoulliwidth Slider(10:10:500, show_value=true)

# ╔═╡ f947a976-8cb6-11eb-2ae7-59eba4c6f40f
url = "https://upload.wikimedia.org/wikipedia/commons/thumb/b/b3/ETH-BIB-Bernoulli%2C_Daniel_%281700-1782%29-Portrait-Portr_10971.tif_%28cropped%29.jpg/440px-ETH-BIB-Bernoulli%2C_Daniel_%281700-1782%29-Portrait-Portr_10971.tif_%28cropped%29.jpg"

# ╔═╡ 5a0b407e-8cb7-11eb-0c0d-c7767a6b0a1d
s = ""

# ╔═╡ fe53ee0c-8cb6-11eb-19bc-2976da1abe16
md"""
Note that we can use *three* sets of double quotes (`"`) to represent a **multi-line** string, or to enclose another string that itself contains quotes.
"""

# ╔═╡ 1894b388-8cac-11eb-2287-97f985df1fbd
HTML(s)

# ╔═╡ 86299112-8cc6-11eb-257a-9d803feac359
html"3
 4"

# ╔═╡ 9eb69e94-8cc6-11eb-323b-5587cc743571
HTML("3
 4")

# ╔═╡ b8d01ae4-8cc6-11eb-1d16-47095532f738
@bind num_breaks Slider(1:10, show_value=true)

# ╔═╡ c5bc97aa-8cc6-11eb-3564-fd96184ff2c3
repeat("
", 5)

# ╔═╡ d5cffd96-8cc6-11eb-2714-975d46d4fa27
"
"^5

# ╔═╡ da6525f0-8cc6-11eb-31ec-cd7d515f054a
"
" * "
" * "
"

# ╔═╡ fc715452-8cc6-11eb-0246-e941f7698cfe
HTML("3 $("
"^num_breaks) 4")

# ╔═╡ 71fbc75e-8bf3-11eb-3ac9-dd5401033c78
md"""
## Math: Bernoulli random variables
"""

# ╔═╡ 7c01f6a6-8bf3-11eb-3c4d-ad7e206a9277
md"""
Recall that a **Bernoulli random variable** models a weighted coin: it takes the value 1, with probability $p$, and 0, with probability $q = (1 - p)$:
"""

# ╔═╡ ba7ffe78-0845-11eb-2847-851a407dd2ec
bernoulli(p) = rand() < p

# ╔═╡ 45bb6df0-8cc7-11eb-100f-37c2a76df464
bernoulli(0.25)

# ╔═╡ dcd279b0-8bf3-11eb-0cb9-95f351626ed1
md"""
Note that `rand() < p` returns a `Bool` (true or false). We are converting to `Int` to get a value 1 or 0.
"""

# ╔═╡ fbada990-8bf3-11eb-2bb7-d786362669e8
md"""
Let's generate (sample) some Bernoulli random variates:
"""

# ╔═╡ b6786ec8-8bf3-11eb-1347-61f231fd3b4c
flips = [Int(bernoulli(0.25)) for i in 1:100]

# ╔═╡ ac98f5da-8bf3-11eb-076f-597ce4455e76
md"""
It is natural to ask what the **mean**, or **expected value**, is:
"""

# ╔═╡ 0e7a04a4-8bf4-11eb-2e9d-fb48c23b8d8c
sample_mean(data) = sum(data) / length(data)

# ╔═╡ 093275e4-8cc8-11eb-136f-3ffe522c4125
sample_mean(flips)

# ╔═╡ 111eccd2-8bf4-11eb-097c-7582f811d146
md"""
If you think about how to calculate the mean (sum up and divide by the total number of flips), it just gives the proportion of 1s, which should be around $p$.

#### Exercise:
Calculate the variance of a Bernoulli random variable.

Hint: What happens when you sum the squares? Remember that you also need to center the data.
"""

# ╔═╡ 4edaec4a-8bf4-11eb-3094-010ebe9b56ab
md"""
## Julia: Make it a type!
"""

# ╔═╡ 9d66e31e-8cad-11eb-3ad0-3980ba66cb0e
md"""

Currently we need one function for sampling from a Bernoulli random variable, a different function to calculate its mean, a different function for its standard deviation, etc.

From a mathematical point of view we have the concept "Bernoulli random variable" and we are calculating properties of that concept. Computationally we can *do the same thing!* by creating a new *object* to represent "a Bernoulli random variable".
"""

# ╔═╡ 8405e310-8bf8-11eb-282b-d93b4fc683aa
struct Bernoulli   # weighted coin flip
	p::Float64
end

# ╔═╡ af2594c4-8cad-11eb-0fff-f59e65102b3f
md"""
We want to be able to sample from it, using `rand`, and take its `mean`.
To do so we will **extend** (sometimes called "overload") the `rand` function from Julia's `Base` library, and the `mean` function from the `Statistics` standard library. Note that we are *adding methods* to these functions; you will do this in the homework.
"""

# ╔═╡ 8aa60da0-8bf8-11eb-0fa2-11aeecb89564
Base.rand(X::Bernoulli) = Int( rand() < X.p )

# ╔═╡ 8d56b6d2-8cc8-11eb-0dbd-b3533dde4aa3
md"""
[Recall: In Julia the convention is that functions and variable names start with lower-case letters; types start with upper case. Maybe with an exception for 1-letter variable names for mathematical objects.]
"""

# ╔═╡ a034c2a6-8bf8-11eb-0f06-0b35a0e8e68d
B = Bernoulli(0.25)

# ╔═╡ c166edb6-8cc8-11eb-0436-f74164ff6ea7
methods(Bernoulli)

# ╔═╡ d25dc130-8cc8-11eb-177f-63a1792494c0
Bernoulli(1//4)

# ╔═╡ e2f45234-8cc8-11eb-2be9-598eb590592a
rand(B)

# ╔═╡ 3ef23da4-8cb4-11eb-0d5f-d5ee8fc56227
md"""
The object `B` really represents "a Bernoulli random variable with probability of success $p$". Since all such random variables are the same, this represents *any* Bernoulli random variable with that probability.

We should use this type any time we need a Bernoulli random variable. If you need this in another notebook you will either need to copy and paste the definition or, better, make your own mini-library. However, note that types like this are already available in the `Distributions.jl` package and the new `MeasureTheory.jl` package.
"""

# ╔═╡ bc5d6fae-8cad-11eb-3351-a734d2366557
rand(B)

# ╔═╡ 2d9c560e-8bf9-11eb-1ac5-f77f7caf776f
Statistics.mean(X::Bernoulli) = X.p

# ╔═╡ ce94541c-8bf9-11eb-1ac9-51e66a017813
mean(B)

# ╔═╡ a057e7ee-8bf9-11eb-2ceb-2dda3718a70a
md"""
## Running the stochastic simulation
"""

# ╔═╡ 4a743662-8cb6-11eb-26a6-d911e60653e4
md"""
Let's take the simulation and run it a few times.
"""

# ╔═╡ e2d764d0-0845-11eb-0031-e74d2f5acaf9
function step!(infectious, p)
	for i in 1:length(infectious)

		if infectious[i] && bernoulli(p)
			infectious[i] = false
		end
	end

	return infectious
end

# ╔═╡ 9282eca0-08db-11eb-2e36-d761594b427c
T = 100

# ╔═╡ fe0aa72c-8b46-11eb-15aa-49ae570e5858
md"""
N = $(@bind N Slider(1:1000, show_value=true, default=70))

p = $(@bind ppp Slider(0:0.01:1, show_value=true, default=0.25))

t = $(@bind t Slider(1:T, show_value=true))
"""

# ╔═╡ 58d8542c-08db-11eb-193a-398ce01b8635
begin
	infected = [true for i in 1:N]

	results = [copy(step!(infected, ppp)) for i in 1:T]
	pushfirst!(results, trues(N))
end

# ╔═╡ 33f9fc36-0846-11eb-18c2-77f92fca3176
function simulate_recovery(p, T)
	infectious = trues(N)
	num_infectious = [N]

	for t in 1:T
		step!(infectious, p)
		push!(num_infectious, count(infectious))
	end

	return num_infectious
end

# ╔═╡ 39a69c2a-0846-11eb-35c1-53c68a9f71e5
p = 0.1

# ╔═╡ cb278624-08dd-11eb-3375-276bfe8d7b3a
begin
	pp = 0.05

	plot(simulate_recovery(pp, T), label="run 1", alpha=0.5, lw=2, m=:o)
	plot!(simulate_recovery(pp, T), label="run 2", alpha=0.5, lw=2, m=:o)

	xlabel!("time t")
	ylabel!("number of light bulbs that are alive")
end

# ╔═╡ f3c85814-0846-11eb-1266-63f31f351a51
all_data = [simulate_recovery(pp, T) for i in 1:30];

# ╔═╡ 01dbe272-0847-11eb-1331-4360a575ff14
begin
	plot(all_data, alpha=0.1, leg=false, m=:o, ms=1,
		size=(500, 400), label="")
	xlabel!("time t")
	ylabel!("number still functioning")
end

# ╔═╡ be8e4ac2-08dd-11eb-2f72-a9da5a750d32
plot!(mean(all_data), leg=true, label="mean",
		lw=3, c=:red, m=:o, alpha=0.5,
		size=(500, 400))

# ╔═╡ 8bc52d58-0848-11eb-3487-ef0d06061042
begin
	plot(replace.(all_data, 0.0 => NaN),
		yscale=:log10, alpha=0.3, leg=false, m=:o, ms=1,
		size=(500, 400))

	plot!(mean(all_data), yscale=:log10, lw=3, c=:red, m=:o, label="mean", alpha=0.5)

	xlabel!("time t")
	ylabel!("number still functioning")
end

# ╔═╡ caa3faa2-08e5-11eb-33fe-cbbc00cfd459
md"""
## Time evolution of the mean: Intuitive derivation
"""

# ╔═╡ 2174aeba-08e6-11eb-09a9-2d6a882a2604
md"""
The mean seems to behave in a rather predictable way over time. Can we derive this?

Let $N_t$ be the number of green light bulbs at time $t$. This decreases because some bulbs fail. Since bulbs fail with probability $p$, the number of bulbs that fail at time $t$ is, on average, $p N_t$. [Note that one time unit corresponds to one *sweep* of the simulation.]

At time $t$ there are $N_t$ green bulbs.
How many decay? Each decays with probability $p$, so *on average* $p N_t$ fail, so are removed from the number of infectious, giving the change

$$\Delta N_t = {\color{lightgreen} N_{t+1}} - {\color{lightgreen}N_t} = -{\color{red} p \, N_t}$$

So

$${\color{lightgreen} N_{t+1}} = {\color{lightgreen}N_t} - {\color{red} p \, N_t}$$

or

$$N_{t+1} = (1 - p) N_t .$$

"""

# ╔═╡ f5756dd6-0847-11eb-0870-fd06ad10b6c7
md"""
We can now take one step backwards:

$$N_{t+1} = (1 - p) (1 - p) N_{t-1} = (1 - p)^2 N_{t-1}$$

and then continue to solve the recurrence:

$$N_t = (1-p)^t \, N_0.$$
"""

# ╔═╡ 113c31b2-08ed-11eb-35ef-6b4726128eff
md"""
Let's compare the exact and numerical results:
"""

# ╔═╡ 6a545268-0846-11eb-3861-c3d5f52c061b
exact = [N * (1-pp)^t for t in 0:T]

# ╔═╡ 4c8827b8-0847-11eb-0fd1-cfbdbdcf392e
begin
	plot(mean(all_data), m=:o, alpha=0.5, label="mean of stochastic simulations",
		size=(500, 400))
	plot!(exact, lw=3, alpha=0.8, label="deterministic model", leg=:right)
	title!("Experiment vs. theory")
	xlabel!("time")
	ylabel!("""number of "greens" """)
end

# ╔═╡ 3cd1ad48-08ed-11eb-294c-f96b0e7c33bb
md"""
They agree well, as they should. The agreement is expected to be better (i.e. the fluctuations smaller) for a larger population.
"""

# ╔═╡ f32b936a-8bf6-11eb-1dd7-fd8c5904bf1f
md"""
## Binomial distribution
"""

# ╔═╡ f8f38028-8bf6-11eb-321b-8f91e38da495
md"""

At time $0$ there are $N_0$ light bulbs. How many will turn ``{\color{red} \text{red}}`` (fail) at the first step? Let's call this $\Delta N_0$.

Intuitively, the mean is $\langle \Delta N_0 \rangle = p N_0$, but in fact $\Delta N_0$ is a random variable! In principle, it could be that no light bulbs fail, or all of them fail, but both of those events have very small probability.

For each of the $N_0$ bulbs, $i=1, \ldots, N_0$, we have a Bernoulli random variable that tells us if bulb $i$ will fail. Let's call them we call $B_0^i$.
Then

$$\Delta N_0 = \sum_{i=1}^{N_0} B_0^i$$
"""

# ╔═╡ 2de1ef6c-8cb1-11eb-3dd9-f3904ec1408b
md"""
Let's make a type to represent the sum of $N$ Bernoullis with probability $p$. This is called a **binomial random variable**. The *only* information that we require is just that, $N$ and $p$.
"""

# ╔═╡ 48fb6ed6-8cb1-11eb-0894-b526e6c43b01
struct Binomial
	N::Int64
	p::Float64
end

# ╔═╡ 713a2644-8cb1-11eb-1904-f301e39d141e
md"""
Note that does not require (or even allow) methods at first, as some other languages would. You can add methods later, and other people can add methods too if they can load your package. (But they do *not* need to modify *your* code.)
"""

# ╔═╡ 511892e0-8cb1-11eb-3814-b98e8e0bbe5c
Base.rand(X::Binomial) = sum(rand(Bernoulli(X.p)) for i in 1:X.N)

# ╔═╡ 1173ebbe-8cb1-11eb-0a21-7d40a2c8a855
rand(Binomial(10, 0.25))

# ╔═╡ dfdaf1dc-8cb1-11eb-0287-f150380d323b
md"""
N = $(@bind binomial_N Slider(1:100, show_value=true, default=1));
p = $(@bind binomial_p Slider(0.0:0.01:1, show_value=true, default=0))

"""

# ╔═╡ ca3db0a8-8cb1-11eb-2f7b-c9343a29ed02
begin
	binomial_data = [rand(Binomial(binomial_N, binomial_p)) for i in 1:10000]

	bar(countmap(binomial_data), alpha=0.5, size=(500, 300), leg=false, bin_width=0.5)
end

# ╔═╡ b3ce9e3a-8c35-11eb-1ad0-81f9b09f963e
md"""
Let's call $q := 1 - p$.
Then for each bulb we are choosing either $p$ (failure) or $q$ (non-failure). (This is the same as flipping $n$ independent, weighted coins.)

The number of ways of choosing such that $k$ bulbs fail is given by the coefficient of $p^k$ in the expansion of $(p + q)^n$, namely the **binomial coefficient**

$$\begin{pmatrix} n \\ k \end{pmatrix} := \frac{n!}{k! \, (n-k)!},$$

where $n! := 1 \times 2 \times \cdots n$ is the factorial of $n$.

"""

# ╔═╡ 2f980870-0848-11eb-3edb-0d4cd1ed5b3d
md"""
## Continuous time

If we look at the graph of the mean as a function of time, it seems to follow a smooth curve. Indeed it makes sense to ask not only how many people have recovered each *day*, but to aim for finer granularity.

Suppose we instead increment time in steps of $\delta t$; the above analysis was for $\delta t = 1$.

Then we will need to adjust the probability of recovery in each time step.
It turns out that to make sense in the limit $\delta t \to 0$, we need to choose the probability $p(\delta t)$ to recover in time $t$ to be proportional to $\delta t$:

$$p(\delta t) \simeq \lambda \, \delta t,$$

where $\lambda$ is the recovery **rate**. Note that a rate is a probability *per unit time*.

We get
"""

# ╔═╡ 6af30142-08b4-11eb-3759-4d2505faf5a0
md"""
$$I(t + \delta t) - I(t) \simeq -\lambda \,\delta t \, I(t)$$
"""

# ╔═╡ c6f9feb6-08f3-11eb-0930-83385ca5f032
md"""
Dividing by $\delta t$ gives

$$\frac{I(t + \delta t) - I(t)}{\delta t} \simeq -\lambda \, I(t)$$

We recognise the left-hand side as the definition of the **derivative** when $\delta t \to 0$. Taking that limit finally gives
"""

# ╔═╡ d8d8e7d8-08b4-11eb-086e-6fdb88511c6a
md"""
$$\frac{dI(t)}{dt} = -\lambda \, I(t)$$

That is, we obtain an **ordinary differential equation** that gives the solution implicitly. Solving this equation with initial condition $I(0) = I_0$ gives
"""

# ╔═╡ 780c483a-08f4-11eb-1205-0b8aaa4b1c2d
md"""
$$I(t) = I_0 \exp(-\lambda \, t).$$
"""

# ╔═╡ a13dd444-08f4-11eb-08f5-df9dd99c8ab5
md"""
Alternatively, we can derive this by recognising the exponential in the limit $\delta t \to 0$ of the following expression, which is basically the expression for compounding interest:
"""

# ╔═╡ cb99fe22-0848-11eb-1f61-5953be879f92
md"""
$$I_{t} = (1 - \lambda \, \delta t)^{(t / \delta t)} I_0$$
"""

# ╔═╡ 8d2858a4-8c38-11eb-0b3b-61a913eed928
md"""
## Discrete to continuous
"""

# ╔═╡ 93da8b36-8c38-11eb-122a-85314d6e1921
function plot_cumulative!(p, N, δ=1; kw...)
    ps = [p * (1 - p)^(n-1) for n in 1:N]
    cumulative = cumsum(ps)

    ys = [0; reduce(vcat, [ [cumulative[n], cumulative[n]] for n in 1:N ])]

    pop!(ys)
    pushfirst!(ys, 0)

    xs = [0; reduce(vcat, [ [n*δ, n*δ] for n in 1:N ])];

    # plot!(xs, ys)
    scatter!([n*δ for n in 1:N], cumulative; kw...)
end

# ╔═╡ f1f0529a-8c39-11eb-372b-95d591a573e2
plotly()

# ╔═╡ 9572eda8-8c38-11eb-258c-739b511de833
begin
	plot(size=(500, 300), leg=false)
	plot_cumulative!(0.1, 30, 1.0, ms=2, c=:red, alpha=1)
#	plot_cumulative!(0.1, 30, 0.5, ms=2, c=:red)
	plot_cumulative!(0.05, 60, 0.5; label="", ms=2, c=:lightgreen, alpha=1)
	plot_cumulative!(0.025, 120, 0.25; label="", ms=1, c=:lightgreen, alpha=1)
	plot_cumulative!(0.0125, 240, 0.125; label="", ms=1, c=:lightgreen, alpha=1)

end

# ╔═╡ 7850b114-8c3b-11eb-276a-df5c332bf6d3
1 - 0.95^2  # almost 10%

# ╔═╡ 9f41d4f2-8c38-11eb-3eae-a1ec0d86d64c
begin
	λ = -log(1 - 0.1)

	plot!(0:0.01:20, t -> 1 - exp(-λ*t), lw=1)
	plot!(0:0.01:20, t -> 1 - exp(-0.1*t), lw=1)

end

# ╔═╡ 148f486c-8c3d-11eb-069f-cd595c5f7177
md"""
What does it mean to talk about a **rate** -- a probability per unit time.
"""

# ╔═╡ 4d61636e-8c3d-11eb-2726-6dc51e8a4f84

# ╔═╡ 3ae9fc0a-8c3d-11eb-09d5-13cefa2d9da5
md"""
How many light bulbs turn red in 1 second, half a second. Looks like 0 / 0. If have billions of light bulbs.
"""

# ╔═╡ c92bf164-8c3d-11eb-128c-7bd2c0ad681e
md"""
People get sick / light bulbs not on discrete time clock. Limit as $\delta t \to 0$

You measure it discretely
"""

# ╔═╡ 1336397c-8c3c-11eb-2ecf-eb017a3a65cd
λ

# ╔═╡ d74bace6-08f4-11eb-2a6b-891e52952f57
md"""
## SIR model
"""

# ╔═╡ dbdf2812-08f4-11eb-25e7-811522b24627
md"""
Now let's extend the procedure to the full SIR model, $S \to I \to R$. Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with probability
"""

# ╔═╡ 238f0716-0903-11eb-1595-df71600f5de7
md"""
Let's denote by $S_t$ and $I_t$ be the number of susceptible and infectious people at time $t$, respectively, and by $N$ the total number of people.

On average, in each sweep each infectious individual has the chance to interact with one other individual. That individual is chosen uniformly at random from the total population of size $N$. But a new infection occurs only if that chosen individual is susceptible, which happens with probability $S_t / N$, and then if the infection is successful, with probability $b$, say.

Hence the change in the number of infectious people after that step is.

The decrease in $S_t$ is also given by $\Delta I_t$.
"""

# ╔═╡ 8e771c8a-0903-11eb-1e34-39de4f45412b
md"""
$$\Delta I_t = I_{t+1} - I_t = b \, I_t \, \left(\frac{S_t}{N} \right)$$
"""

# ╔═╡ e83fc5b8-0904-11eb-096b-8da3a1acba12
md"""
It is useful to normalize by $N$, so we define

$$s_t := \frac{S_t}{N}; \quad i_t := \frac{I_t}{N}; \quad r_t := \frac{R_t}{N}$$
"""

# ╔═╡ d1fbea7a-0904-11eb-377d-690d7a16aa7b
md"""
Including recovery with probability $c$ we obtain the **discrete-time SIR model**:
"""

# ╔═╡ dba896a4-0904-11eb-3c47-cbbf6c01e830
md"""
$$\begin{align}
s_{t+1} &= s_t - b \, s_t \, i_t \\
i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\
r_{t+1} &= r_t + c \, i_t
\end{align}$$
"""

# ╔═╡ 267cd19e-090d-11eb-0676-0f88b57da937
md"""
Again we can obtain this from the stochastic process by taking expectations (exercise!). [Hint: Ignore recovery to start with and take variables $Y_t^i$ that are $0$ if the person is susceptible and 1 if it is infected.]
"""

# ╔═╡ 4e3c7e62-090d-11eb-3d16-e921405a6b16
md"""
And again we can allow the processes to occur in steps of length $\delta t$ and take the limit $\delta t \to 0$. With rates $\beta$ and $\gamma$ we obtain the standard (continuous-time) **SIR model**:
"""

# ╔═╡ 72061c66-090d-11eb-14c0-df619958e2b6
md"""
$$\begin{align}
\textstyle \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\
\textstyle \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) &- \gamma \, i(t)\\
\textstyle \frac{dr(t)}{dt} &= &+ \gamma \, i(t)
\end{align}$$
"""

# ╔═╡ c07367be-0987-11eb-0680-0bebd894e1be
md"""
We can think of this as a model of a chemical reaction with species S, I and R. The term $s(t) i(t)$ is known as the [**mass action**](https://en.wikipedia.org/wiki/Law_of_mass_action) form of interaction.

Note that no analytical solutions of these (simple) nonlinear ODEs are known as a function of time! (However, [parametric solutions are known](https://arxiv.org/abs/1403.2160).)
"""

# ╔═╡ f8a28ba0-0915-11eb-12d1-336f291e1d84
md"""
Below is a simulation of the discrete-time model. Note that the simplest numerical method to solve (approximately) the system of ODEs, the **Euler method**, basically reduces to solving the discrete-time model!  A whole suite of more advanced ODE solvers is provided in the [Julia `DiffEq` ecosystem](https://diffeq.sciml.ai/dev/).
"""

# ╔═╡ d994e972-090d-11eb-1b77-6d5ddb5daeab
begin
	NN = 100

	SS = NN - 1
	II = 1
	RR = 0
end

# ╔═╡ 050bffbc-0915-11eb-2925-ad11b3f67030
ss, ii, rr = SS/NN, II/NN, RR/NN

# ╔═╡ 1d0baf98-0915-11eb-2f1e-8176d14c06ad
p_infection, p_recovery = 0.1, 0.01

# ╔═╡ 28e1ec24-0915-11eb-228c-4daf9abe189b
TT = 1000

# ╔═╡ 349eb1b6-0915-11eb-36e3-1b9459c38a95
function discrete_SIR(s0, i0, r0, T=1000)

	s, i, r = s0, i0, r0

	results = [(s=s, i=i, r=r)]

	for t in 1:T

		Δi = p_infection * s * i
		Δr = p_recovery * i

		s_new = s - Δi
		i_new = i + Δi - Δr
		r_new = r      + Δr

		push!(results, (s=s_new, i=i_new, r=r_new))

		s, i, r = s_new, i_new, r_new
	end

	return results
end

# ╔═╡ 39c24ef0-0915-11eb-1a0e-c56f7dd01235
SIR = discrete_SIR(ss, ii, rr)

# ╔═╡ 442035a6-0915-11eb-21de-e11cf950f230
begin
	ts = 1:length(SIR)
	discrete_time_SIR_plot = plot(ts, [x.s for x in SIR],
		m=:o, label="S", alpha=0.2, linecolor=:blue, leg=:right, size=(400, 300))
	plot!(ts, [x.i for x in SIR], m=:o, label="I", alpha=0.2)
	plot!(ts, [x.r for x in SIR], m=:o, label="R", alpha=0.2)

	xlims!(0, 500)
end

# ╔═╡ Cell order:
# ╠═9a0cec14-08db-11eb-3cfa-4d1c327c63f1
# ╠═fb6cdc08-8b44-11eb-09f5-43c167aa53fd
# ╟─f3aad4f0-8cc2-11eb-1a25-535297327c65
# ╟─b6b055b6-8cae-11eb-29e5-b507c1a2b9bf
# ╟─bcfaedfa-8cae-11eb-10a1-cb7be7dc2e6b
# ╟─d2c19564-8b44-11eb-1077-ddf6d1395b59
# ╟─ac4f6944-08d9-11eb-0b3c-f5e0a8b8c17e
# ╟─4ca399f4-8b45-11eb-2d2b-8189e04fc804
# ╟─57080632-8b45-11eb-1003-05afb2331b25
# ╟─139ecfec-8b46-11eb-2649-2f77833d749a
# ╟─f9a75ac4-08d9-11eb-3167-011eb698a32c
# ╟─17812c7c-8cac-11eb-1d0a-6512415f6938
# ╠═178631ec-8cac-11eb-1117-5d872ba7f66e
# ╟─179a4db2-8cac-11eb-374f-0f24dc81ebeb
# ╠═17bbf532-8cac-11eb-1e3f-c54072021208
# ╟─17cf895a-8cac-11eb-017e-c79ffcab60b1
# ╠═a38fe2b2-8cae-11eb-19e8-d563e82855d3
# ╠═17e0d142-8cac-11eb-2d6a-fdf175f5d419
# ╠═18da7920-8cac-11eb-07f4-e109298fd5f1
# ╟─17fe87a0-8cac-11eb-2938-2d9cd19ecc0f
# ╟─1829091c-8cac-11eb-1b77-c5ed7dd1261b
# ╟─1851dd6a-8cac-11eb-18e4-87dbe1714be0
# ╟─a9447530-8cb6-11eb-38f7-ff69a640e3c4
# ╟─c1fde6ba-8cb6-11eb-2170-af6bc84c01a7
# ╠═18755e3e-8cac-11eb-37bf-1dfa5fbe730a
# ╠═f947a976-8cb6-11eb-2ae7-59eba4c6f40f
# ╠═5a0b407e-8cb7-11eb-0c0d-c7767a6b0a1d
# ╟─fe53ee0c-8cb6-11eb-19bc-2976da1abe16
# ╠═1894b388-8cac-11eb-2287-97f985df1fbd
# ╠═86299112-8cc6-11eb-257a-9d803feac359
# ╠═9eb69e94-8cc6-11eb-323b-5587cc743571
# ╠═b8d01ae4-8cc6-11eb-1d16-47095532f738
# ╠═c5bc97aa-8cc6-11eb-3564-fd96184ff2c3
# ╠═d5cffd96-8cc6-11eb-2714-975d46d4fa27
# ╠═da6525f0-8cc6-11eb-31ec-cd7d515f054a
# ╠═fc715452-8cc6-11eb-0246-e941f7698cfe
# ╟─71fbc75e-8bf3-11eb-3ac9-dd5401033c78
# ╟─7c01f6a6-8bf3-11eb-3c4d-ad7e206a9277
# ╠═ba7ffe78-0845-11eb-2847-851a407dd2ec
# ╠═45bb6df0-8cc7-11eb-100f-37c2a76df464
# ╟─dcd279b0-8bf3-11eb-0cb9-95f351626ed1
# ╟─fbada990-8bf3-11eb-2bb7-d786362669e8
# ╠═b6786ec8-8bf3-11eb-1347-61f231fd3b4c
# ╟─ac98f5da-8bf3-11eb-076f-597ce4455e76
# ╠═0e7a04a4-8bf4-11eb-2e9d-fb48c23b8d8c
# ╠═093275e4-8cc8-11eb-136f-3ffe522c4125
# ╟─111eccd2-8bf4-11eb-097c-7582f811d146
# ╟─4edaec4a-8bf4-11eb-3094-010ebe9b56ab
# ╟─9d66e31e-8cad-11eb-3ad0-3980ba66cb0e
# ╠═8405e310-8bf8-11eb-282b-d93b4fc683aa
# ╟─af2594c4-8cad-11eb-0fff-f59e65102b3f
# ╠═8aa60da0-8bf8-11eb-0fa2-11aeecb89564
# ╟─8d56b6d2-8cc8-11eb-0dbd-b3533dde4aa3
# ╠═a034c2a6-8bf8-11eb-0f06-0b35a0e8e68d
# ╠═c166edb6-8cc8-11eb-0436-f74164ff6ea7
# ╠═d25dc130-8cc8-11eb-177f-63a1792494c0
# ╠═e2f45234-8cc8-11eb-2be9-598eb590592a
# ╟─3ef23da4-8cb4-11eb-0d5f-d5ee8fc56227
# ╠═bc5d6fae-8cad-11eb-3351-a734d2366557
# ╠═2d9c560e-8bf9-11eb-1ac5-f77f7caf776f
# ╠═ce94541c-8bf9-11eb-1ac9-51e66a017813
# ╟─a057e7ee-8bf9-11eb-2ceb-2dda3718a70a
# ╟─4a743662-8cb6-11eb-26a6-d911e60653e4
# ╠═e2d764d0-0845-11eb-0031-e74d2f5acaf9
# ╠═9282eca0-08db-11eb-2e36-d761594b427c
# ╠═58d8542c-08db-11eb-193a-398ce01b8635
# ╟─fe0aa72c-8b46-11eb-15aa-49ae570e5858
# ╠═33f9fc36-0846-11eb-18c2-77f92fca3176
# ╠═39a69c2a-0846-11eb-35c1-53c68a9f71e5
# ╠═cb278624-08dd-11eb-3375-276bfe8d7b3a
# ╠═f3c85814-0846-11eb-1266-63f31f351a51
# ╠═01dbe272-0847-11eb-1331-4360a575ff14
# ╟─be8e4ac2-08dd-11eb-2f72-a9da5a750d32
# ╟─8bc52d58-0848-11eb-3487-ef0d06061042
# ╟─caa3faa2-08e5-11eb-33fe-cbbc00cfd459
# ╟─2174aeba-08e6-11eb-09a9-2d6a882a2604
# ╟─f5756dd6-0847-11eb-0870-fd06ad10b6c7
# ╟─113c31b2-08ed-11eb-35ef-6b4726128eff
# ╟─6a545268-0846-11eb-3861-c3d5f52c061b
# ╟─4c8827b8-0847-11eb-0fd1-cfbdbdcf392e
# ╟─3cd1ad48-08ed-11eb-294c-f96b0e7c33bb
# ╟─f32b936a-8bf6-11eb-1dd7-fd8c5904bf1f
# ╟─f8f38028-8bf6-11eb-321b-8f91e38da495
# ╟─2de1ef6c-8cb1-11eb-3dd9-f3904ec1408b
# ╠═48fb6ed6-8cb1-11eb-0894-b526e6c43b01
# ╟─713a2644-8cb1-11eb-1904-f301e39d141e
# ╠═511892e0-8cb1-11eb-3814-b98e8e0bbe5c
# ╠═1173ebbe-8cb1-11eb-0a21-7d40a2c8a855
# ╟─dfdaf1dc-8cb1-11eb-0287-f150380d323b
# ╠═ca3db0a8-8cb1-11eb-2f7b-c9343a29ed02
# ╟─b3ce9e3a-8c35-11eb-1ad0-81f9b09f963e
# ╟─2f980870-0848-11eb-3edb-0d4cd1ed5b3d
# ╟─6af30142-08b4-11eb-3759-4d2505faf5a0
# ╟─c6f9feb6-08f3-11eb-0930-83385ca5f032
# ╟─d8d8e7d8-08b4-11eb-086e-6fdb88511c6a
# ╟─780c483a-08f4-11eb-1205-0b8aaa4b1c2d
# ╟─a13dd444-08f4-11eb-08f5-df9dd99c8ab5
# ╟─cb99fe22-0848-11eb-1f61-5953be879f92
# ╟─8d2858a4-8c38-11eb-0b3b-61a913eed928
# ╠═93da8b36-8c38-11eb-122a-85314d6e1921
# ╟─f1f0529a-8c39-11eb-372b-95d591a573e2
# ╟─9572eda8-8c38-11eb-258c-739b511de833
# ╟─7850b114-8c3b-11eb-276a-df5c332bf6d3
# ╟─9f41d4f2-8c38-11eb-3eae-a1ec0d86d64c
# ╟─148f486c-8c3d-11eb-069f-cd595c5f7177
# ╟─4d61636e-8c3d-11eb-2726-6dc51e8a4f84
# ╟─3ae9fc0a-8c3d-11eb-09d5-13cefa2d9da5
# ╟─c92bf164-8c3d-11eb-128c-7bd2c0ad681e
# ╟─1336397c-8c3c-11eb-2ecf-eb017a3a65cd
# ╟─d74bace6-08f4-11eb-2a6b-891e52952f57
# ╟─dbdf2812-08f4-11eb-25e7-811522b24627
# ╟─238f0716-0903-11eb-1595-df71600f5de7
# ╟─8e771c8a-0903-11eb-1e34-39de4f45412b
# ╟─e83fc5b8-0904-11eb-096b-8da3a1acba12
# ╟─d1fbea7a-0904-11eb-377d-690d7a16aa7b
# ╟─dba896a4-0904-11eb-3c47-cbbf6c01e830
# ╟─267cd19e-090d-11eb-0676-0f88b57da937
# ╟─4e3c7e62-090d-11eb-3d16-e921405a6b16
# ╟─72061c66-090d-11eb-14c0-df619958e2b6
# ╟─c07367be-0987-11eb-0680-0bebd894e1be
# ╟─f8a28ba0-0915-11eb-12d1-336f291e1d84
# ╠═442035a6-0915-11eb-21de-e11cf950f230
# ╠═d994e972-090d-11eb-1b77-6d5ddb5daeab
# ╠═050bffbc-0915-11eb-2925-ad11b3f67030
# ╠═1d0baf98-0915-11eb-2f1e-8176d14c06ad
# ╠═28e1ec24-0915-11eb-228c-4daf9abe189b
# ╠═349eb1b6-0915-11eb-36e3-1b9459c38a95
# ╠═39c24ef0-0915-11eb-1a0e-c56f7dd01235


text:
# Handers for execute_request and related messages, which are
# the core of the IPython protocol: execution of Julia code and
# returning results.

#######################################################################
const text_plain = MIME("text/plain")
const image_svg = MIME("image/svg+xml")
const image_png = MIME("image/png")
const image_jpeg = MIME("image/jpeg")
const text_html = MIME("text/html")
const text_latex = MIME("text/latex") # IPython expects this
const text_latex2 = MIME("application/x-latex") # but this is more standard?

# backwards compatibility with old mimewritable API
if method_exists(mimewritable, (String,Any))
    mimewritable_(mime, x) = mimewritable(mime, x)
else
    mimewritable_(mime, x) = mimewritable(mime, typeof(x))
end

# return a String=>Any dictionary to attach as metadata
# in IPython display_data and pyout messages
metadata(x) = Dict()

# return a String=>String dictionary of mimetype=>data for passing to
# IPython display_data and pyout messages.
function display_dict(x)
    data = (ASCIIString=>ByteString)[ "text/plain" =>
                                          sprint(writemime, "text/plain", x) ]
    if mimewritable_(image_svg, x)
        data[string(image_svg)] = stringmime(image_svg, x)
    end
    if mimewritable_(image_png, x)
        data[string(image_png)] = stringmime(image_png, x)
    elseif mimewritable_(image_jpeg, x) # don't send jpeg if we have png
        data[string(image_jpeg)] = stringmime(image_jpeg, x)
    end
    if mimewritable_(text_html, x)
        data[string(text_html)] = stringmime(text_html, x)
    end
    if mimewritable_(text_latex, x)
        data[string(text_latex)] = stringmime(text_latex, x)
    elseif mimewritable_(text_latex2, x)
        data[string(text_latex)] = stringmime(text_latex2, x)
    end
    return data
end

# queue of objects to display at end of cell execution
const displayqueue = Any[]

# remove x from the display queue
function undisplay(x)
    i = findfirst(displayqueue, x)
    if i > 0
        splice!(displayqueue, i)
    end
    return x
end

#######################################################################

# avoid deprecation warning for change in split API
if VERSION >= v"0.4-"
    splitkeep(str, splitter) = split(str, splitter, keep=true)
else
    splitkeep(str, splitter) = split(str, splitter, true)
end

# return the content of a pyerr message for exception e
function pyerr_content(e, msg::String="")
    tb = map(utf8, splitkeep(sprint(Base.show_backtrace, :execute_request_0x535c5df2,
                                    catch_backtrace(), 1:typemax(Int)), "\n"))
    if !isempty(tb) && ismatch(r"^\s*in\s+include_string\s+", tb[end])
        pop!(tb) # don't include include_string in backtrace
    end
    ename = string(typeof(e))
    evalue = sprint(showerror, e)
    unshift!(tb, evalue) # fperez says this needs to be in traceback too
    if !isempty(msg)
        unshift!(tb, msg)
    end
    ["execution_count" => _n,
     "ename" => ename, "evalue" => evalue,
     "traceback" => tb]
end

#######################################################################
# Similar to the ipython kernel, we provide a mechanism by
# which modules can register thunk functions to be called after
# executing an input cell, e.g. to "close" the current plot in Pylab.
# Modules should only use these if isdefined(Main, IJulia) is true.

const postexecute_hooks = Function[]
push_postexecute_hook(f::Function) = push!(postexecute_hooks, f)
pop_postexecute_hook(f::Function) = splice!(postexecute_hooks, findfirst(postexecute_hooks, f))

const preexecute_hooks = Function[]
push_preexecute_hook(f::Function) = push!(preexecute_hooks, f)
pop_preexecute_hook(f::Function) = splice!(preexecute_hooks, findfirst(pretexecute_hooks, f))

# similar, but called after an error (e.g. to reset plotting state)
const posterror_hooks = Function[]
push_posterror_hook(f::Function) = push!(posterror_hooks, f)
pop_posterror_hook(f::Function) = splice!(posterror_hooks, findfirst(posterror_hooks, f))

#######################################################################

# global variable so that display can be done in the correct Msg context
execute_msg = Msg(["julia"], ["username"=>"julia", "session"=>"????"], Dict())

# note: 0x535c5df2 is a random integer to make name collisions in
# backtrace analysis less likely.
function execute_request_0x535c5df2(socket, msg)
    code = msg.content["code"]
    @vprintln("EXECUTING ", code)
    global execute_msg = msg
    global _n, In, Out, _, __, ___, ans
    silent = msg.content["silent"] || ismatch(r";\s*$", code)

    # present in spec but missing from notebook's messages:
    store_history = get(msg.content, "store_history", !silent)

    _n += 1
    if store_history
        In[_n] = code
    end
    send_ipython(publish,
                 msg_pub(msg, "pyin",
                         ["execution_count" => _n,
                          "code" => code]))

    # "; ..." cells are interpreted as shell commands for run
    code = replace(code, r"^\s*;.*$",
                   m -> string(replace(m, r"^\s*;", "Base.repl_cmd(`"),
                               "`)"), 0)

    # a cell beginning with "? ..." is interpreted as a help request
    code = replace(code, r"^\s*\?", "Base.@help ")

    try
        for hook in preexecute_hooks
            hook()
        end
        ans = result = include_string(code, "In[$_n]")
        if silent
            result = nothing
        elseif result != nothing
            ___ = __ # 3rd result from last
            __ = _ # 2nd result from last
            _ = result
            if store_history
                if result != Out # workaround for Julia #3066
                    Out[_n] = result
                    # note: don't use result in eval, since things
                    # get confused if result is a symbol:
                    eval(Main, :($(symbol(string("_",_n))) = Out[$_n]))
                end
            end
        end

        user_variables = Dict()
        user_expressions = Dict()
        for v in get(msg.content, "user_variables", String[]) # gone in IPy3
            user_variables[v] = eval(Main,parse(v))
        end
        for (v,ex) in msg.content["user_expressions"]
            user_expressions[v] = eval(Main,parse(ex))
        end

        for hook in postexecute_hooks
            hook()
        end

	# flush pending stdio
        flush_cstdio() # flush writes to stdout/stderr by external C code
	send_stream(takebuf_string(read_stdout.buffer), "stdout")
	send_stream(takebuf_string(read_stderr.buffer), "stderr")

        undisplay(result) # dequeue if needed, since we display result in pyout
        display() # flush pending display requests

        if result != nothing

            # Work around for Julia issue #265 (see # #7884 for context)
            # We have to explicitly invoke the correct metadata method.
            result_metadata = invoke(metadata, (typeof(result),), result)

            send_ipython(publish,
                         msg_pub(msg, "pyout",
                                 ["execution_count" => _n,
                                 "metadata" => result_metadata,
                                 "data" => display_dict(result) ]))
        end

        send_ipython(requests,
                     msg_reply(msg, "execute_reply",
                               ["status" => "ok", "execution_count" => _n,
                               "payload" => [],
                               "user_variables" => user_variables,
                                "user_expressions" => user_expressions]))
    catch e
        try
            # flush pending stdio
            flush_cstdio() # flush writes to stdout/stderr by external C code
            send_stream(takebuf_string(read_stdout.buffer), "stdout")
            send_stream(takebuf_string(read_stderr.buffer), "stderr")
            for hook in posterror_hooks
                hook()
            end
        catch
        end
        empty!(displayqueue) # discard pending display requests on an error
        content = pyerr_content(e)
        send_ipython(publish, msg_pub(msg, "pyerr", content))
        content["status"] = "error"
        send_ipython(requests, msg_reply(msg, "execute_reply", content))
    end
end

#######################################################################


text:
function client(master, port, log_dir)
    jobid = ENV["SLURM_JOB_ID"]
    nodeid = parse(Int,ENV["SLURM_NODEID"])
    localid = parse(Int,ENV["SLURM_LOCALID"])

    open("$(log_dir)/client-$(jobid)-$(nodeid)-$(localid).log", "w+") do log_file
        client_logger = SimpleLogger(log_file, Logging.Debug)
        @debug "Client $(nodeid)#$(localid) starting"
        @debug "Master at $(master):$(port)"
        @debug "What is thy bidding, my master?"
        run_id = -1
        while true
            try
                cl = connect(master, port)
                write(cl, "$(run_id)\n")
                run_id = parse(Int, strip(readline(cl)))
                run_id == -1 && break
                cmd = replace(strip(readline(cl)), "[local-id]"=>"$(localid)")
                @debug "Supposed to do $(cmd), why not."
                for c in split(cmd, ";")
                    cc = split(c)
                    run(`$(cc)`)
                end
                @debug "I'm all done, gimme more."
            catch e
                @error "$(e)"
                break
            end
        end
    end
end


text:
export Message

"""
A [`Message`](@ref)'s type.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-types).
"""
@enum MessageType begin
    MT_DEFAULT
    MT_RECIPIENT_ADD
    MT_RECIPIENT_REMOVE
    MT_CALL
    MT_CHANNEL_NAME_CHANGE
    MT_CHANNEL_ICON_CHANGE
    MT_CHANNEL_PINNED_MESSAGE
    MT_GUILD_MEMBER_JOIN
    MT_USER_PREMIUM_GUILD_SUBSCRIPTION
    MT_USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_1
    MT_USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_2
    MT_USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_3
    MT_CHANNEL_FOLLOW_ADD
end
@boilerplate MessageType :export :lower

"""
A [`Message`](@ref)'s activity type.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-activity-types).
"""
@enum MessageActivityType MAT_JOIN MAT_SPECTATE MAT_LISTEN MAT_JOIN_REQUEST
@boilerplate MessageActivityType :export :lower

"""
A [`Message`](@ref) activity.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-activity-structure).
"""
struct MessageActivity
    type::MessageActivityType
    party_id::Optional{String}
end
@boilerplate MessageActivity :constructors :docs :lower :merge :mock

"""
A Rich Presence [`Message`](@ref)'s application information.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-application-structure).
"""
struct MessageApplication
    id::Snowflake
    cover_image::Optional{String}
    description::String
    icon::String
    name::String
end
@boilerplate MessageApplication :constructors :docs :lower :merge :mock

"""
A message sent to a [`DiscordChannel`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object).
"""
struct Message
    id::Snowflake
    channel_id::Snowflake
    # MessageUpdate only requires the ID and channel ID.
    guild_id::Optional{Snowflake}
    author::Optional{User}
    member::Optional{Member}
    content::Optional{String}
    timestamp::Optional{DateTime}
    edited_timestamp::OptionalNullable{DateTime}
    tts::Optional{Bool}
    mention_everyone::Optional{Bool}
    mentions::Optional{Vector{User}}
    mention_roles::Optional{Vector{Snowflake}}
    attachments::Optional{Vector{Attachment}}
    embeds::Optional{Vector{Embed}}
    reactions::Optional{Vector{Reaction}}
    nonce::OptionalNullable{Snowflake}
    pinned::Optional{Bool}
    webhook_id::Optional{Snowflake}
    type::Optional{MessageType}
    activity::Optional{MessageActivity}
    application::Optional{MessageApplication}
end
@boilerplate Message :constructors :docs :lower :merge :mock


text:
function ground(video, action, framenumber)
    background("white")
    sethue("blue")
    return framenumber
end

function ground_color(color_bg, color_pen, framenumber)
    background(color_bg)
    sethue(color_pen)
    return framenumber
end

function latex_title(args...)
    fontsize(20)
    latex(L"E=mc^2", 0, -200)
end

function circ(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return p
end

function circ_ret_trans(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return Transformation(p, 0.0)
end

function path!(points, pos, color)
    sethue(color)
    push!(points, pos)
    circle.(points, 2, :fill)
end

function rad(p1, p2, color)
    sethue(color)
    line(p1, p2, :stroke)
end

@testset "Dancing circles (gif)" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:25,
                ground,
                Rotation(0.0),
                Translation(Point(25, 25), Point(25, 25)),
            ),
            Action(latex_title),
            Action(
                Rel(-24:0),
                :red_ball,
                (args...) -> circ(p1, "red"),
                Rotation(from_rot, to_rot),
            ),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:25,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(:same, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "dancing.gif",
    )

    @test_reference "refs/dancing_circles_16.png" load("images/0000000016.png")
    @test isfile("dancing.gif")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
    rm("images/palette.bmp")
    rm("dancing.gif")
end

@testset "Dancing circles (mp4)" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:2,
                ground,
                Rotation(0.0),
                Translation(Point(25, 25), Point(25, 25)),
            ),
            Action(latex_title),
            Action(
                Rel(-1:0),
                :red_ball,
                (args...) -> circ(p1, "red"),
                Rotation(from_rot, to_rot),
            ),
            Action(
                1:2,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:2,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(:same, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:2, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "dancing.mp4",
    )

    # The `y` for isapprox was determined experimentally on a Fedora 32 OS.
    # On that machine, the time duration was found to be `0.067` seconds.
    # The `atol` was also experimentally determined based upon VideoIO's
    # `get_duration` function.
    @test isapprox(VideoIO.get_duration("dancing.mp4"), 0.07, atol = 0.01)
    rm("dancing.mp4")
end

@testset "Dancing circles layered" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            Action(
                1:25,
                ground,
                Rotation(π / 2, π / 2, O),
                Translation(Point(25, 25), Point(25, 25));
                in_global_layer = true,
            ),
            Action(1:25, latex_title),
            Action(1:25, :red_ball, (args...) -> circ(p1, "red"), Rotation(to_rot)),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:25,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(1:25, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Dancing circles layered return Transformation" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            Action(
                1:25,
                ground,
                Rotation(π / 2, π / 2, O),
                Translation(Point(25, 25), Point(25, 25));
                in_global_layer = true,
            ),
            Action(1:25, latex_title),
            Action(
                1:25,
                :red_ball,
                (args...) -> circ_ret_trans(p1, "red"),
                Rotation(to_rot),
            ),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ_ret_trans(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(1:25, (video, args...) -> path!(path_of_red, pos(:red_ball), "red")),
            Action(1:25, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Drawing grid" begin
    video = Video(500, 500)
    javis(
        video,
        [
            Action(1:40, ground),
            Action(1:10, draw_grid(direction = "BL", line_gap = 25)),
            Action(zero_lines(direction = "BL", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "BR", line_gap = 25)),
            Action(zero_lines(direction = "BR", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "TL", line_gap = 25)),
            Action(zero_lines(direction = "TL", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "TR", line_gap = 25)),
            Action(zero_lines(direction = "TR", line_thickness = 10)),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/grid_drawing_bl.png" load("images/0000000008.png")
    @test_reference "refs/grid_drawing_br.png" load("images/0000000018.png")
    @test_reference "refs/grid_drawing_tl.png" load("images/0000000028.png")
    @test_reference "refs/grid_drawing_tr.png" load("images/0000000038.png")
    for i in 1:40
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Use returned angle" begin
    video = Video(500, 500)
    p = Point(100, 0)
    javis(
        video,
        [
            Action(1:10, ground),
            Action(:circ, (args...) -> circ_ret_trans(p), Rotation(2π)),
            Action(
                (args...) -> line(Point(-200, 0), Point(-200, -10 * ang(:circ)), :stroke),
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/circle_angle.png" load("images/0000000008.png")
    for i in 1:10
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

astar(args...) = star(Point(-100, -100), 30)
acirc(args...) = circle(Point(100, 100), 30)

@testset "morphing star2circle and back" begin
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:20,
                :framenumber,
                (args...) -> ground_color("white", "black", args[3]),
            ),
            Action(1:10, (args...) -> circle(Point(-100, 0), val(:framenumber), :fill)),
            Action(1:10, morph(astar, acirc)),
            Action(11:20, morph(acirc, astar)),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/star2circle5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle15.png" load("images/0000000015.png")
    for i in 1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_nicholas(args...)
    background("white")
    sethue("blue")
    setline(3)
end

function house_of_nicholas(; p1 = O, width = 100, color = "black")
    # sethue(color)
    #     .p1
    # .p2   .p3
    #
    # .p4   .p5
    width2 = div(width, 2)
    p2 = p1 + Point(-width2, width2)
    p3 = p1 + Point(width2, width2)
    p4 = p2 + Point(0, width)
    p5 = p3 + Point(0, width)
    line(p4, p5, :stroke)
    line(p5, p2, :stroke)
    line(p2, p4, :stroke)
    line(p4, p3, :stroke)
    line(p3, p1, :stroke)
    line(p1, p2, :stroke)
    line(p2, p3, :stroke)
    setline(8)
    line(p3, p5, :stroke)
end

@testset "House of Nicholas line_width" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:50, ground_nicholas),
            Action(
                (args...) -> house_of_nicholas();
                subactions = [
                    SubAction(1:25, appear(:fade_line_width)),
                    SubAction(Rel(25), disappear(:fade_line_width)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/nicholas15.png" load("images/0000000015.png")
    @test_reference "refs/nicholas25.png" load("images/0000000025.png")
    @test_reference "refs/nicholas35.png" load("images/0000000035.png")
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Squeezing a circle using scale" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:100, ground),
            Action(:start_scale, (args...) -> (1.0, 1.0)),
            Action(
                (args...) -> circ();
                subactions = [
                    SubAction(1:25, Scaling((1.0, 1.5))),
                    SubAction(Rel(25), Scaling((2.0, 1.0))),
                    SubAction(Rel(25), Scaling(:start_scale)),
                    SubAction(Rel(25), Scaling(2.0)),
                ],
            ),
            Action((args...) -> circ(Point(-100, 0))),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/squeeze15.png" load("images/0000000015.png")
    @test_reference "refs/squeeze25.png" load("images/0000000025.png")
    @test_reference "refs/squeeze35.png" load("images/0000000035.png")
    @test_reference "refs/squeeze65.png" load("images/0000000065.png")
    @test_reference "refs/squeeze85.png" load("images/0000000085.png")
    for i in 1:100
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_opacity(args...)
    background("white")
    sethue("black")
    setopacity(0.5)
end

function square_opacity(p1, w)
    setopacity(1.0)
    sethue("red")
    rect(p1, w, w, :fill)
end

@testset "Circle/square appear opacity" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:50, ground_opacity),
            Action(
                (args...) -> circ();
                subactions = [
                    SubAction(1:25, appear(:fade)),
                    SubAction(26:50, disappear(:fade)),
                ],
            ),
            Action(
                5:50,
                (args...) -> square_opacity(Point(-100, 0), 60);
                subactions = [
                    SubAction(1:15, appear(:fade)),
                    SubAction(Rel(20), Translation(100, 50)),
                    SubAction(Rel(5), disappear(:fade)),
                    # for global frames 46-50 it should still be disappeared
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/circlerSquare07opacity.png" load("images/0000000007.png")
    @test_reference "refs/circlerSquare25opacity.png" load("images/0000000025.png")
    @test_reference "refs/circlerSquare42opacity.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Scaling circle" begin
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(1:50, ground),
            Action(:set_scale, (args...) -> 2),
            Action(
                (args...) -> circ(),
                subactions = [
                    SubAction(1:15, Scaling(0.0, :set_scale)),
                    SubAction(36:50, Scaling(0.0)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")

    # test using appear and disappear
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(1:50, ground),
            BackgroundAction(1:50, (args...) -> scale(2)),
            Action(
                (args...) -> circ(),
                subactions = [
                    SubAction(1:15, appear(:scale)),
                    SubAction(36:50, disappear(:scale)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "test default kwargs" begin
    video = Video(500, 500)
    pathname = javis(video, [Action(1:10, ground), Action(1:10, morph(astar, acirc))])
    path, ext = splitext(pathname)
    @test ext == ".gif"
    @test isfile(pathname)
    rm(pathname)
end

@testset "test @error .mp3" begin
    video = Video(500, 500)
    @test_logs (:error,) javis(
        video,
        [Action(1:10, ground), Action(1:10, morph(astar, acirc))];
        pathname = "test.mp3",
    )
end


text:
#!/usr/bin/env julia
#=
This is to use a fixed set of ID as reference, and then sample more and more ID
from the rest to mask, and then to impute them back. The purpose is to test what
happes if a lot more ID are masked.
=#

using DelimitedFiles, Plots, Statistics
dat = readdlm("num.txt")
dat = reshape(dat, 2, :)'
n = length(dat[:, 1])

x=50:100:1950
plot(x, dat[:, 1],
     xlabel="Masked set size",
     label="Error rates",
     dpi=300,
     ylabel="Error rates",
     leg=false)

savefig("fix-2-less2more.png")

max_stars_repo_name: xijiang/Refactered-NSG

max_issues_repo_name: xijiang/Refactered-NSG

max_forks_repo_name: xijiang/Refactered-NSG


text:
function mean_Lx_plot(sim; species="electron", save=true)
    Lx(file) = compute_L(:x, file, species) |> unit_L
    @time m_Lx = mean(Lx, sim)
    @time m_Lx_plus = mean(Lx, sim, cond=lx->lx>zero(lx))
    @time m_Lx_minus = mean(Lx, sim, cond=lx->lx<zero(lx))
    l = get_parameter(sim, :constant, :m)

    ts = get_time.(sim) .|> unit_t

    plt = Plots.plot(ts, m_Lx,
        title="Average Lx for l = $l",
        xlabel="t",
        ylabel="Lx",
        label="Total Lx")
    Plots.plot!(plt, ts, m_Lx_plus, label="Lx > 0")
    Plots.plot!(plt, ts, m_Lx_minus, label="Lx < 0")

    save && Plots.savefig(plt, "mean_Lx.png")
    plt
end

"""
    Lx_section_plot(file, slice_location; ϵ, species="electron", colormap=:RdYlBu_11, max_size=10^4, label="", save=false)

Plot a section through the Lx of the given `species`, along the direction given by `slice_location`.
The slice has thicknes `ϵ` and the plot is downsampled to `max_size`.
"""
function Lx_section_plot(file, slice_location;
                         ϵ=2e-3unit_l,
                         species="electron",
                         colormap=:RdYlBu_11,
                         max_size=10^4,
                         label="",
                         save=false)

    Lx = uconvert(unit_L, unit_l, compute_L(:x, file, species))
    if all(iszero.(Lx))
        @warn "Lx is 0"
        return nothing
    end
    Lx_section = downsample(slice(Lx, :x, slice_location, ϵ), max_size)
    if isempty(Lx_section)
        @warn "Empty slice"
        return nothing
    end

    cl = max(abs.(extrema(Lx))...)

    bg = Plots.cgrad(colormap, categorical=true)[6]
    t = sprint(show, get_time(file), context=:compact => true)

    plt = Plots.plot(Lx_section;
        xlabel = "y",
        ylabel = "z",
        clims = (-cl,cl),
        cb_title = "Lx",
        title = "Lx at x = $slice_location and t = $t",
        seriescolor = colormap,
        background_inside = bg,
        framestyle = :box,
        label)

    save && Plots.savefig(plt, "Lx_x$(slice_location)_eps$(ϵ)_t$t.png")

    return plt
end


text:
# Translation from character to DP A B C D E F G pins

const CHAR_TRANSLATOR = Dict(
    '0' => 0b01111110,
    '1' => 0b00110000,
    '2' => 0b01101101,
    '3' => 0b01111001,
    '4' => 0b00110011,
    '5' => 0b01011011,
    '6' => 0b01011111,
    '7' => 0b01110000,
    '8' => 0b01111111,
    '9' => 0b01111011,

    'a' => 0b01110111, # A
    'A' => 0b01110111,
    'b' => 0b01111011,
    'B' => 0b01111111,
    'c' => 0b00001101,
    'C' => 0b01111011,
    'd' => 0b01111011,
    'D' => 0b01111011, # d
    'e' => 0b01101111,
    'E' => 0b01001111,
    'f' => 0b01000111, # F
    'F' => 0b01000111,
    'g' => 0b01111011,
    'G' => 0b01111011, # g
    'h' => 0b00010111,
    'H' => 0b00110111,
    'i' => 0b00110000, # I
    'I' => 0b00110000,
    'j' => 0b01111100, # J
    'J' => 0b01111100,
    # 'k'
    'l' => 0b00000110,
    'L' => 0b00001110,
    # 'm'
    # 'n'
    'o' => 0b00011101,
    'O' => 0b01111110,
    'p' => 0b01100111, # P
    'P' => 0b01100111,
    'q' => 0b01110011,
    'Q' => 0b01110011, # q
    # 'r'
    's' => 0b01011011, # S
    'S' => 0b01011011,
    # 't'
    'u' => 0b00111110, # U
    'U' => 0b00111110,
    # 'v'
    # 'w'
    # 'x'
    'y' => 0b00111011,
    'Y' => 0b00111011, # y
    # 'z'

    '-' => 0b00000001,
    '_' => 0b00001000,
    ' ' => 0b00000000
)

const NUM_TRANSLATOR = Dict(
    UInt8(0) => 0b01111110,
    UInt8(1) => 0b00110000,
    UInt8(2) => 0b01101101,
    UInt8(3) => 0b01111001,
    UInt8(4) => 0b00110011,
    UInt8(5) => 0b01011011,
    UInt8(6) => 0b01011111,
    UInt8(7) => 0b01110000,
    UInt8(8) => 0b01111111,
    UInt8(9) => 0b01111011,

    UInt8(10) => 0b00000001, # -
    UInt8(11) => 0b01001111, # E
    UInt8(12) => 0b00110111, # H
    UInt8(13) => 0b00001110, # L
    UInt8(14) => 0b01100111, # P
    UInt8(15) => 0b00000000  # blank
)


text:
using Test
using WikstromTerelius: ElGamal, PrimeGenerator, prove, verify, Simulator, Enc, Dec, gen_shuffle, Verifier, PoSChallenge, Shuffle

import WikstromTerelius: step, challenge, PoSChallenge

@enum VState Config Init PermCommit PoSCommit

###
struct HonestVerifier{T} <: Verifier
    challenge::PoSChallenge
end

HonestVerifier(challenge::PoSChallenge) = HonestVerifier{Config}(challenge)
HonestVerifier{T}(verifier::HonestVerifier) where T = HonestVerifier{T}(verifier.challenge)

PoSChallenge(verifier::HonestVerifier{PoSCommit}) = verifier.challenge

step(verifier::HonestVerifier{Config}, proposition::Shuffle) = HonestVerifier{Init}(verifier)
step(verifier::HonestVerifier{Init}, 𝐜) = HonestVerifier{PermCommit}(verifier)
step(verifier::HonestVerifier{PermCommit}, 𝐜̂, t) = HonestVerifier{PoSCommit}(verifier)
#step(verifier::HonestVerifier{PoSCommit}, s) = HonestVerifier{End}(verifier)

challenge(verifier::HonestVerifier{Init}) = (verifier.challenge.𝐡, verifier.challenge.𝐡[1])
challenge(verifier::HonestVerifier{PermCommit}) = verifier.challenge.𝐮
challenge(verifier::HonestVerifier{PoSCommit}) = verifier.challenge.c

p = 23
g = PrimeGenerator(3, p)

sk = 5
pk = g^sk

enc = Enc(pk, g)
dec = Dec(sk)

𝐦 = [g, g^2, g^3]
𝐞 = enc(𝐦, [2, 3, 4])

N = length(𝐞)

𝐡 = [g^i for i in 2:N+1]

𝐫′ = [4, 2, 3]
proposition, secret = gen_shuffle(enc, 𝐞, 𝐫′) # In practice total of random factors can't match as it reveals

(; 𝛙) = secret
(; 𝐞, 𝐞′) = proposition
@test dec(𝐞)[𝛙] == dec(𝐞′)

𝐡 = [g^i for i in 2:N+1]
𝐮 = [3, 4, 5]
c = 9

chg = PoSChallenge(𝐡, 𝐮, c)

verifier = HonestVerifier(chg)

simulator = prove(proposition, secret, verifier)

@test verify(simulator)


text:
using Test

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

import MathOptFormat
const SDPA = MathOptFormat.SDPA

const SDPA_TEST_FILE = "test.sdpa"
const SDPA_MODELS_DIR = joinpath(@__DIR__, "models")

function set_var_and_con_names(model::MOI.ModelLike)
    variable_names = String[]
    for j in MOI.get(model, MOI.ListOfVariableIndices())
        var_name_j = "v" * string(j.value)
        push!(variable_names, var_name_j)
        MOI.set(model, MOI.VariableName(), j, var_name_j)
    end

    idx = 0
    constraint_names = String[]
    for i in Iterators.flatten((
        MOI.get(model, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}()),
        MOI.get(model, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle}())))
        idx += 1
        con_name_i = "c" * string(idx)
        push!(constraint_names, con_name_i)
        MOI.set(model, MOI.ConstraintName(), i, con_name_i)
    end

    return (variable_names, constraint_names)
end

function test_write_then_read(model_string::String)
    model1 = SDPA.Model()
    MOIU.loadfromstring!(model1, model_string)
    (variable_names, constraint_names) = set_var_and_con_names(model1)

    MOI.write_to_file(model1, SDPA_TEST_FILE)
    model2 = SDPA.Model()
    MOI.read_from_file(model2, SDPA_TEST_FILE)
    set_var_and_con_names(model2)
    if MOI.get(model1, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        MOI.set(model2, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        attr = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
        obj = MOI.get(model2, attr)
        MOI.set(model2, attr, MOIU.operate(-, Float64, obj))
    end

    MOIU.test_models_equal(model1, model2, variable_names, constraint_names)
end

function test_read(filename::String, model_string::String)
    model1 = SDPA.Model()
    MOIU.loadfromstring!(model1, model_string)
    (variable_names, constraint_names) = set_var_and_con_names(model1)

    model2 = SDPA.Model()
    MOI.read_from_file(model2, filename)
    set_var_and_con_names(model2)

    MOIU.test_models_equal(model1, model2, variable_names, constraint_names)
end

@test sprint(show, SDPA.Model()) == "A SemiDefinite Programming Algorithm Format (SDPA) model"

@testset "Support errors" begin
    @testset "$set variable bound" for set in [
            MOI.EqualTo(1.0),
            MOI.LessThan(1.0),
            MOI.GreaterThan(1.0),
            MOI.Interval(1.0, 2.0),
            MOI.Semiinteger(1.0, 2.0),
            MOI.Semicontinuous(1.0, 2.0),
            MOI.Integer(),
            MOI.ZeroOne()
        ]
        model_string = """
        variables: x
        minobjective: 1x
        c: x in $set
        """
        model = SDPA.Model()
        @test !MOI.supports_constraint(model, MOI.SingleVariable, typeof(set))
        err = MOI.UnsupportedConstraint{MOI.SingleVariable, typeof(set)}
        @test_throws err MOIU.loadfromstring!(model, model_string)
    end
end

@testset "Deleted variables error with $T" for T in [Int, Float64]
    model = SDPA.Model(; number_type = T)
    x = MOI.add_variable(model)
    MOI.delete(model, x)
    y = MOI.add_variable(model)
    fy = MOI.SingleVariable(y)
    MOI.add_constraint(model, MOIU.vectorize([one(T) * fy]), MOI.Nonnegatives(1))
    err = ErrorException("Non-contiguous variable indices not supported. This might be due to deleted variables.")
    @test_throws err MOI.write_to_file(model, SDPA_TEST_FILE)
end

@testset "Objective function with $T" for T in [Int, Float64]
    model = SDPA.Model(; number_type = T)
    @test !MOI.supports(model, MOI.ObjectiveFunction{MOI.SingleVariable}())
    @test !MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}())
end

@testset "Read errors" begin
    @testset "Non-empty model" begin
        model = SDPA.Model()
        MOI.add_variable(model)
        err = ErrorException("Cannot read in file because model is not empty.")
        @test_throws err MOI.read_from_file(model,
            joinpath(SDPA_MODELS_DIR, "example_A.sdpa"))
    end

    @testset "Bad number of blocks" begin
        model = SDPA.Model()
        err = ErrorException("The number of blocks (3) does not match the length of the list of blocks dimensions (2).")
        @test_throws err MOI.read_from_file(model,
            joinpath(SDPA_MODELS_DIR, "bad_blocks.sdpa"))
    end

    @testset "Bad number of variables" begin
        model = SDPA.Model()
        err = ErrorException("The number of variables (3) does not match the length of the list of coefficients for the objective function vector of coefficients (2).")
        @test_throws err MOI.read_from_file(model,
            joinpath(SDPA_MODELS_DIR, "bad_vars.sdpa"))
    end

    @testset "Wrong number of values in entry" begin
        model = SDPA.Model()
        err = ErrorException("Invalid line specifying entry: 0 1 2 2. There are 4 values instead of 5.")
        @test_throws err MOI.read_from_file(model,
            joinpath(SDPA_MODELS_DIR, "bad_entry.sdpa"))
    end

    @testset "Non-diagonal entry in diagonal block" begin
        model = SDPA.Model()
        err = ErrorException("Invalid line specifying entry: 0 1 1 2 1.0. `1 != 2` while block 1 has dimension 2 so it is a diagonal block.")
        @test_throws err MOI.read_from_file(model,
            joinpath(SDPA_MODELS_DIR, "bad_diag.sdpa"))
    end
end

@testset "Write errors" begin
    @testset "Nonzero constant in objective" begin
        model = SDPA.Model()
        MOIU.loadfromstring!(model, """
            variables: x
            minobjective: x + 1
        """)
        err = ErrorException(
            "Nonzero constant in objective function not supported. Note that " *
            "the constant may be added by the substitution of a bridged variable."
        )
        @test_throws err MOI.write_to_file(model, SDPA_TEST_FILE)
    end

    # TODO NLP not supported test.
end

@testset "Model name" begin
    model = SDPA.Model()
    MOI.set(model, MOI.Name(), "FooBar")
    MOI.write_to_file(model, SDPA_TEST_FILE)
    @test readlines(SDPA_TEST_FILE) == ["\"FooBar", "0", "0", "", ""]
end

write_read_models = [
    ("min ScalarAffine", """
        variables: x, y
        minobjective: 1.2x + -1y
    """),
    ("max ScalarAffine", """
        variables: x, y
        maxobjective: 1.2x + -1y
    """),
    ("VectorAffineFunction in Nonnegatives", """
        variables: x, y
        minobjective: 1.2x
        c1: [1.1 * x, y + 1] in Nonnegatives(2)
    """),
    ("VectorAffineFunction in PositiveSemidefiniteConeTriangle", """
        variables: x, y, z
        minobjective: 1.2x
        c1: [1.1x, y + 1, 2x + z] in PositiveSemidefiniteConeTriangle(2)
    """),
]
@testset "Write/read $model_name" for (model_name, model_string) in
    write_read_models
    test_write_then_read(model_string)
end

example_models = [
    ("example_A.sdpa", """
        variables: x, y
        minobjective: 10x + 20y
        c1: [x + 1, 0, x + 2] in PositiveSemidefiniteConeTriangle(2)
        c2: [5y + 3, 4y, 6y + 4] in PositiveSemidefiniteConeTriangle(2)
    """),
    ("example_B.sdpa", """
        variables: x
        minobjective: 1x
        c1: [0, 2x + 2, 0] in PositiveSemidefiniteConeTriangle(2)
    """),
]
@testset "Read and write/read $model_name" for (model_name, model_string) in example_models
    test_read(joinpath(SDPA_MODELS_DIR, model_name), model_string)
    test_write_then_read(model_string)
end

# Clean up.
sleep(1.0)  # Allow time for unlink to happen.
rm(SDPA_TEST_FILE, force = true)