text:
const DiagnosticSeverity = Dict("Error" => 1, "Warning" => 2, "Information" => 3, "Hint" => 4)
const TextDocumentSyncKind = Dict("None" => 0, "Full" => 1, "Incremental" => 2)
const WatchKind = Dict("Create" => 1,
                       "Change" => 2,
                       "Delete" => 3)

const TextDocumentReason = Dict("Manual" => 1,
                                "AfterDelay" => 2,
                                "FocusOut" => 3)

const FileChangeType = Dict("Created" => 1, "Changed" => 2, "Deleted" => 3)
const FileChangeType_Created = 1
const FileChangeType_Changed = 2
const FileChangeType_Deleted = 3
const DocumentHighlightKind = Dict("Text" => 1, "Read" => 2, "Write" => 3)

haskeynotnull(d::Dict, k) = haskey(d, k) && d[k] != nothing

macro json_read(arg)
    ex = quote
        $((arg))

        function $((arg.args[2]))(d::Dict)
        end
    end
    fex = :($((arg.args[2]))())
    for a in arg.args[3].args
        if !(a isa LineNumberNode)
            fn = string(a.args[1])
            if a.args[2] isa Expr && length(a.args[2].args) > 1 && a.args[2].head == :curly && a.args[2].args[1] == :Union
                isnullable = true
                t = a.args[2].args[3]
            else
                isnullable = false
                t = a.args[2]
            end
            if t isa Expr && t.head == :curly && t.args[2] != :Any
                f = :($(t.args[2]).(d[$fn]))
            elseif t != :Any
                f = :($(t)(d[$fn]))
            else
                f = :(d[$fn])
            end
            if isnullable
                f = :(haskeynotnull(d,$fn) ? $f : nothing)
            end
            push!(fex.args, f)
        end
    end
    push!(ex.args[end].args[2].args, fex)
    return esc(ex)
end

include("basic.jl")
include("configuration.jl")
include("document.jl")
include("providers.jl")

const MessageType = Dict("Error" => 1, "Warning" => 2, "Info" => 3, "Log" => 4)

function Message(t::Int, text::AbstractString)
    Dict("jsonrpc" => "2.0", "method" => "window/showMessage", "params" => Dict("type" => t, "message" => text))
end

mutable struct ShowMessageParams
    typ::Integer
    message::String
end

JSON.lower(a::ShowMessageParams) = Dict("type" => a.typ, "message" => a.message)

mutable struct MessageActionItem
    title::String
end

mutable struct ShowMessageRequestParams
    typ::Integer
    message::String
    actions::Union{Nothing,Vector{MessageActionItem}}
end

function JSON.lower(a::ShowMessageRequestParams)
    d = Dict("type" => a.typ, "message" => a.message)
    if a.actions isa Nothing
        d["actions"] = JSON.lower(a.actions)
    end
    return d
end

mutable struct LogMessageParams
    typ::Integer
    message::String
end

JSON.lower(a::LogMessageParams) = Dict("type" => a.typ, "message" => a.message)

mutable struct CancelParams
    id::Union{String,Int64}
end
CancelParams(d::Dict) = CancelParams(d["id"])


text:
struct FreeMemoryBlock
    start_address::UInt16
    length::UInt32
end

ext: jl

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

using ValidatedNumerics
using FactCheck

setprecision(Interval, Float64)

facts("setdiff") do
    x = 2..4
    y = 3..5

    d = setdiff(x, y)

    @fact typeof(d) --> Vector{Interval{Float64}}
    @fact length(d) --> 1
    @fact d --> [2..3]
    @fact setdiff(y, x) --> [4..5]

    x = 2..4
    y = 2..5

    @fact typeof(d) --> Vector{Interval{Float64}}
    @fact length(setdiff(x, y)) --> 0
    @fact setdiff(y, x) --> [4..5]

    x = 2..5
    y = 3..4
    @fact setdiff(x, y) --> [2..3, 4..5]

#    X = IntervalBox(2..4, 3..5)
#    Y = IntervalBox(3..5, 4..6)

    #@fact setdiff(X, Y) --> IntervalBox(2..3, 3..4)
end


text:
module VoronoiCells

using Compat

using VoronoiDelaunay
import VoronoiDelaunay: getx, gety, getgena, getgenb

export
	# Types
	IndexablePoint2D,
	Tessellation,

	# Functions
	voronoicells,
	voronoiarea,
	polyarea,
	clip,
	density,
	large2small,
	small2large

include("Types.jl")
include("Area.jl")
include("Clipping.jl")
include("Corners.jl")
include("Density.jl")
include("Misc.jl")

end # module


text:
# TODO: This runs but is outdated and isn't a good reference anymore.

module Wager

using Dates
using Discord

# Guild ID -> username -> token count.
const TOKENS = Dict{Discord.Snowflake, Dict{String, UInt}}()
const TOKEN_START = 100
const TOKEN_INTERVAL = Minute(30)
const TOKEN_INCREMENT = 100

"""
Add users to the token cache via a `GuildCreate` event.
"""
function add_users(c::Client, e::GuildCreate)
    if !haskey(TOKENS, e.guild.id)
        TOKENS[e.guild.id] = Dict()
    end

    guild = TOKENS[e.guild.id]

    for m in e.guild.members
        if !haskey(guild, m.user.username)
            guild[m.user.username] = TOKEN_START
        end
    end
end

"""
Give out tokens to all users on an interval.
"""
function distribute_tokens(c::Client)
    while isopen(c)
        for g in keys(TOKENS)
            for u in keys(g)
                g[u] += TOKEN_INCREMENT
            end
        end
        sleep(TOKEN_INTERVAL)
    end
end

"""
Insert a guild and/or user from a message into the token cache if they don't exist.
"""
function ensure_updated(m::Message)
    if !haskey(TOKENS, m.guild_id)
        TOKENS[m.guild_id] = Dict()
    end
    if !haskey(TOKENS[m.guild_id], m.author.username)
        TOKENS[m.guild_id][m.author.username] = TOKEN_START
    end
end

"""
Get the token count for the user who sent a message.
"""
token_count(m::Message) = get(get(TOKENS, m.guild_id, Dict()), m.author.username, 0)

"""
Reply to a message with the author's token count via a `MessageCreate` event.
"""
function reply_token_count(c::Client, m::Message)
    ensure_updated(m)
    reply(c, m, "You have $(token_count(m)) tokens.")
end

"""
Reply to a message with the guild's token leaderboard via a `MessageCreate` event.
"""
function reply_token_leaderboard(c::Client, m::Message)
    ensure_updated(m)

    # Get user => token count pairs by token count in descending order.
    sorted = sort(collect(TOKENS[m.guild_id]); by=p -> p.second, rev=true)

    entries = String[]
    for i in 1:min(10, length(sorted))  # Display at most 10.
        user, tokens = sorted[i]
        push!(entries, "$user: $tokens")
    end

    reply(c, m, join(entries, "\n"))
end

"""
Transfer tokens from one user to another via a `MessageCreate` event.
"""
function send_tokens(c::Client, m::Message)
    ensure_updated(m)

    words = split(m.content)
    if length(words) < 3
        return reply(c, m, "Invalid !send command.")
    end

    tokens = tryparse(UInt, words[2])
    tokens === nothing && return reply(c, m, "'$(words[2])' is not a valid number.")

    recipient = words[3]
    if !haskey(TOKENS[m.guild_id], recipient)
        return reply(c, m, "Couldn't find user '$recipient' in this guild.")
    end
    if token_count(m) < tokens
        return reply(c, m, "You don't have enough tokens to give.")
    end

    TOKENS[m.guild_id][m.author.username] -= tokens
    TOKENS[m.guild_id][recipient] += tokens
    reply(c, m, "You sent $tokens tokens to $recipient.")
end

"""
Wager a user's tokens via a `MessageCreate` event. The user either doubles their bet or
loses it.
"""
function wager_tokens(c::Client, m::Message)
    ensure_updated(m)

    words = split(m.content)
    if length(words) < 2
        return reply(c, m, "Invalid !wager command.")
    end

    tokens = try
        parse(UInt, words[2])
    catch
        return reply(c, m, "'$(words[2])' is not a valid number of tokens.")
    end
    if token_count(m) < tokens
        return reply(c, m, "You don't have enough tokens to wager.")
    end

    if rand() > 0.5
        TOKENS[m.guild_id][m.author.username] += tokens
        reply(c, m, "You win!")
    else
        TOKENS[m.guild_id][m.author.username] -= tokens
        reply(c, m, "You lose.")
    end
end

function main()
    c = Client(ENV["DISCORD_TOKEN"]; prefix='!')
    add_handler!(c, GuildCreate, add_users)
    add_command!(c, :count, reply_token_count)
    add_command!(c, :leaderboard, reply_token_leaderboard)
    add_command!(c, :send, send_tokens)
    add_command!(c, :wager, wager_tokens)
    open(c)
    @async distribute_tokens(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Wager.main()
    wait(c)
end


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

mutable struct SubnetAssociation <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    securityRules::Any # spec type: Union{ Nothing, Vector{SecurityRule} } # spec name: securityRules

    function SubnetAssociation(;id=nothing, securityRules=nothing)
        o = new()
        validate_property(SubnetAssociation, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(SubnetAssociation, Symbol("securityRules"), securityRules)
        setfield!(o, Symbol("securityRules"), securityRules)
        o
    end
end # type SubnetAssociation

const _property_map_SubnetAssociation = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("securityRules")=>Symbol("securityRules"))
const _property_types_SubnetAssociation = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("securityRules")=>"Vector{SecurityRule}")
Base.propertynames(::Type{ SubnetAssociation }) = collect(keys(_property_map_SubnetAssociation))
Swagger.property_type(::Type{ SubnetAssociation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SubnetAssociation[name]))}
Swagger.field_name(::Type{ SubnetAssociation }, property_name::Symbol) =  _property_map_SubnetAssociation[property_name]

function check_required(o::SubnetAssociation)
    true
end

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


text:
# Test the finite differencing technique
@testset "finite_difference" begin

function f!(fvec, x)
    fvec[1] = (x[1]+3)*(x[2]^3-7)+18
    fvec[2] = sin(x[2]*exp(x[1])-1)
end

r = nlsolve(f!, [ 0.1; 1.2])

@test norm(r.zero - [ 0; 1]) < 1e-8

end


text:
import Base: getproperty, setproperty!, setindex!
import Sockets

struct PDict
    dict::Dict{Symbol, Any}
    proto::Dict{Symbol, Function}
    PDict(x...) = new(Dict{Symbol, Any}(x...), Dict{Symbol, Function}())
end

function getproperty(d::PDict, field::Symbol)
    dict = getfield(d, :dict)
    field in keys(dict) && return dict[field]

    proto = getfield(d, :proto)
    field in keys(proto) && return proto[field](d)

    # throw(KeyError(field))
    nothing
end

function setproperty!(d::PDict, field::Symbol, value)
    getfield(d, :dict)[field] = value
end

function setindex!(d::PDict, f::Function, field::Symbol)
    getfield(d, :proto)[field] = f
end

function context(http::HTTP.Stream)
    req = PDict(:_http=>http, :method=>http.message.method, :uri=>HTTP.URIs.URI(http.message.target))
    req[:header] = this -> this.header = Dict(this._http.message.headers)

    res = PDict(:headers=>Pair{String, String}[], :status=>500)

    req, res
end

function app()
    routing_rules, server = [], nothing

    function _add_handler(method)
        function handle(f::Function, path, fs...)
            push!(routing_rules, (method, path, fs..., f))
        end
        function handle(path, fs...)
            push!(routing_rules, (method, path, fs...))
        end
        handle
    end

    function _listen(address="127.0.0.1", port=3001)
        routing_tree = build_routing_tree(routing_rules)
        server = Sockets.listen(parse(Sockets.IPAddr, address), port)

        HTTP.listen(address, port, server=server) do http::HTTP.Stream
            req, res = context(http)
            routing_tree(req, res)
            HTTP.setstatus(http, res.status)
            for (k, v) in res.headers
                HTTP.setheader(http, k => v)
            end
            HTTP.startwrite(http)
            res.body != nothing && write(http, res.body)
        end
    end

    (
        route = _add_handler("HOOK"),
        [method => _add_handler(uppercase(string(method))) for method in methods]...,
        listen = _listen,
        close = () -> server != nothing && close(server)
    )
end


text:
@struct immutable cpu_sample
    cpu::UInt32
    res::UInt32
end

function dump_sample_dict(io, sample_dict)
    for (k,v) in sample_dict
        println(io, PERF_SAMPLE_TYPE[k]," = ", v)
    end
end

immutable counter_value
  value::UInt64
  id::UInt64
end

immutable time_value
  time_enabled::UInt64
  time_running::UInt64
