text:
using DecFP
is_ch_type(::Val{:Decimal32})  = true
is_ch_type(::Val{:Decimal64})  = true
is_ch_type(::Val{:Decimal128})  = true
is_ch_type(::Val{:Decimal})  = true

result_type(::Val{:Decimal32}, scale)  = Vector{Dec32}
result_type(::Val{:Decimal64}, scale)  = Vector{Dec64}
result_type(::Val{:Decimal128}, scale)  = Vector{Dec128}

function dec_type_by_precision(precision_str::String)
    precision = parse(Int, precision_str)

    precision in 1:9 && return :Decimal32
    precision in 10:18 && return :Decimal64
    precision in 19:38 && return :Decimal128
    error("Decimal error: unsupported precision $(precision)")
end
function result_type(::Val{:Decimal}, precision_str, scale)
    return result_type(
            Val(dec_type_by_precision(precision_str)),
            scale
        )
end

function read_decimal(::Type{DecT}, ::Type{IntT},
                     sock::ClickHouseSock, num_rows, scale_str) where {DecT, IntT}
    scale = parse(Int, scale_str)
    data = chread(sock, Vector{IntT}, num_rows)
    return DecT.(data, -scale)
end

function write_decimal(::Type{DecT}, ::Type{IntT},
    sock::ClickHouseSock, data, scale_str) where {DecT, IntT}
    scale = parse(Int, scale_str)
    tmp = Vector{IntT}(undef, length(data))
    for (i,v) in enumerate(data)
        (sign, value, exp) = sigexp(convert(DecT, v))

        (exp != -scale) &&
            error("Decimal:Wrong exponent in input data, expected $(scale) got $(exp)")

        tmp[i] = IntT(sign * value)
    end
    chwrite(sock, tmp)
end

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
                                            ::Val{:Decimal32}, scale_str)
    return read_decimal(Dec32, Int32, sock, num_rows, scale_str)
end

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
    ::Val{:Decimal64}, scale_str)
    return read_decimal(Dec64, Int64, sock, num_rows, scale_str)
end

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
    ::Val{:Decimal128}, scale_str)
    return read_decimal(Dec128, Int128, sock, num_rows, scale_str)
end

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
    ::Val{:Decimal}, precision_str,  scale_str)
    return read_col_data(sock, num_rows,
    Val(dec_type_by_precision(precision_str)) ,scale_str)
end

function write_col_data(sock::ClickHouseSock,
                    data, ::Val{:Decimal32}, scale_str)

    return write_decimal(Dec32, Int32, sock, data, scale_str)
end

function write_col_data(sock::ClickHouseSock,
    data, ::Val{:Decimal64}, scale_str)

    return write_decimal(Dec64, Int64, sock, data, scale_str)
end

function write_col_data(sock::ClickHouseSock,
    data, ::Val{:Decimal128}, scale_str)

    return write_decimal(Dec128, Int128, sock, data, scale_str)
end

function write_col_data(sock::ClickHouseSock, data,
    ::Val{:Decimal}, precision_str,  scale_str)
    return write_col_data(sock, data,
    Val(dec_type_by_precision(precision_str)) ,scale_str)
end


text:
# # Statistics

## load packages
using PlotPlants
#------------------------------------------------------------------------------

# ## Linear regression
# ### General linear model
# [`line_regress`](@ref) returns a [`LinearRegressionResult`](@ref) struct,
#     inside which are `lr` (GLM linear regression result), `inter`
#     (intercept), `slope`, `inter_p` (P value for intercept), `slope_p` (P
#     value for slope), `inter_ci` (intercept confidence internal), `slope_ci`
#     (slope confidence interval), and `df` (dataframe for plotting purpose).
#     By default, `df` was sorted from low X to high X. The `lowerY` and
#     `upperY` are useful to plot the confidence intervals as in
#     [`plot_line_regress!`](@ref).
xs = rand(10) .* 10;
ys = xs .* 0.05 .+ rand(10) .+ 1;
lf = line_regress(xs, ys);
@show lf.df;
lf = line_regress(xs, ys; sorting=false);
@show lf.df;
#------------------------------------------------------------------------------

# ### Set intercept = 0
lf = line_regress(xs, ys; intercept=false);
@show lf.inter;
@show lf.inter_p;
@show lf.inter_ci;
#------------------------------------------------------------------------------

# ### Is slope different from a given value
# By default [`test_slope`](@ref) returns the p value of whether the slope is
#     different from 0. Yet, you may also customize whether you want to compare
#     the slope to a given value, e.g., 1. Also, you may set the intercept as
#     0.
@show test_slope(xs, ys);
@show test_slope(xs, ys; slope=1);
@show test_slope(xs, ys; intercept=false);
#------------------------------------------------------------------------------

# ## Calculate the scatter density
# ### General calculation
# Calculating the density of large dataset can be extremely slow, and thus
#     avoiding repeating the calculation is important. Function
#     [`calculate_density`](@ref) calculates the density and returns a
#     dataframe that store the density information (sorted from low to high).
#     If total number of values exceeds 1000, a progress bar will be shown. Use
#     this function with [`plot_density!`](@ref).
N  = 100;
xs = [rand(N); rand(N).* 0.5; rand(N).* 0.2; rand(N).* 0.1];
ys = [rand(N); rand(N).* 0.5; rand(N).* 0.2; rand(N).* 0.1];
calculate_density(xs,ys);
N  = 1001;
xs = [rand(N); rand(N).* 0.5; rand(N).* 0.2; rand(N).* 0.1];
ys = [rand(N); rand(N).* 0.5; rand(N).* 0.2; rand(N).* 0.1];
calculate_density(xs,ys);
#------------------------------------------------------------------------------

# ### Multi-threading calculation
## the density will be calculated using 4 threads
calculate_density(xs,ys,4);
#------------------------------------------------------------------------------

ext: jl


text:
using LinearAlgebra

# For now just centered at 0
freezing_time(x, θ, μ, F) = freezing_time(x, θ, F)
function freezing_time(x, θ, μ, F::Union{Boomerang, Boomerang1d})
    if μ == 0
        if θ*x >= 0.0
            return π - atan(x/θ)
        else
            return atan(-x/θ)
        end
    else
        u = x^2 - 2μ*x + θ^2
        u < 0 && return Inf
        t1 = mod(2atan((sqrt(u) - θ)/(2μ - x)), 2pi)
        t2 = mod(-2atan((sqrt(u) + θ)/(2μ - x)), 2pi)
        x == 0 && return max(t1, t2)
        return min(t1, t2)
    end
end

function freezing_time!(tfrez, t, x, θ, f, Z::Union{BouncyParticle, Boomerang})
    for i in eachindex(x)
        if f[i]
            tfrez[i] = t + freezing_time(x[i], θ[i], Z.μ[i], Z)
        end
    end
    tfrez
end

function refresh_sticky_vel!(θ, θf, f, F::Union{BouncyParticle, Boomerang})
    ρ̄ = sqrt(1-F.ρ^2)
    for i in eachindex(θ)
        if f[i]
            θ[i] = F.ρ*θ[i] +  ρ̄*randn()
        else
            θf[i] = abs( F.ρ*θf[i] +  ρ̄*randn())*sign(θf[i])
        end
    end
    θ, θf
end

function subnormsq(∇ϕx, θ)
    res = 0.0
    for i in eachindex(θ)
        if θ[i] == 0
            continue
        else
            res += ∇ϕx[i]^2
        end
    end
    res
end

function sdot(x,y,θ)
    res = 0.0
    @inbounds for i in eachindex(θ)
        if θ[i] == 0.0
            continue
        else
            res += x[i]*y[i]
        end
    end
    res
end

function reflect_sticky!(∇ϕx, x, θ, f, Flow::Union{BouncyParticle, Boomerang})
    c = (2*sdot(∇ϕx, θ, θ)/subnormsq(∇ϕx, θ))
    for i in eachindex(θ)
        if f[i]
            θ[i] -= c*∇ϕx[i]
        end
    end
    θ
end

function smove_forward!(τ, t, x, θ, f, Z::Union{BouncyParticle, ZigZag})
    t += τ
    for i in eachindex(x)
        if f[i]
            x[i] += θ[i]*τ
        end
    end
    t, x, θ
end

function smove_forward!(τ, t, x, θ, f, B::Union{Boomerang, FactBoomerang})
    s, c = sincos(τ)
    for i in eachindex(x)
        if f[i]
            x[i], θ[i] = (x[i] - B.μ[i])*c + θ[i]*s + B.μ[i],
                        -(x[i] - B.μ[i])*s + θ[i]*c
        end
    end
    t + τ, x, θ
end

function sevent(t, x, θ, f, Z::Union{BouncyParticle,Boomerang})
    t, copy(x), copy(θ), copy(f)
end

function sticky_pdmp_inner!(Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, b, t′, f, θf, tfrez, tref, told, (acc, num),
        Flow::Union{BouncyParticle, Boomerang}, κ, args...; strong_upperbounds = false, factor=1.5, adapt=false)

    # f[i] is true if i is free
    # frez[i] is the time to freeze, if free, the time to unfreeze, if frozen
    while true
        tᶠ, i = findmin(tfrez) # could be implemented with a queue
        tt, j = findmin([tref, tᶠ, t′])
        τ = tt - t
        t, x, θ = smove_forward!(τ, t, x, θ, f, Flow)
        # move forward
        if j == 1 # refreshments of velocities
            θ, θf = refresh_sticky_vel!(θ, θf, f, Flow)
            tref = t + waiting_time_ref(Flow) # regenerate refreshment time
            b = ab(x, θ, c, Flow) # regenerate reflection time
            told = t
            t′, _ = next_time(t, b, rand())
            tfrez = freezing_time!(tfrez, t, x, θ, f, Flow)
            for i in eachindex(f) # make function later...
                if !f[i]
                    tfrez[i] = t - log(rand())/(κ[i]*abs(θf[i]))
                end
            end
        elseif j == 2 # get frozen or unfrozen in i
            if f[i] # if free
                if abs(x[i]) > 1e-8
                    tfrez[i] = t + freezing_time(x[i], θ[i], Flow.μ[i], Flow) # wrong zero of curve
                    error("x[i] = $(x[i]) !≈ 0 at $(tfrez[i])")
                end
                x[i] = -0*θ[i]
                θf[i], θ[i] = θ[i], 0.0 # stop and save speed
                f[i] = false # change tag
                tfrez[i] = t - log(rand())/(κ[i]*abs(θf[i])) # sticky time
                # tfrez[i] = t - log(rand()) # option 2
                if !(strong_upperbounds) #not strong upperbounds, draw new waiting time
                    b = ab(x, θ, c, Flow) # regenerate reflection time
                    told = t
                    t′, _ = next_time(t, b, rand())
                end
            else # is frozen ->  unfreeze
                @assert x[i] == 0 && θ[i] == 0
                θ[i], θf[i] = θf[i], 0.0 # restore speed
                f[i] = true # change tag
                tfrez[i] = t + freezing_time(x[i], θ[i], Flow.μ[i], Flow)
                b = ab(x, θ, c, Flow) # regenerate reflection time
                told = t
                t′, _ = next_time(t, b, rand())
            end
        else #   t′ usual bouncy particle / boomerang step
            ∇ϕx = ∇ϕ!(∇ϕx, x, args...)
            ∇ϕx = grad_correct!(∇ϕx, x, Flow)
            l, lb = λ(∇ϕx, θ, Flow), sλ̄(b, t - told) # CHECK if depends on f
            num += 1
            if rand()*lb <= l # reflect!
                acc += 1
                if l > lb
                    !adapt && error("Tuning parameter `c` too small.")
                    c *= factor
                end
                θ = reflect_sticky!(∇ϕx, x, θ, f, Flow)
                b = ab(x, θ, c, Flow) # regenerate reflection time
                told = t
                t′, _ = next_time(t, b, rand())
                tfrez = freezing_time!(tfrez, t, x, θ, f, Flow)
            else # nothing happened
                b = ab(x, θ, c, Flow)
                told = t
                t′, _ = next_time(t, b, rand())
                continue
            end
        end
        push!(Ξ, sevent(t, x, θ, f, Flow))
        return t, x, θ, t′, tref, tfrez, told, f, θf, (acc, num), c, b
    end