end

function extract_read_format(buf, attr)
    values = counter_value[]
    value = id = time_enabled = time_running = 0
    nr = 1
    if (attr.read_type & PERF_FORMAT_GROUP) != 0
      nr = read(buf, UInt64)
    else
      value = read(buf, UInt64)
    end
    if (attr.read_type & PERF_FORMAT_TOTAL_TIME_ENABLED) != 0
      time_enabled = read(buf, UInt64)
    end
    if (attr.read_type & PERF_FORMAT_TOTAL_TIME_RUNNING) != 0
      time_running = read(buf, UInt64)
    end
    if (attr.read_type & PERF_FORMAT_GROUP) != 0
      for i = 1:nr
        value = read(buf, UInt64)
        if (attr.read_type & PERF_FORMAT_ID) != 0
          id = read(buf, UInt64)
        end
        push!(values, counter_value(value, id))
      end
    else
      if (attr.read_type & PERF_FORMAT_ID) != 0
        id = read(buf, UInt64)
      end
      push!(values, counter_value(value, id))
    end
    (time_value(time_enabled,time_running), values)
end

extract_sample_kinds(data, attr::perf_file_attr) = extract_sample_kinds(data, attr.attr)
function extract_sample_kinds(data, attr)
    sample_type = attr.sample_type
    sample_data = Dict{UInt64, Any}()
    buf = IOBuffer(data)
    for (kind, T) in ((PERF_SAMPLE_IDENTIFIER, UInt64),
                      (PERF_SAMPLE_IP, UInt64),
                      (PERF_SAMPLE_TID, tid_record),
                      (PERF_SAMPLE_TIME, UInt64),
                      (PERF_SAMPLE_ADDR, UInt64),
                      (PERF_SAMPLE_ID, UInt64),
                      (PERF_SAMPLE_STREAM_ID, UInt64),
                      (PERF_SAMPLE_CPU, cpu_sample),
                      (PERF_SAMPLE_PERIOD, UInt64))
      if (sample_type & kind) != 0
          sample_data[kind] = unpack(buf, T)
      end
    end
    if (sample_type & PERF_SAMPLE_READ) != 0
      sample_data[PERF_SAMPLE_READ] = extract_read_format(buf, attr)
    end
    if (sample_type & PERF_SAMPLE_CALLCHAIN) != 0
      nr = read(buf, UInt64)
      sample_data[PERF_SAMPLE_CALLCHAIN] = read(buf, UInt64, nr)
    end
    if (sample_type & PERF_SAMPLE_RAW) != 0
      size = read(buf, UInt32)
      sample_data[PERF_SAMPLE_RAW] = read(buf, UInt8, size)
    end
    if (sample_type & PERF_SAMPLE_BRANCH_STACK) != 0
      error("Unsupported")
    end
    if (sample_type & PERF_SAMPLE_REGS_USER) != 0
      abi = read(buf, UInt64)
      sample_data[PERF_SAMPLE_REGS_USER] = abi == 0 ? UInt64[] :
        read(buf, UInt64, count_ones(attr.sample_regs_user))
    end
    if (sample_type & PERF_SAMPLE_STACK_USER) != 0
      size = read(buf, UInt64)
      if size == 0
        sample_data[PERF_SAMPLE_STACK_USER] = UInt8[]
      else
        data = read(buf, UInt8, size)
        dyn_size = read(buf, UInt64)
        resize!(data, dyn_size)
        sample_data[PERF_SAMPLE_STACK_USER] = data
      end
    end
    for kind in (PERF_SAMPLE_WEIGHT,
                      PERF_SAMPLE_DATA_SRC,
                      PERF_SAMPLE_TRANSACTION)
      if (sample_type & kind) != 0
          sample_data[kind] = unpack(buf, UInt64)
      end
    end
    if (sample_type & PERF_SAMPLE_REGS_INTR) != 0
      abi = read(buf, UInt64)
      sample_data[PERF_SAMPLE_REGS_INTR] = (abi, read(buf, count_ones(attr.sample_regs_intr)))
    end
    sample_data
end

function is_perf_context_ip(ip)
  UInt64(ip) in (PERF_CONTEXT_HV,
    PERF_CONTEXT_KERNEL, PERF_CONTEXT_USER, PERF_CONTEXT_GUEST,
    PERF_CONTEXT_GUEST_KERNEL, PERF_CONTEXT_GUEST_USER)
end


text:
function learnable_normal(x::Float64, y::Float64)
    l = learnable(:l)
    m = learnable(:m)
    q = rand(:q, Normal(l, 3 + m^2))
    return q
end

@time @testset "Convergence for learning - MLE 1" begin
    sel = target([(:q, ) => 6.0])
    params = learnables([(:l, ) => 3.0,
                         (:m, ) => 10.0])
    ret, cl, w = generate(sel, params, learnable_normal, 5.0, 3.0)
    for i in 1 : 40
        δ, params, acc = mle(params, cl)
        δ < 1e-12 && acc && break
    end
    @test params[:l] ≈ 6.0 atol = 1e-2
    @test params[:m] ≈ 0.0 atol = 1e-2
end

@time @testset "Convergence for learning - MAP 1" begin
end

function model()
    slope = rand(:slope, Normal(-1.0, exp(0.5)))
    intercept = rand(:intercept, Normal(1.0, exp(2.0)))
end

function var()
    slope_mu = learnable(:slope_mu)
    slope_log_std = learnable(:slope_log_std)
    intercept_mu = learnable(:intercept_mu)
    intercept_log_std = learnable(:intercept_log_std)
    slope = rand(:slope, Normal(slope_mu, exp(slope_log_std)))
    intercept = rand(:intercept, Normal(intercept_mu, exp(intercept_log_std)))
end

@time @testset "Convergence for learning - ADVI 1" begin
    sel = target()
    params = learnables([(:slope_mu, ) => 0.0,
                         (:slope_log_std, ) => 0.0,
                         (:intercept_mu, ) => 0.0,
                         (:intercept_log_std, ) => 0.0])
    opt = Momentum(0.25, 0.9)
    for i in 1 : 1000
        params, _ = advi(opt, sel, params, var, (), model, ();
                         gs_samples = 100)
    end
    @test params[:slope_mu] ≈ -1.0 atol = 7e-2
    @test params[:slope_log_std] ≈ 0.5 atol = 7e-2
    @test params[:intercept_mu] ≈ 1.0 atol = 7e-2
    @test params[:intercept_log_std] ≈ 2.0 atol = 7e-2
end

#@testset "Convergence for learning - geo VIMCO 1" begin
#    sel = target()
#    params = learnables([(:slope_mu, ) => 0.0,
#                         (:slope_log_std, ) => 0.0,
#                         (:intercept_mu, ) => 0.0,
#                         (:intercept_log_std, ) => 0.0])
#    opt = Momentum(0.3, 0.9)
#    for i in 1 : 50
#        params, _ = adgv(opt, sel, params, var, (), model, ();
#                         est_samples = 20, gs_samples = 100)
#        display(params)
#    end
#    @test params[:slope_mu] ≈ -1.0 atol = 7e-2
#    @test params[:slope_log_std] ≈ 0.5 atol = 7e-2
#    @test params[:intercept_mu] ≈ 1.0 atol = 7e-2
#    @test params[:intercept_log_std] ≈ 2.0 atol = 7e-2
#end


text:
function getproposal(p::Prior, nparams)

  newparams = zeros(Float64, nparams)
  update_newparams!(newparams, p)

  return newparams
end

update_newparams!(newparams,p::Prior) = update_newparams!(newparams, 1, p.distribution...)

@inline function update_newparams!(newparams, i, x, y...)
  newparams[i] = rand(x)
  update_newparams!(newparams, i + 1, y...)
end
@inline function update_newparams!(newparams,i,x)
  newparams[i] = rand(x)
end

particleperturbationkernel(x0, scale) = rand(Uniform(x0 - scale, x0 + scale))

function copyparticle(particle::ParticleSMC)
  return ParticleSMC(copy(particle.params),
copy(particle.weight), copy(particle.scales), copy(particle.distance), [1])
end

function copyparticle(particle::ParticleSMCModel)
  return ParticleSMCModel(copy(particle.params),
copy(particle.weight), copy(particle.scales), copy(particle.model), copy(particle.distance), [1])
end

function copyparticle(particle::ParticleRejection)
  return ParticleRejection(copy(particle.params), copy(particle.distance),
  [1])
end

function copyparticle(particle::ParticleRejectionModel)
  return ParticleRejectionModel(copy(particle.params), copy(particle.model),
  copy(particle.distance), [1])
end

function perturbparticle(particle)

  newparticle = copyparticle(particle)
  newparams = zeros(Float64, length(newparticle.params))

  for i in 1:length(newparams)
    newparams[i] = particleperturbationkernel(newparticle.params[i], newparticle.scales[i])
  end

  newparticle.params = newparams

  return newparticle
end

function kernel_prob(p1, p2)

    prob = 1
    for i in 1:length(p1.params)
      prob = prob * pdf(Uniform(p2.params[i] - p2.scales[i], p2.params[i] + p2.scales[i]), p1.params[i])
    end

    return prob
end

function perturbmodel(ABCsetup, mstar, modelprob)

    prob = ABCsetup.modelkern
    mprob = ones(Float64, length(modelprob))
    mprob[modelprob.==0.0] = 0.0
    nsurvivingmodels = sum(mprob)
    mprob[mprob.> 0.0] = (1 - prob) / (nsurvivingmodels - 1)
    mprob[mstar] = prob

    wsample(1:ABCsetup.nmodels, mprob)
end

function getmodelfreq(particles, ABCsetup)

  freq = zeros(Int64, ABCsetup.nmodels)
  models = map(x -> x.model, particles)
  for i in 1:ABCsetup.nmodels
    freq[i] = sum(models.==i)
  end

  return freq
end

function getmodelprob(currmodel, prevmodel, modelprob, ABCsetup)

  prob = ABCsetup.modelkern

  if currmodel == prevmodel
    return prob
  elseif sum(modelprob.>0.0) > 1
    return (1 - prob) / (sum(modelprob.>0.0) - 1)
  else
    return prob
  end

end

function smcweights(particles, oldparticles, prior)

  weights = zeros(Float64, length(particles))

  for i in 1:length(particles)
    numerator = priorprob(particles[i].params, prior)
    denom = 0.0
    for j in 1:length(particles)
      denom = denom + kernel_prob(particles[i], oldparticles[j])
    end
    weights[i] = numerator / (oldparticles[i].weight * denom)
  end

  weights = weights ./ sum(weights)
  for i in 1:length(particles)
    particles[i].weight = weights[i]
  end

  return particles, weights
end

function denom_modelfunc(p, prevmodelprob, ABCsetup)

  denom_m = 0.0
  for i in 1:ABCsetup.nmodels
    denom_m = denom_m + prevmodelprob[i] * getmodelprob(p.model, i, prevmodelprob, ABCsetup)
  end

  return denom_m
end

function denom_particlefunc(p, particles, prevmodelprob)

  denom_p = 0.0
  for i in particles
    if p.model == i.model
      denom_p = denom_p + ((i.weight * kernel_prob(p, i)) / prevmodelprob[p.model])
    end
  end

  return denom_p
end

function smcweightsmodel(particles, oldparticles, ABCsetup, prevmodelprob)

  numerator = zeros(Float64, length(particles))
  modeldenominator = zeros(Float64, length(particles))
  particledenominator = zeros(Float64, length(particles))

  #calculate numerator
  for i in 1:length(particles)
    numerator[i] = priorprob(particles[i].params, ABCsetup.Models[particles[i].model].prior)
    particledenominator[i] = denom_particlefunc(particles[i], oldparticles, prevmodelprob)
    modeldenominator[i] = denom_modelfunc(particles[i], prevmodelprob, ABCsetup)
  end

  weights = numerator ./ (modeldenominator .* particledenominator)
  weights = weights ./ sum(weights)

  for i in 1:length(particles)
    particles[i].weight = weights[i]
  end

  return particles, weights
end


text:
export
    term_pl02,
    termpt,
    timdef,
    timout,
    tipbod,
    tisbod,
    tkvrsn,
    tparse,
    tpictr,
    trace,
    tsetyr,
    twopi,
    twovec,
    tyear

"""

Find terminator points on a target body. The caller specifies
half-planes, bounded by the illumination source center-target center
vector, in which to search for terminator points.

The terminator can be either umbral or penumbral. The umbral
terminator is the boundary of the region on the target surface
where no light from the source is visible. The penumbral
terminator is the boundary of the region on the target surface
where none of the light from the source is blocked by the target
itself.

The surface of the target body may be represented either by a
triaxial ellipsoid or by topographic data.

### Arguments ###

- `method`: Computation method
- `ilusrc`: Illumination source
- `target`: Name of target body
- `et    `: Epoch in ephemeris seconds past J2000 TDB
- `fixref`: Body-fixed, body-centered target body frame
- `abcorr`: Aberration correction
- `corloc`: Aberration correction locus
- `obsrvr`: Name of observing body
- `refvec`: Reference vector for cutting half-planes
- `rolstp`: Roll angular step for cutting half-planes
- `ncuts `: Number of cutting planes
- `schstp`: Angular step size for searching
- `soltol`: Solution convergence tolerance
- `maxn  `: Maximum number of entries in output arrays

### Output ###

- `npts  `: Counts of terminator points corresponding to cuts
- `points`: Terminator points
- `epochs`: Times associated with terminator points
- `trmvcs`: Terminator vectors emanating from the observer

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/termpt_c.html)
"""
function termpt(method, ilusrc, target, et, fixref, abcorr, corloc, obsrvr, refvec, rolstp,
                ncuts, schstp, soltol, maxn)
    @checkdims 3 refvec
    npts = zeros(SpiceInt, maxn)
    points = Array{SpiceDouble}(undef, 3, maxn)
    epochs = Array{SpiceDouble}(undef, maxn)
    trmvcs = Array{SpiceDouble}(undef, 3, maxn)

    ccall((:termpt_c, libcspice), Cvoid,
          (Cstring, Cstring, Cstring, SpiceDouble, Cstring, Cstring, Cstring, Cstring,
           Ref{SpiceDouble}, SpiceDouble, SpiceInt, SpiceDouble, SpiceDouble, SpiceInt,
           Ref{SpiceInt}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          method, ilusrc, target, et, fixref, abcorr, corloc, obsrvr, refvec, rolstp,
          ncuts, schstp, soltol, maxn, npts, points, epochs, trmvcs)
    handleerror()
    valid_points = npts .>= 1
    npts[valid_points], cmatrix_to_array(points)[valid_points],
        epochs, cmatrix_to_array(trmvcs)[valid_points]
end

@deprecate term_pl02 termpt

"""
    timdef(action, item, value="")

Set and retrieve the defaults associated with calendar input strings.

### Arguments ###

- `action`: The kind of action to take, either `:SET` or `:GET`
- `item`: The default item of interest. The items that may be requested are:
    - `:CALENDAR` with allowed values:
        - `"GREGORIAN"`
        - `"JULIAN"`
        - `"MIXED"`
    - `:SYSTEM` with allowed values:
        - `"TDB"`
        - `"TDT"`
        - `"UTC"`
    - `:ZONE` with allowed values (`0 <= HR < 13` and `0 <= MN < 60`):
        - `"EST"`
        - `"EDT"`
        - `"CST"`
        - `"CDT"`
        - `"MST"`
        - `"MDT"`
        - `"PST"`
        - `"PDT"`
        - `"UTC+\$HR"`
        - `"UTC-\$HR"`
        - `"UTC+\$HR:\$MN"`
        - `"UTC-\$HR:\$MN"`

### Output ###

Returns the value associated with the default item.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/timdef_c.html)
"""
function timdef(action, item, value="")
    lenout = 10
    val = fill(SpiceChar(0), lenout)
    val[1:length(value)] .= collect(value)
    ccall((:timdef_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceInt, Ref{SpiceChar}),
          string(action), string(item), lenout, val)
    handleerror()
    chararray_to_string(val)
end

"""
    timout(et, pictur, lenout=128)

This routine converts an input epoch represented in TDB seconds past the TDB
epoch of J2000 to a character string formatted to the specifications of a
user's format picture.

### Arguments ###

- `et`: An epoch in seconds past the ephemeris epoch J2000
- `pictur`: A format specification for the output string
- `lenout`: The length of the output string plus 1 (default: 128)

### Output ###

Returns a string representation of the input epoch.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/timout_c.html)
"""
function timout(et, pictur, lenout=128)
    string = Array{SpiceChar}(undef, lenout)
    ccall((:timout_c, libcspice), Cvoid,
          (Cdouble, Cstring, SpiceInt, Ref{SpiceChar}),
          et, pictur, lenout, string)
    handleerror()
    chararray_to_string(string)
end

"""
    tipbod(ref, body, et)

Return a 3x3 matrix that transforms positions in inertial coordinates to
positions in body-equator-and-prime-meridian coordinates.

### Arguments ###

- `ref`: Name of inertial reference frame to transform from
- `body`: ID code of body
- `et`: Epoch of transformation

### Output ###

Returns transformation matrix from intertial position to prime meridian.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tipbod_c.html)
"""
function tipbod(ref, body, et)
    tipm = Array{SpiceDouble}(undef, 3, 3)
    ccall((:tipbod_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceDouble, Ref{SpiceDouble}),
          ref, body, et, tipm)
    handleerror()
    permutedims(tipm)
end

"""
    tisbod(ref, body, et)

Return a 6x6 matrix that transforms states in inertial coordinates to
states in body-equator-and-prime-meridian coordinates.

### Arguments ###

- `ref`: Name of inertial reference frame to transform from
- `body`: ID code of body
- `et`: Epoch of transformation

### Output ###

Returns transformation matrix from intertial state to prime meridian.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tisbod_c.html)
"""
function tisbod(ref, body, et)
    tsipm = Array{SpiceDouble}(undef, 6, 6)
    ccall((:tisbod_c, libcspice), Cvoid, (Cstring, SpiceInt, SpiceDouble, Ref{SpiceDouble}),
        ref, body, et, tsipm)
    handleerror()
    permutedims(tsipm)
end

"""
    tkvrsn(item=:TOOLKIT)

Given an item such as the Toolkit or an entry point name, return the latest
version string.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tkvrsn_c.html)
"""
function tkvrsn(item=:TOOLKIT)
    out = ccall((:tkvrsn_c, libcspice), Cstring, (Cstring,), string(item))
    return unsafe_string(out)
end

"""
    tparse(string)

Parse a time string and return seconds past the J2000 epoch on a formal calendar.

### Arguments ###

- `string`: Input time string in UTC

### Output ###

Returns UTC expressed in seconds since J2000.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tparse_c.html)
"""
function tparse(string)
    lenout = 128
    sp2000 = Ref{SpiceDouble}()
    errmsg = Array{SpiceChar}(undef, lenout)
    ccall((:tparse_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceDouble}, Ref{SpiceChar}),
          string, lenout, sp2000, errmsg)
    msg = chararray_to_string(errmsg)
    if !isempty(msg)
        throw(SpiceError(msg))
    end
    sp2000[]
end

"""
    tpictr(sample, lenout=80)

Given a sample time string, create a time format picture suitable for use by
the routine [`timout`](@ref).

### Arguments ###

- `sample`: A sample time string
- `lenout`: The length for the output picture string (default: 80)

### Output ###

Returns a format picture that describes sample.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tpictr_c.html)
"""
function tpictr(sample, lenout=80)
    lenerr = 128
    pictur = Array{SpiceChar}(undef, lenout)
    ok = Ref{SpiceBoolean}()
    errmsg = Array{SpiceChar}(undef, lenerr)
    ccall((:tpictr_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceBoolean}, Ref{SpiceChar}),
          sample, lenout, lenerr, pictur, ok, errmsg)
    if !Bool(ok[])
        throw(SpiceError(chararray_to_string(errmsg)))
    end
    chararray_to_string(pictur)
end

function _trace(matrix)
    ccall((:trace_c, libcspice), SpiceDouble, (Ref{SpiceDouble},), matrix)
end

"""
    trace(matrix)

!!! warning "Deprecated"
    Use `LinearAlgebra.tr(matrix)` instead.
"""
trace

@deprecate trace LinearAlgebra.tr

"""
    tsetyr(year)

Set the lower bound on the 100 year range.

### Arguments ###

-`year`: Lower bound on the 100 year interval of expansion

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tsetyr_c.html)
"""
function tsetyr(year)
    ccall((:tsetyr_c, libcspice), Cvoid, (SpiceInt,), year)
end

function _twopi()
    ccall((:twopi_c, libcspice), SpiceDouble, ())
end

"""
    twopi()

!!! warning "Deprecated"
    Use `2π` instead.
"""
twopi

@deprecate twopi() 2π

"""
    twovec(axdef, indexa, plndef, indexp)

Find the transformation to the right-handed frame having a given vector as a
specified axis and having a second given vector lying in a specified coordinate
plane.

### Arguments ###

- `axdef`: Vector defining a principal axis
- `indexa`: Principal axis number of axdef (X=1, Y=2, Z=3)
- `plndef`: Vector defining (with axdef) a principal plane
- `indexp`: Second axis number (with indexa) of principal plane

### Output ###

Returns output rotation matrix.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/twovec_c.html)
"""
function twovec(axdef, indexa, plndef, indexp)
    @checkdims 3 axdef plndef
    mout = Array{SpiceDouble}(undef, 3, 3)
    ccall((:twovec_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, Ref{SpiceDouble}, SpiceInt, Ref{SpiceDouble}),
          axdef, indexa, plndef, indexp, mout)
    handleerror()
    permutedims(mout)
end

"""
    tyear()

Returns the number of seconds per tropical year.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/tyear_c.html)
"""
function tyear()
    ccall((:tyear_c, libcspice), SpiceDouble, ())
end


text:
### MCSamplerState represents the Monte Carlo samplers' internal state ("local variables")

abstract MCSamplerState{F<:VariateForm}

abstract MHSamplerState{F<:VariateForm} <: MCSamplerState{F}

abstract HMCSamplerState{F<:VariateForm} <: MCSamplerState{F}

abstract LMCSamplerState{F<:VariateForm} <: MCSamplerState{F}

### Abstract Monte Carlo sampler typesystem

## Root Monte Carlo sampler

abstract MCSampler

## Family of samplers based on Metropolis-Hastings

abstract MHSampler <: MCSampler

## Family of Hamiltonian Monte Carlo samplers

abstract HMCSampler <: MCSampler

## Family of Langevin Monte Carlo samplers

abstract LMCSampler <: MCSampler

function initialize_step!{F<:VariateForm}(
  sstate::HMCSamplerState{F},
  parameter::Parameter{Continuous, F},
  sampler::HMCSampler,
  tuner::DualAveragingMCTuner
)
  sstate.oldhamiltonian = hamiltonian(sstate.pstate.logtarget, sstate.momentum)

  leapfrog!(sstate, parameter.gradlogtarget!)
  sstate.newhamiltonian = hamiltonian(sstate.pstate.logtarget, sstate.momentum)

  sstate.ratio = sstate.newhamiltonian-sstate.oldhamiltonian
  sstate.a = 2*(exp(sstate.ratio) > 0.5)-1

  while exp(sstate.ratio)^sstate.a > 2^(-sstate.a)
    sstate.tune.step = (2^sstate.a)*sstate.tune.step
    leapfrog!(sstate, parameter.gradlogtarget!)
    sstate.newhamiltonian = hamiltonian(sstate.pstate.logtarget, sstate.momentum)

    sstate.ratio = sstate.newhamiltonian-sstate.oldhamiltonian
  end
end