end

function sspdmp(∇ϕ!, t0, x0, θ0, T, c, Flow::Union{BouncyParticle, Boomerang},
        κ, args...;  strong_upperbounds = false, adapt=false, factor=2.0)
    t, x, θ, ∇ϕx = t0, deepcopy(x0), deepcopy(θ0), deepcopy(θ0)
    told = t0
    θf = 0*θ # tags
    f = [true for _ in eachindex(x)]
    Ξ = Trace(t0, x0, θ0, f, Flow)
    push!(Ξ, sevent(t, x0, θ0, f, Flow))
    tref = waiting_time_ref(Flow) #refreshment times
    tfrez = zero(x)
    tfrez = freezing_time!(tfrez, t0, x0, θ0, f, Flow) #freexing times
    num = acc = 0
    b = ab(x, θ, c, Flow)
    t′, _ = next_time(t, b, rand()) # reflection time
    while t < T
        t, x, θ, t′, tref, tfrez, told, f, θf, (acc, num), c, b = sticky_pdmp_inner!(Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, b, t′, f, θf, tfrez, tref, told, (acc, num),
                Flow, κ, args...; strong_upperbounds = strong_upperbounds, factor = factor, adapt = adapt)
    end
    return Ξ, (t, x, θ), (acc, num), c
end


text:
using MCMCBenchmarks, CSV

ProjDir = @__DIR__
cd(ProjDir)

# Load results
folder = "results/2019_10_14T13_55_00/"
results = CSV.read(folder * "results.csv", copycols=true)
allowmissing!(results)

# Make plots
include("primary_plots.jl")


text:
#!/usr/bin/env julia

# A Julia script that reproduces the second part of expected results of
# the test module 'test/intexpTest.cpp'.
# The other parts are performed in 'scripts/test/intexp.m' and
# 'scripts/test/intexp.mac'.
#
# From a shell, run the script as:
#   julia /path/to/intexp.jl
#
# From Julia, the script can be run as:
#   include("path/to/intexp.jl")

# The test cases are performed in three scripts because each mathematical
# structure is (natively) supported by one programming language/mathematical
# software better than by the others.

# continued from test cases in 'intexp.m'

function intexpTest()
    f = -3//4
    println("f = ", f)
    println("f^0 = ", f^0)
    println("f^7 = ", f^7)
end

intexpTest()

# test cases continue in 'intexp.mac'.

max_stars_repo_name: jkovacic/math

max_issues_repo_name: jkovacic/math

max_forks_repo_name: jkovacic/math


text:
julia> RpcApi.SimpleOrganism[RpcApi.SimpleOrganism(0x00000000000395fc,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_52::String = "TPIvuGqq"
        local var_41::Int8 = -113
        local var_40::Bool = false
        local var_30::Bool = false
        local var_29::Int16 = 20638
        function func_24(var_22::Int16=-19845,var_23::Bool=false)
            Creature.stepDown(o)
            o.mem[var_56] = var_56
            Creature.stepLeft(o)
            Creature.stepUp(o)
            Creature.stepRight(o)
            Creature.stepRight(o)
            Creature.stepRight(o)
            Creature.stepRight(o)
            if var_22 < var_22
                Creature.stepLeft(o)
            end
            return var_22
        end
        local var_19::Int16 = -25546
        function func_7(var_5::Int8=28,var_6::Bool=false)
            local var_31::Int64 = -2617427854565851962
            local var_60::Int8 = 120
            Creature.eatLeft(c,o,Int(var_5))
            var_5 = var_5 - var_5
            Creature.eatUp(c,o,Int(var_5))
            Creature.stepLeft(o)
            Creature.eatUp(c,o,Int(var_5))
            Creature.stepRight(o)
            if var_5 <= var_5
                Creature.stepRight(o)
            end
            return var_5
        end
        local var_1::Int8 = 79
        o.mem[var_29] = var_19
        var_1 = var_1 * var_1
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        if var_40 !== var_40
            local var_57::Int8 = -13
            Creature.eatUp(c,o,Int(var_57))
            Creature.stepLeft(o)
            var_57 = var_57 / var_57
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        if var_29 <= var_19
            local var_58::Bool = false
            Creature.stepDown(o)
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        o.mem[var_29] = var_19
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_59::String = "IXldA7Sx"
            local var_55::Int64 = 553205973710463888
            Creature.stepUp(o)
            Creature.stepRight(o)
            Creature.stepDown(o)
            Creature.stepUp(o)
            var_55 = var_55 * var_55
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        return true
    end),[100,300,50,0,1,1,1],2,1000,1,196827254,0x000003a0,Dict(-25546=>-25546,20638=>-25546,14444=>-25546),[238,118],70),RpcApi.SimpleOrganism(0x0000000000044766,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_52::String = "TPIvuGqq"
        local var_41::Int8 = -113
        local var_40::Bool = false
        local var_29::Int16 = 20638
        function func_24(var_22::Int16=-19845,var_23::Bool=false)
            Creature.stepUp(o)
            Creature.stepDown(o)
            if var_23 > var_23
                Creature.stepLeft(o)
                Creature.stepRight(o)
                Creature.stepRight(o)
            end
            Creature.stepDown(o)
            Creature.stepRight(o)
            var_23 = Bool(abs(var_23 - var_23))
            return var_22
        end
        local var_19::Int16 = -25546
        function func_7(var_5::Int8=28,var_6::Bool=false)
            local var_59::Int8 = 61
            local var_31::Int64 = -2617427854565851962
            Creature.eatLeft(c,o,Int(var_5))
            var_31 = var_31 - var_31
            Creature.stepDown(o)
            Creature.stepLeft(o)
            var_59 = var_59 - var_5
            Creature.eatUp(c,o,Int(var_5))
            Creature.eatUp(c,o,Int(var_59))
            Creature.stepUp(o)
            Creature.eatRight(c,o,Int(var_59))
            Creature.stepDown(o)
            Creature.eatDown(c,o,Int(var_59))
            var_6 = var_6 * var_6
            Creature.stepRight(o)
            Creature.stepDown(o)
            Creature.eatDown(c,o,Int(var_59))
            return var_5
        end
        local var_1::Int8 = 79
        o.mem[var_29] = var_19
        var_1 = var_1 * var_1
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        var_52 = var_52[1:if length(var_52) > length(var_52) > 0
                        0
                    else
                        div(length(var_52),length(var_52))
                    end]
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_55::Int64 = 553205973710463888
            Creature.stepRight(o)
            var_55 = var_55 / var_55
            Creature.stepDown(o)
            Creature.stepUp(o)
            var_55 = var_55 * var_55
        end
        var_41 = var_41 - var_41
        return true
    end),[100,300,50,0,1,1,1],23,1000,1,110674095,0x000003ad,Dict(-25546=>-25546,20638=>-25546,14444=>-25546),[241,121],68),RpcApi.SimpleOrganism(0x000000000003c1c8,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_61(var_59::Int8=19,var_60::Int8=-88)
            local var_62::Int16 = 32507
            if var_59 === var_60
            end
            Creature.eatLeft(c,o,Int(var_59))
            return var_59
        end
        local var_52::String = "TPIvuGqq"
        local var_41::Int8 = -113
        local var_40::Bool = false
        local var_30::Bool = false
        local var_29::Int16 = 20638
        function func_24(var_22::Int16=-19845,var_23::Bool=false)
            local var_63::Int8 = 127
            local var_56::Int16 = -25543
            Creature.stepUp(o)
            var_22 = var_56 / var_22
            var_22 = var_56 - var_22
            var_23 = var_23 & var_23
            Creature.stepLeft(o)
            return var_22
        end
        local var_19::Int16 = -25546
        function func_7(var_5::Int8=28,var_6::Bool=false)
            local var_31::Int64 = -2617427854565851962
            Creature.eatLeft(c,o,Int(var_5))
            Creature.eatDown(c,o,Int(var_5))
            var_5 = var_5 - var_5
            Creature.stepDown(o)
            Creature.stepRight(o)
            Creature.stepUp(o)
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
                Creature.stepRight(o)
                Creature.stepDown(o)
            end
            Creature.stepUp(o)
            Creature.stepUp(o)
            var_5 = var_5 * var_5
            Creature.stepUp(o)
            Creature.stepDown(o)
            return var_5
        end
        local var_1::Int8 = 79
        o.mem[var_29] = var_19
        var_1 = var_1 * var_1
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        if var_40 !== var_40
            local var_57::Int8 = -13
            var_57 = var_57 * var_57
            Creature.eatLeft(c,o,Int(var_57))
            Creature.eatUp(c,o,Int(var_57))
            Creature.eatLeft(c,o,Int(var_57))
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        o.mem[var_29] = var_19
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_55::Int64 = 553205973710463888
            Creature.stepRight(o)
            Creature.stepLeft(o)
            Creature.stepLeft(o)
            var_55 = var_55 * var_55
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        return true
    end),[100,300,50,0,1,1,1],40,1000,1,74445826,0x000003a2,Dict(-25546=>-25546,20638=>-25546,14444=>-25546),[241,119],74),RpcApi.SimpleOrganism(0x00000000000385b9,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_52::String = "TPIvuGqq"
        local var_41::Int8 = -113
        local var_40::Bool = false
        function func_39(var_37::Int64=-7847296428555035582,var_38::Bool=true)
            local var_61::String = "zN9a26yX"
            Creature.stepUp(o)
            return var_37
        end
        local var_30::Bool = false
        function func_24(var_22::Int16=-19845,var_23::Bool=false)
            local var_44::Int64 = -2474270459538129053
            Creature.eatRight(c,o,Int(var_48))
            var_23 = var_23 * var_23
            var_44 = var_44 * var_44
            Creature.stepLeft(o)
            o.mem[var_22] = var_22
            Creature.stepDown(o)
            var_22 = var_22 - var_22
            Creature.stepDown(o)
            Creature.stepRight(o)
            var_23 = var_23 | var_23
            return var_22
        end
        local var_19::Int16 = -25546
        function func_7(var_5::Int8=28,var_6::Bool=false)
            local var_59::String = "iqDQDiUN"
            local var_57::String = "rgvYYKtL"
            Creature.eatLeft(c,o,Int(var_5))
            Creature.eatLeft(c,o,Int(var_5))
            var_6 = var_6 & var_6
            Creature.stepLeft(o)
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
            end
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
            end
            Creature.eatUp(c,o,Int(var_5))
            var_57 = var_57[1:if length(var_57) > length(var_57) > 0
                            0
                        else
                            div(length(var_57),length(var_57))
                        end]
            var_59 = var_59 * var_59
            Creature.eatLeft(c,o,Int(var_5))
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
            end
            Creature.stepDown(o)
            return var_5
        end
        local var_1::Int8 = 79
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = var_1 * var_1
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        if var_40 !== var_40
            local var_56::Bool = true
            Creature.stepLeft(o)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        if var_29 <= var_19
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        o.mem[var_19] = var_19
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_55::Int8 = -87
            local var_60::Int8 = 109
            Creature.stepRight(o)
            Creature.stepRight(o)
            var_55 = var_55 * var_55
            Creature.eatRight(c,o,Int(var_55))
            Creature.stepUp(o)
        end
        var_41 = var_1 + var_1
        return true
    end),[100,300,50,0,1,1,1],44,1000,1,54904827,0x000003a6,Dict(-25546=>-25546,20638=>-25546),[238,119],73),RpcApi.SimpleOrganism(0x000000000003c41f,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_52::String = "TPIvuGqq"
        local var_41::Int8 = -113
        local var_40::Bool = false
        local var_30::Bool = false
        local var_29::Int16 = 20638
        function func_24(var_22::Int16=-19845,var_23::Bool=false)
            Creature.stepRight(o)
            var_22 = var_22 * var_22
            var_22 = var_22 / var_22
            Creature.stepRight(o)
            Creature.stepUp(o)
            Creature.stepRight(o)
            Creature.stepLeft(o)
            Creature.stepUp(o)
            Creature.stepLeft(o)
            return var_22
        end
        local var_19::Int16 = -25546
        function func_7(var_5::Int8=28,var_6::Bool=false)
            local var_31::Int64 = -2617427854565851962
            Creature.eatDown(c,o,Int(var_5))
            Creature.stepDown(o)
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
                Creature.stepRight(o)
            end
            var_31 = var_31 * var_31
            if var_5 < var_5
            end
            Creature.stepDown(o)
            var_5 = var_5 - var_5
            Creature.eatRight(c,o,Int(var_5))
            Creature.eatUp(c,o,Int(var_5))
            for i::Int8 = 1:div(var_5,_LOOP_STEPS_DIVIDER)
                Creature.stepRight(o)
                Creature.stepRight(o)
            end
            Creature.stepUp(o)
            Creature.stepDown(o)
            Creature.stepDown(o)
            Creature.stepUp(o)
            Creature.stepDown(o)
            return var_5
        end
        local var_1::Int8 = 79
        var_1 = var_1 * var_1
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        if var_40 !== var_40
            local var_57::Int8 = -13
            local var_59::Int64 = 7137138221165316371
            Creature.eatUp(c,o,Int(var_57))
            var_57 = var_57 / var_57
            Creature.stepLeft(o)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        if var_29 <= var_19
            local var_58::Bool = false
            Creature.stepDown(o)
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        o.mem[var_29] = var_19
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        var_29 = var_19 + var_19
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            Creature.stepUp(o)
            Creature.stepRight(o)
            Creature.stepUp(o)
            Creature.stepLeft(o)
            var_55 = var_55 * var_55
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_60::String = "Kq9rAKfT"
        end
        return true
    end),[100,300,50,0,1,1,1],40,1000,1,54292311,0x000003a5,Dict(-25546=>-25546,20638=>-25546,14444=>-25546),[234,117],79)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


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