function initialize_step!(
  sstate::HMCSamplerState{Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::HMCSampler,
  tuner::DualAveragingMCTuner,
  momentum::Real
)
  sstate.momentum = momentum
  initialize_step!(sstate, parameter, sampler, tuner)
end

function initialize_step!(
  sstate::HMCSamplerState{Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::HMCSampler,
  tuner::DualAveragingMCTuner,
  momentum::RealVector
)
  sstate.momentum[:] = momentum
  initialize_step!(sstate, parameter, sampler, tuner)
end

tuner_state(parameter::Parameter, sampler::MCSampler, tuner::VanillaMCTuner) = BasicMCTune(NaN, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::MHSampler, tuner::VanillaMCTuner) = BasicMCTune(1., 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::HMCSampler, tuner::VanillaMCTuner) =
  BasicMCTune(sampler.leapstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::LMCSampler, tuner::VanillaMCTuner) =
  BasicMCTune(sampler.driftstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::MHSampler, tuner::AcceptanceRateMCTuner) = BasicMCTune(1., 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::HMCSampler, tuner::AcceptanceRateMCTuner) =
  BasicMCTune(sampler.leapstep, 0, 0, tuner.period)

tuner_state(parameter::Parameter, sampler::LMCSampler, tuner::AcceptanceRateMCTuner) =
  BasicMCTune(sampler.driftstep, 0, 0, tuner.period)

reset!(tune::BasicMCTune, sampler::MCSampler, tuner::VanillaMCTuner) =
  ((tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN))

function reset!(tune::BasicMCTune, sampler::HMCSampler, tuner::AcceptanceRateMCTuner)
  tune.step = sampler.leapstep
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

function reset!(tune::BasicMCTune, sampler::LMCSampler, tuner::AcceptanceRateMCTuner)
  tune.step = sampler.driftstep
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

reset!{S<:ValueSupport, F<:VariateForm}(
  sstate::MCSamplerState{F},
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MCSampler,
  tuner::MCTuner
) =
  reset!(sstate.tune, sampler, tuner)

hamiltonian(logtarget::Real, momentum::Real) = logtarget-0.5*abs2(momentum)

hamiltonian(logtarget::Real, momentum::RealVector) = logtarget-0.5*dot(momentum, momentum)

function leapfrog!(
  pstate::ParameterState{Continuous, Univariate},
  pstate0::ParameterState{Continuous, Univariate},
  momentum0::Real,
  step::Real,
  gradlogtarget!::Function
)
  local momentum::Real

  momentum = momentum0+0.5*step*pstate0.gradlogtarget
  pstate.value = pstate0.value+step*momentum
  gradlogtarget!(pstate)
  momentum = momentum+0.5*step*pstate.gradlogtarget

  return momentum
end

function leapfrog!(
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  pstate0::ParameterState{Continuous, Multivariate},
  momentum0::RealVector,
  step::Real,
  gradlogtarget!::Function
)
  momentum[:] = momentum0+0.5*step*pstate0.gradlogtarget
  pstate.value[:] = pstate0.value+step*momentum
  gradlogtarget!(pstate)
  momentum[:] = momentum+0.5*step*pstate.gradlogtarget
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/PCRE2_jll.jl
baremodule PCRE2_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libpcre2_8

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libpcre2_8_handle = C_NULL
libpcre2_8_path = ""

if Sys.iswindows()
    const libpcre2_8 = "libpcre2-8-0.dll"
elseif Sys.isapple()
    const libpcre2_8 = "@rpath/libpcre2-8.0.dylib"
else
    const libpcre2_8 = "libpcre2-8.so.0"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libpcre2_8_handle = dlopen(libpcre2_8)
    global libpcre2_8_path = dlpath(libpcre2_8_handle)
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libpcre2_8_path() = libpcre2_8_path

end  # module PCRE2_jll

max_stars_repo_name: dmageeLANL/julia

max_issues_repo_name: dmageeLANL/julia

max_forks_repo_name: dmageeLANL/julia


text:
module BlockOpt

import Base: copy

import LinearAlgebra:
    diagm,
    BLAS,
    Symmetric,
    mul!,
    norm,
    dot,
    issymmetric,
    rank,
    qr,
    I,
    pinv,
    eigvals

import Statistics: mean

import Pkg: activate

import ForwardDiff: ForwardDiff, Dual, jacobian

import TRS: trs_small, trs

import Weave

import RecipesBase: RecipesBase, @series, @recipe, @userplot

import Printf: @sprintf, @printf

import Dates: Date, Second, now, Minute

export optimize

# Types
export Model, Driver

# See model.jl
export model,
    name,
    obj,
    grad!,
    grad,
    objective,
    gradient,
    initial_iterate,
    formula,
    dimension,
    directory,
    objective!,
    gradient!,
    initial_iterate!,
    formula!,
    final,
    hess_sample,
    hessAD

# See options.jl
export DriverOptions,
    samples,
    Δ_max,
    δ_tol,
    ϵ_tol,
    max_iterations,
    weave_level,
    log_level,
    samples!,
    Δ_max!,
    δ_tol!,
    ϵ_tol!,
    max_iterations!,
    weave_level!,
    log_level!,
    LogLevel,
    INFO,
    DEBUG,
    WARN,
    ERROR,
    WeaveLevel,
    NONE,
    ALL

# See driver.jl
export driver,
    orth,
    pflag,
    QN_update,
    SR1,
    PSB,
    S_update,
    S_update_a,
    S_update_b,
    S_update_c,
    S_update_d,
    S_update_e,
    S_update_f

# See trace.jl
export trace,
    Δt,
    evaluations,
    trs_timer,
    trs_counter,
    ghs_timer,
    ghs_counter,
    f_vals,
    ∇f_norms,
    Δ_vals,
    p_norms,
    ρ_vals

# See recipes.jl
export rhotrace,
    rhotrace!,
    steptrace,
    steptrace!,
    radiustrace,
    radiustrace!,
    objtrace,
    objtrace!,
    gradtrace,
    gradtrace!

# See simulation.jl
export weave

include("util.jl")

include("model.jl")

include("options.jl")

include("driver.jl")

include("lib/trace.jl")

include("lib/backend.jl")

include("simulation.jl")

"""
    optimize(model::Model, driver::Driver)

An entry-point into the minimization iteration with the given `model` subject to the specified `driver`.
"""
optimize(model::Model, driver::Driver) = optimize!(Simulation(model, driver))

"""
    optimize(f::Function, ∇f!::Function, x₀::AbstractArray)

Attempts to determine the unconstrained minimum of `f` via a first-order method with the initial iterate
given by `x₀`. The gradient ∇f! must be specified as an inplace operation.
"""
optimize(f::Function, ∇f!::Function, x₀::AbstractArray) = optimize(
    Model("Missing"; objective = f, gradient = ∇f!, initial_iterate = x₀),
    Driver(),
)

include("lib/show.jl")

include("lib/recipes.jl")

const BLAS_THREADS = Ref{Int}(1)

blas_threads!(threads::Int) = (BLAS_THREADS[] = threads) > 1 && BLAS.set_num_threads(threads)

function __init__()
    blas_threads!(Threads.nthreads())
end

end # module


text:
using DifferentialEquations
using DiffEqCallbacks
using Distributions
using StatsBase
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

@enum InfectionStatus Susceptible Infected Recovered

Base.zero(::Type{InfectionStatus}) = Infected

function sir_abm!(du,u,p,t)
    (β,c,γ,δt) = p
    N = length(u)
    # Initialize du to u
    for i in 1:N
        du[i] = u[i]
    end
    for i in 1:N # loop through agents
        # If recovered
        if u[i]==Recovered
            continue
        # If susceptible
        elseif u[i]==Susceptible
            ncontacts = rand(Poisson(c*δt))
            while ncontacts > 0
                j = sample(1:N)
                if j==i
                    continue
                end
                a = u[j]
                if a==Infected && rand() < β
                    du[i] = Infected
                    break
                end
                ncontacts -= 1
            end
        # If infected
        else u[i]==Infected
            if rand() < γ
                du[i] = Recovered
            end
        end
    end
    nothing
end;

δt = 0.1
tf = 40.0
tspan = (0.0,tf);

β = 0.05
c = 10.0
γ = rate_to_proportion(0.25,δt)
p = [β,c,γ,δt];

N = 1000
I0 = 10
u0 = Array{InfectionStatus}(undef,N)
for i in 1:N
    if i <= I0
        s = Infected
    else
        s = Susceptible
    end
    u0[i] = s
end

Random.seed!(1234);

susceptible(u) = count(i == Susceptible for i in u)
infected(u) = count(i == Infected for i in u)
recovered(u) = count(i == Recovered for i in u);

saved_values = SavedValues(Float64, Tuple{Int64,Int64,Int64})
cb = SavingCallback((u,t,integrator)->(susceptible(u),infected(u),recovered(u)),
    saved_values,
    saveat=0:δt:tf);

prob_abm = DiscreteProblem(sir_abm!,u0,tspan,p);

sol_abm = solve(prob_abm,
    solver = FunctionMap(),
    dt = δt,
    callback = cb,
    dense = false,
    save_on = false);

df_abm = DataFrame(saved_values.saveval)
rename!(df_abm,[:S,:I,:R])
df_abm[!,:t] = saved_values.t;

@df df_abm plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_abm,
    solver=FunctionMap,
    dt=δt,
    callback=cb,
    dense=false,
    save_on=false)


text:
#=

This file specifies how to pretty-print Scattering module types

=#

# Overload the show method for ScatteringMatrix
function Base.show(io::IO, x::ScatteringMatrix)

    println(io, "  f₁₁| $(length(x.f₁₁))-element vector")
    println(io, "  f₁₂| $(length(x.f₁₂))-element vector")
    println(io, "  f₂₂| $(length(x.f₂₂))-element vector")
    println(io, "  f₃₃| $(length(x.f₃₃))-element vector")
    println(io, "  f₃₄| $(length(x.f₃₄))-element vector")
    print(io, "  f₄₄| $(length(x.f₄₄))-element vector")
end


text:
include("../src/MsgPackRpcClient.jl")

using MsgPackRpcClient

session = MsgPackRpcClientSession.create()
MsgPackRpcClientSession.create_sock(session, "localhost", 5000)
MsgPackRpcClientSession.create_sock(session, "localhost", 5001)
MsgPackRpcClientSession.create_sock(session, "localhost", 5002)

the_number_of_call = 10000
results = {}

@sync begin
@async for i in 1:the_number_of_call
  future  = call(session, "hello"; sync = false, sock = session.socks.pool[session.ptr])
  push!(results, get(future))
  MsgPackRpcClientSession.rotate(session)
end
end

@sync begin
@async for r in results
  println(r)
end
end

session.destroy()

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
rndict = Dict(1 => 3, 2 => 3,
              3 => 1, 4 => 1,
              5 => 2, 6 => 2,
              7 => 3, 8 => 3,
              9 => 1,
              99 => 0)

struct OpCode
    instr::Int
    mode::Array{Int, 1}
    N::Int
end

function OpCode(c::Int)
    digs = digits(c)
    oc = digs[1] + (length(digs) > 1 ? 10*digs[2] : 0)
    N = rndict[oc]
    mode = zeros(Int, N)
    mode[1:length(digs[3:end])] = digs[3:end]
    return OpCode(oc, mode, N)
end

function intcode_computer(program::Array{Int,1}, input::Channel{Int}, output::Channel{Union{Int, Symbol}}; pc::Int=1)

    Np = length(program)
    Np_new = 10*Np
    for j=1:(Np_new-Np)
        push!(program, 0)
    end

    jump_flag = false
    rel_base = 0

    function load(addr::Int, mode::Int)::Int
        #println("Load: $(addr) -- $(mode)")
        if mode == 0
            return program[addr+1]
        elseif mode == 1
            return addr
        elseif mode == 2
            return program[addr + rel_base + 1]
        else
            error("Unknown mode $(mode).")
        end
    end

    function store(data::Int, addr::Int, mode::Int)
        #println("Store: $(data) --> $(addr) -- $(mode)")
        if mode == 0
            program[addr+1] = data
        elseif mode == 1
            error("Invalid mode 1 for store instruction.")
        elseif mode == 2
            program[addr + rel_base + 1] = data
        else
            error("Unknown mode $(mode).")
        end
    end

    function load_input(dst::Int, mode::Int)
        #println("Waiting for next input: store to $(dst) -- $(mode)")
        inp = take!(input)
        store(inp, dst, mode)
    end

    function store_output(data::Int)
        #println("Tried to store output $(data)")
        put!(output, data)
    end

    function jtrue(x::Int, y::Int)
        if x != 0
            pc = y+1
            jump_flag = true
        end
    end

    function jfalse(x::Int, y::Int)
        if x == 0
            pc = y+1
            jump_flag = true
        end
    end

    lt(x::Int, y::Int)::Int = Int(x < y)

    eq(x::Int, y::Int)::Int = Int(x == y)

    while program[pc] != 99
        jump_flag = false
        oc = OpCode(program[pc])
        #println(oc)
        if oc.instr == 1
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            store(x+y, program[pc+3], oc.mode[3])

        elseif oc.instr == 2
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            store(x*y, program[pc+3], oc.mode[3])

        elseif oc.instr == 3
            load_input(program[pc+1], oc.mode[1])

        elseif oc.instr == 4
            x = load(program[pc+1], oc.mode[1])
            store_output(x)

        elseif oc.instr == 5
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            jtrue(x, y)
        elseif oc.instr == 6
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            jfalse(x, y)

        elseif oc.instr == 7
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            store(lt(x,y), program[pc+3], oc.mode[3])

        elseif oc.instr == 8
            x = load(program[pc+1], oc.mode[1])
            y = load(program[pc+2], oc.mode[2])
            store(eq(x,y), program[pc+3], oc.mode[3])

        elseif oc.instr == 9
            rel_base += load(program[pc+1], oc.mode[1])

        else
            @error "Something went wrong! Got code $(oc.instr)."
        end

        pc += (jump_flag ? 0 : oc.N+1)
    end
    println("Computer finished!")
    put!(output, :done);
end


text:
# scatter plotting with pyplot.
using PyPlot
using Random
x,y = randn(100),randn(100)
scatter(x,y)
display(gcf())


text:
"""
```julia
HMCSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=10)
```

Draw samples from the true posterior via Hamiltonian Monte Carlo.

**Keywords arguments**
    - `ϵ::T` : convergence criteria
    - `nBurnin::Int` : Number of samples discarded before starting to save samples
    - `samplefrequency::Int` : Frequency of sampling
"""
mutable struct HMCSampling{T<:Real} <: SamplingInference{T}
    nBurnin::Integer # Number of burnin samples
    samplefrequency::Integer # Frequency at which samples are saved
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 # Number of samples
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    sample_store::AbstractArray{T,3}
    x::SubArray{T,2,Matrix{T}}#,Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
    y::LatentArray{SubArray}
    function HMCSampling{T}(nBurnin::Int,samplefrequency::Int,ϵ::T,nIter::Integer,Stochastic::Bool,nSamples::Integer,nSamplesUsed::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool) where T
        return new{T}(nBurnin,samplefrequency,ϵ,nIter,Stochastic,nSamples,nSamplesUsed,MBIndices,ρ,flag)
    end
end

function HMCSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=10) where {T<:Real}
    HMCSampling{Float64}(nBurnin,samplefrequency,ϵ,0,false,1,1,[1],1.0,true)
end

function Base.show(io::IO,inference::HMCSampling{T}) where {T<:Real}
    print(io,"Gibbs Sampler")
end

function grad_log_joint_pdf(model::MCGP{T,L,HMCSampling{T}},x) where {T,L}
    grad_log_likelihood(model.likelihood,get_y(model),x) + sum(grad_log_gp_prior.(model.f,x))
end

function init_sampler(inference::HMCSampling{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer) where {T<:Real}
    inference.nSamples =
    inference.MBIndices = collect(1:nSamples)
    inference.sample_store = zeros(T,nLatent,nFeatures,n)
    return inference
end

function sample_parameters!(model::MCGP{T,L,HMCSampling{T}},nSamples::Int,callback) where {T,L}
    f_init = vcat(get_f(model)...)
    logpdf(x) = log_joint_model(model,x)
    gradlogpdf(x) = (log_joint_model(model,x),grad_log_joint_model(model,x))
    metric = i.metric
    h = Hamiltonian(metric, logpdf, gradlogpdf)
    int = Leapfrog(find_good_eps(h, f_init))
    prop = NUTS{MultinomialTS,GeneralisedNoUTurn}(int)
    adaptor = StanHMCAdaptor(n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, int.ϵ))
    samples, stats = sample(h, prop, f_init, nSamples, adaptor, n_adapts; progress=true)
end


text:
"""
    itemtype(ci::CodeInfo, item)

Infer a concrete DataType from item. Item might be slotnumber or ssavalue etc.
"""
itemtype(ci::CodeInfo, item) = typeof(item)
itemtype(ci::CodeInfo, item::DataType) = item
itemtype(ci::CodeInfo, item::Number) = typeof(item)
itemtype(ci::CodeInfo, item::SlotNumber) = ci.slottypes[item.id]
itemtype(ci::CodeInfo, item::SSAValue) = ci.ssavaluetypes[item.id]
itemtype(ci::CodeInfo, item::TypedSlot) = itemtype(ci, SlotNumber(item.id))
itemtype(ci::CodeInfo, item::Compiler.Const) = itemtype(ci, item.val)

hasitemtype(ci::CodeInfo, item, type::DataType) = itemtype(ci, item) <: type
hasitemtype(ci::CodeInfo, item, types::Array{DataType,1}) = any([itemtype(ci, item) <: type for type in types])
function hasitemtype(ci::CodeInfo, items::Array, type::DataType)
    if hasname(items[1], :(getindex)) && isa(items[2], SlotNumber)
        datatype = ci.slottypes[items[2].id]
        return datatype.parameters[1] <: type
    else
        # return any([itemtype(ci, item) <: type for item in items])
        return any([hasitemtype(ci, item, type) for item in items])
    end
end

"""
    argtypes!(ci::CodeInfo, argtypes::Dict, funcs::Dict, items::Array)

Infer argtypes and update argtypes if items[1] is in funcs.
"""
argtypes!(ci::CodeInfo, argtypes::Dict, funcs::Dict, item) = nothing
function argtypes!(ci::CodeInfo, argtypes::Dict, funcs::Dict, items::Array)
    if isa(items[1], GlobalRef) && items[1].name in keys(funcs)
        argtypes[Symbol(items[1].name)] = itemtype.((ci,), items[2:end])
    end
    argtypes!.((ci,), (argtypes,), (funcs,), items)
    return nothing
end


text:
# Pie charts

function PlotPieChart(label_ids, values::AbstractArray{T}, count::Integer, x::Real, y::Real, radius::Real, normalize::Bool = false, label_fmt::String = "%.1f", angle::Real = 90.) where {T<:ImPlotData}
    ImPlot_PlotPieChart(label_ids, values, count, x, y, radius, normalize, label_fmt, angle)
end

function PlotPieChart(label_ids, values::AbstractArray{T}, count::Integer, x::Real, y::Real, radius::Real, normalize::Bool = false, label_fmt::String = "%.1f", angle::Real = 90.) where {T<:Real}
    ImPlot_PlotPieChart(label_ids, Float64.(values), count, x, y, radius, normalize, label_fmt, angle)
end

function PlotPieChart(values::AbstractArray{T}, count::Integer, x, y, radius;
                      normalize::Bool = false, label_fmt::String = "%.1f", angle0 = 90.0,
                      label_ids::Vector{String} = ["" for _ in 1:length(values)]) where {T<:ImPlotData}

    ImPlot_PlotPieChart(label_ids, values, count, x, y, radius, normalize, label_fmt, angle0)
end

function PlotPieChart(values::AbstractArray{T}, count::Integer, x, y, radius;
                      normalize::Bool = false, label_fmt::String = "%.1f", angle0 = 90.0,
                      label_ids::Vector{String} = ["" for _ in 1:length(values)]) where {T<:Real}

    ImPlot_PlotPieChart(label_ids, Float64.(values), count, x, y, radius, normalize,
                            label_fmt, angle0)
end


text:
function englishiness(ptarray::Array{UInt8, 1})
    EnglishFrequencies = Dict( # from http://scottbryce.com/cryptograms/stats.htm
      'e'=>12.51,
      't'=>9.25,
      'a'=>8.04,
      'o'=>7.6,
      'i'=>7.26,
      'n'=>7.09,
      's'=>6.54,
      'r'=>6.12,
      'h'=>5.49,
      'l'=>4.14,
      'd'=>3.99,
      'c'=>3.06,
      'u'=>2.71,
      'm'=>2.53,
      'f'=>2.3,
      'p'=>2.0,
      'g'=>1.96,
      'w'=>1.92,
      'y'=>1.73,
      'b'=>1.54,
      'v'=>0.99,
      'k'=>0.67,
      'x'=>0.19,
      'j'=>0.16,
      'q'=>0.11,
      'z'=>0.09,
    )
    s=0
    for ll in ptarray
        llchar = Char(ll)
        if isalpha(llchar)
            s+=get(EnglishFrequencies, lowercase(llchar), -5)
        else
            s+=-30
        end
    end
    return s
end

function genBestSingleByteXORDecrypt(ctbytes::Array{UInt8, 1})
    bestdecrypt = []
    bestscore = -Inf
    bestkey = 0x00
    for xx in 0x00:0xff
        if isascii(Char(xx))
            ptbytes = ctbytes $ xx
            ptscore = englishiness(ptbytes)
            if ptscore > bestscore
                bestscore = ptscore
                bestdecrypt = ptbytes
                bestkey = xx
            end
        end
    end
    return (bestscore, bestdecrypt, bestkey)
end

# from challenge 5
function stringToBytes(X::String)
    return [UInt8(X[ii]) for ii in 1:length(X)]
end;

function repeatToLength(a::Array{UInt8,1}, len::Int64)
    if len < length(a)
        return a[1:len]
    else
        b = repeat(a, outer=convert(Int64, ceil(len / length(a))))
        return b[1:len]
    end
end

function repeatingKeyXOR(pt::Array{UInt8,1}, key::Array{UInt8,1})
    return pt $ repeatToLength(key, length(pt))
end

function repeatingKeyXOR(pt::String, key::String)
    return repeatingKeyXOR(stringToBytes(pt), stringToBytes(key))
end

# from challenge 6
function hammingDistance(a::Array{UInt8,1}, b::Array{UInt8,1})
    differentbits = a $ b
    d = 0
    for ii in 1:length(differentbits)
        c = bits(differentbits[ii])
        for jj in 1:8
            if c[jj] == '1'
                d+=1
            end
        end
    end
    return d
end

function hammingDistance(a::String, b::String)
    return hammingDistance(stringToBytes(a), stringToBytes(b))
end;

# from challenge 7
function decryptAES128ECBFileWithKey(fname, key)
    decrypt = chomp(readstring(`openssl enc -d -a -aes-128-ecb -K $key -in $fname`))
    return decrypt
end

# from challenge 9
  function padToNBytes(a::Array{UInt8,1}, nbytes)
      newa = Array{UInt8,1}(nbytes)
      ntopad = nbytes - length(a)
      newa[1:length(a)] = a
      newa[(length(a)+1):nbytes] = UInt8(ntopad)
      return newa
  end

  function padToNBytes(a::String, nbytes)
      return String(padToNBytes(stringToBytes(a), nbytes))
  end;

# assumes ciphertxt is base64encoded, and key is in hex
function decryptAES128ECBWithKey(ciphertext::String, key::String)
    decrypt = chomp(
        readstring(
            pipeline(`echo -n $ciphertext`,
                     `base64 --decode`,
                     `openssl enc -nopad -d -aes-128-ecb -K $key`,
                     `base64`)))
    return decrypt
end

# takes bytes and returns bytes
function decryptAES128ECBWithKey(cipherbytes::Array{UInt8,1}, keybytes::Array{UInt8,1})
    ciphertext = base64encode(cipherbytes)
    keytext = bytes2hex(keybytes)
    return base64decode(decryptAES128ECBWithKey(ciphertext, keytext))
end

function encryptAES128ECBWithKey(plaintext::String, key::String)
    encrypt = chomp(
        readstring(
            pipeline(`echo -n $plaintext`,
                     `base64 --decode`,
                     `openssl enc -nopad -e -aes-128-ecb -K $key`,
                     `base64`)))
    return encrypt
end

function encryptAES128ECBWithKey(plaintextbytes::Array{UInt8,1}, keybytes::Array{UInt8,1})
    plaintext = base64encode(plaintextbytes)
    keytext = bytes2hex(keybytes)
    return base64decode(encryptAES128ECBWithKey(plaintext, keytext))
end

function encryptCBCByHand(ptbytes::Array{UInt8,1}, blocksize, ivbytes::Array{UInt8,1}, keybytes::Array{UInt8,1})
    cipherbytes = Array{UInt8,1}()
    nblocks = convert(Int64, ceil(length(ptbytes)/blocksize))
    ptbytes = padToNBytes(ptbytes, blocksize*nblocks)
    lastblock = ivbytes
    for ii in 1:nblocks
        bytestoencrypt = ptbytes[((ii-1)*blocksize+1):(ii*blocksize)] $ lastblock
        encryptedbytes = encryptAES128ECBWithKey(bytestoencrypt, keybytes)
        append!(cipherbytes, encryptedbytes)
        lastblock = encryptedbytes
    end
    return cipherbytes
end;

function decryptCBCByHand(ctbytes::Array{UInt8,1}, blocksize, ivbytes::Array{UInt8,1}, keybytes::Array{UInt8,1})
    ptbytes = Array{UInt8,1}()
    nblocks = convert(Int64, length(ctbytes)/blocksize)
    lastblock = ivbytes
    for ii in 1:nblocks
        bytestodecrypt = ctbytes[((ii-1)*blocksize+1):(ii*blocksize)]
        decryptedbytes = decryptAES128ECBWithKey(bytestodecrypt, keybytes) $ lastblock
        append!(ptbytes, decryptedbytes)
        lastblock = bytestodecrypt
    end
    return ptbytes
end;

function isThisECBOrCBC(ctbytes::Array{UInt8,1}, blocksize=16)
    nblocks = convert(Int64, length(ctbytes)/blocksize)
    interblocksimilarity = []
    for ii in 1:(nblocks-1)
        for jj in (ii+1):nblocks
            blockii = ctbytes[((ii-1)*blocksize+1):(ii*blocksize)]
            blockjj = ctbytes[((jj-1)*blocksize+1):(jj*blocksize)]
            samebits = ~(blockii $ blockjj)
            similarity = 0
            for cc in samebits
                ccbits = bits(cc)
                for dd in 1:length(ccbits)
                    if ccbits[dd]=='1'
                        similarity+=1
                    end
                end
            end
            push!(interblocksimilarity, similarity/(8*blocksize))
        end
    end
    pexactsame = sum(interblocksimilarity.==1)/length(interblocksimilarity)
    cipher = "cbc"
    if pexactsame > 0.5
        cipher = "ecb"
    end
    return (cipher, pexactsame, interblocksimilarity)