mutable struct VirtualNetworkGatewayIPConfigurationPropertiesFormat <: SwaggerModel
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, String } # spec name: privateIPAllocationMethod
    subnet::Any # spec type: Union{ Nothing, SubResource } # spec name: subnet
    publicIPAddress::Any # spec type: Union{ Nothing, SubResource } # spec name: publicIPAddress
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function VirtualNetworkGatewayIPConfigurationPropertiesFormat(;privateIPAllocationMethod=nothing, subnet=nothing, publicIPAddress=nothing, provisioningState=nothing)
        o = new()
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(VirtualNetworkGatewayIPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type VirtualNetworkGatewayIPConfigurationPropertiesFormat

const _property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("subnet")=>Symbol("subnet"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_VirtualNetworkGatewayIPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("privateIPAllocationMethod")=>"String", Symbol("subnet")=>"SubResource", Symbol("publicIPAddress")=>"SubResource", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualNetworkGatewayIPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayIPConfigurationPropertiesFormat[property_name]

const _allowed_VirtualNetworkGatewayIPConfigurationPropertiesFormat_privateIPAllocationMethod = ["Static", "Dynamic"]

function check_required(o::VirtualNetworkGatewayIPConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkGatewayIPConfigurationPropertiesFormat }, name::Symbol, val)
    if name === Symbol("privateIPAllocationMethod")
        Swagger.validate_param(name, "VirtualNetworkGatewayIPConfigurationPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayIPConfigurationPropertiesFormat_privateIPAllocationMethod)
    end
end


text:
@doc """

    ar2(y, x, p, βbar, Bbarinv, ϕbar, Vbarinv, υbar, δbar, βold, ϕold, σ2old, varind, lev2ind, factors, numobs, varassign)

Description:
Estimate observed variable hyperparameters conditional on the factors.
Refer to pg. 1002 of Otrok and Whiteman (1998), or Chib and Greenberg (1994).

Inputs:
- y             = observed series vector.
- x             = data matrix of regressors in the observation equation (intercept + factors).
- p             = number of idiosyncratic error lags.
- βbar          = prior mean of factor loading (factor coefficient in measurement equation).
- Bbarinv       = prior variance of factor loading (factor coefficient in measurement equation).
- ϕbar          = prior mean of phi (idiosyncratic AR polynomial).
- Vbarinv       = prior precision of phi (idiosyncratic AR polynomial).
- υbar          = inverted gamma parameter of innovation variances.
- δbar          = inverted gamma parameter of innovation variances.
- βold          = old draw of factor loading.
- ϕold          = old draw of phi (idiosyncratic AR polynomial).
- σ2old         = old draw of innovation variances.
- varind        = measurement eq. regressor index.
- lev2ind       = level-2 factor index.
- factors       = factor draw.
- numobs        = total number of time periods / observations in the sample.
- varassign     = list containing vector of observed series indeces assigned to each level-2 factor.

Outputs:
- βnew      = β coefficient hyperparameter estimates.
- ϕnew      = ϕ coefficient hyperparameter estimates.
- σ2new     = σ^2 coefficient hyperparameter estimate.
- factors   = factor estimates.
"""

function ar2(y, x, p, βbar, Bbarinv, ϕbar, Vbarinv, υbar, δbar, βold, ϕold, σ2old, varind, lev2ind, factors, numobs, varassign)

    # Guarantee that y is a vector
    y = vec(y)

    # Define ỹ = first p observations of series vector
    ỹ = y[1:p]

    # Define x̃ = first p observations of regressor matrix
    x̃ = x[1:p, :]

    # Create Φ matrix (companion matrix of idiosyncratic error autoregression)
    Φ = [ϕold'; I(p - 1) zeros(p - 1)]

    # Covariance matrix component of the first p errors (Σ in σ2Σ)
    vecΣ = inv(I - Φ ⊗ Φ) * vec([1; zeros(p - 1)] * [1; zeros(p - 1)]')
    Σ = reshape(vecΣ, p, p)
    Σ = Hermitian(Σ)

    # Compute the Cholesky factor of Σ
    # notation is such that Q Q' = Σ
    # so that Q is the lower factor of Σ
    Q = cholesky(Σ).L

    # Define ỹ∗1 = Q^{-1} ỹ and x̃∗1 = Q^{-1} x̃
    ỹstar1 = inv(Q) * ỹ
    x̃star1 = inv(Q) * x̃

    # Define ỹ∗2 as a (T-p × 1) with t-th row ϕ(L)y
    # and x̃∗2 as a (T-p × 3) with t-th row ϕ(1)
    ỹstar2 = gendiff(y, ϕold)
    x̃star2 = [gendiff(ones(numobs), ϕold) gendiff(factors[:, 1], ϕold) gendiff(factors[:, lev2ind], ϕold)]

    # Define e = last T-p "observed" idiosyncratic errors
    e = y - x * βold
    e = e[(p+1):numobs]

    # Define E = [ lag(e,1) lag(e,2) ... lag(e,p) ]
    E = zeros(numobs, p)
    for j = 1:p
        E[:, j] = lag(y - x * βold, j, default = 0.0)
    end
    E = E[(p+1):numobs, :]

    # Define x̃∗ = [x̃∗1 ; x̃∗2] and ỹ∗ = [ỹ∗1 ; ỹ∗2]
    x̃star = [x̃star1; x̃star2]
    ỹstar = [ỹstar1; ỹstar2]

    #############################################
    #############################################
    ## Create objects for posterior densities

    # Define V̅ and B̅
    Vbar = inv(Vbarinv)
    Bbar = inv(Bbarinv)

    # Define V = V̅ + σ²E'E
    V = Vbar + σ2old * E' * E

    # Define B = B̅ + ̅σ²x̃∗'x̃∗
    B = Bbar + σ2old * x̃star' * x̃star

    # Define ̂ϕ = inv(V)(V̅ ̅ϕ + σ²E'e)
    ϕhat = inv(V) * (Vbar * ϕbar + σ2old * E' * e)

    #############################################
    #############################################
    ## Draw β

    # Initialize sign of level-1 factor loading
    # of first series in dataset (must be positive)
    β2sign = false

    # Initialize sign of level-2 factor loading
    # of first series in corresponding level-2 group
    # (must be positive)
    β3sign = false

    # Initialize β
    βnew = βold

    # Keep track of wrong sign persistence
    β2signmax = 0
    β3signmax = 0

    whileind = 0
    while β2sign == false || β3sign == false

        whileind += 1
        #=println("Series index: $varind")
        println("While loop index: $whileind")
        println("β2sigmax: $β2signmax")
        println("β3signmax: $β3signmax")
        =#

        # Draw new β
        #βnew = sim_MvNormal(inv(B) * (Bbar * βbar + inv(σ2old) * x̃star' * ỹstar), inv(B))
        β0 = zeros(size(x)[2])
        Σ0 = Matrix(I, size(β0)[1], size(β0)[1]) .* 1000.0
        ## Posterior parameters in N(β1,Σ1)
        β1 = inv(inv(Σ0) + inv(σ2old) * transpose(x) * x) * (inv(Σ0) * β0 + inv(σ2old) * transpose(x) * y)
        β1 = vec(β1)
        Σ1 = inv(inv(Σ0) + inv(σ2old) * transpose(x) * x)
        Σ1 = Hermitian(Σ1)

        ## Generate new β
        βnew = sim_MvNormal(β1, Σ1)

        # Check if new β draw satisfies
        # sign identification restrictions
        if varind == 1
            β2sign = (βnew[2] > 0)
            if β2sign == false
                β2signmax += 1
            end
        else
            β2sign = true
        end
        if varind == varassign[lev2ind][1]
            β3sign = (βnew[3] > 0)
            if β3sign == false
                β3signmax += 1
            end
        else
            β3sign = true
        end

        if β2signmax > 100

            # Reflect level-1 factor over x-axis
            factors[:, 1] = -factors[:, 1]

            # Redefine x
            x = [ones(numobs) factors[:, 1] factors[:, 1+lev2ind]]

            # Define x̃ = first p observations of regressor matrix
            x̃ = x[1:p, :]

            # Create Φ matrix (companion matrix of idiosyncratic error autoregression)
            Φ = [ϕold'; I(p - 1) zeros(p - 1)]

            # Covariance matrix component of the first p errors (Σ in σ2Σ)
            vecΣ = inv(I - Φ ⊗ Φ) * vec([1; zeros(p - 1)] * [1; zeros(p - 1)]')
            Σ = reshape(vecΣ, p, p)
            Σ = Hermitian(Σ)

            # Compute the Cholesky factor of Σ
            # notation is such that Q Q' = Σ
            # so that Q is the lower factor of Σ
            Q = cholesky(Σ).L

            # Define ỹ∗1 = Q^{-1} ỹ and x̃∗1 = Q^{-1} x̃
            ỹstar1 = inv(Q) * ỹ
            x̃star1 = inv(Q) * x̃

            # Define ỹ∗2 as a (T-p × 1) with t-th row ϕ(L)y
            # and x̃∗2 as a (T-p × 3) with t-th row ϕ(1)
            ỹstar2 = gendiff(y, ϕold)
            x̃star2 = [gendiff(ones(numobs), ϕold) gendiff(factors[:, 1], ϕold) gendiff(factors[:, lev2ind], ϕold)]

            # Define e = last T-p "observed" idiosyncratic errors
            e = y - x * βold
            e = e[(p+1):numobs]

            # Define E = [ lag(e,1) lag(e,2) ... lag(e,p) ]
            E = zeros(numobs, p)
            for j = 1:p
                E[:, j] = lag(y - x * βold, j, default = 0.0)
            end
            E = E[(p+1):numobs, :]

            # Define x̃∗ = [x̃∗1 ; x̃∗2] and ỹ∗ = [ỹ∗1 ; ỹ∗2]
            x̃star = [x̃star1; x̃star2]
            ỹstar = [ỹstar1; ỹstar2]

            # Reset `β2signmax`
            β2signmax = 0
        end

        if β3signmax > 1000

            # Reflect level-1 factor over x-axis
            factors[:, 1+lev2ind] = -factors[:, 1+lev2ind]

            # Redefine x
            x = [ones(numobs) factors[:, 1] factors[:, 1+lev2ind]]

            # Define x̃ = first p observations of regressor matrix
            x̃ = x[1:p, :]

            # Create Φ matrix (companion matrix of idiosyncratic error autoregression)
            Φ = [ϕold'; I(p - 1) zeros(p - 1)]

            # Covariance matrix component of the first p errors (Σ in σ2Σ)
            vecΣ = inv(I - Φ ⊗ Φ) * vec([1; zeros(p - 1)] * [1; zeros(p - 1)]')
            Σ = reshape(vecΣ, p, p)
            Σ = Hermitian(Σ)

            # Compute the Cholesky factor of Σ
            # notation is such that Q Q' = Σ
            # so that Q is the lower factor of Σ
            Q = cholesky(Σ).L

            # Define ỹ∗1 = Q^{-1} ỹ and x̃∗1 = Q^{-1} x̃
            ỹstar1 = inv(Q) * ỹ
            x̃star1 = inv(Q) * x̃

            # Define ỹ∗2 as a (T-p × 1) with t-th row ϕ(L)y
            # and x̃∗2 as a (T-p × 3) with t-th row ϕ(1)
            ỹstar2 = gendiff(y, ϕold)
            x̃star2 = [gendiff(ones(numobs), ϕold) gendiff(factors[:, 1], ϕold) gendiff(factors[:, lev2ind], ϕold)]

            # Define e = last T-p "observed" idiosyncratic errors
            e = y - x * βold
            e = e[(p+1):numobs]

            # Define E = [ lag(e,1) lag(e,2) ... lag(e,p) ]
            E = zeros(numobs, p)
            for j = 1:p
                E[:, j] = lag(y - x * βold, j, default = 0.0)
            end
            E = E[(p+1):numobs, :]

            # Define x̃∗ = [x̃∗1 ; x̃∗2] and ỹ∗ = [ỹ∗1 ; ỹ∗2]
            x̃star = [x̃star1; x̃star2]
            ỹstar = [ỹstar1; ỹstar2]

            # Reset `β2signmax`
            β3signmax = 0
        end

    end

    #############################################
    #############################################
    ## Draw ϕ

    #=
    # Draw new ϕ candidate
    ϕnew = sim_MvNormal(ϕhat, inv(V))

    # Check for stationary of new ϕ draw
    coef = [-reverse(ϕnew); 1]
    coef = vec(coef)
    root = roots(Polynomial(reverse(coef)))
    rootmod = abs.(root)
    accept = min(rootmod...) >= 1.001

    # Complete
    if accept == 0
        ϕnew = ϕold
    else
        # Define Ψ(ϕold) (pg. 1002)
        Ψold = det(sigmat(ϕold, p))^(-1 / 2) * exp(-(1 / (2 * σ2old)) * (ỹ - x̃ * βnew)' * inv(sigmat(ϕold, p)) * (ỹ - x̃ * βnew))
        Ψold = first(Ψold)

        # Define Ψ(ϕnew) (pg. 1002)
        Ψnew = det(sigmat(ϕnew, p))^(-1 / 2) * exp(-(1 / (2 * σ2old)) * (ỹ - x̃ * βnew)' * inv(sigmat(ϕnew, p)) * (ỹ - x̃ * βnew))
        Ψnew = first(Ψnew)

        # Determine acceptance probability
        if Ψold == 0
            accept = 1
        else
            u = rand()
            accept = (u <= Ψnew / Ψold)
        end

        # Pick new ϕ
        ϕnew = ϕnew * accept + ϕold * (1 - accept)
    end
    =#

    # Generate ϕ^j
    ## Prior parameters in N(c0,B0)
    c0 = zeros(size(ϕold)[1])
    B0 = Matrix(I, size(c0)[1], size(c0)[1]) .* 1000.0
    ## Generate e^⋆
    e_star = y - x * βnew
    ## Generate E^⋆
    E_star = zeros(numobs, length(ϕold))
    E_star = E_star[(1+length(ϕold)):end, :]
    for i = 1:(size(E_star)[2]) # iterate over variables in X
        e_temp = lag(e_star, i)
        E_star[:, i] = e_temp[(1+length(ϕold)):end]
    end
    e_star = e_star[(1+length(ϕold)):end]
    ## Posterior parameters in N(c1,B1)
    c1 = inv(inv(B0) + inv(σ2old) * transpose(E_star) * E_star) * (inv(B0) * c0 + inv(σ2old) * transpose(E_star) * e_star)
    B1 = inv(inv(B0) + inv(σ2old) * transpose(E_star) * E_star)
    B1 = Hermitian(B1)
    ## Generate new ϕ
    ϕnew = sim_MvNormal(c1, B1)

    #############################################
    #############################################
    ## Draw σ2

    # Draw new σ2 from inverse gamma distribution
    #σ2new = rand(InverseGamma((υbar + numobs) / 2, (δbar + norm(ỹstar - x̃star * βnew)^2) / 2))
    ν0 = 0.002
    δ0 = 0.002
    ## Posterior parameters in IG(ν1/2, δ1/2)
    ν1 = ν0 + numobs
    δ1 = δ0 + norm(y - x * βnew)^2
    ## Generate new σ2
    σ2new = rand(InverseGamma(ν1 / 2, δ1 / 2))

    #############################################
    #############################################
    ## Return new hyperparameter draws
    return βnew, σ2new, ϕnew, factors
end


text:
using PackageCompiler

if !isdefined(Main, :builddir)
    builddir = "build"
end
create_app(".", builddir, force=true)


text:
using Libtask

@testset "tapedfunction" begin
    # Test case 1: stack allocated objects are deep copied.
    @testset "Instruction{typeof(__new__)}" begin
        mutable struct S
            i::Int
            S(x, y) = new(x + y)
        end

        tf = Libtask.TapedFunction(S, 1, 2)
        s1 = tf(1, 2)
        @test s1.i == 3
        newins = findall(x -> isa(x, Libtask.Instruction{typeof(Libtask.__new__)}), tf.tape)
        @test length(newins) == 1
    end

    @testset "Compiled Tape" begin
        function g(x, y)
            if x>y
                r= string(sin(x))
            else
                r= sin(x) * cos(y)
            end
            return r
        end

        tf = Libtask.TapedFunction(g, 1., 2.)
        ctf = Libtask.compile(tf)
        r = ctf(1., 2.)

        @test typeof(r) === Float64
    end
end

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
# String matching functions for Aaron
using StringDistances

# Fairly arbitrarily chosen list of stop words
const stopwords = [
    "the",
    "is",
    "at",
    "which",
    "on",
    "in",
    "for",
    "with"
]

function clean_string(x)
    xlwr = lowercase(x)
    xcln = replace(xlwr, r"[^-a-z]", " ")     # keep only letters
    return xcln
end

function cut_stopwords(x_arr, stopwords)
    res = [z for z in x_arr if z ∉ stopwords]
    return res
end

function is_anagram(a, b)
    if length(a) ≠ length(b)
        res = false
    else
        a_sorted = sort(collect(a))
        b_sorted = sort(collect(b))
        res = a_sorted == b_sorted
    end
    return res
end

# Given two strings, this function returns a boolean indicating
# whether or not the first letters of all words in the phrases
# are the same.
function same_first_letters(x, y)
    x_arr = split(x)
    y_arr = split(y)
    n_words = length(x_arr)
    res = true
    if n_words ≠ length(y_arr)
        warn("Attempted to compare phrases with differing number of words")
        res = false
    else
        for i = 1:n_words
            if first(x_arr[i]) ≠ first(y_arr[i])
                res = false
                break
            end
        end
    end
    return res
end

## testing
same_first_letters("this is a phrase", "this is also phrase")
same_first_letters("this is a phrase", "this also phrase")

# This function uses moving window to compare words
# combination from two strings (strings are arrays of words).
function found_match(x_arr, y_arr, thresh)
    # By convention, y_arr will be longer
    len_y = length(y_arr)
    len_x = length(x_arr)
    x_str = join(x_arr, " ")
    windowsize = len_x

    for i = 1:(len_y - windowsize + 1)
        y_str = join(y_arr[i:(i + windowsize - 1)], " ")
        n_char = length(y_str)

        if !same_first_letters(x_str, y_str)
            continue
        end
        println("Comparing: \'$y_str\' and \'$x_str\'")
        similarity = compare(DamerauLevenshtein(), y_str, x_str)

        if n_char ≤ 8 && similarity ≥ thresh
            println("Matched \'$y_str\' and \'$x_str\'")
            return true
        elseif n_char ≤ 12 && similarity ≥ thresh - 0.05
            println("Matched \'$y_str\' and \'$x_str\' with similarity $similarity")
            return true
        elseif n_char ≤ 18 && similarity ≥ thresh - 0.1
            println("Matched \'$y_str\' and \'$x_str\' with similarity $similarity")
            return true
        elseif similarity ≥ thresh - 0.15
            println("Matched \'$y_str\' and \'$x_str\' with similarity $similarity")
            return true
        end
    end
    return false
end

function is_string_match(x, y, thresh, stopwords)
    x_cln = clean_string(x)
    y_cln = clean_string(y)

    x_arr = cut_stopwords(split(x_cln), stopwords)
    y_arr = cut_stopwords(split(y_cln), stopwords)

    len_x = length(x_arr)
    len_y = length(y_arr)

    if len_x < len_y
        res = found_match(x_arr, y_arr, thresh)
    else
        # switching position of x and y
        res = found_match(y_arr, x_arr, thresh)
    end
    return res
end

# testing our functions
target1 = "pericardial effusion"
candidate1 = "Metastatic non-small cell lung cancer.  Dyspnea. No pulsus on exam. Evaluate for perracardial perracardial effusion/interval change, signs of tamponade."

threshold1 = 0.75

is_string_match(target1, candidate1, threshold1, stopwords)
is_string_match(candidate1, target1, threshold1, stopwords)

x1 = "left ventricular left ventricular hypertophy"
y1 = "left ventricular fxn"

is_string_match(x1, y1, 0.95, stopwords)
is_string_match(y1, x1, 0.95, stopwords)


text:
module Coefplots
    using StatsModels
    using PGFPlotsX
    using Printf
    using DataStructures
    import FixedEffectModels: FixedEffectModel
    import GLFixedEffectModels: GLFixedEffectModel
    import Base: issingletontype, parse

    include("color.jl")
    include("preamble.jl")
    include("options.jl")
    include("pgfplot.jl")
    include("singlecoefplot.jl")
    include("coefplot.jl")
    include("multicoefplot.jl")
    include("parse.jl")
    include("plot.jl")
    include("esplot.jl")

    export setxtitle!, setytitle!, setname!, includenote!, equidist!, setlegends!, addcomponent!, clearcomponents!, coefrename, coefrename!, transpose!
    export plot, parse, esplot, esparse
end


text:
# Trigger several successive GC sweeps. This is more comprehensive than running just a
# single sweep, since freeable objects may need more than one sweep to be appropriately
# marked and freed.
gcscrub() = (GC.gc(); GC.gc(); GC.gc(); GC.gc())

#############
# Benchmark #
#############

mutable struct Benchmark
    samplefunc
    params::Parameters
end

params(b::Benchmark) = b.params

function loadparams!(b::Benchmark, params::Parameters, fields...)
    loadparams!(b.params, params, fields...)
    return b
end

function Base.show(io::IO, b::Benchmark)
    str = string("Benchmark(evals=", params(b).evals,
                 ", seconds=", params(b).seconds,
                 ", samples=", params(b).samples, ")")
    print(io, str)
end

######################
# compatiblity hacks #
######################

run_result(b::Benchmark, p::Parameters = b.params; kwargs...) = Base.invokelatest(_run, b, p; kwargs...)
lineartrial(b::Benchmark, p::Parameters = b.params; kwargs...) = Base.invokelatest(_lineartrial, b, p; kwargs...)

##############################
# progress logging utilities #
##############################

# As used in ProgressLogging.jl
# https://github.com/JunoLab/ProgressLogging.jl/blob/v0.1.0/src/ProgressLogging.jl#L11
const ProgressLevel = LogLevel(-1)

"""
    _withprogress(
        name::AbstractString,
        group::BenchmarkGroup;
        kwargs...,
    ) do progressid, nleaves, ndone
        ...
    end

Execute do block with following arguments:

* `progressid`: logging ID to be used for `@logmsg`.
* `nleaves`: total number of benchmarks counted at the root benchmark group.
* `ndone`: number of completed benchmarks

They are either extracted from `kwargs` (for sub-groups) or newly created
(for root benchmark group).
"""
function _withprogress(
    f,
    name::AbstractString,
    group::BenchmarkGroup;
    progressid = nothing,
    nleaves = NaN,
    ndone = NaN,
    _...,
)
    if progressid !== nothing
        return f(progressid, nleaves, ndone)
    end
    progressid = uuid4()
    nleaves = length(leaves(group))
    @logmsg(ProgressLevel, name, progress = NaN, _id = progressid)
    try
        return f(progressid, nleaves, 0)
    finally
        @logmsg(ProgressLevel, name, progress = "done", _id = progressid)
    end
end

#############
# execution #
#############

# Note that trials executed via `run` and `lineartrial` are always executed at top-level
# scope, in order to allow transfer of locally-scoped variables into benchmark scope.

function _run(b::Benchmark, p::Parameters; verbose = false, pad = "", kwargs...)
    params = Parameters(p; kwargs...)
    @assert params.seconds > 0.0 "time limit must be greater than 0.0"
    params.gctrial && gcscrub()
    start_time = Base.time()
    trial = Trial(params)
    params.gcsample && gcscrub()
    s = b.samplefunc(params)
    push!(trial, s[1:end-1]...)
    return_val = s[end]
    iters = 2
    while (Base.time() - start_time) < params.seconds && iters ≤ params.samples
         params.gcsample && gcscrub()
         push!(trial, b.samplefunc(params)[1:end-1]...)
         iters += 1
    end
    return sort!(trial), return_val
end

"""
    run(b::Benchmark[, p::Parameters = b.params]; kwargs...)

Run the benchmark defined by [`@benchmarkable`](@ref).
"""
Base.run(b::Benchmark, p::Parameters = b.params; progressid=nothing, nleaves=NaN, ndone=NaN, kwargs...) =
    run_result(b, p; kwargs...)[1]

"""
    run(group::BenchmarkGroup[, args...]; verbose::Bool = false, pad = "", kwargs...)

Run the benchmark group, with benchmark parameters set to `group`'s by default.
"""
Base.run(group::BenchmarkGroup, args...; verbose::Bool = false, pad = "", kwargs...) =
    _withprogress("Benchmarking", group; kwargs...) do progressid, nleaves, ndone
        result = similar(group)
        gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
        i = 1
        for id in keys(group)
            @logmsg(ProgressLevel, "Benchmarking", progress = ndone / nleaves, _id = progressid)
            verbose &&
                println(pad, "($(i)/$(length(group))) benchmarking ", repr(id), "...")
            took_seconds = @elapsed begin
                result[id] = run(
                    group[id],
                    args...;
                    verbose = verbose,
                    pad = pad * "  ",
                    kwargs...,
                    progressid = progressid,
                    nleaves = nleaves,
                    ndone = ndone,
                )
            end
            ndone += group[id] isa BenchmarkGroup ? length(leaves(group[id])) : 1
            verbose && (println(pad, "done (took ", took_seconds, " seconds)"); i += 1)
        end
        return result
    end

function _lineartrial(b::Benchmark, p::Parameters = b.params; maxevals = RESOLUTION, kwargs...)
    params = Parameters(p; kwargs...)
    estimates = zeros(maxevals)
    completed = 0
    params.gctrial && gcscrub()
    start_time = time()
    for evals in eachindex(estimates)
        params.gcsample && gcscrub()
        params.evals = evals
        estimates[evals] = first(b.samplefunc(params))
        completed += 1
        ((time() - start_time) > params.seconds) && break
    end
    return estimates[1:completed]
end

function warmup(item; verbose::Bool = true)
    return run(item;
               verbose = verbose, samples = 1,
               evals = 1, gctrial = false,
               gcsample = false)
end

####################
# parameter tuning #
####################

# The tuning process is as follows:
#
#   1. Using `lineartrial`, take one sample of the benchmark for each `evals` in `1:RESOLUTION`.
#
#   2. Extract the minimum sample found in this trial. Hopefully, this value will be
#      reasonably close to the true benchmark time. At the very least, we can be certain
#      that the minimum sample overcomes discretization noise, assuming a sufficient time
#      budget for the trial (which technically should never need to be more than a couple of
#      seconds if `RESOLUTION` ∼ 1μs).
#
#   3. Using the sample from step 2 as a reasonable estimate `t`, calculate `evals`. This
#      is done by indexing into a pregenerated table that maps potential `t` values to the
#      appropriate `evals` values. The actual generation of this table is accomplished via
#      heuristic, driven mainly by empircal insight. It's built on the following postulates:
#          - Small variations in `t` should not cause large variations in `evals`
#          - `RESOLUTION` <= 1000ns
#          - `evals` <= `RESOLUTION` if `RESOLUTION` is in nanoseconds and `t` >= 1ns
#          - Increasing `t` should generally decrease `evals`
#
#   4. By default, tune the `samples` parameter as a function of `t`, `evals`, and the time
#      budget for the benchmark.

# the logistic function is useful for determining `evals` for `1 < t < RESOLUTION`
logistic(u, l, k, t, t0) = round(Int, ((u - l) / (1 + exp(-k * (t - t0)))) + l)

const EVALS = Vector{Int}(undef, 9000) # any `t > length(EVALS)` should get an `evals` of 1
for t in 1:400    (EVALS[t] = logistic(1006, 195, -0.025, t, 200)) end # EVALS[1] == 1000, EVALS[400] == 200
for t in 401:1000 (EVALS[t] = logistic(204, -16, -0.01, t, 800))   end # EVALS[401] == 200, EVALS[1000] == 10
for i in 1:8      (EVALS[((i*1000)+1):((i+1)*1000)] .= 11 - i)     end # linearly decrease from EVALS[1000]

guessevals(t) = t <= length(EVALS) ? EVALS[t] : 1

"""
    tune!(group::BenchmarkGroup; verbose::Bool = false, pad = "", kwargs...)

Tune a `BenchmarkGroup` instance. For most benchmarks, `tune!` needs to perform many
evaluations to determine the proper parameters for any given benchmark - often more
evaluations than are performed when running a trial. In fact, the majority of total
benchmarking time is usually spent tuning parameters, rather than actually running
trials.
"""
tune!(group::BenchmarkGroup; verbose::Bool = false, pad = "", kwargs...) =
    _withprogress("Tuning", group; kwargs...) do progressid, nleaves, ndone
        gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
        i = 1
        for id in keys(group)
            @logmsg(ProgressLevel, "Tuning", progress = ndone / nleaves, _id = progressid)
            verbose && println(pad, "($(i)/$(length(group))) tuning ", repr(id), "...")
            took_seconds = @elapsed tune!(
                group[id];
                verbose = verbose,
                pad = pad * "  ",
                kwargs...,
                progressid = progressid,
                nleaves = nleaves,
                ndone = ndone,
            )
            ndone += group[id] isa BenchmarkGroup ? length(leaves(group[id])) : 1
            verbose && (println(pad, "done (took ", took_seconds, " seconds)"); i += 1)
        end
        return group
    end

"""
    tune!(b::Benchmark, p::Parameters = b.params; verbose::Bool = false, pad = "", kwargs...)

Tune a `Benchmark` instance.
"""
function tune!(b::Benchmark, p::Parameters = b.params;
               progressid=nothing, nleaves=NaN, ndone=NaN,  # ignored
               verbose::Bool = false, pad = "", kwargs...)
    warmup(b, verbose = false)
    estimate = ceil(Int, minimum(lineartrial(b, p; kwargs...)))
    b.params.evals = guessevals(estimate)
    return b
end

#############################
# @benchmark/@benchmarkable #
#############################

function prunekwargs(args...)
    @nospecialize
    firstarg = first(args)
    if isa(firstarg, Expr) && firstarg.head == :parameters
        return prunekwargs(drop(args, 1)..., firstarg.args...)
    else
        core = firstarg
        params = collect(drop(args, 1))
        for ex in params
            if isa(ex, Expr) && ex.head == :(=)
                ex.head = :kw
            end
        end
        if isa(core, Expr) && core.head == :kw
            core.head = :(=)
        end
        return core, params
    end
end

function hasevals(params)
    for p in params
        if isa(p, Expr) && p.head == :kw && first(p.args) == :evals
            return true
        end
    end
    return false
end

function collectvars(ex::Expr, vars::Vector{Symbol} = Symbol[])
    if ex.head == :(=)
        lhs = first(ex.args)
        if isa(lhs, Symbol)
            push!(vars, lhs)
        elseif isa(lhs, Expr) && lhs.head == :tuple
            append!(vars, lhs.args)
        end
    elseif (ex.head == :comprehension || ex.head == :generator)
        arg = ex.args[1]
        isa(arg, Expr) && collectvars(arg, vars)
    else
        for arg in ex.args
            isa(arg, Expr) && collectvars(arg, vars)
        end
    end
    return vars
end

function quasiquote!(ex::Expr, vars::Vector{Expr})
    if ex.head === :($)
        lhs = ex.args[1]
        rhs = isa(lhs, Symbol) ? gensym(lhs) : gensym()
        push!(vars, Expr(:(=), rhs, ex))
        return rhs
    elseif ex.head !== :quote
        for i in 1:length(ex.args)
            arg = ex.args[i]
            if isa(arg, Expr)
                ex.args[i] = quasiquote!(arg, vars)
            end
        end
    end
    return ex
end

@doc raw"""
    @benchmark <expr to benchmark> [setup=<setup expr>]

Run benchmark on a given expression.

# Example

The simplest usage of this macro is to put it in front of what you want
to benchmark.

```julia-repl
julia> @benchmark sin(1)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     13.610 ns (0.00% GC)
  median time:      13.622 ns (0.00% GC)
  mean time:        13.638 ns (0.00% GC)
  maximum time:     21.084 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     998
```

You can interpolate values into `@benchmark` expressions:

```julia
# rand(1000) is executed for each evaluation
julia> @benchmark sum(rand(1000))
BenchmarkTools.Trial:
  memory estimate:  7.94 KiB
  allocs estimate:  1
  --------------
  minimum time:     1.566 μs (0.00% GC)
  median time:      2.135 μs (0.00% GC)
  mean time:        3.071 μs (25.06% GC)
  maximum time:     296.818 μs (95.91% GC)
  --------------
  samples:          10000
  evals/sample:     10

# rand(1000) is evaluated at definition time, and the resulting
# value is interpolated into the benchmark expression
julia> @benchmark sum($(rand(1000)))
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     101.627 ns (0.00% GC)
  median time:      101.909 ns (0.00% GC)
  mean time:        103.834 ns (0.00% GC)
  maximum time:     276.033 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     935
```
"""
macro benchmark(args...)
    _, params = prunekwargs(args...)
    tmp = gensym()
    return esc(quote
        local $tmp = $BenchmarkTools.@benchmarkable $(args...)
        $BenchmarkTools.warmup($tmp)
        $(hasevals(params) ? :() : :($BenchmarkTools.tune!($tmp)))
        $BenchmarkTools.run($tmp)
    end)
end

function benchmarkable_parts(args)
    @nospecialize
    core, params = prunekwargs(args...)

    # extract setup/teardown if present, removing them from the original expression
    setup, teardown = nothing, nothing
    delinds = Int[]
    for i in eachindex(params)
        ex = params[i]
        if ex.args[1] == :setup
            setup = ex.args[2]
            push!(delinds, i)
        elseif ex.args[1] == :teardown
            teardown = ex.args[2]
            push!(delinds, i)
        end
    end
    deleteat!(params, delinds)

    if isa(core, Expr)
        quote_vars = Expr[]
        core = quasiquote!(core, quote_vars)
        if !isempty(quote_vars)
            setup = Expr(:block, setup, quote_vars...)
        end
    end

    return core, setup, teardown, params
end

"""
    @benchmarkable <expr to benchmark> [setup=<setup expr>]

Create a `Benchmark` instance for the given expression. `@benchmarkable`
has similar syntax with `@benchmark`. See also [`@benchmark`](@ref).
"""
macro benchmarkable(args...)
    core, setup, teardown, params = benchmarkable_parts(args)
    map!(esc, params, params)

    # extract any variable bindings shared between the core and setup expressions
    setup_vars = isa(setup, Expr) ? collectvars(setup) : []
    core_vars = isa(core, Expr) ? collectvars(core) : []
    out_vars = filter(var -> var in setup_vars, core_vars)

    # generate the benchmark definition
    return quote
        generate_benchmark_definition($__module__,
                                      $(Expr(:quote, out_vars)),
                                      $(Expr(:quote, setup_vars)),
                                      $(esc(Expr(:quote, core))),
                                      $(esc(Expr(:quote, setup))),
                                      $(esc(Expr(:quote, teardown))),
                                      Parameters($(params...)))
    end
end

# `eval` an expression that forcibly defines the specified benchmark at
# top-level in order to allow transfer of locally-scoped variables into
# benchmark scope.
#
# The double-underscore-prefixed variable names are not particularly hygienic - it's
# possible for them to conflict with names used in the setup or teardown expressions.
# A more robust solution would be preferable.
function generate_benchmark_definition(eval_module, out_vars, setup_vars, core, setup, teardown, params)
    @nospecialize
    corefunc = gensym("core")
    samplefunc = gensym("sample")
    type_vars = [gensym() for i in 1:length(setup_vars)]
    signature = Expr(:call, corefunc, setup_vars...)
    signature_def = Expr(:where, Expr(:call, corefunc,
                                  [Expr(:(::), setup_var, type_var) for (setup_var, type_var) in zip(setup_vars, type_vars)]...)
                    , type_vars...)
    if length(out_vars) == 0
        invocation = signature
        core_body = core
    elseif length(out_vars) == 1
        returns = :(return $(out_vars[1]))
        invocation = :($(out_vars[1]) = $(signature))
        core_body = :($(core); $(returns))
    else
        returns = :(return $(Expr(:tuple, out_vars...)))
        invocation = :($(Expr(:tuple, out_vars...)) = $(signature))
        core_body = :($(core); $(returns))
    end
    return Core.eval(eval_module, quote
        @noinline $(signature_def) = begin $(core_body) end
        @noinline function $(samplefunc)(__params::$BenchmarkTools.Parameters)
            $(setup)
            __evals = __params.evals
            __gc_start = Base.gc_num()
            __start_time = time_ns()
            __return_val = $(invocation)
            for __iter in 2:__evals
                $(invocation)
            end
            __sample_time = time_ns() - __start_time
            __gcdiff = Base.GC_Diff(Base.gc_num(), __gc_start)
            $(teardown)
            __time = max((__sample_time / __evals) - __params.overhead, 0.001)
            __gctime = max((__gcdiff.total_time / __evals) - __params.overhead, 0.0)
            __memory = Int(Base.fld(__gcdiff.allocd, __evals))
            __allocs = Int(Base.fld(__gcdiff.malloc + __gcdiff.realloc +
                               __gcdiff.poolalloc + __gcdiff.bigalloc,
                               __evals))
            return __time, __gctime, __memory, __allocs, __return_val
        end
        $BenchmarkTools.Benchmark($(samplefunc), $(params))
    end)
end

######################
# convenience macros #
######################

# These macros provide drop-in replacements for the
# Base.@time, Base.@elapsed macros and Base.@allocated, which use
# @benchmark but yield only the minimum time.

"""
    @belapsed expression [other parameters...]

Similar to the `@elapsed` macro included with Julia,
this returns the elapsed time (in seconds) to
execute a given expression.   It uses the `@benchmark`
macro, however, and accepts all of the same additional
parameters as `@benchmark`.  The returned time
is the *minimum* elapsed time measured during the benchmark.
"""
macro belapsed(args...)
    return esc(quote
        $BenchmarkTools.time($BenchmarkTools.minimum($BenchmarkTools.@benchmark $(args...)))/1e9
    end)
end

"""
    @ballocated expression [other parameters...]

Similar to the `@allocated` macro included with Julia,
this returns the number of bytes allocated when executing
a given expression.   It uses the `@benchmark`
macro, however, and accepts all of the same additional
parameters as `@benchmark`.  The returned allocations
correspond to the trial with the *minimum* elapsed time measured
during the benchmark.
"""
macro ballocated(args...)
    return esc(quote
        $BenchmarkTools.memory($BenchmarkTools.minimum($BenchmarkTools.@benchmark $(args...)))
    end)
end

"""
    @btime expression [other parameters...]

Similar to the `@time` macro included with Julia,
this executes an expression, printing the time
it took to execute and the memory allocated before
returning the value of the expression.

Unlike `@time`, it uses the `@benchmark`
macro, and accepts all of the same additional
parameters as `@benchmark`.  The printed time
is the *minimum* elapsed time measured during the benchmark.
"""
macro btime(args...)
    _, params = prunekwargs(args...)
    bench, trial, result = gensym(), gensym(), gensym()
    trialmin, trialallocs = gensym(), gensym()
    tune_phase = hasevals(params) ? :() : :($BenchmarkTools.tune!($bench))
    return esc(quote
        local $bench = $BenchmarkTools.@benchmarkable $(args...)
        $BenchmarkTools.warmup($bench)
        $tune_phase
        local $trial, $result = $BenchmarkTools.run_result($bench)
        local $trialmin = $BenchmarkTools.minimum($trial)
        local $trialallocs = $BenchmarkTools.allocs($trialmin)
        println("  ",
                $BenchmarkTools.prettytime($BenchmarkTools.time($trialmin)),
                " (", $trialallocs , " allocation",
                $trialallocs == 1 ? "" : "s", ": ",
                $BenchmarkTools.prettymemory($BenchmarkTools.memory($trialmin)), ")")
        $result
    end)
end


text:
import Base: view

mutable struct ViewableFunction
    callable::Function
end

function view(vf::ViewableFunction, idx...)
    return vf
end

function (vf::ViewableFunction)(args...)
    return vf.callable(args...)
end

# We only want to make *Functions* viewable.
# Everything else should either be
# (1) a mutable struct with `view` defined; or
# (2) something else that causes a problem.
function make_viewable(obj)
    return obj
end

function make_viewable(obj::Function)
    return ViewableFunction(obj)
end

# Function should be idempotent
function make_viewable(obj::ViewableFunction)
    return obj
end

max_stars_repo_name: dpmerrell/BatchMatFac

max_issues_repo_name: dpmerrell/BatchMatFac

max_forks_repo_name: dpmerrell/BatchMatFac


text:
function blockMatching(X::Array{Float32,3}, Par::PAR, Neighbor_arr::Array{Int32,2}, Num_arr::Array{Int32,1}, SelfIndex_arr::Array{Int32,1})::Array{Int32,2}

    L = length(Num_arr)

    Init_Index = zeros(Int32, Par.patnum[], L)

    for  i in 1 : L
        Patch = view(X, :, :, SelfIndex_arr[i])
        Nj = Num_arr[i]

        Dist = zeros(Float32,Nj)
        for j in 1:Nj
            @inbounds Neighbors = view(X,:,:,Neighbor_arr[j, i])
            Dist[j] = distance(Neighbors, Patch)
        end

        index = Int32.(sortperm(Dist, alg = QuickSort))

        Init_Index[:, i] .= Neighbor_arr[view(index,1:Par.patnum[]), i]
    end

    Init_Index
end

function distance(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
    @assert size(A) == size(B) "size of matrices does not match, got $(size(A)) and $(size(B))"
    out = zero(T)
    for k=1:size(B, 2)
        for l=1:size(B, 1)
            @inbounds x = A[l,k] - B[l,k]
            out += abs2(x)
        end
    end
    return out
end


text:
# Julian C bindings for Joyent's http-parser library.
# see: https://github.com/joyent/http-parser
#
isdefined(Base, :__precompile__) && __precompile__()

module HttpParser

depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(depsjl) ? include(depsjl) : error("HttpParser not properly ",
    "installed. Please run\nPkg.build(\"HttpParser\")")

using HttpCommon
using Compat
import Compat: String

import Base.show

# Export the structs and the C calls.
export Parser,
       ParserSettings,
       http_parser_init,
       http_parser_execute,
       http_method_str,
       http_should_keep_alive,
       upgrade,
       parse_url

# The id pool is used to keep track of incoming requests.
id_pool = 0

# A composite type that matches bit-for-bit a C struct.
type Parser
    # parser + flag = Single byte
    type_and_flags::Cuchar

    state::Cuchar
    header_state::Cuchar
    index::Cuchar

    nread::UInt32
    content_length::UInt64

    http_major::Cushort
    http_minor::Cushort
    status_code::Cushort
    method::Cuchar

    # http_errno + upgrade = Single byte
    errno_and_upgrade::Cuchar

    data::Any
    id::Int
end
Parser() = Parser(
    convert(Cuchar, 0),
    convert(Cuchar, 0),
    convert(Cuchar, 0),
    convert(Cuchar, 0),

    convert(UInt32, 0),
    convert(UInt64, 0),

    convert(Cushort, 0),
    convert(Cushort, 0),
    convert(Cushort, 0),
    convert(Cuchar, 0),

    convert(Cuchar, 0),

    convert(Ptr{UInt8}, C_NULL),
    (global id_pool += 1)
)

# Datatype Tuples for the different `cfunction` signatures for callback functions
const HTTP_CB      = (Int, (Ptr{Parser},))
const HTTP_DATA_CB = (Int, (Ptr{Parser}, Ptr{Cchar}, Csize_t,))

# A composite type that is expecting C functions to be run as callbacks.
type ParserSettings
    on_message_begin_cb::Ptr{Void}
    on_url_cb::Ptr{Void}
    on_status_complete_cb::Ptr{Void}
    on_header_field_cb::Ptr{Void}
    on_header_value_cb::Ptr{Void}
    on_headers_complete_cb::Ptr{Void}
    on_body_cb::Ptr{Void}
    on_message_complete_cb::Ptr{Void}
    on_chunk_header::Ptr{Void}
    on_chunk_complete::Ptr{Void}
end

ParserSettings(on_message_begin_cb, on_url_cb, on_status_complete_cb, on_header_field_cb, on_header_value_cb, on_headers_complete_cb, on_body_cb, on_message_complete_cb) = ParserSettings(on_message_begin_cb, on_url_cb, on_status_complete_cb, on_header_field_cb, on_header_value_cb, on_headers_complete_cb, on_body_cb, on_message_complete_cb, C_NULL, C_NULL)

function show(io::IO,p::Parser)
    print(io,"libhttp-parser: v$(version()), ")
    print(io,"HTTP/$(p.http_major).$(p.http_minor), ")
    print(io,"Content-Length: $(p.content_length)")
end

function version()
    ver = ccall((:http_parser_version, lib), Culong, ())
    major = (ver >> 16) & 255
    minor = (ver >> 8) & 255
    patch = ver & 255
    return VersionNumber(major, minor, patch)
end

"Intializes the Parser object with the correct memory."
function http_parser_init(parser::Parser,isserver=true)
    ccall((:http_parser_init, lib), Void, (Ptr{Parser}, Cint), &parser, !isserver)
end

"Run a request through a parser with specific callbacks on the settings instance."
function http_parser_execute(parser::Parser, settings::ParserSettings, request)
    ccall((:http_parser_execute, lib), Csize_t,
           (Ptr{Parser}, Ptr{ParserSettings}, Cstring, Csize_t,),
            Ref(parser), Ref(settings), convert(Cstring, pointer(request)), sizeof(request))
    if errno(parser) != 0
        throw(HttpParserError(errno(parser)))
    end
end

"Returns a string version of the HTTP method."
function http_method_str(method::Int)
    val = ccall((:http_method_str, lib), Cstring, (Int,), method)
    return unsafe_string(val)
end

# Is the request a keep-alive request?
function http_should_keep_alive(parser::Ptr{Parser})
    ccall((:http_should_keep_alive, lib), Int, (Ptr{Parser},), Ref(parser))
end

"Pauses the parser."
pause(parser::Parser) = ccall((:http_parser_pause,lib), Void, (Ptr{Parser}, Cint), Ref(parser), one(Cint))
"Resumes the parser."
resume(parser::Parser) = ccall((:http_parser_pause,lib), Void,(Ptr{Parser}, Cint), Ref(parser), zero(Cint))
"Checks if this is the final chunk of the body."
isfinalchunk(parser::Parser) = ccall((:http_parser_pause,lib), Cint, (Ptr{Parser},), Ref(parser)) == 1

upgrade(parser::Parser) = (parser.errno_and_upgrade & 0b10000000)>0
errno(parser::Parser) = parser.errno_and_upgrade & 0b01111111
errno_name(errno::Integer) = unsafe_string(ccall((:http_errno_name,lib),Cstring,(Int32,),errno))
errno_description(errno::Integer) = unsafe_string(ccall((:http_errno_description,lib),Cstring,(Int32,),errno))

immutable HttpParserError <: Exception
    errno::Int32
    HttpParserError(errno::Integer) = new(Int32(errno))
end

show(io::IO, err::HttpParserError) = print(io,"HTTP Parser Exception: ",errno_name(err.errno),"(",string(err.errno),"):",errno_description(err.errno))

@enum(UrlFields, UF_SCHEMA           = Cint(0),
                 UF_HOST             = Cint(1),
                 UF_PORT             = Cint(2),
                 UF_PATH             = Cint(3),
                 UF_QUERY            = Cint(4),
                 UF_FRAGMENT         = Cint(5),
                 UF_USERINFO         = Cint(6),
                 UF_MAX              = Cint(7))

immutable ParserUrl
    field_set::UInt16 # Bitmask of (1 << UF_*) values
    port::UInt16      # Converted UF_PORT string
    field_data::NTuple{Cint(UF_MAX)*2, UInt16}
    ParserUrl() = new(zero(UInt16), zero(UInt16), ntuple(i->zero(UInt16), Cint(UF_MAX)*2))
end

"Parse a URL"
function parse_url(url::AbstractString; isconnect::Bool = false)
    parsed = Dict{Symbol,AbstractString}()
    purl_ref = Ref(ParserUrl())
    res = ccall((:http_parser_parse_url, lib), Cint,
                 (Cstring, Csize_t, Cint, Ptr{ParserUrl}),
                  url, sizeof(url), Cint(isconnect), purl_ref)
    res > 0 && return parsed, -1
    purl = purl_ref[]

    for (i,uf) in enumerate(instances(UrlFields))
        !((purl.field_set & (1 << Cint(uf)) > 0) && uf != UF_MAX) && continue
        off = purl.field_data[2*(i-1)+1]
        len = purl.field_data[2*(i-1)+2]
        parsed[Symbol(uf)] = url[(off+1):(off+len)]
    end
    return parsed, purl.port
end

end # module HttpParser


text:
# Fallback functions for conjugates

posterior_canon(pri, G::IncompleteFormulation, x) = posterior_canon(pri, suffstats(G, x))
posterior_canon(pri, G::IncompleteFormulation, x, w) = posterior_canon(pri, suffstats(G, x, w))

posterior(pri::P, ss::SufficientStats) where {P<:Distribution} = Base.convert(P, posterior_canon(pri, ss))
posterior(pri::P, G::IncompleteFormulation, x) where {P<:Distribution} = Base.convert(P, posterior_canon(pri, G, x))
posterior(pri::P, G::IncompleteFormulation, x, w) where {P<:Distribution} = Base.convert(P, posterior_canon(pri, G, x, w))

posterior_rand(pri, ss::SufficientStats) = Base.rand(posterior_canon(pri, ss))
posterior_rand(pri, G::IncompleteFormulation, x) = Base.rand(posterior_canon(pri, G, x))
posterior_rand(pri, G::IncompleteFormulation, x, w) = Base.rand(posterior_canon(pri, G, x, w))

posterior_rand!(r::Array, pri, ss::SufficientStats) = Base.rand!(posterior_canon(pri, ss), r)
posterior_rand!(r::Array, pri, G::IncompleteFormulation, x) = Base.rand!(posterior_canon(pri, G, x), r)
posterior_rand!(r::Array, pri, G::IncompleteFormulation, x, w) = Base.rand!(posterior_canon(pri, G, x, w), r)

posterior_mode(pri, ss::SufficientStats) = Distributions.mode(posterior_canon(pri, ss))
posterior_mode(pri, G::IncompleteFormulation, x) = Distributions.mode(posterior_canon(pri, G, x))
posterior_mode(pri, G::IncompleteFormulation, x, w) = Distributions.mode(posterior_canon(pri, G, x, w))

fit_map(pri, G::IncompleteFormulation, x) = complete(G, pri, posterior_mode(pri, G, x))
fit_map(pri, G::IncompleteFormulation, x, w) = complete(G, pri, posterior_mode(pri, G, x, w))

posterior_randmodel(pri, G::IncompleteFormulation, x) = complete(G, pri, posterior_rand(pri, G, x))
posterior_randmodel(pri, G::IncompleteFormulation, x, w) = complete(G, pri, posterior_rand(pri, G, x, w))


text:
## R code 11.4
m11.1 <- quap(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a ,
        a ~ dnorm( 0 , 10 )
    ) , data=d )

## R code 11.7
m11.2 <- quap(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a + b[treatment] ,
        a ~ dnorm( 0 , 1.5 ),
        b[treatment] ~ dnorm( 0 , 10 )
    ) , data=d )
prior <- extract.prior( m11.2 , n=1e4 )
p <- sapply( 1:4 , function(k) inv_logit( prior$a + prior$b[,k] ) )

## R code 11.9
m11.3 <- quap(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a + b[treatment] ,
        a ~ dnorm( 0 , 1.5 ),
        b[treatment] ~ dnorm( 0 , 0.5 )
    ) , data=d )

## R code 11.11
m11.4 <- ulam(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a[actor] + b[treatment] ,
        a[actor] ~ dnorm( 0 , 1.5 ),
        b[treatment] ~ dnorm( 0 , 0.5 )
    ) , data=dat_list , chains=4 , log_lik=TRUE )

## R code 11.19
dat_list2 <- list(
    pulled_left = d$pulled_left,
    actor = d$actor,
    side = d$side,
    cond = d$cond )
m11.5 <- ulam(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a[actor] + bs[side] + bc[cond] ,
        a[actor] ~ dnorm( 0 , 1.5 ),
        bs[side] ~ dnorm( 0 , 0.5 ),
        bc[cond] ~ dnorm( 0 , 0.5 )
    ) , data=dat_list2 , chains=4 , log_lik=TRUE )

## R code 11.25
dat <- with( d_aggregated , list(
    left_pulls = left_pulls,
    treatment = treatment,
    actor = actor,
    side = side,
    cond = cond ) )

m11.6 <- ulam(
    alist(
        left_pulls ~ dbinom( 18 , p ) ,
        logit(p) <- a[actor] + b[treatment] ,
        a[actor] ~ dnorm( 0 , 1.5 ) ,
        b[treatment] ~ dnorm( 0 , 0.5 )
    ) , data=dat , chains=4 , log_lik=TRUE )

## R code 11.29
dat_list <- list(
    admit = d$admit,
    applications = d$applications,
    gid = ifelse( d$applicant.gender=="male" , 1 , 2 )
)
m11.7 <- ulam(
    alist(
        admit ~ dbinom( applications , p ) ,
        logit(p) <- a[gid] ,
        a[gid] ~ dnorm( 0 , 1.5 )
    ) , data=dat_list , chains=4 )

## R code 11.32
dat_list$dept_id <- rep(1:6,each=2)
m11.8 <- ulam(
    alist(
        admit ~ dbinom( applications , p ) ,
        logit(p) <- a[gid] + delta[dept_id] ,
        a[gid] ~ dnorm( 0 , 1.5 ) ,
        delta[dept_id] ~ dnorm( 0 , 1.5 )
    ) , data=dat_list , chains=4 , iter=4000 )

## R code 11.45
dat <- list(
    T = d$total_tools ,
    P = d$P ,
    cid = d$contact_id )

# intercept only
m11.9 <- ulam(
    alist(
        T ~ dpois( lambda ),
        log(lambda) <- a,
        a ~ dnorm(3,0.5)
    ), data=dat , chains=4 , log_lik=TRUE )

# interaction model
m11.10 <- ulam(
    alist(
        T ~ dpois( lambda ),
        log(lambda) <- a[cid] + b[cid]*P,
        a[cid] ~ dnorm( 3 , 0.5 ),
        b[cid] ~ dnorm( 0 , 0.2 )
    ), data=dat , chains=4 , log_lik=TRUE )

## R code 11.49
dat2 <- list( T=d$total_tools, P=d$population, cid=d$contact_id )
m11.11 <- ulam(
    alist(
        T ~ dpois( lambda ),
        lambda <- exp(a[cid])*P^b[cid]/g,
        a[cid] ~ dnorm(1,1),
        b[cid] ~ dexp(1),
        g ~ dexp(1)
    ), data=dat2 , chains=4 , log_lik=TRUE )

## R code 11.53
# compute the offset
d$log_days <- log( d$days )

# fit the model
m11.12 <- quap(
    alist(
        y ~ dpois( lambda ),
        log(lambda) <- log_days + a + b*monastery,
        a ~ dnorm( 0 , 1 ),
        b ~ dnorm( 0 , 1 )
    ), data=d )

## R code 11.56
code_m11.13 <- "
data{
    int N; // number of individuals
    int K; // number of possible careers
    int career[N]; // outcome
    vector[K] career_income;
}
parameters{
    vector[K-1] a; // intercepts
    real<lower=0> b; // association of income with choice
}
model{
    vector[K] p;
    vector[K] s;
    a ~ normal( 0 , 1 );
    b ~ normal( 0 , 0.5 );
    s[1] = a[1] + b*career_income[1];
    s[2] = a[2] + b*career_income[2];
    s[3] = 0; // pivot
    p = softmax( s );
    career ~ categorical( p );
}
"

## R code 11.59
code_m11.14 <- "
data{
    int N; // number of observations
    int K; // number of outcome values
    int career[N]; // outcome
    real family_income[N];
}
parameters{
    vector[K-1] a; // intercepts
    vector[K-1] b; // coefficients on family income
}
model{
    vector[K] p;
    vector[K] s;
    a ~ normal(0,1.5);
    b ~ normal(0,1);
    for ( i in 1:N ) {
        for ( j in 1:(K-1) ) s[j] = a[j] + b[j]*family_income[i];
        s[K] = 0; // the pivot
        p = softmax( s );
        career[i] ~ categorical( p );
    }
}
"

dat_list <- list( N=N , K=3 , career=career , family_income=family_income )
m11.14 <- stan( model_code=code_m11.14 , data=dat_list , chains=4 )

## R code 11.61
# binomial model of overall admission probability
m_binom <- quap(
    alist(
        admit ~ dbinom(applications,p),
        logit(p) <- a,
        a ~ dnorm( 0 , 1.5 )
    ), data=d )

# Poisson model of overall admission rate and rejection rate
# 'reject' is a reserved word in Stan, cannot use as variable name
dat <- list( admit=d$admit , rej=d$reject )
m_pois <- ulam(
    alist(
        admit ~ dpois(lambda1),
        rej ~ dpois(lambda2),
        log(lambda1) <- a1,
        log(lambda2) <- a2,
        c(a1,a2) ~ dnorm(0,1.5)
    ), data=dat , chains=3 , cores=3 )

m11.15 <- ulam(
    alist(
        days_to_event|adopted==1 ~ exponential( lambda ),
        days_to_event|adopted==0 ~ custom(exponential_lccdf( !Y | lambda )),
        lambda <- 1.0/mu,
        log(mu) <- a[color_id],
        a[color_id] ~ normal(0,1)
    ), data=dat , chains=4 , cores=4 )


text:
@with_kw struct model

    foodweb::Matrix{Int64}
    S::Int64
    basal::Int64
    nonbasal::Int64
    top::Int64
    nontop::Int64
    reorder::Vector{Int64}
end

function model(foodweb::Matrix{Int64})

    # Prune prey with only specialist prey
    foodweb = prune(foodweb)

    # Identify basal species
    nprey = sum(foodweb, dims = 1)[1, :]
    basal = findall(nprey .== 0)

    # Identify top predators
    npred = sum(foodweb, dims = 2)[:, 1]
    top = findall(npred .== 0)

    # Identify everyone else
    middle = findall(x -> !(x in basal) && !(x in top), 1:S)

    # Sort the species
    order = vcat(basal, middle, top)

    model(foodweb = foodweb[order, order],
          S = size(foodweb, 1),
          basal = length(basal),
          nonbasal = size(foodweb, 1) - length(basal),
          top = length(top),
          nontop = size(foodweb, 1) - length(top),
          reorder = order)

end

function computeLL(θ, m)

    loglik = eltype(θ[1])(0.0)

    @inbounds for pred in 1:m.nonbasal
        for prey in 1:m.nontop
            ϕ = 0.99 * exp(-((θ.n[prey] - θ.c[pred] * θ.n[m.basal + pred]) / (θ.r[pred] / 2)) ^ 2)
            loglik += logpdf(Bernoulli(ϕ), m.foodweb[prey, m.basal + pred])
        end
    end

    return loglik
end

function (m::model)(θ)

    logpost = computeLL(θ, m)

    logpost += logpdf.(Beta(2, 8), θ.r) |> sum

    isnan(logpost) && return -Inf

    # @show logpost.value
    # @show ForwardDiff.value.(θ.c)

    return logpost
end

max_stars_repo_name: clint-leach/niche-model

max_issues_repo_name: clint-leach/niche-model

max_forks_repo_name: clint-leach/niche-model


text:
using ExaTron

datafile = "./data/"*ARGS[1]
pq_val = parse(Float64, ARGS[2])
va_val = parse(Float64, ARGS[3])
max_iter = parse(Int, ARGS[4])
gpu = parse(Bool, ARGS[5])

println("Data: ", datafile, " rho_pq = ", pq_val, " rho_va = ", va_val, " max_iter = ", max_iter, " gpu = ", gpu)

# Warm-up
ExaTron.admm_rect_gpu(datafile;
                      iterlim=1, rho_pq=pq_val, rho_va=va_val, scale=1e-4, use_polar=true, use_gpu=gpu)

# Run
ExaTron.admm_rect_gpu(datafile;
                      iterlim=max_iter, rho_pq=pq_val, rho_va=va_val, scale=1e-4, use_polar=true, use_gpu=gpu)