end;


text:
using JuMP
import Ipopt
import NLopt

function P1(z)
    model = Model(Ipopt.Optimizer)
    JuMP.set_silent(model)
    JuMP.@variable(model, x1)
    JuMP.@variable(model, x2)
    JuMP.@NLexpression(model, nl_expr, z[1] * ((x1 - 3)^2 + (x2 - 1)^2) + z[2] * ((x1 - 1)^2 + (x2 - 1)^2))

    JuMP.@NLobjective(model, Min, nl_expr)
    JuMP.@constraint(model, c1, x1 + 2 * x2 <= 2)
    JuMP.@constraint(model, c2, x1 - 2 * x2 <= 2)
    JuMP.@constraint(model, c3, -x1 + 2 * x2 <= 2)
    JuMP.@constraint(model, c4, -x1 - 2 * x2 <= 2)
    JuMP.optimize!(model)
    return value.([x1, x2]), JuMP.objective_value(model)
end
function P2(v)
    model = Model(Ipopt.Optimizer)
    JuMP.set_silent(model)

    JuMP.@variable(model, z)
    JuMP.@variable(model, x1)
    JuMP.@variable(model, x2)
    JuMP.@objective(model, Min, z)
    JuMP.@constraint(model, c1, x1 + 2 * x2 <= 2)
    JuMP.@constraint(model, c2, x1 - 2 * x2 <= 2)
    JuMP.@constraint(model, c3, -x1 + 2 * x2 <= 2)
    JuMP.@constraint(model, c4, -x1 - 2 * x2 <= 2)
    # Gammax = Γ([x1, x2])
    # z[1] * ((x1 - 3)^2 + (x2 - 1)^2) + z[2] * ((x1 - 1)^2 + (x2 - 1)^2)
    JuMP.@NLconstraint(model, c5, (x1 - 3)^2 + (x2 - 1)^2 - z - v[1] <= 0)
    JuMP.@NLconstraint(model, c6, (x1 - 1)^2 + (x2 - 1)^2 - z - v[2] <= 0)
    JuMP.optimize!(model)
    return (value.([x1, x2]), value(z))
end
function dual_objective(u1, u2, u3, u4, w1, w2)
    model2 = Model(Ipopt.Optimizer)
    JuMP.set_silent(model2)
    JuMP.@variable(model2, x1)
    JuMP.@variable(model2, x2)
    JuMP.@NLexpression(model2, my_expr, w1 * (x1 - 3)^2 + (x2 - 1)^2 + w2 * (x1 - 1)^2 + (x2 - 1)^2)
    JuMP.@NLobjective(model2, Min,
        u1 * (x1 + 2 * x2 - 2) +
        u2 * (x1 - 2 * x2 - 2) +
        u3 * (-x1 + 2 * x2 - 2) +
        u4 * (-x1 - 2 * x2 - 2) +
        my_expr
    )
    JuMP.optimize!(model2)
    return JuMP.objective_value(model2)
end
# function D2(v)
#     model = Model(Ipopt.Optimizer)
#     @variable(model, u >= 0)
#     @variable(model, w1 >= 0)
#     @variable(model, w2 >= 0)
#     @NLobjective(
#         model,
#         Max,
#         (-(w1^2 + w2^2) / (4 * u) - u) + w1 * (1 - v[1]) + w2 * (1 - v[2])
#     )
#     @constraint(model, c1, w1 + w2 == 1)
#     optimize!(model)
#     return (value.([w1, w2]), value(u))
# end

# function D2(v)
#     model = Model(Ipopt.Optimizer)
#     @variable(model, u1 >= 0)
#     @variable(model, u2 >= 0)
#     @variable(model, u3 >= 0)
#     @variable(model, u4 >= 0)
#     @variable(model, w1 >= 0)
#     @variable(model, w2 >= 0)

#     @NLexpression(model, my_expr,
#         u3 * (u2 + u4 - u1 - u3 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
#         u1 * (u2 + u4 - u1 - u3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
#         u4 * (-4 + u1 + u3 - u2 - u4 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
#         u2 * (-4 + u1 + u3 - u2 - u4 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
#         w1 * (-3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 +
#         (1 / 2) * (-u1 + u2 - u3 + u4)^2 +
#         w2 * (-1 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2)

#     @NLobjective(
#         model,
#         Max,
#         my_expr - w1 * v[1] - w2 * v[2]
#     )
#     @constraint(model, c1, w1 + w2 == 1)
#     optimize!(model)
#     return (value.([w1, w2]), value.([u1, u2, u3, u4]))
# end

function D2(v)
    model = Model(Ipopt.Optimizer)
    @variable(model, u1 >= 0)
    @variable(model, u2 >= 0)
    @variable(model, u3 >= 0)
    @variable(model, u4 >= 0)
    @variable(model, w1 >= 0)
    @variable(model, w2 >= 0)

    @NLexpression(model, my_expr,
        u3 * (u2 + u4 - u1 - u3 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
        u1 * (u2 + u4 - u1 - u3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
        u4 * (-4 + u1 + u3 - u2 - u4 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
        u2 * (-4 + u1 + u3 - u2 - u4 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
        w1 * (-3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 +
        (1 / 2) * (-u1 + u2 - u3 + u4)^2 +
        w2 * (-1 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2)

    @NLobjective(
        model,
        Max,
        my_expr - w1 * v[1] - w2 * v[2]
    )
    @constraint(model, c1, w1 + w2 == 1)
    optimize!(model)
    return (value.([w1, w2]), value.([u1, u2, u3, u4]))
end
# function D2para(v)
#     model = Model(Ipopt.Optimizer)
#     @variable(model, u >= 0)
#     @variable(model, w1 >= 0)
#     @variable(model, w2 >= 0)
#     exprrr = ParameterD2()
#     @NLexpression(model,my_expr_obj,exprrr)
#     @NLobjective(
#         model,
#         Max,
#         exprrr - w1*v[1]-w2*v[2]
#     )
#     @constraint(model, c1, w1 + w2 == 1)
#     optimize!(model)
#     return (value.([w1, w2]), value(u))
# end

# function D2(v0::Vector)

#     function ExplicitFunction(uw::Vector, grad::Vector, v::Vector)
#         println("called my own OFunction")
#         if length(grad) > 0
#             grad[1] = 0
#             grad[2] = 0.5 / sqrt(x[2])
#         end
#         w1, w2, u1, u2, u3, u4 = uw
#         return u3 * (u2 + u4 - u1 - u3 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
#                u1 * (u2 + u4 - u1 - u3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
#                u4 * (-4 + u1 + u3 - u2 - u4 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
#                u2 * (-4 + u1 + u3 - u2 - u4 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
#                w1 * (-3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 +
#                (1 / 2) * (-u1 + u2 - u3 + u4)^2 +
#                w2 * (-1 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 - w1 * v[1] - w2 * v[2]
#     end
#     function ImplicitFunction(uw::Vector, grad::Vector, v::Vector)

#         if length(grad) > 0
#             grad[1] = 0
#             grad[2] = 0.5 / sqrt(x[2])
#         end
#         model2 = Model(Ipopt.Optimizer)
#         JuMP.set_silent(model2)
#         JuMP.@variable(model2, x1)
#         JuMP.@variable(model2, x2)
#         JuMP.@NLexpression(model2, my_expr, uw[1] * (x1 - 3)^2 + (x2 - 1)^2 + uw[2] * (x1 - 1)^2 + (x2 - 1)^2)
#         JuMP.@NLobjective(model2, Min,
#             uw[3] * (x1 + 2 * x2 - 2) +
#             uw[4] * (x1 - 2 * x2 - 2) +
#             uw[5] * (-x1 + 2 * x2 - 2) +
#             uw[6] * (-x1 - 2 * x2 - 2) +
#             my_expr
#         )
#         JuMP.optimize!(model2)
#         return JuMP.objective_value(model2) - uw[1] * v[1] - uw[2] * v[2]
#     end

#     function constraint(uw::Vector, _grad::Vector)
#         uw[1] + uw[2] - 1
#     end

#     opt = NLopt.Opt(:LN_COBYLA, 6) # LN_COBYLA LD_MMA
#     opt.lower_bounds = repeat([0], 6)
#     # opt.xtol_rel = 1e-2

#     opt.max_objective = (x, g) -> ImplicitFunction(x, g, v0)
#     NLopt.equality_constraint!(opt, constraint, 1e-8)

#     (minf, minx, ret) = NLopt.optimize(opt, repeat([0.9], 6))
#     numevals = opt.numevals # the number of function evaluations
#     w1, w2, u1, u2, u3, u4 = minx
#     return (value.([w1, w2]), value.([u1, u2, u3, u4]))
# end

function Γ(x)
    [(x[1] - 3)^2 + (x[2] - 1)^2; (x[1] - 1)^2 + (x[2] - 1)^2]
end


text:
@testset "I393" begin
    s = """
    type MyType{Int64}
    end
    """
    msgs = lintstr(s)
    @test msgs[1].code == :I393
    @test msgs[1].variable == "Int64"
    @test contains(msgs[1].message, "using an existing type as type parameter name is probably a typo")

    @test messageset(lintstr("""
    type MyType{Int64} <: Float64
    end
    """)) == Set([:I393])
end

@testset "E513" begin
    msgs = lintstr("""
    type MyType{T<:Int}
    end
    """)
    @test messageset(msgs) == Set([:E513])
    @test msgs[1].variable == "T <: Int"
    @test contains(msgs[1].message, "leaf type as a type constraint makes no sense")

    msgs = lintstr("""
    type MyType{T<:Int, Int<:Real}
    end
    """)
    @test messageset(msgs) == Set([:E513, :E538])
    @test msgs[1].variable == "T <: Int"
    @test contains(msgs[1].message, "leaf type as a type constraint makes no sense")
    @test msgs[2].variable == "Int"
    @test contains(msgs[2].message, "known type in parametric data type, use {T<:...}")
end

@testset "E538" begin
    msgs = lintstr("""
    type MyType{Int<:Real}
    end
    """)
    @test messageset(msgs) == Set([:E538])
    @test contains(msgs[1].message, "known type in parametric data type, use {T<:...}")

    msgs = lintstr("""
    type SomeType
    end
    type MyType{SomeType<:Real}
    end
    """)
    @test messageset(msgs) == Set([:E538])
    @test contains(msgs[1].message, "known type in parametric data type, use {T<:...}")
end

# TODO: this inner constructor syntax is deprecated
@testset "Types" begin
    @test_broken isempty(lintstr("""
    type MyType{T}
        t::T
        MyType(x::T) = new(x)
    end
    """))

    @test_broken isempty(lintstr("""
    type MyType{T<:Integer}
        t::T
        MyType(x) = new(convert(T, x))
    end
    """))

    s = """
    type MyType <: Integer
        t::Int
        function MyTypo()
            new(1)
        end
    end
    """
    msgs = lintstr(s)
    @test_broken msgs[1].code == :E517
    @test_broken msgs[1].variable == "MyTypo"
    @test_broken contains(msgs[1].message, "constructor-like function name doesn't match type MyType")

    @test_broken isempty(lintstr("""
    type MyType{T}
        b::T
        MyType{S}(y::S) = new(convert(T,y))
    end
    """))
end

s = """
type MyType
    t::Int
    NotATypo() = 1
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
@compat abstract type SomeAbsType end
@compat abstract type SomeAbsNum <: Number end
@compat abstract type SomeAbsVec1{T} <: Array{T,1} end
@compat abstract type SomeAbsVec2{T} end
@compat primitive type MyBitsType 8 end

type MyType{T<:SomeAbsType}
    t::T
    MyType(x) = new(convert(T, x))
end
"""
@test_broken isempty(lintstr(s))

s = """
type MyType{T<:Integer}
    t::T
    function MyType(x, someFunc::Function)
        o = new(convert(T, x))
        finalizer(o, someFunc) # forgot to return o
    end
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :E611
@test_broken contains(msgs[1].message, "constructor doesn't seem to return the constructed object")

s = """
type MyType{T<:Integer}
    t::T
    function MyType(x, someFunc::Function)
        o = new(convert(T, x))
        finalizer(o, someFunc) # didn't forget to return o
        return o
    end
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)

s = """
type MyType
    a
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I691
@test msgs[1].variable == "a"
@test contains(msgs[1].message, "a type is not given to the field which can be slow")

s = """
type MyType
    @lintpragma("Ignore untyped field a")
    a
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
type MyType
    a::Array{Float64}
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I692
@test msgs[1].variable == "a"
@test contains(msgs[1].message, "array field has no dimension which can be slow")

s = """
type MyType
    @lintpragma("Ignore dimensionless array field a")
    a::Array{Float64}
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
type MyType
    lintpragma("Ignore dimensionless array field a")
    a::Array{Float64}
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E425
@test contains(msgs[1].message, "use lintpragma macro inside type declaration")

s = """
@compat primitive type 8 a end
"""
msgs = lintstr(s)
@test msgs[1].code == :E101
@test contains(msgs[1].message, "this expression must be a Symbol")

s = """
type MyType
    a::Int
    b::Int
    MyType(x::Int,y::Int) = new(x,y)
    MyType(x::Int) = MyType(x, 0)
    function MyType()
        v = new(0, 0)
        v[2] = convert(Int, rand()* 5.0) # assume we'd define setindex! somewhere
        v
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
type MyType
    a::Int
    b::Int
    function(x)
        new(x)
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E417
@test contains(msgs[1].message, "anonymous function inside type definition")
@test_broken msgs[2].code == :E517
@test_broken msgs[2].variable == ""
@test_broken contains(msgs[2].message, "constructor-like function name doesn't match type MyType")

s = """
type MyType
    a::Int
    b::Int
    MyType(x) = new(x)
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :I671
@test_broken contains(msgs[1].message, "new is provided with fewer arguments than fields")

@testset "Inner Constructors" begin
    @test isempty(lintstr("""
    type MyType{T}
        b::T
        (::Type{MyType}){T}(x::T) = new{T}(x)
    end
    """))

    @test isempty(lintstr("""
    type MyType{T}
        b::T
        (::Type{MyType}){T<:Integer}(x::T) = new{T}(x)
    end
    """))
end

s = """
type MyType
    a::Int
    b::Int
    MyType(x) = new(x, 0, 0)
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :E435
@test_broken contains(msgs[1].message, "new is provided with more arguments than fields")

s = """
type MyType
    a::Int
    b::Int
    MyType() = new()
end
"""
msgs = lintstr(s)
@test isempty(msgs) # ok

s = """
type MyType
    a::Int
    b::Int
    function MyType(x)
        @lintpragma("Ignore short new argument")
        new(x)
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

# TODO: this inner constructor syntax is deprecated
s = """
type MyType{T}
    b::T
    MyType(x::T) = new(x)
    MyType(x::Int) = MyType{T}(convert(T, x))
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)

s = """
type myType{T}
    b::T
    myType(x::T) = new(x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I771
@test contains(msgs[1].message, "type names should start with an upper case")


text:
struct ScoreParams
    word_groups::Vector{Vector{String}}
    group_scores::Vector{Float64}
    required_groups::Vector{Int}
    gap_penalty::Float64
    normalize::Bool
    word_to_group::Dict{String, Int}

    function ScoreParams(word_groups, group_scores, required_groups, gap_penalty; normalize=false)
        word_to_group = Dict(word => i for (i, words) in enumerate(word_groups) for word in words)
        return new(word_groups, group_scores, required_groups, gap_penalty, normalize, word_to_group)
    end
end

function ScoreParams(pairs::AbstractVector{<:Pair}, gap_penalty::Real; normalize=false)
    word_groups = Vector{String}[]
    group_scores = Float64[]
    required_groups = Int[]

    for (i, (words, score)) in enumerate(pairs)
        push!(word_groups, words)
        if score == Inf
            push!(group_scores, 0.)
            push!(required_groups, i)
        else
            push!(group_scores, score)
        end
    end

    return ScoreParams(word_groups, group_scores, required_groups, gap_penalty; normalize=normalize)
end

function score(components, params::ScoreParams)
    group_seen = falses(length(params.word_groups))
    s = 0.

    for word in components
        g = get(params.word_to_group, word, 0)
        if g > 0
            if !group_seen[g]
                s += params.group_scores[g]
                group_seen[g] = true
            end
        else
            s -= params.gap_penalty
        end
    end

    params.normalize && (s /= length(components))

    return all(group_seen[i] for i in params.required_groups) ? s : -Inf
end

score(a::Acronym, params::ScoreParams) = score(a.components, params)

function score_word(word::String, params::ScoreParams)
    g = params.word_to_group[word]
    return g == 0 ? -params.gap_penalty : params.group_scores[g]
end

struct ScoredAcronymFamily
    family::AcronymFamily
    params::ScoreParams
    best::Acronym
    best_score::Float64

    ScoredAcronymFamily(family, params, best) = new(family, params, best, score(best, params))
    ScoredAcronymFamily(family, params) = new(family, params, best_acronym(family, params)...)
end


text:
function linearindex(shape, idxs)
   stride = 1
   idx = 1
   @inbounds for d in 1:min(length(shape), length(idxs))
      i = idxs[d]
      n = shape[d]
      @assert 1<=i<=n "1<=$i<=$n   shape=$shape idxs=$idxs"
      idx += stride*(i-1)
      stride *= n
   end
   return idx
end


text:
@inline function (-){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C})
   Flex{negate(S), C}(-a.hi, -a.lo)
end

for T in (:CLCL, :OPOP)
   @eval begin
       function (-){Q<:Qualia,C<:Clay}(a::Flex{$T,Q,C}, b::Flex{$T,Q,C})
           lo = (-)(a.lo, b.hi, RoundDown)
           hi = (-)(a.hi, b.lo, RoundUp)
           Flex{$T,Q,C}(lo,hi)
       end
   end
end

function (-){S<:Sculpt, T<:Sculpt,Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C})
    aLoIsOpen, aHiIsOpen = boundries(S)
    bLoIsOpen, bHiIsOpen = boundries(T)
    sculpting = boundries(aLoIsOpen|bLoIsOpen, aHiIsOpen|bHiIsOpen)

    lo = (-)(a.lo, b.hi, RoundDown)
    hi = (-)(a.hi, b.lo, RoundUp)

    Flex{sculpting,Q,C}(lo,hi)
end

(-){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::C) = (-)(a, Flex{S,Q,C}(b))
(-){S<:Sculpt,Q<:Qualia,C<:Clay}(a::C, b::Flex{S,Q,C}) = (-)(Flex{S,Q,C}(a), b)


text:
using Random
using Distributions
using FunctionWrappers

const Age=UInt8

include("params/households.jl")
include("params/friendship.jl")
include("params/progression.jl")
include("params/hospital.jl")
include("params/phonetracking.jl")
include("params/spreading.jl")

abstract type InfectionModulation end

struct SimParams
  household_ptrs::Vector{Tuple{PersonIdx,PersonIdx}}  # (i1,i2) where i1 and i2 are the indices of first and last member of the household

  ages::Vector{Age}
  genders::BitVector

  progressions::Vector{Progression} # not sure if progressions should be there

  hospital_kernel_params::Union{Nothing, HospitalInfectionParams}  # nothing if hospital kernel not active
  friendship_kernel_params::Union{Nothing, FriendshipKernelParams}

  constant_kernel_param::Float64
  household_kernel_param::Float64

  hospital_detections::Bool
  mild_detection_prob::Float64

  backward_tracking_prob::Float32
  backward_detection_delay::TimeDiff

  forward_tracking_prob::Float32
  forward_detection_delay::TimeDiff

  quarantine_length::Float32
  testing_time::TimeDiff

  phone_tracking_params::Union{Nothing, PhoneTrackingParams}

  infection_modulation_function::Union{Nothing,InfectionModulation}
  spreading_params::Union{Nothing, SpreadingParams}
end

include("params/modulations.jl")

numindividuals(params::SimParams) = length(params.household_ptrs)
progressionof(params::SimParams, person_id::Integer) = params.progressions[person_id]
severityof(params::SimParams, person_id::Integer) = progressionof(params, person_id).severity
householdof(params::SimParams, person_id::Integer) = UnitRange(params.household_ptrs[person_id]...)
age(params::SimParams, person_id::Integer) = params.ages[person_id]
gender(params::SimParams, person_id::Integer) = params.genders[person_id]

socialcompetence(params::SimParams, person_id::Integer) =
  nothing!=params.friendship_kernel_params && socialcompetence(params.friendship_kernel_params, person_id)

ishealthcare(params::SimParams, person_id::Integer) =
  nothing!=params.hospital_kernel_params && ishealthcare(params.hospital_kernel_params, person_id)
uses_phone_tracking(params::SimParams, person_id::Integer) =
  nothing!=params.phone_tracking_params && uses_phone_tracking(params.phone_tracking_params, person_id)

spreading(params::SimParams, person_id::Integer) = isnothing(params.spreading_params) ? 1.0 : spreading(params.spreading_params, person_id)

function load_params(rng=MersenneTwister(0);
        population::Union{AbstractString,DataFrame},
        infection_modulation_name::Union{Nothing,AbstractString}=nothing,
        infection_modulation_params::NamedTuple=NamedTuple{}(),
        kwargs...
        )

  individuals_df::DataFrame = isa(population, AbstractString) ? load_individuals(population) : population

  num_individuals = individuals_df |> nrow

  dist_incubation_time = LogNormal(1.3669786931887833, 0.5045104580676582)
  dist_symptom_onset_time = Gamma(0.8738003969079596, 2.9148873266517685)
  dist_hospitalization_time = Gamma(1.1765988120148885, 2.6664347368236787)
  dist_mild_recovery_time = Uniform(11, 17)
  dist_severe_recovery_time = Uniform(4*7, 8*7)
  dist_death_time = LogNormal(2.610727719719777, 0.44476420066780653)

  progressions = Vector{Simulation.Progression}(undef, num_individuals);
  resample!(rng, progressions, individuals_df.age,
    dist_incubation_time,
    dist_symptom_onset_time,
    dist_hospitalization_time,
    dist_mild_recovery_time,
    dist_severe_recovery_time,
    dist_death_time
  )

  infection_modulation_function = isnothing(infection_modulation_name) ? nothing : make_infection_modulation(infection_modulation_name; infection_modulation_params...)

  make_params(
    rng,
    individuals_df=individuals_df,
    progressions=progressions,
    infection_modulation_function=infection_modulation_function;
    kwargs...
  )
end

function make_params(
  rng::AbstractRNG=MersenneTwister(0);
  individuals_df::DataFrame,
  progressions::AbstractArray{Progression},

  infection_modulation_function=nothing,

  hospital_kernel_param::Float64=0.0,
  healthcare_detection_prob::Float64=0.8,
  healthcare_detection_delay::Float64=1.0,

  constant_kernel_param::Float64=1.0,
  household_kernel_param::Float64=1.0,
  friendship_kernel_param::Float64=0.0,

  hospital_detections::Bool=true,
  mild_detection_prob::Float64=0.0,

  backward_tracking_prob::Float64=1.0,
  backward_detection_delay::Float64=1.0,

  forward_tracking_prob::Float64=1.0,
  forward_detection_delay::Float64=1.0,

  quarantine_length::Float64=14.0,
  testing_time::Float64=1.0,

  phone_tracking_usage::Real=0.0,
  phone_detection_delay::Real=0.25,
  phone_tracking_usage_by_household::Bool=false,

  spreading_alpha::Union{Nothing,Real}=nothing,
  spreading_x0::Real=1,
  spreading_truncation::Real=Inf
)
  sort!(individuals_df, :household_index)

  num_individuals = individuals_df |> nrow

  @assert num_individuals == length(progressions)

  household_ptrs = make_household_ptrs(individuals_df.household_index)

  friendship_kernel_params =  if 0 == friendship_kernel_param; nothing
                        elseif 0.0 < friendship_kernel_param
                          FriendshipKernelParams(
                            friendship_kernel_param,
                            Age.(individuals_df.age),
                            individuals_df.gender,
                            individuals_df.social_competence)
                        else error("bad condition for friendship kernel")
                        end

  hospital_kernel_params =  if 0 == hospital_kernel_param; nothing
                            elseif 0.0 < hospital_kernel_param;
                              HospitalInfectionParams(
                                individuals_df.ishealthcare,
                                (UInt32(1):UInt32(num_individuals))[individuals_df.ishealthcare],
                                hospital_kernel_param,
                                healthcare_detection_prob,
                                healthcare_detection_delay
                              )
                            else error("hospital_kernel_param must be postive or 0, got $hospital_kernel_param")
                            end

  phone_tracking_params = if 0 == phone_tracking_usage; nothing
                      elseif 0.0 < phone_tracking_usage <= 1.0
                        phone_tracking_usage_by_household ?
                          PhoneTrackingParams(rng, num_individuals, phone_tracking_usage, phone_detection_delay, 1, household_ptrs) :
                          PhoneTrackingParams(rng, num_individuals, phone_tracking_usage, phone_detection_delay, 1)
                      else error("tracking_app_usage must be nonnegative, got $phone_tracking_usage")
                      end

  spreading_params = if nothing == spreading_alpha; nothing
                  elseif 0.0 < spreading_alpha
                    SpreadingParams(rng, num_individuals, alpha=spreading_alpha, x0=spreading_x0, truncation=spreading_truncation)
                  else error("spreading_alpha must be larger than 0, got $spreading_alpha")
                  end

  params = SimParams(
    household_ptrs,
    individuals_df.age,
    individuals_df.gender,

    progressions,

    hospital_kernel_params,
    friendship_kernel_params,

    constant_kernel_param,
    household_kernel_param,

    hospital_detections,
    mild_detection_prob,

    backward_tracking_prob,
    backward_detection_delay,

    forward_tracking_prob,
    forward_detection_delay,

    quarantine_length, # quarantine length
    testing_time, # testing time

    phone_tracking_params,
    infection_modulation_function,
    spreading_params
  )
  params
end

function saveparams(dict, p::SimParams)
  dict["num_individuals"] = numindividuals(p)
  dict["constant/kernel_param"] = p.constant_kernel_param
  dict["household/kernel_param"] = p.household_kernel_param
  dict["quarantine/duration"] = p.quarantine_length
  dict["detections/hospital"] = p.hospital_detections
  dict["detections/mild"] = p.mild_detection_prob

  dict["infection_modulation"] = p.infection_modulation_function

  dict["tracking/backward_prob"] = p.backward_tracking_prob
  dict["tracking/backward_detection_delay"] = p.backward_detection_delay
  dict["tracking/forward_tracking_prob"] = p.forward_tracking_prob
  dict["tracking/forward_detection_delay"] = p.forward_detection_delay
  dict["tracking/testing_time"] = p.testing_time

  saveparams(dict, p.progressions, "progressions/")
  nothing==p.hospital_kernel_params || saveparams(dict, p.hospital_kernel_params, "hospital/")
  nothing==p.phone_tracking_params || saveparams(dict, p.phone_tracking_params, "phone_tracking/")
  nothing==p.spreading_params || saveparams(dict, p.spreading_params, "spreading")
  nothing
end


text:
# Extended Gtk functions

function gtk_widget_add_tick_callback(widget::GtkWidget, callback, userdata, notify)
    @ccall libgtk.gtk_widget_add_tick_callback(
        widget::Ptr{GObject},
        callback::Ptr{Nothing},
        userdata::Ptr{Nothing},
        notify::Ptr{Nothing}
    )::Cuint
end

function gtk_widget_remove_tick_callback(widget::GtkWidget, id::Cuint)
    @ccall libgtk.gtk_widget_remove_tick_callback(widget::Ptr{GObject}, id::Cuint)::Nothing
end

function gtk_button_new_from_icon_name(name::String, size::Cuint)
    @ccall libgtk.gtk_button_new_from_icon_name(name::Ptr{Cchar}, size::Cuint)::Ptr{GObject}
end

function icon_button(name::String, size::Cuint = GtkIconSize.BUTTON)
    return Gtk.GtkButtonLeaf(gtk_button_new_from_icon_name(name, size))
end

gdk_keyval_name(keyval::Cuint) = @ccall libgdk.gdk_keyval_name(keyval::Cuint)::Ptr{UInt8}
keyval_name(keyval::Cuint) = Gtk.bytestring(gdk_keyval_name(keyval))

function gtk_grid_attach_next_to(grid::GtkGrid, child::GtkWidget, sibling::GtkWidget, side::Symbol, width::Int, height::Int)
    side = getfield(Gtk.GtkPositionType, side)
    @ccall libgtk.gtk_grid_attach_next_to(
        grid::Ptr{GObject},
        child::Ptr{GObject},
        sibling::Ptr{GObject},
        side::Cint,
        width::Cint,
        height::Cint
    )::Cvoid
end


text:
module ContextExprRules

using ExprRules
using StatsBase

include("utils.jl")
include("rulenode_operators.jl")
include("context.jl")
include("constraints.jl")
include("propagators.jl")
include("grammar.jl")
include("sampling.jl")

export	containedin,

	change_expr,
	swap_node,
	get_rulesequence,
	rulesoftype,
	rulesonleft,

	GrammarContext,
	addparent!,

	Constraint,
	ValidatorConstraint,
	PropagatorConstraint,
	ComesAfter,
	Ordered,

	propagate,

	ContextGrammar,
	@cgrammar,
	addconstraint!,
	ContextExpressionIterator,
	propagate_contraints

end # module


text:
function work_on_subproblems!(subworker::SubWorker{T,A,S},
                              work::Work,
                              progress::ProgressQueue{T},
                              x̄::RunningAverage{A},
                              δ::RunningAverage{T},
                              decisions::Decisions{A},
                              r::IteratedValue{T}) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{T,A,S}} = fetch(subworker)
    if isempty(subproblems)
       # Workers has nothing do to, return.
       return
    end
    while true
        t::Int = try
            wait(work)
            take!(work)
        catch err
            if err isa InvalidStateException
                # Master closed the work channel. Worker finished
                return
            end
        end
        if t == -1
            # Worker finished
            return
        end
        ξ::A = fetch(decisions, t)
        if t > 1
            update_subproblems!(subproblems, ξ, fetch(r,t-1))
        end
        @sync for (i,subproblem) ∈ enumerate(subproblems)
            @async begin
                take!(δ, i)
                take!(x̄, i)
                put!(δ, i, norm(subproblem.x - ξ, 2)^2, subproblem.π)
                reformulate_subproblem!(subproblem, ξ, fetch(r,t))
                Q::T = subproblem()
                put!(x̄, i, subproblem.π)
                put!(progress, (t,subproblem.id,Q))
            end
        end
    end
end


text:
prog = "1 + 1"

ex1 = Meta.parse(prog)
ex1

typeof(ex1)

ex1.head

ex1.args

ex2 = Expr(:call, :+, 1, 1)

ex1 == ex2

dump(ex2)

ex3 = Meta.parse("(4 + 4) / 2")

Meta.show_sexpr(ex3)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
abstract type AbstractSolver end
abstract type AbstractLinearSolver <: AbstractSolver end
abstract type AbstractNonLinearSolver <: AbstractSolver end

abstract type AbstractDirectSolver <: AbstractLinearSolver end
abstract type AbstractIterativeSolver <: AbstractLinearSolver end

##### PRECONDITIONERS #####
# abstract type AbstractPreconditioner end
#
# struct ILU <: AbstractPreconditioner
# 	drop_tol::Float64
# end
# (p::ILU)(J) = ilu(J, τ = p.drop_tol)

# Default backslash solver #
struct BackslashSolver <: AbstractDirectSolver end # single core direct solver

(l::BackslashSolver)(K,RHS,model) = K\RHS

# MUMPS multithreaded sparse direct solver #
struct MUMPS <: AbstractDirectSolver end  # TODO

# GMRES from IterativeSolvers package

@kwdef mutable struct GMRESIterativeSolver{T<:Real, Tl, Tr} <: AbstractIterativeSolver
	"Tolerance for solver"
	tol::T = 1e-4

	"Number of restarts"
	restart::Int64 = 200

	"Maximum number of iterations"
	maxiter::Int64 = 100

	"Display information during iterations"
	verbose::Bool = false

	"Record information"
	log::Bool = true

	"Start with zero guess"
	initially_zero::Bool = true

	"Left preconditioner"
	Pl::Tl = IterativeSolvers.Identity()

	"Right preconditioner"
	Pr::Tr = IterativeSolvers.Identity()
end

function (l::GMRESIterativeSolver{T, Tl, Tr})(J, rhs, model) where {T, Tl, Tr}
	#J_map = v -> apply(J, v)
	#dim = ndofs(model.dofhandler)
	#Jmap = LinearMaps.LinearMap{T}(J_map, dim, dim ; ismutating = false)
	(l.Pl isa AbstractPreconditioner) && (l.Pl = l.Pl(J))
	res = IterativeSolvers.gmres(J, rhs; tol = l.tol, log = l.log, verbose = l.verbose, restart = l.restart, maxiter = l.maxiter, initially_zero = l.initially_zero, Pl = l.Pl, Pr = l.Pr)
	(res[2].iters >= l.maxiter) && (@warn "IterativeSolvers.gmres iterated maxIter = $(res[2].iters) times without achieving the desired tolerance.\n")
	return res[1]
end

@kwdef struct ConjugateGradients <: AbstractIterativeSolver
    package::Symbol
    kwargs::NamedTuple = NamedTuple()
end

@kwdef struct MINRES <: AbstractIterativeSolver
    package::Symbol
    kwargs::NamedTuple = NamedTuple()
end

@kwdef struct NewtonRaphson{T<:AbstractLinearSolver,R1<:Real,R2<:Real,R3<:Real} <: AbstractNonLinearSolver
	alpha::R1 = 1.0
    atol::R2 = 1e-4 # absolute residual norm tolerance
    max_iter_number::Int = 10 # maximum non linear iterations
	max_iter_atol::R3 = 0.0 # absolute residual norm tolerance at max iteration
    max_div_iter::Int = 4 # maximum diverging iterations before restart
    autodiff::Bool = false # Has the jacobian to be computed numerically
    iter_on_nonassociated_plasticity::Bool = false
    linear_solver::T = BackslashSolver()
end

struct NLsolveNewton{LS,R1<:Real,R2<:Real} <: AbstractNonLinearSolver
    xtol::R1
	ftol::R2
	iterations::Int
	store_trace::Bool
	show_trace::Bool
	extended_trace::Bool
	linesearch::LS

	# Custom inner constructor used to import the relevant packages and perform some checks
	function NLsolveNewton{LS,R1,R2}(xtol,ftol,iterations,store_trace,show_trace,extended_trace,linesearch) where {LS,R1<:Real,R2<:Real}
		# load modules
		@info("loading NLsolve package")
        @eval(Rheologies, import NLsolve)
		@info("loading LineSearches package")
        @eval(Rheologies, import LineSearches)

		# Some checks
		@assert xtol >= 0
		@assert ftol >= 0
		return new{LS,R1,R2}(xtol,ftol,iterations,store_trace,show_trace,extended_trace,linesearch)
	end
end

# general outer constructor without type parameters
NLsolveNewton(xtol::R1,ftol::R2,iterations,store_trace,show_trace,extended_trace,linesearch::LS) where {LS,R1<:Real,R2<:Real} = NLsolveNewton{LS,R1,R2}(xtol,ftol,iterations,store_trace,show_trace,extended_trace,linesearch)
# kwargs method
NLsolveNewton(;xtol = 0.0,ftol = 1e-8,iterations = 1000,store_trace = false,show_trace = false,extended_trace = false,linesearch = LineSearches.Static) = NLsolveNewton(xtol,ftol,iterations,store_trace,show_trace,extended_trace,linesearch)


text:
#!/usr/bin/env julia

data = ARGS
@assert length(data) >= 1

for row in data
    println(row,"\t",length(row))
end


text:
using Distributions, GLM, DataFrames, Random, LinearAlgebra
Random.seed!(0)

n = 100
beta0, beta1, beta2, beta3 = 10, 30, 60, 90
sig, sigX = 25, 5
etaVals = [200.0, 100.0 , 50.0, 10.0, 1.0, 0.1, 0.01, 0.0]

function createDataFrame(eta)
    Random.seed!(1)
    x1 = round.(collect(1:n) + sigX*rand(Normal(),n),digits = 5)
    x2 = round.(collect(1:n) + sigX*rand(Normal(),n),digits = 5)
    x3 = round.(x1 +2*x2 + eta*rand(Normal(),n),digits = 5)
    y = beta0 .+ beta1*x1 + beta2*x2 + beta3*x3 + sig*rand(Normal(),n)
    return DataFrame(Y = y, X1 = x1, X2 = x2, X3 = x3)
end

VIF3() = 1/(1-r2(lm(@formula(X3 ~ X1 + X2),df)))

for eta in etaVals
    print("eta = $(eta): ")
    df = createDataFrame(eta)
    glmOK = true
    try
        global model = lm(@formula(Y ~ X1 + X2 + X3),df)
    catch err
        println("Exception with GLM: ", err)
        glmOK = false
    end

    if glmOK
        covMat = vcov(model)
        sigVec = sqrt.(diag(covMat))
        corrmat = round.(covMat ./ (sigVec*sigVec'),digits=5)
        println("Corr(X1,X3) = ", corrmat[2,4],
                ",\t Corr(X2,X3) = ",corrmat[3,4],
                ",\t VIF3 = ", round.(VIF3(),digits=2) )
    else
         println("\t In this case we may use SVD or ridge regression if needed.")
    end
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook