text:
using Statistics
using StatsBase: sample, quantile
using EvoTrees
using EvoTrees: sigmoid, logit
using Plots
using GraphRecipes
using Random: seed!
using StaticArrays

# prepare a dataset
seed!(123)
X = rand(1_000, 2) .* 2
Y = sin.(X[:,1] .* π) .+ X[:,2]
Y = Y .+ randn(size(Y)) .* 0.1 #logit(Y)
# Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# make a grid
grid_size = 101
range = 2
X_grid = zeros(grid_size^2,2)
for j in 1:grid_size
    for i in 1:grid_size
        X_grid[grid_size*(j-1) + i,:] .= [(i-1) / (grid_size-1) * range,  (j-1) / (grid_size-1) * range]
    end
end
Y_grid = sin.(X_grid[:,1] .* π) .+ X_grid[:,2]

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(T=Float64,
    loss=:linear, metric=:mse,
    nrounds=100, nbins = 16,
    λ=0.0, γ=0.0, η=0.05,
    max_depth = 3, min_weight = 1.0,
    rowsample=0.8, colsample=1.0)

edges = EvoTrees.get_edges(X_train, params1.nbins)
X_bin = EvoTrees.binarize(X_train, edges)
@time model = model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)

# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 25, metric=:mae)
@time pred_train_linear = predict(model, X_train)
# @time pred_eval_linear = predict(model, X_eval)
# mean(abs.(pred_train_linear .- Y_train))
# sqrt(mean((pred_train_linear .- Y_train) .^ 2))

μ = mean(Y_train)

x_perm_1 = sortperm(X_train[:,1])
x_perm_2 = sortperm(X_train[:,2])
p1 = plot(X_train[:,1], Y_train, ms = 3, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var1"), yaxis = ("target"), leg = false, cbar=true)
p1_bin = plot(X_bin[:,1], Y_train, ms = 3, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msw=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("var1"), yaxis = ("target"), legend = false, cbar=true, label = "")
plot!(fill(μ, params1.nbins), lw=3, color="#66ffcc", background_color = "white", seriestype=:line, label="predict", leg=true)
# savefig(p1, "var1.svg")
# plot!(X_train[:,1][x_perm_1], pred_train_linear[x_perm_1], color = "red", mswidth=0, msize=3, label = "Linear", st=:scatter, leg=false)

p2 = plot(X_train[:,2], Y_train, ms = 3, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("target"), leg = false, cbar=true)
p2_bin = plot(X_bin[:,2], Y_train, ms = 3, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("target"), leg = false, cbar=true, label="")
plot!(fill(μ, params1.nbins), lw=3, color="#66ffcc", background_color = "white", seriestype=:line, label="predict", leg = true)
# savefig(p2, "var2.svg")
# plot!(X_train[:,2][x_perm_2], pred_train_linear[x_perm_2], color = "red", mswidth=0, msize=3, st=:scatter, label = "Predict")

p = plot(p1,p2, layout=(2,1))
savefig(p, "blog/raw_one_ways.svg")
savefig(p, "blog/raw_one_ways.png")
p = plot(p1_bin, p2_bin, layout=(2,1))
savefig(p, "blog/bin_one_ways.svg")
savefig(p, "blog/bin_one_ways.png")

# train iteration
# plot left vs right points
left_id = X_bin[:,2] .== 1
p = plot(X_bin[left_id,2], Y_train[left_id], ms = 3, color="darkred", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("residual"), leg = true, cbar=true, label="left")
plot!(X_bin[.!left_id,2], Y_train[.!left_id], ms = 3, color="#33ccff", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("residual"), leg = true, cbar=true, label="right")
plot!(fill(μ, params1.nbins), lw=3, color="#66ffcc", background_color = "white", seriestype=:line, label="predict", leg = true)

# residuals
residuals = Y_train .- μ
left_id = X_bin[:,2] .== 1
p = plot(X_bin[left_id,2], residuals[left_id], ms = 3, color="darkred", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("residual"), leg = true, cbar=true, label="left")
plot!(X_bin[.!left_id,2], residuals[.!left_id], ms = 3, color="#33ccff", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("var2"), yaxis = ("residual"), leg = true, cbar=true, label="right")
plot!(fill(0, params1.nbins), lw=3, color="#66ffcc", background_color = "white", seriestype=:line, label="", leg = true)
savefig(p, "blog/first_split.svg")
savefig(p, "blog/first_split.png")

left_res = residuals[left_id]
function loss(x::Vector, val)
    sum((x .- val).^2)
end
eval_pts = -2.1:0.01:0.5
left_loss = loss.(Ref(left_res), eval_pts)
p = plot(left_res, left_res .* 0, ms = 5, color="darkred", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("predict"), yaxis = ("loss"), leg = true, cbar=true, label="observation")
# plot!(left_res, left_loss, ms = 3, color="#33ccff", msw=0, background_color = "white", seriestype=:scatter, xaxis = ("residual"), yaxis = ("loss"), leg = true, cbar=true, label="")
plot!(eval_pts, left_loss, lw=3, color="#66ffcc", background_color = "white", seriestype=:line, label="loss", leg = true)
savefig(p, "blog/left_parabole.svg")
savefig(p, "blog/left_parabole.png")
loss(left_res, 0.0) - loss(left_res, mean(left_res))

# add single pt parabol
left_res_1 = sort(left_res)[35]
left_loss_1 = loss.(Ref([left_res_1]), eval_pts)
plot!(eval_pts, left_loss_1, lw=0.5, color="#33ccff", background_color = "white", seriestype=:line, label="", leg = true)

###########################
# raw compute
###########################
𝑖 = collect(1:size(X_train,1))
δ, δ² = zeros(SVector{model.K, Float64}, size(X_train, 1)), zeros(SVector{model.K, Float64}, size(X_train, 1))
𝑤 = zeros(SVector{1, Float64}, size(X_train, 1)) .+ 1
pred = zeros(SVector{model.K,Float64}, size(X_train,1)) .+= μ
EvoTrees.update_grads!(params1.loss, params1.α, pred, Y_train, δ, δ², 𝑤)
∑δ, ∑δ², ∑𝑤 = sum(δ[left_id]), sum(δ²[left_id]), sum(𝑤[left_id])
gain = EvoTrees.get_gain(params1.loss, ∑δ, ∑δ², ∑𝑤, params1.λ)

#####################################
# 3D visualisation of data
#####################################
# p = plot(X_train[:,1], X_train[:,2], Y_train, zcolor=Y_train, color=cgrad(["red","#3399ff"]), msize=5, markerstrokewidth=0, leg=false, cbar=true, w=1, st=:scatter)
p = plot(X_grid[:,1], X_grid[:,2], Y_grid, zcolor=Y_grid, color=cgrad(["#555555", "#eeeeee"]), msize=5, markerstrokewidth=0, leg=false, cbar=true, st=:scatter, xaxis="var1", yaxis="var2")
plot!(X_train[:,1], X_train[:,2], Y_train, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msize=4, markerstrokewidth=0, st=:scatter)
savefig(p, "blog/data_3D.svg")
savefig(p, "blog/data_3D.png")
plot(X_train[:,1], X_train[:,2], pred_train_linear, zcolor=Y_train, m=(5, 0.9, :rainbow, Plots.stroke(0)), leg=false, cbar=true, w=1, st=:scatter)
plot(X_train[:,1], X_train[:,2], pred_train_linear, zcolor=Y_train, st=[:surface], leg=false, cbar=true, fillcolor=:rainbow, markeralpha=1.0)

p_bin = plot(X_bin[:,1], X_bin[:,2], Y_train, zcolor=Y_train, color=cgrad(["darkred", "#33ccff"]), msize=4, markerstrokewidth=0, st=:scatter, leg=false, cbar=true)

gr()

params1 = EvoTreeRegressor(
    loss=:linear, metric=:mse,
    nrounds=100, nbins = 100,
    λ = 0.0, γ=0.0, η=0.5,
    max_depth = 2, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

anim = @animate for i=1:20
    params1.nrounds = (i-1)*5+1
    model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = Inf)
    pred_train_linear = predict(model, X_train)
    x_perm = sortperm(X_train[:,1])
    plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
    plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
end
gif(anim, "blog/anim_fps1.gif", fps = 1)

# tree vec
function treevec(tree)
    source, target = zeros(Int, max(1, length(tree.nodes)-1)), zeros(Int, max(1, length(tree.nodes)-1))
    count_s, count_t = 1, 1
    for i in 1:length(tree.nodes)
        if tree.nodes[i].split
            source[count_s] = i
            source[count_s+1] = i
            target[count_t] = tree.nodes[i].left
            target[count_t+1] = tree.nodes[i].right
            count_s += 2
            count_t += 2
        elseif i ==1
            source[i] = i
            target[i] = i
        end
    end
    return source, target
end

# plot tree
function nodenames(tree)
    names = []
    for i in 1:length(tree.nodes)
        if tree.nodes[i].split
            push!(names, "feat: " * string(tree.nodes[i].feat) * "\n< " * string(round(tree.nodes[i].cond, sigdigits=3)))
        else
            push!(names, "pred:\n" * string(round(tree.nodes[i].pred[1], sigdigits=3)))
        end
    end
    return names
end

tree1 = model.trees[2]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p1 = graphplot(source, target, method=:tree, names = nodes, linecolor=:brown, nodeshape=:hexagon, fontsize=8, fillcolor="#66ffcc")

tree1 = model.trees[3]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p2 = graphplot(source, target, method=:tree, names = nodes, linecolor=:brown, nodeshape=:hexagon, fontsize=8, fillcolor="#66ffcc")

tree1 = model.trees[50]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p3 = graphplot(source, target, method=:tree, names = nodes, linecolor=:brown, nodeshape=:hexagon, fontsize=8, fillcolor="#66ffcc")

tree1 = model.trees[90]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p4 = graphplot(source, target, method=:tree, names = nodes, linecolor=:brown, nodeshape=:hexagon, fontsize=8, fillcolor="#66ffcc")

p = plot(p1,p2,p3,p4)
savefig(p, "blog/tree_group.svg")
savefig(p1, "blog/tree_1.svg")
savefig(p, "blog/tree_group.png")
savefig(p1, "blog/tree_1.png")

################################
# animation on one way plot
################################
# prepare a dataset
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(
    loss=:linear, metric=:mse,
    nrounds=10, nbins = 32,
    λ = 0.0, γ=0.0, η=0.5,
    max_depth = 3, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)

anim = @animate for i=0:4
    params1.nrounds = i
    # i == 1 ? params1.η = 0.001 : params1.η = 0.5
    model = fit_evotree(params1, X_train, Y_train, X_eval = X_eval, Y_eval = Y_eval, print_every_n = Inf)
    pred_train_linear = predict(model, X_train)
    x_perm = sortperm(X_train[:,1])
    plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
    plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
end
gif(anim, "blog/anim_fps1.gif", fps = 1)

ps = [p1,p2,p3,p4]
anim = @animate for i=0:4
    # i == 1 ? params1.η = 0.001 : params1.η = 0.5
    # if i == 1
    #     p = plot(p1)
    # elseif i ==2
    #     p = plot(p1,p2)
    # end
    if i==0
        plot(foreground_color_subplot=:white)
    else
        plot(ps[1:i]...)
    end
end
gif(anim, "blog/anim_tree.gif", fps = 1)

tree1 = model.trees[2]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p1 = graphplot(source, target, method=:tree, names = nodes, edgecolor=:brown, nodeshape=:hexagon, fontsize=8, nodecolor="#66ffcc", nodestrokewidth=0)

tree1 = model.trees[3]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p2 = graphplot(source, target, method=:tree, names = nodes, edgecolor=:brown, nodeshape=:hexagon, fontsize=8, nodecolor="#66ffcc", nodestrokewidth=0)

tree1 = model.trees[4]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p3 = graphplot(source, target, method=:tree, names = nodes, edgecolor=:brown, nodeshape=:hexagon, fontsize=8, nodecolor="#66ffcc", nodestrokewidth=0)

tree1 = model.trees[5]
source, target = treevec(tree1)
nodes = nodenames(tree1)
seed!(1)
p4 = graphplot(source, target, method=:tree, names = nodes, edgecolor=:brown, nodeshape=:hexagon, fontsize=8, nodecolor="#66ffcc", nodestrokewidth=0)

p = plot(p1,p2,p3,p4)
savefig(p, "blog/tree_group.svg")
savefig(p1, "blog/tree_1.svg")
savefig(p, "blog/tree_group.png")
savefig(p1, "blog/tree_1.png")


text:
using CPLEX

struct MO_GAP_Inst
    m::Model
    a::Matrix{Int}
    b::Vector{Int}
    c1::Matrix{Int}
    c2::Matrix{Int}
end

function parseGAP(fname, solver=CplexSolver())
    f = open(fname)
    n,m = parse.(Int, split(readline(f)))

    a = zeros(Int, n, m)
    c = zeros(Int, n, m)

    for i = 1:n
        cij = Int[]
        while length(cij) < m
            append!(cij, parse.(Int, split(readline(f))))
        end
        @assert length(cij) == m
        c[i, :] .= cij
    end

    for i = 1:n
        aij = Int[]
        while length(aij) < m
            append!(aij, parse.(Int, split(readline(f))))
        end
        @assert length(aij) == m
        a[i, :] .= aij
    end

    b = Int[]
    while length(b) < n
        append!(b, parse.(Int, split(readline(f))))
    end

    c2 = reshape(shuffle(MersenneTwister(0), vec(c)), n, m)

    #m = m÷4

    mod = vModel(solver=solver)
    @variable(mod, x[1:n, 1:m], Bin)
    @addobjective(mod, Min, sum(c[i,j]*x[i,j] for i = 1:n, j = 1:m))
    @addobjective(mod, Min, sum(c2[i,j]*x[i,j] for i = 1:n, j = 1:m))
    @constraint(mod, [i=1:n], sum(a[i,j]*x[i,j] for j = 1:m) <= b[i])
    @constraint(mod, [j=1:m], sum(x[i,j] for i = 1:n) == 1)

    mod,a,b,c,c2
end

function parseGAP_orlib(fname, solver=CplexSolver(CPX_PARAM_SCRIND = 0))
    f = open(fname)
    nbinst = parse(Int, chomp(readline(f)))
    res = MO_GAP_Inst[]

    for _ = 1:nbinst

        n,m = parse.(Int, split(readline(f)))
        a = zeros(Int, n, m)
        c = zeros(Int, n, m)

        for i = 1:n
            cij = Int[]
            while length(cij) < m
                append!(cij, parse.(Int, split(readline(f))))
            end
            @assert length(cij) == m
            c[i, :] .= cij
        end

        for i = 1:n
            aij = Int[]
            while length(aij) < m
                append!(aij, parse.(Int, split(readline(f))))
            end
            @assert length(aij) == m
            a[i, :] .= aij
        end

        b = Int[]
        while length(b) < n
            append!(b, parse.(Int, split(readline(f))))
        end

        c2 = reshape(shuffle(MersenneTwister(0), vec(c)), n, m)
        #m = m÷4

        mod = vModel(solver=solver)
        @variable(mod, x[1:n, 1:m], Bin)
        @addobjective(mod, Min, sum(c[i,j]*x[i,j] for i = 1:n, j = 1:m))
        @addobjective(mod, Min, sum(c2[i,j]*x[i,j] for i = 1:n, j = 1:m))
        @constraint(mod, [i=1:n], sum(a[i,j]*x[i,j] for j = 1:m) <= b[i])
        @constraint(mod, [j=1:m], sum(x[i,j] for i = 1:n) == 1)

        push!(res, MO_GAP_Inst(mod,a,b,c,c2))
    end
    close(f)
    res
end

ext: jl


text:
module Mutation
    using ..Representation
    using Distributions, Random

    # Uncorrelated mutation one step size
    function uncorr_mutation_onestepsize(population::Representation.Population)
        t_value = generate_t_value(length(population.members))
        new_members = Organism[]
        for member in population.members
            σ′ = generate_new_sigma(member.σ, t_value)
            new_gauss_kernels = Tuple{Float64, Float64, Float64}[]
            for kernel in member.gauss_kernels
                new_kernel = tuple([mutate_param(param, σ′[1]) for param  in kernel]...)
                push!(new_gauss_kernels, new_kernel)
            end
            mutated_organism = Representation.Organism(new_gauss_kernels,
                                                        σ′,
                                                        member.α,
                                                        missing,
                                                        missing)
            push!(new_members, mutated_organism)
        end
        return Representation.Population(new_members)
    end
    # Generate tau for 1 step
    function generate_t_value(population_length)
        constant = 1 / population_length^(1/2)
        return rand(Normal(0.0, constant), 1)[1]
    end

    function mutate_param(param, σ′)
        return param + σ′ * rand(Normal(0.0, 1), 1)[1]
    end
    # Generate sigma one step
    function generate_new_sigma(σ, t_value)
        return σ * exp(t_value*rand(Normal(0.0, 1), 1)[1])
    end
    # Generate taus for n steps
    function generate_tau_values(population_length)
        τ = 1 / ((2 * population_length)^(1/2))
        τ′ = 1 / ((2 * (population_length^(1/2)))^(1/2))
        return τ, τ′
    end
    # Generate sigma n step
    function generate_new_sigma_nsteps(σ, τ, τ′)
        return σ * exp(τ * rand(Normal(0.0, 1), 1)[1] + τ′ * rand(Normal(0.0, 1), 1)[1])
    end
    # Uncorrelated mutation n step sizes
    function uncorr_mutation_n_stepsize(population::Representation.Population)
        τ, τ′ = generate_tau_values(length(population.members))
        new_members = Organism[]
        for member in population.members
            new_gauss_kernels = Tuple{Float64, Float64, Float64}[]
            new_sigmas = Float64[]
            for (index, kernel) in enumerate(member.gauss_kernels)
                new_kernel = Float64[]
                for (param, σ) in zip(kernel, member.σ[((index - 1) * 3) + 1:(index-1) * 3 + 3])
                    σ′ = generate_new_sigma_nsteps(σ, τ, τ′)
                    push!(new_kernel, mutate_param(param, σ′))
                    push!(new_sigmas, σ′)
                end
                push!(new_gauss_kernels, tuple(new_kernel...))
            end
            mutated_organism = Representation.Organism(new_gauss_kernels,
                                                        new_sigmas,
                                                        member.α,
                                                        missing,
                                                        missing)
            push!(new_members, mutated_organism)
        end
        return Representation.Population(new_members)
    end

end  # module Mutation


text:
using  FortranFiles: FString
import Random

struct CodegenTask
   jtype :: DataType
   ftype :: String
   sz    :: Dims
   var   :: String
end

Base.size(task::CodegenTask) = sizeof(task.jtype) * prod(task.sz)

function fdecl(task::CodegenTask)
   if task.sz == (1,)
      code = "$(task.ftype) :: $(task.var)"
   else
      code = "$(task.ftype), allocatable :: $(task.var)("
      code *= join( (":" for _ in task.sz), "," ) * ")"
   end
   return code
end

function fmkdata(task::CodegenTask)
   codes = String[]
   if task.sz != (1,)
      dims = join( ("$dim" for dim in task.sz), "," )
      push!(codes, "allocate($(task.var)($dims))")
   end
   push!(codes, "call random($(task.var))")
   return codes
end

function fwrdata(tasks::Vector{CodegenTask})
   vars = join((task.var for task in tasks), ",")
   code = "write(lun) $(vars)"
   return code
end

function jwrdata(tasks::Vector{CodegenTask})
   vars   = join((task.var for task in tasks), ", ")
   nbytes = sum(size(task) for task in tasks)
   codes  = [ "$(vars) = popfirst!(data)",
              "nwritten = write(f, $(vars))",
              "@test nwritten == $(nbytes)" ]
   return codes
end

function jrddata(tasks::Vector{CodegenTask})
   if isempty(tasks)
      return [ "read(f)" ]
   end
   specs = String[]
   tests = String[]
   for task in tasks
      spec = "$(task.jtype)"
      typ  = "$(task.jtype)"
      if task.sz != (1,)
         dims = join( ("$dim" for dim in task.sz), "," )
         mixdims = join( ( rand(["$dim", "Int64($dim)", "Int32($dim)", "Int16($dim)"]) for dim in task.sz), "," )
         if length(task.sz) == 1
            spec = rand([ "Array{$(typ)}(undef, $(dims))", "($spec, $(mixdims))" ])
         else
            spec = rand([ "Array{$(typ)}(undef, $(dims))", "($spec, $(mixdims))", "($spec, ($(mixdims)))" ])
         end
         typ = "Array{$(typ),$(length(task.sz))}"
      end
      push!(specs, spec)
      push!(tests, "@test typeof($(task.var)) == $(typ)")
      push!(tests, "@test sizeof($(task.var)) == $(size(task))")
   end
   vars = join((task.var for task in tasks), ", ")
   vartup = (length(tasks)==1) ? vars : "($(vars))"
   specstr = join(specs, ", ")
   codes = [ "$(vars) = read(f, $(specstr))",
             tests...,
             "push!(data, $(vartup))" ]
   return codes
end

function jfrdata(tasks::Vector{CodegenTask})
   if isempty(tasks)
      return [ "@fread f" ]
   end
   specs = String[]
   tests = String[]
   decls = String[]
   for task in tasks
      typ = "$(task.jtype)"
      if task.sz == (1,)
         spec = "$(task.var)::$(typ)"
      else
         dims = join( ("$dim" for dim in task.sz), "," )
         choice = rand(1:3)
         if choice==1
            # read into preallocated array
            decl = "$(task.var) = Array{$(typ)}(undef, $(dims))"
            push!(decls, decl)
            spec = "$(task.var)"
         elseif choice==2
            # read into freshly constructed array
            spec = "$(task.var)::Array{$(typ)}(undef, $(dims))"
         else
            # ditto, using short syntax (T,n,m) or (T,(n,m))
            spec = ( (length(task.sz)==1) || rand(Bool) ) ?
                   "$(task.var)::($(typ),$(dims))" :
                   "$(task.var)::($(typ),($(dims)))"
         end
         typ = "Array{$(typ),$(length(task.sz))}"
      end
      push!(specs, spec)
      push!(tests, "@test typeof($(task.var)) == $(typ)")
      push!(tests, "@test sizeof($(task.var)) == $(size(task))")
   end
   vars = join((task.var for task in tasks), ", ")
   vartup = (length(tasks)==1) ? vars : "($(vars))"
   specstr = join(specs, " ")
   codes = [ decls...,
             "@fread f $(specstr)",
             tests...,
             "push!(data, $(vartup))" ]
   return codes
end

function gencode(nscalar=5, narray=3, nstrlen=3; seed=1)
# will generate (8+nstrlen)*(nscalar+3*narray) CodegenTasks
   @static if VERSION >= v"0.7.0-beta2.171"
      Random.seed!(seed)
   else
      Random.srand(seed)
   end

   numtypes = [
      ( Int8 ,      "integer(kind=int8)"   ),
      ( Int16,      "integer(kind=int16)"  ),
      ( Int32,      "integer(kind=int32)"  ),
      ( Int64,      "integer(kind=int64)"  ),
      ( Float32,    "real(kind=real32)"    ),
      ( Float64,    "real(kind=real64)"    ),
      ( ComplexF32, "complex(kind=real32)" ),
      ( ComplexF64, "complex(kind=real64)" ) ]
   strtypes = [
      ( FString{n}, "character(len=$n)" ) for n in rand(1:200,nstrlen) ]
   types = vcat(numtypes, strtypes)

   dims0 = [ ( 1, ) for i=1:nscalar ]
   dims1 = [ ( rand(1:10_000), ) for i=1:narray ]
   dims2 = [ ( rand(1:100), rand(1:100) ) for i=1:narray ]
   dims3 = [ ( rand(1:10), rand(1:10), rand(1:10) ) for i=1:narray ]
   dims  = vcat(dims0, dims1, dims2, dims3)

   fdeclf = open("fdecl.f90", "w")
   fwrf = open("fwrite.f90", "w")
   jwrf = open("jwrite.jl", "w")
   jrdf = open("jread.jl", "w")
   jskf = open("jskip.jl", "w")
   jfrf = open("jfread.jl", "w")
   print(jwrf, "function writedata(f::FortranFile, data)\n")
   print(jrdf, "function readdata(f::FortranFile)\n   data = Any[]\n   @test eof(f) == false\n")
   print(jskf, "function skipdata(f::FortranFile)\n   data = Any[]\n")
   print(jfrf, "function freaddata(f::FortranFile)\n   data = Any[]\n")

   tasks = CodegenTask[]
   itask = 1
   for (jt,ft) in types, sz in dims
      var = "var$(itask)"
      push!( tasks, CodegenTask(jt, ft, sz, var) )
      itask += 1
   end

   for task in tasks
      println(fdeclf, "   ", fdecl(task))
      for code in fmkdata(task)
         println(fwrf, "   ", code)
      end
   end

   Random.shuffle!(tasks)
   taskgroups = Vector{CodegenTask}[]
   itask = 1
   while itask <= length(tasks)
      jtask = min(itask + rand(0:3), length(tasks))
      push!( taskgroups, tasks[itask:jtask] )
      itask = jtask + 1
   end

   for tg in taskgroups
      println(fwrf, "   ", fwrdata(tg))
      for line in jwrdata(tg); println(jwrf, "   ", line); end
      for line in jrddata(tg); println(jrdf, "   ", line); end
      for line in jrddata(tg[1:rand(0:end-1)]); println(jskf, "   ", line); end
      for line in jfrdata(tg); println(jfrf, "   ", line); end
   end

   print(jwrf, "end\n")
   print(jrdf, "   @test eof(f) == true\n   return data\nend\n")
   print(jskf, "   return data\nend\n")
   print(jfrf, "   return data\nend\n")
   close(fwrf)
   close(jwrf)
   close(jrdf)
   close(jskf)
   close(jfrf)
end

gencode(1,1,2,seed=12345678)


text:
module output
    using DataFrames, DataStructures

    function to_string(instance)
        sok_dict, sok_arr = DataStructures.OrderedDict(key=> value.data for (key, value) in instance._sok_dict), [item.data for item in instance._sok_arr]
        max_lengths = Dict(key => ceil(Int, maximum(insert!([length(string(item)) for item in sok_dict[key]], 1, length(key))) * 1.5) for key in keys(sok_dict))
        max_lengths["index"] = maximum([length(string(length(sok_arr))), length("index")])

        print_string = string(
            "[ ", instance.name, " | ", length(sok_dict), " × ", length(sok_arr), " ]\n",
            "| ", lpad("", max_lengths["index"], " "), " | ",
            join([string(lpad(key, max_lengths[key], " "), " | ") for key in Base.keys(sok_dict)]),
            "\n"
        )

        print_arr, arr_string = length(sok_arr) > 20 ? sok_arr[1:20] : sok_arr, ""
        for i in eachindex(print_arr)
            arr = print_arr[i]
            arr_string = string(arr_string, "| ", lpad(i, max_lengths["index"], " "), " |")
            for key in keys(arr)
                arr_string = string(arr_string, lpad(arr[key], max_lengths[key] + 1, " "), " | ")
            end
            arr_string = string(arr_string, "\n")
        end
        print_string = string(print_string, arr_string)[1:end - 1]

        return print_string
    end

    function to_dataframe(instance)
        return DataFrames.DataFrame(Dict(key=> value.data for (key, value) in instance._sok_dict))
    end

    function to_dict(instance, orientation)
        sok_dict, sok_arr = DataStructures.OrderedDict(key=> value.data for (key, value) in instance._sok_dict), [item.data for item in instance._sok_arr]

        if orientation == "Dict"
            result = sok_dict
        elseif orientation == "Array"
            result = sok_arr
        end

        return result
    end

    function to_csv(instance, csv_path)
        sok_dict, sok_arr = DataStructures.OrderedDict(key=> value.data for (key, value) in instance._sok_dict), [item.data for item in instance._sok_arr]

        try
            open(csv_path, "w") do f_csv
                csv_string = string(join(keys(sok_dict), ", "), "\n")
                for sok in sok_arr
                    csv_string = string(csv_string, join(values(sok), ", "), "\n")
                end
                write(f_csv, csv_string[1:end - 1])
            end
            result = true
        catch error
            result = false
        end

        return result
    end
end


text:
#=
conjugates.jl

25/08/2015
Adham Beyki, odinay@gmail.com

TODO:
    add @check_args
    remove the check in delitem!
=#

abstract Conjugate

#=
##################################
###### Gaussian1DGaussian1D ######
##################################
Gaussian1DGaussian1D constructs a Bayesian likelihood-prior pair with
a Gaussian1D likelihood with known variance, where its mean is endowed
with another Gaussian1D.

Parameters:
    m0      (1x1) mean hyper-parameter of the Gaussian1D prior
    v0      (1X1) variance hyper-parameter of the Gaussian1D prior
    vv      (1x1) fixed (known) variance of the likelihood

MODEL:
    mu | m0, v0         ~ Gaussian1D(m0, v0)
    x_1,...,x_n | mu,vv ~ Gaussian1D(mu, vv)
=#
type Gaussian1DGaussian1D <: Conjugate
    m0::Float64             # mean hyper-parameter of the prior
    v0::Float64             # variance hyper-parameter of the prior

    mu::Float64             # mean of the likelihood
    vv::Float64             # fixed variance of the likelihood

    nn::Int                 # number of data points the component contains

    Gaussian1DGaussian1D(m0::Float64, v0::Float64, mu::Float64, vv::Float64, nn::Int) = new(m0, v0, mu, vv, nn)
end

Gaussian1DGaussian1D(m0::Real, v0::Real, vv::Real) = Gaussian1DGaussian1D(
        convert(Float64, m0),
        convert(Float64, v0),
        convert(Float64, m0),
        convert(Float64, vv),
        0)
Base.deepcopy(me::Gaussian1DGaussian1D) = Gaussian1DGaussian1D(me.m0, me.v0, me.mu, me.vv, me.nn)

function Base.show(io::IO, me::Gaussian1DGaussian1D)
    println(io, "Gaussian1DGaussian1D conjugate component")
    println(io, "likelihood parameters: mu=", me.mu, ", vv=", me.vv)
    println(io, "prior parameters     : m0=", me.m0, ", v0=", me.v0)
    println(io, "number of data points: nn=", me.nn)
end

# adding/removing data points to/from the conjugate component
function additem!(me::Gaussian1DGaussian1D, xx::Float64)
    me.nn += 1
    me.mu = ((me.nn-1)/me.nn)*me.mu + xx/me.nn
end
function delitem!(me::Gaussian1DGaussian1D, xx::Float64)
    if me.nn < 1
        throw("there is no more datapoints in the conjugate component")
    elseif me.nn == 1
        me.mu = me.m0
        me.nn = 0
    else
        me.mu = me.mu*(me.nn/(me.nn-1)) - xx/(me.nn-1)
        me.nn -= 1
    end
end

function posterior(me::Gaussian1DGaussian1D)

    if me.nn>0
        denom = me.vv + me.nn*me.v0
        mu = (me.nn * me.v0 * me.mu + me.vv * me.m0) / denom
        vv = me.v0 * me.vv / denom
    else
        mu = me.m0
        vv = me.v0
    end

    Gaussian1D(mu, vv)
end

function logpredictive(me::Gaussian1DGaussian1D, xx::Float64)
    # log p(xx|x1,...,xn)

    post = posterior(me)
    mu = post.mu
    vv = post.vv

    dummy = vv + me.vv
    val = -(xx - mu)^2 / (2 * dummy) - 0.5 * log(2 * pi * dummy)

    val
end

loglikelihood(me::Gaussian1DGaussian1D, xx::Float64) = logpdf(posterior(me), xx)

##################################
###### MultinomialDirichlet ######
##################################
#=
MultinomialDirichlet constructs a Bayesian likelihood-prior pair with a Multinomial likelihood
where its probability vector is drawn from a symmetric Dirichlet prior.

Parameters:
    dd: (1x1)  Cardinality, i.e. length of the multinomial probability vector
    pp: (ddx1) Multinomial probability vector
    aa: (1x1)  Concentration parameter of the Dirichlet prior

    mm: (1x1)  total count in xx
    mi: (ddx1) individual counts
    nn: (1x1)  number of data points

Model:

    pp | aa                 ~ Dirichlet(aa/dd,...,aa/dd)
    x_1,...,x_n | mm, mi, pp ~ Multinomial(mm, mi, pp)
=#

type MultinomialDirichlet <: Conjugate
    dd::Int                               # Cardinality
    aa::Float64                           # Dirichlet concentration hyper-parameter
    pp::Vector{Float64}                   # Multinomial probability vector

    mm::Int                               # total count in data
    mi::Vector{Int}                       # individual count in data
    nn::Int                               # number of data points

    MultinomialDirichlet(dd::Int, aa::Float64, pp::Vector{Float64}, mm::Int, mi::Vector{Int}, nn::Int) = new(dd, aa, pp, mm, mi, nn)
end
MultinomialDirichlet(dd::Int, aa::Float64) = MultinomialDirichlet(dd, aa/dd, zeros(Float64, dd), 0, zeros(Int, dd), 0)
Base.deepcopy(me::MultinomialDirichlet) = MultinomialDirichlet(me.dd, me.aa, deepcopy(me.pp), me.mm, deepcopy(me.mi), me.nn)

# displaying the MultinomialDirichlet component
function Base.show(io::IO, me::MultinomialDirichlet)
    println(io, "MultinomialDirichlet component")
    println(io, "cardinality: dd=$(me.dd), Dirichlet prior parameter: aa=$(me.aa)")
    println(io, "data: mm=$(me.mm), nn=$(me.nn)")
end

# adding/removing data points to/from the component
function additem!(me::MultinomialDirichlet, xx::Int)
    me.mi[xx] += 1
    me.mm += 1
    me.nn += 1
end
function delitem!(me::MultinomialDirichlet, xx::Int)
    me.mi[xx] -= 1
    me.mm -= 1
    me.nn -= 1
end
function additem!(me::MultinomialDirichlet, xx::Sent)
    xx_len = length(xx.w)
    @inbounds for ii = 1:xx_len
        me.mi[xx.w[ii]] += xx.c[ii]
        me.mm += xx.c[ii]
    end
    me.nn += 1
end
function delitem!(me::MultinomialDirichlet, xx::Sent)
    xx_len = length(xx.w)
    @inbounds for ii = 1:xx_len
        me.mi[xx.w[ii]] -= xx.c[ii]
        me.mm -= xx.c[ii]
    end
    me.nn -= 1
end

logpredictive(me::MultinomialDirichlet, xx::Int) = log((me.aa + me.mi[xx]) / (me.aa*me.dd + me.mm))
function logpredictive(me::MultinomialDirichlet, xx::Sent)
    xx_len = length(xx.w)
    lll = me.aa + me.mi
    mm = 0
    mi = zeros(Int, xx_len)
    @inbounds for ii = 1:xx_len
        lll[xx.w[ii]] += xx.c[ii]

        mm += xx.c[ii]
        mi[ii] = xx.c[ii]
    end

    ll = lgamma(mm+1) - sum(lgamma(mi+1)) + lgamma(me.aa*me.dd + me.mm) -
            lgamma(me.aa*me.dd + me.mm + mm) +
            sum(lgamma(lll) - lgamma(me.aa+me.mi))

    ll
end

function posterior(me::MultinomialDirichlet)
    aa = me.mi + me.aa
    return Dirichlet(aa)
end

loglikelihood(me::MultinomialDirichlet, xx::Int) = mean(posterior(me))[xx]
function loglikelihood(me::MultinomialDirichlet, xx::Sent)
    ll = 0.0
    alpha = mean(posterior(me))
    @inbounds for i=1:length(xx)
        ll += xx.c[i] * alpha[xx.w[i]]
    end
    ll
end


text:
module TestStaticArrays

include("preamble.jl")
using BangBang: implements
using InitialValues: Init
using StaticArrays:
    MArray, MMatrix, MVector, SArray, SMatrix, SVector, StaticArrays, StaticVector

@testset "implements(push!, _)" begin
    @test !implements(push!, SVector(0))
    @test !implements(push!, SMatrix{1, 1}(0))
    @test !implements(push!, MVector(0))
    @test !implements(push!, MMatrix{1, 1}(0))
end

@testset "implements(setindex!, _)" begin
    @test !implements(setindex!, SVector(0))
    @test !implements(setindex!, SMatrix{1, 1}(0))
    @test implements(setindex!, MVector(0))
    @test implements(setindex!, MMatrix{1, 1}(0))
end

@testset "push!!" begin
    @test push!!(SVector(1, 2, 3), 4) === SVector(1, 2, 3, 4)
    if versionof(StaticArrays) >= v"0.12"  # not sure when it was fixed
        @test push!!(MVector(1, 2, 3), 4) ==ₜ MVector(1, 2, 3, 4)
    else
        # For StaticArrays 0.8
        @test_broken push!!(MVector(1, 2, 3), 4) ==ₜ MVector(1, 2, 3, 4)
        @test push!!(MVector(1, 2, 3), 4)::StaticVector == [1, 2, 3, 4]
    end
end

@testset "append!! with Init(append!!)" begin
    @test append!!(Init(append!!), SVector(1, 2)) === SVector(1, 2)
    @test append!!(SVector(1, 2), Init(append!!)) === SVector(1, 2)
    @test append!!(Init(append!!), MVector(1, 2)) ==ₜ MVector(1, 2)
    xs = MVector(1, 2)
    @test append!!(xs, Init(append!!)) === xs
end

@testset "setindex!!(SArray, ...)" begin
    @test setindex!!(SVector(1, 2, 3), 200, 2) === SVector(1, 200, 3)
    @test_broken setindex!!(SVector(1, 2, 3), 0.2, 2) === SVector(1.0, 0.2, 3.0)
end

@testset "setindex!!(MArray, ...)" begin
    xs = MVector(1, 2, 3)
    @test setindex!!(xs, 200, 2) === xs ==ₜ MVector(1, 200, 3)
    @test_broken setindex!!(xs, 0.2, 2) === xs ==ₜ MVector(1.0, 0.2, 3.0)
end

@testset "empty!!" begin
    @test empty!!(SVector(1)) === SVector{0,Int}()
    @test isempty(empty!!(MVector(1)))
    if versionof(StaticArrays) >= v"0.12"  # not sure when it was fixed
        @test empty!!(MVector(1)) ==ₜ MVector{0,Int}()
    else
        # For StaticArrays 0.8
        @test_broken empty!!(MVector(1)) ==ₜ MVector{0,Int}()
    end
end

@testset "Empty" begin
    @test push!!(Empty(SVector), 1) === SVector(1)
    @test push!!(Empty(SArray), 1) === SVector(1)
    @test push!!(Empty(MVector), 1) ==ₜ MVector(1)
    @test push!!(Empty(MArray), 1) ==ₜ MVector(1)
end

end  # module


text:
function drectangle(
    p,
    x1::T,
    x2::T,
    y1::T,
    y2::T,
)::T where {T<:Float64}
    return -min(min(min(-y1 + p[2], y2 - p[2], -x1 + p[1], x2 - p[1])))
end

function dcircle(
    p,
    xc::T,
    yc::T,
    r::T
)::T where {T<:Float64}
    return sqrt((p[1] - xc) .^ 2 + (p[2] - yc) .^ 2) - r
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
# functions related to noncentral F distribution

import .RFunctions:
    nfdistpdf,
    nfdistlogpdf,
    nfdistcdf,
    nfdistccdf,
    nfdistlogcdf,
    nfdistlogccdf,
    nfdistinvcdf,
    nfdistinvccdf,
    nfdistinvlogcdf,
    nfdistinvlogccdf


text:
using FooBar
using Test

@testset "FooBar.jl" begin
    # Write your tests here.
end


text:
function u_fcm()
end

function c_fcm(memberships::Vector{<:Real}, patterns::Vector{<:Vector{<:Real}}; m::Real)
    N = length(patterns)
    d = length(patterns[1])

    numerator = zeros(d)
    denominator = 0.0
    for j = 1:N
        w = (memberships[j] ^ m)
        numerator += w * patterns[j]
        denominator += w
    end
    c = numerator / denominator
    c
end

function euclidean_distance(a, b; squared::Bool=false)
    d_squared = sum((a - b).^2)
    return squared ? d_squared : sqrt(d_squared)
end


text:
function process_data(path)
    if typeof(path) != String
    	(network, bus_shift)=process_raw(true, path[1]);
	network=process_rop(true, path[2], network, bus_shift);
	return network
    elseif typeof(path) == String
    	data = PowerModels.parse_file(path);
    	(network, bus_shift)=process_raw(data);
    	return network
    else
    	println("The data format is not recognized or mismatched")
    end
end

function update_format(network)
    bus_index=Dict{String,Any}();
    series = 0
    for i in keys(network["bus"]) series += 1
        bus_index=merge(bus_index, Dict(i=>series));
        network["bus"][i]=merge(network["bus"][i], Dict("index"=>series));
        network["bus"][i]=merge(network["bus"][i], Dict("gen"=>false));
        network["bus"][i]=merge(network["bus"][i], Dict("Pd"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("Qd"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("GL"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("BL"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("sh_max"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("sh_min"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("bcs"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("gfM"=>0.0));
        network["bus"][i]=merge(network["bus"][i], Dict("bfM"=>0.0));
        if network["bus"][i]["bus_type"] == 3
        	network=merge(network, Dict("slack"=>parse(Int,i)))
        end
    end
    for i in keys(network["load"])
        I_index=network["load"][i]["load_bus"];
        I=bus_index[string(I_index)];
        ST=network["load"][i]["status"];
        PL=ST*(network["load"][string(i)]["pd"])+network["bus"][string(I_index)]["Pd"];
        QL=ST*(network["load"][string(i)]["qd"])+network["bus"][string(I_index)]["Qd"];
        network["bus"][string(I_index)]=merge(network["bus"][string(I_index)], Dict("Pd"=>PL, "Qd"=>QL));
        #network["bus"][string(I)]=merge(net_bus[string(I)], Dict("Qd"=>QL));
    end
    for i in keys(network["shunt"])
        I_index=network["shunt"][i]["shunt_bus"];
        I=bus_index[string(I_index)];
        ST=network["shunt"][i]["status"];
        GL=ST*(network["shunt"][string(i)]["gs"])+network["bus"][string(I_index)]["GL"];
        BL=ST*(network["shunt"][string(i)]["bs"])+network["bus"][string(I_index)]["BL"];
        network["bus"][string(I_index)]=merge(network["bus"][string(I_index)], Dict("GL"=>GL, "BL"=>BL));
    end
    for i in keys(network["gen"])
        I_index=network["gen"][string(i)]["gen_bus"];
        I=bus_index[string(I_index)];
        network["gen"][i]=merge(network["gen"][i], Dict("gen_bus_index"=>I));
        network["gen"][i]=merge(network["gen"][i], Dict("gen_bus"=>I_index));
        network["gen"][i]=merge(network["gen"][i], Dict("gov_droop"=>0.0));
        network["gen"][i]=merge(network["gen"][i], Dict("total_apf"=>0.0));
        network["gen_cost_order"] = max(length(network["gen"][string(i)]["cost"]),network["gen_cost_order"]);
    end
    for i in keys(network["branch"])
	I_index=network["branch"][string(i)]["f_bus"];
	I=bus_index[string(I_index)];
	J_index=network["branch"][string(i)]["t_bus"];
	J=bus_index[string(J_index)];
	network["branch"][i]=merge(network["branch"][i], Dict("f_bus_index"=>I));
	network["branch"][i]=merge(network["branch"][i], Dict("t_bus_index"=>J));
	network["branch"][i]=merge(network["branch"][i], Dict("gfM"=>0.0));
	network["branch"][i]=merge(network["branch"][i], Dict("bfM"=>0.0));
    end
    return(network,bus_index)
end

function process_raw(network)
    network=merge(network, Dict("nbal"=>Dict{String,Any}(), "sshunt"=>Dict{String,Any}(), "gen_cost_order"=>0));
    network["sshunt"]=merge(network["sshunt"], Dict("bus_index"=>zeros(Int32,0)));

    (network, bus_shift) = update_format(network);

    for i=1:length(network["bus"])
        network["nbal"]=merge(network["nbal"], Dict(string(i)=>Dict()));
        network["nbal"][string(i)]=merge(network["nbal"][string(i)], Dict("gen"=>zeros(Int32,0)), Dict("Sij"=>zeros(Int32,0)), Dict("Sji"=>zeros(Int32,0)), Dict("shunt"=>zeros(Int32,0)));
    end

    for g=1:length(network["gen"]) gen=network["gen"][string(g)];
        if (gen["gen_bus_index"]*gen["gen_status"]>0)
        push!(network["nbal"][string(gen["gen_bus_index"])]["gen"],g);
        end
    end

    for b=1:length(network["branch"]) br=network["branch"][string(b)];
        if (br["f_bus_index"]*br["br_status"]>0) push!(network["nbal"][string(br["f_bus_index"])]["Sij"], b); end
        if (br["t_bus_index"]*br["br_status"]>0) push!(network["nbal"][string(br["t_bus_index"])]["Sji"], b); end
    end

    return(network, bus_shift)
end


text:
import Measurements
@testset "Measurements+Optim" begin
    #example problem in #823
    f(x)=(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    xmes=zeros(Measurements.Measurement{Float64},2)
    xfloat = zeros(2)
    resmes = optimize(f,xmes)
    resfloat = optimize(f,xfloat)
    #given an initial value, they should give the exact same answer
    @test all(Optim.minimizer(resmes) .|> Measurements.value .== Optim.minimizer(resfloat))
    @test Optim.minimum(resmes) .|> Measurements.value .== Optim.minimum(resfloat)

end


text:
using LinearAlgebra: diagm, lmul!

export raising, lowering, number, X, Y, X_Y, XY, flip_flop,
       decay, dephasing, dipole_drive, parametric_drive, rwa_dipole

######################################################
# Primitives
######################################################

# ladder operators
"""
    raising(q::QSystem, ϕ::Real=0.0)

raising/creation/a ladder operator. Given an eignenstate of the number operator `|n⟩` on a
`QSystem q`, `raising(q)|n⟩ = √(n+1)|n+1⟩`. Optionally an additional phase exp(2πiϕ) is applied.

# Examples
```jldoctest
julia> q = DuffingTransmon("q", 3, DuffingSpec(5, -0.2))
DuffingTransmon("q", 3, DuffingSpec(5, -0.2))

julia> raising(q)
3×3 Array{Float64,2}:
 0.0  0.0      0.0
 1.0  0.0      0.0
 0.0  1.41421  0.0
```
"""
raising(q::QSystem, ϕ::Real=0.0) = raising(dimension(q), ϕ)

# custom diagm because it's twice as fast for small matrices
function raising(dim::Integer,  ϕ::Real=0.0)
    m = zeros(typeof(complex(ϕ)), dim, dim)
    fac = exp(1im * 2π * ϕ)
    for i in 1:(dim -1 )
        m[i + 1, i] = sqrt(i) * fac
    end
    return m
end

"""
    lowering(q::QSystem, ϕ::Real=0.0, factor::Real=1))

lowering/destruction/a† ladder operator. Given an eignenstate of the number operator `|n⟩` on a
`QSystem q`, `lowering(q)|n⟩ = √(n+1)|n+1⟩`. Optionally apply an additional phase ϕ and scaling `factor`.

# Examples
```jldoctest
julia> q = DuffingTransmon("q", 3, DuffingSpec(5, -0.2))
DuffingTransmon("q", 3, DuffingSpec(5, -0.2))

julia> lowering(q)
3×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.41421
 0.0  0.0  0.0
```
"""
lowering(q::QSystem, ϕ::Real=0, factor::Real=1) = lowering(dimension(q), ϕ, factor)

# custom diagm because it's twice as fast for small matrices
function lowering(dim::Integer,  ϕ::Real=0, factor::Real=1)
    m = zeros(typeof(complex(float(ϕ))), dim, dim)
    efac = exp(-1im * 2π * ϕ) * factor
    for i in 1:(dim - 1)
        m[i, i + 1] = sqrt(i) * efac
    end
    return m
end

"""
    number(q::QSystem, factor::Real=1)

Number operator (a†a) on a QSystem `q` with an optional scaling factor which gives a dephasing
matrix if `factor != 1`.
"""
number(q::QSystem, factor::Real=1) = number(dimension(q), factor)

# custom diagm because it's twice as fast for small matrices
function number(dim::Integer, factor::T=1) where {T<:Real}
    m = zeros(complex(float(T)), dim, dim)
    for i in 2:dim
        m[i, i] = (i - 1) * factor
    end
    return m
end

"""
    X(q::QSystem)

X (a + a†) operator on a QSystem `q`.
"""
function X(q::QSystem)
    # Rather than just returning raising(q) + lowering(q) we specialize this implementation because
    # it more than twice as fast with less than half the allocations. See a4add5ca0 for details.
    diag_elements = sqrt.(1:(dimension(q)-1))
    diagm(-1 => diag_elements, +1 => diag_elements)
end

"""
Apply a phase ϕ (units of τ) to an X operator

# Examples
```jldoctest
julia> q = DuffingTransmon("q", 2, DuffingSpec(5, -0.2))
DuffingTransmon("q", 2, DuffingSpec(5, -0.2))

julia> X(q, 0.25)
2×2 Array{Complex{Float64},2}:
         0.0+0.0im  6.12323e-17-1.0im
 6.12323e-17+1.0im          0.0+0.0im

julia> X(q, 0.25) ≈ Y(q)
true
```
"""
X(q::QSystem, ϕ::Real) = raising(q, ϕ) + lowering(q, ϕ)

X(qs::Vector{<:QSystem}, ϕs::Vector{<:Real}) = reduce(⊗, [X(q, ϕ) for (q, ϕ) in zip(qs, ϕs)])
X(qs::Vector{<:QSystem}) = reduce(⊗, [X(q) for q in qs])

"""
    Y(q::QSystem)

Y (-ia + ia†) operator on a QSystem `q`.
"""
function Y(q::QSystem)
    diag_elements = sqrt.(1:(dimension(q)-1))
    diagm(-1 => 1im*diag_elements, +1 => -1im*diag_elements)
end

Y(q::QSystem, ϕ::Real) = 1im*(raising(q, ϕ) - lowering(q, ϕ))

Y(qs::Vector{<:QSystem}, ϕs::Vector{<:Real}) = reduce(⊗, [Y(q, ϕ) for (q, ϕ) in zip(qs, ϕs)])
Y(qs::Vector{<:QSystem}) = reduce(⊗, [Y(q) for q in qs])

X_Y(qs::Vector{<:QSystem}, ϕs::Vector{<:Real}) = X(qs, ϕs) + Y(qs, ϕs)
X_Y(qs::Vector{<:QSystem}) = X(qs) + Y(qs)

"""
    flip_flop(a::QSystem, b::QSystem)

Bilinear photon exchange Hamiltonian on `a` and `b`: `ab† + a†b`. For qubits corresponds to Pauli operator
`σ⁺σ⁻ + σ⁻σ⁺`
"""
flip_flop(a::QSystem, b::QSystem) = raising(a)⊗lowering(b) + lowering(a)⊗raising(b)

"""
    flip_flop(a::QSystem, b::QSystem, ϕ::Real)

Apply an additional relative phase ϕ: `exp(2πiϕ)ab† + exp(-2πiϕ)a†b`
"""
flip_flop(a::QSystem, b::QSystem, ϕ::Real) = exp(-1im*2π*ϕ)*raising(a)⊗lowering(b) + exp(1im*2π*ϕ)*lowering(b)⊗raising(b)

"""
    XY(a::QSystem, b::QSystem)

Bilinear XY Hamiltonian on `a` and `b` which is (up to a scale) equivalent to a "flip-flop"
Hamiltonian. XY(a,b) = XᵃXᵇ + YᵃYᵇ = 2*(ab† + a†b)
"""
XY(a::QSystem, b::QSystem) = lmul!(2.0, flip_flop(a,b))

"""
    XY(a::QSystem, b::QSystem, ϕ::Real)

Apply an additional phase rotation to the XY Hamiltonian 2*(exp(2πiϕ)ab† + exp(-2πiϕ)a†b)
"""
XY(a::QSystem, b::QSystem, ϕ::Real) =  lmul!(2.0, flip_flop(a,b,ϕ))

"""
    decay(qs::QSystem, γ:Real)

T1 decay for a QSystem.

## args
* `qs`: a QSystem.
* `γ`: a decay rate in frequency units. Note T1 = 1/(2πγ).

## returns
The lindblad operator for decay.
"""
function decay(qs::QSystem, γ::Real)
    return lowering(qs, 0, sqrt(γ))
end

"""
    dephasing(qs::QSystem, γ::Real)

Dephasing for a QSystem.

## args
* `qs`: a QSystem.
* `γ`: a decay rate in frequency units. Note Tϕ = 1/(2πγ).

## returns
The lindblad operator for dephasing.
"""
function dephasing(qs::QSystem, γ::Real)
    return number(qs, sqrt(2γ))
end

"""
    dipole_drive(qs::QSystem, drive::Function, rotation_rate::Real=0.0)

Given some function of time, return a function applying a time dependent
dipole Hamiltonian. Note that this does not use the rotating wave approximation
and therefore requires a real valued drive. See also `rwa_dipole`

## args
* `qs`: a QSystem.
* `drive`: a function of time returning a real value.
* `rotation_rate`: the rotation rate of a rotating frame.

## returns
A function of time.
"""
function dipole_drive(qs::QSystem, drive::Function, rotation_rate::Real=0.0)
    function ham(t)
        pulse::Real = drive(t)
        return pulse * X(qs, rotation_rate * t)
    end
    return ham
end

"""
    rwa_dipole(qs::QSystem, drive::Function)

Given some function of time, return a function applying a time dependent
dipole Hamiltonian under the rotating wave approximation (RWA).

## args
* `qs`: a QSystem.
* `drive`: a function of time returning a real or complex value. The real
    part couples to X and the imaginary part couples to Y.

## returns
A function of time.
"""
function rwa_dipole(qs::QSystem, drive::Function)
    x_ham = X(qs)
    y_ham = Y(qs)
    function ham(t)
        pulse = drive(t)
        return real(pulse) * x_ham + imag(pulse) * y_ham
    end
    return ham
end

"""
    parametric_drive(qs::QSystem, drive::Function)

Given some function of time, return a function applying a
time dependent Hamiltonian.

## args
* `qs`: a QSystem with a method of `hamiltonian` accepting a function of time.
* `drive`: a function of time returning a real value.

## returns
A function of time.
"""
function parametric_drive(qs::QSystem, drive::Function)
    ham(t) = hamiltonian(qs, drive(t))
    return ham
end


text:
push!(LOAD_PATH,"/home/michel/git/PIPS.jl/src")
using Revise
using PIPS
using Ipopt
using Test
using MPI

MPI.Init()

@testset "C API PIPS" begin
    # First of all, test that hs071 example works
    include("hs071_test.jl")
end

@testset "C API Ipopt" begin
    # First of all, test that hs071 example works
    include("hs071_test_ipopt.jl")
end

@testset "MathOptInterface Ipopt" begin
    include("MOI_Ipopt_wrapper.jl")
end

@testset "MathOptInterface PIPS" begin
    include("MOI_wrapper.jl")
end
MPI.Finalize()


text:
using HCIToolbox
using StableRNGs
using Statistics
using Test

rng = StableRNG(8799)

@testset "HCIToolbox.jl" begin
    include("morphology.jl")
    include("masking.jl")
    include("angles.jl")
    include("inject.jl")
    include("scaling.jl")
    include("geometry.jl")
    include("profiles.jl")
end


text:
# ---------------- #
# Helper functions #
# ---------------- #
import Base: *

# ckronx.m -- DONE
function ckronx(b::AbstractMatrix{TM}, c::Array,
                ind::AbstractArray{Int}=1:length(b)) where TM<:AbstractMatrix
    d = length(ind)  # 26
    n = Array{Int}(d)  # 27
    for i=1:d  # 28
        n[i] = size(b[ind[i]], 2)
    end

    if prod(n) != size(c, 1)  # 29-31
        msg = "b and c are not conformable (b suggests size(c, 1) should be $(prod(n)))"
        error(msg)
    end

    z = c'  # 32
    mm = 1  # 33
    for i=1:d
        m = Int(length(z) / n[i])  # 35
        z = reshape(z, m, n[i])  # 36
        z = b[ind[i]] * z'  # 37
        mm = mm * size(z, 1)  # 38
    end
    reshape(z, mm, size(c, 2))  # 40
end

# dprod.m  - DONE
function row_kron(A::AbstractMatrix{S}, B::AbstractMatrix{T}) where {S,T}
    out = _allocate_row_kron_out(A, B)
    row_kron!(out, A, B)
end

function row_kron!(out::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
    # get input dimensions
    nobsa, na = size(A)
    nobsb, nb = size(B)

    nobsa != nobsb && error("A and B must have same number of rows")

    # fill in each element. To do this we make sure we access each array
    # consistent with its column major memory layout.
    @inbounds for ia=1:na, ib=1:nb, t=1:nobsa
        out[t, nb*(ia-1) + ib] = A[t, ia] * B[t, ib]
    end
    out
end

function row_kron!(out::SparseMatrixCSC, A::SparseMatrixCSC, B::SparseMatrixCSC)
    colptr = out.colptr
    rowval = out.rowval
    nzval = out.nzval

    na = size(A, 2)
    nb = size(B, 2)
    av = A.nzval
    bv = B.nzval

    colptr[1] = 1
    ck = 1  # running total for colptr
    rv_ix = 0  # which entry of rowval and nzval should be filled
    col_ix = 1 # which entry of colptr should be filled

    @inbounds for i in 1:na  # columns of a
        start_a = A.colptr[i]
        len_a = A.colptr[i+1] - start_a
        b_end = 1
        for j in 2:nb+1  # columns of b
            ix_b = b_end
            b_end = B.colptr[j]

            ix_a = 0

            while ix_a < len_a && ix_b < b_end
                if A.rowval[start_a+ix_a] == B.rowval[ix_b]
                    # found one!
                    nzval[rv_ix+=1] = av[start_a + ix_a] * bv[ix_b]
                    rowval[rv_ix] = A.rowval[start_a+ix_a]
                    ix_a += 1
                    ck += 1
                elseif A.rowval[start_a+ix_a] < B.rowval[ix_b]
                    ix_a += 1
                else
                    ix_b += 1
                end
            end

            colptr[col_ix+=1] = ck
        end

    end

    out
end

function row_kron!(out::SparseMatrixCSC, A::AbstractMatrix{T}, B::SparseMatrixCSC) where T
    colptr = out.colptr
    rowval = out.rowval
    nzval = out.nzval
    my_zero = zero(T)

    na = size(A, 2)
    nb = size(B, 2)
    br = B.rowval
    bv = B.nzval

    colptr[1] = 1
    ck = 1  # running total for colptr
    rv_ix = 0  # which entry of rowval and nzval should be filled
    c_ix = 1 # which entry of colptr should be filled

    @inbounds for i in 1:na  # columns of a
        ix_b = 1
        for j in 2:nb+1  # columns of b
            b_end = B.colptr[j]

            while ix_b < b_end  # rows
                row = br[ix_b]
                if A[row, i] != my_zero
                    nzval[rv_ix+=1] = A[row, i] * bv[ix_b]
                    rowval[rv_ix] = row
                    ck += 1
                end
                ix_b += 1
            end
            colptr[c_ix+=1] = ck
        end
    end
    out
end

function row_kron!(out::SparseMatrixCSC, A::SparseMatrixCSC, B::AbstractMatrix{T}) where T
    colptr = out.colptr
    rowval = out.rowval
    nzval = out.nzval
    my_zero = zero(T)

    na = size(A, 2)
    nb = size(B, 2)
    ar = A.rowval
    av = A.nzval

    colptr[1] = 1
    ck = 1  # running total for colptr
    rv_ix = 0  # which entry of rowval and nzval should be filled
    c_ix = 1 # which entry of colptr should be filled

    @inbounds for i in 1:na  # columns of a
        for j in 1:nb  # columns of b
            for ptr in A.colptr[i]:(A.colptr[i+1]-1)  # rows
                row = ar[ptr]
                if B[row, j] != my_zero
                    nzval[rv_ix+=1] = av[ptr] * B[row, j]
                    rowval[rv_ix] = row
                    ck += 1
                end
            end
            colptr[c_ix+=1] = ck
        end
    end
    out
end

# ckronxi.m -- DONE
ckronxi(b::Matrix{T}, c, ind=1:length(b)) where {T<:Number} = b \ c  # 23

function ckronxi(b::Array, c, ind=1:length(b))
    d = length(ind)  # 25
    n = Int[size(b[ind[i]], 2) for i=1:d]  #26-27
    prod(n) != size(c, 1) && error("b and c are not conformable")  # 28-30

    z = c'  # 31
    mm = 1  # 32
    for i=1:d  # 33
        m = round(Int, prod(size(z)) / n[i])  # 34
        z = reshape(z, m, n[i])  # 35
        z = b[ind[i]] \ z'  # 36
        mm *= size(z, 1)  # 37
    end  # 38
    reshape(z, mm, size(c, 2))  # 39
end

struct RowKron{T<:Tuple{Vararg{<:AbstractMatrix}}}
    B::T
end

function RowKron(B::AbstractMatrix...)
    nrow = map(x -> size(x, 1), B)
    if any(x -> x != nrow[1], nrow)
        msg = "All matrices must have the same number of rows"
        throw(DimensionMismatch(msg))
    end

    RowKron(B)
end

Base.length(rk::RowKron) = length(rk.B)
@generated Base.eltype(::RowKron{T}) where {T} =
    reduce(promote_type, map(eltype, T.parameters))

Base.issparse(rk::RowKron) = map(issparse, rk.B)

function Base.size(rk::RowKron, i::Int)
    if i == 1
        size(rk.B[1], 1)
    elseif i == 2
        prod(map(A -> size(A, 2), rk.B))
    else
        1
    end
end

Base.size(rk::RowKron) = (size(rk, 1), size(rk, 2))

sizes(rk::RowKron, i::Integer) = collect(map(A -> size(A, i), rk.B))

for (f, op, transp) in ((:A_mul_B!, :identity, false),
                        (:Ac_mul_B!, :ctranspose, true),
                        (:At_mul_B!, :transpose, true))

    # code to do type checks for input args
    checks = transp ? quote
        size(out, 1) == size(rk, 2) || throw(DimensionMismatch())
        size(out, 2) == size(c, 2) || throw(DimensionMismatch())
        size(c, 1) == size(rk, 1) || throw(DimensionMismatch())
        fill!(out, 0.0)
    end : quote
        size(out, 1) == size(rk, 1) || throw(DimensionMismatch())
        size(out, 2) == size(c, 2) || throw(DimensionMismatch())
        size(c, 1) == size(rk, 2) || throw(DimensionMismatch())
        fill!(out, 0.0)
    end

    # the outer loop has a different range depending on transpose
    outer_loop_range = transp ? :(1:size(out, 1)) : :(1:size(c, 1))

    # How we fill `out` is also different...
    fill_out = transp ? quote
        if _last_sparse
            for ptr in (Bend.colptr[ccB]):(Bend.colptr[ccB+1]-1)
                _r = Bend.rowval[ptr]
                _v = $(op)(Bend.nzval[ptr])
                full_B = d[_r, end] * _v
                for _c in 1:ccol
                    out[i, _c] += full_B * c[_r, _c]
                end
            end
        else
            for _r in 1:nrow
                full_B = d[_r, end] * Bend[_r, ccB]
                for _c in 1:ccol
                    out[i, _c] += full_B * c[_r, _c]
                end
            end
        end
    end : quote
        if _last_sparse
            for ptr in (Bend.colptr[ccB]):(Bend.colptr[ccB+1]-1)
                _r = Bend.rowval[ptr]
                _v = $(op)(Bend.nzval[ptr])
                full_B = d[_r, end] * _v

                for _c in 1:ccol
                    out[_r, _c] += full_B * c[i, _c]
                end
            end
        else
            for _r in 1:nrow
                full_B = d[_r, end] * Bend[_r, ccB]
                for _c in 1:ccol
                    out[_r, _c] += full_B * c[i, _c]
                end
            end
        end
    end

    @eval begin
    function Base.$(f)(out::StridedVecOrMat, rk::RowKron, c::StridedVecOrMat)
        $(checks)

        _is_sparse = map(x -> isa(x, SparseMatrixCSC), rk.B)
        _last_sparse = _is_sparse[end]

        nrow = size(rk, 1)
        Nb = length(rk)
        d = Array{eltype(rk)}(nrow, Nb)
        d[:, 1] = one(eltype(rk))

        n_col_B = sizes(rk, 2)
        cur_col_B = copy(n_col_B) + 1

        # simplify notation
        ccol = size(c, 2)
        Bend = rk.B[end]

        for i in $(outer_loop_range)  # loop over rows of out
            if cur_col_B[end] > n_col_B[end]
                cur_col_B[end] = 1
                j = Nb
                while j > 1
                    j -= 1
                    cur_col_B[j] += 1
                    cur_col_B[j] <= n_col_B[j] && break
                    cur_col_B[j] = 1
                end

                for col_d in j+1:Nb
                    _j = col_d - 1
                    _col_j = cur_col_B[_j]
                    this_B = rk.B[_j]

                    if _is_sparse[_j]
                        d[:, col_d] = 0.0

                        # fill in non-zero rows
                        for ptr in this_B.colptr[_col_j]:(this_B.colptr[_col_j+1]-1)
                            _r = this_B.rowval[ptr]
                            _v = $(op)(this_B.nzval[ptr])
                            d[_r, col_d] = d[_r, col_d-1] * _v
                        end

                    else
                        @simd for row_d in 1:nrow
                            d[row_d, col_d] = d[row_d, col_d-1] * this_B[row_d, _col_j]
                        end
                    end
                end

            end

            # now we use d[:, end] .* B[end][:, cur_col_B[end]] .* c[i]
            ccB = cur_col_B[end]

            $(fill_out)

            cur_col_B[end] += 1
        end
        out
    end
    end  # @eval begin
end

function *(rk::RowKron, c::StridedVector)
    out = zeros(promote_type(eltype(rk), eltype(c)), size(rk, 1))
    A_mul_B!(out, rk, c)
    out
end

function *(rk::RowKron, c::StridedMatrix)
    out = zeros(promote_type(eltype(rk), eltype(c)), size(rk, 1), size(c, 2))
    A_mul_B!(out, rk, c)
    out
end

function Base.At_mul_B(rk::RowKron, c::StridedVector)
    out = zeros(promote_type(eltype(rk), eltype(c)), size(rk, 2))
    At_mul_B!(out, rk, c)
    out
end

function Base.At_mul_B(rk::RowKron, c::StridedMatrix)
    out = zeros(promote_type(eltype(rk), eltype(c)), size(rk, 2), size(c, 2))
    At_mul_B!(out, rk, c)
    out
end

# cdprodx.m -- DONE
cdprodx(b::Matrix{T}, c, ind=1:prod(size(b))) where {T<:Number} = b*c  # 39

function cdprodx(b::AbstractArray{T}, c::StridedVecOrMat,
ind::AbstractArray{Int}=1:prod(size(b))) where T<:AbstractMatrix
    _check_cdprodx(b, c, ind)
    rk = RowKron(b[ind]...)
    rk*c
end

# nodeunif.m -- DONE
function nodeunif(n::Int, a::Int, b::Int)
    x = linspace(a, b, n)
    return x, x
end

function nodeunif(n::AbstractArray{T}, a::AbstractArray, b::AbstractArray) where T<:Integer
    d = length(n)
    xcoord = Array{AbstractVector}(d)
    for k=1:d
        xcoord[k] = linspace(a[k], b[k], n[k])
    end
    return gridmake(xcoord...), xcoord
end

function lookup(table::AbstractVector, x::Real, p::Int=0)
    ind = searchsortedfirst(table, x) - 1
    m = length(table)

    if ind == m && p >= 2
        tosub = 1
        @inbounds for i in m-1:-1:1
            if table[i] == table[end]
                tosub += 1
            else
                break
            end
        end
        return m - tosub
    end

    if ind == 0 && (p == 1 || p == 3)
        ix = 0
        @inbounds for i in 1:m
            if table[1] == table[i]
                ix += 1
            else
                break
            end
        end
        return ix
    end
    ind
end

# lookup.c -- DONE
function lookup(table::AbstractVector, x::AbstractVector, p::Int=0)
    n = length(table)
    m = length(x)
    out = Array{Int}(m)

    # lower enbound adjustment
    numfirst = 1
    t1 = table[1]
    for i in 2:n
        if table[i] == t1
            numfirst += 1
        else
            break
        end
    end

    # upper endpoint adjustment
    tn = table[n]
    if p >= 2
        n -= 1
        for i in n:-1:1
            if table[i] == tn
                n -= 1
            else
                break
            end
        end
    end

    n1 = n - 1
    n2 = n - 1

    if n - numfirst < 1  # only one unique value in table
        if p == 1 || p == 3
            for i=1:m
                out[i] = numfirst
            end
        else
            for i in 1:m
                if table[1] <= x[i]
                    out[i] = numfirst
                else
                    out[i] = 0
                end
            end
        end
        return out
    end

    jlo = 1

    for i in 1:m
        inc = 1
        xi = x[i]
        if xi >= table[jlo]
            jhi = jlo + 1
            while jhi <= n && xi >= table[jhi]
                jlo = jhi
                jhi += inc
                if jhi > n
                    jhi = n+1
                end
            end
        else
            jhi = jlo
            jlo -= 1
            while jlo > 0 && xi < table[jlo]
                jhi = jlo
                jlo -= inc
                if jlo < 1
                    jlo = 0
                    break
                else
                    inc += inc
                end
            end
        end

        while jhi - jlo > 1
            j = (jhi + jlo) >> 1
            if xi >= table[j]
                jlo = j
            else
                jhi = j
            end
        end

        out[i] = jlo

        if jlo < 1
            jlo = 1
            if p == 1 || p == 3
                out[i] = numfirst
            end
        end

        if jlo == n1
            jlo = n2
        end
    end

    out
end

"""
    lookup(table::AbstractVector, x, p::Int=0)

For a sorted vector `table`, return the index of the _last_ point smaller than
`x`. `p` specifies endpoint adjustments and takes one of 4 values:

1. `p=0`: no adjustment made
2. `p=1`: if `x < minimum(table)`, then return `sum(table .== table[1])`
3. `p=2`: if `x > maximum(table)`, then return `length(table)-sum(table .== table[end])`
4. `p=3`: do both adjustments 1 and 2

If `p > 3`, the `p=3` case is applied

When `x` is a vector, a vector of integers is returned
"""
lookup

# utility to expand the order input if needed
# used in basis_structure.jl and interp.jl
_check_order(N::Int, order::Int) = fill(order, 1, N)
_check_order(N::Int, order::Vector) = reshape(order, 1, N)

function _check_order(N::Int, order::Matrix)
    if size(order, 2) == N
        return order
    end

    if size(order, 1) == N
        m = size(order, 2)
        return reshape(order, m, N)
    end

    error("invalid order argument. Should have $N columns")
end

function _check_cdprodx(b::Array, c, ind::AbstractArray{Int})
    _ind_min, _ind_max = extrema(ind)
    @assert _ind_min > 0 && _ind_max <= length(b) "ind not conformable"
end

function _nnz_per_row(A::Matrix{T}) where T
    counts = zeros(Int, size(A, 1))
    my_zero = zero(T)
    @inbounds for col in 1:size(A, 2), row in 1:size(A, 1)
        if A[row, col] != my_zero
            counts[row] += 1
        end
    end
    counts
end

function _nnz_per_row(A::SparseMatrixCSC{T}) where T
    counts = zeros(Int, size(A, 1))
    @inbounds for col in 1:size(A, 2), ptr in A.colptr[col]:(A.colptr[col+1]-1)
        counts[A.rowval[ptr]] += 1
    end
    counts
end

function _row_kron_sparse_out_nnz(A, B)
    # get number of non-zeros in out
    acounts = _nnz_per_row(A)
    bcounts = _nnz_per_row(B)
    k = 0
    @inbounds @simd for _r in 1:size(A, 1)
        k += acounts[_r]*bcounts[_r]
    end
    k
end

function _allocate_row_kron_out(::Type{SparseMatrixCSC},
                                A::AbstractMatrix{T},
                                B::AbstractMatrix{S}) where {T,S}
    nobsa, na = size(A)
    nobsb, nb = size(B)
    k = _row_kron_sparse_out_nnz(A, B)
    SparseMatrixCSC(nobsa, na*nb,
        ones(Int, na*nb+1),          # colptr
        Array{Int}(k),                # rowval
        Array{promote_type(S,T)}(k)  # nzval
    )
end

function _allocate_row_kron_out(::Type{Matrix},
                                A::AbstractMatrix{T},
                                B::AbstractMatrix{S}) where {T,S}
    nobsa, na = size(A)
    nobsb, nb = size(B)
    Array{promote_type(S,T)}(nobsa, na*nb)
end

_allocate_row_kron_out(A::SparseMatrixCSC, B::SparseMatrixCSC) =
    _allocate_row_kron_out(SparseMatrixCSC, A, B)

_allocate_row_kron_out(A::AbstractMatrix, B::SparseMatrixCSC) =
    _allocate_row_kron_out(SparseMatrixCSC, A, B)

_allocate_row_kron_out(A::SparseMatrixCSC, B::AbstractMatrix) =
    _allocate_row_kron_out(SparseMatrixCSC, A, B)

_allocate_row_kron_out(A::AbstractMatrix, B::AbstractMatrix) =
    _allocate_row_kron_out(Matrix, A, B)


text:
"""
    cudnnPoolingForward(x; mode, nanOpt, window, padding, stride, alpha)
    cudnnPoolingForward(x, d::cudnnPoolingDescriptor; alpha)
    cudnnPoolingForward!(y, x; mode, nanOpt, window, padding, stride, alpha, beta)
    cudnnPoolingForward!(y, x, d::cudnnPoolingDescriptor; alpha, beta)

Return pooled `x`, overwriting `y` if provided, according to keyword arguments or the
pooling descriptor `d`. Please see the [cuDNN
docs](https://docs.nvidia.com/deeplearning/cudnn/api/index.html#cudnnPoolingForward) for
details.

The dimensions of `x,y` tensors that are less than 4-D are assumed to be padded on the left
with 1's. The first `n-2` are spatial dimensions, the last two are always assumed to be
channel and batch.

The arguments `window`, `padding`, and `stride` can be specified as `n-2` dimensional
vectors, tuples or a single integer which is assumed to be repeated `n-2` times. If any of
the entries is larger than the corresponding `x` dimension, the `x` dimension is used
instead.

Arguments:
* `mode = CUDNN_POOLING_MAX`: Pooling method, other options are `CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING`, `CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING`, `CUDNN_POOLING_MAX_DETERMINISTIC`
* `nanOpt = CUDNN_NOT_PROPAGATE_NAN`: NAN propagation policy, the other option is `CUDNN_PROPAGATE_NAN`
* `window = 2`: Pooling window size
* `padding = 0`: Padding assumed around `x`
* `stride = window`: How far to shift pooling window at each step
* `alpha=1, beta=0` can be used for scaling, i.e. `y .= alpha*op(x1) .+ beta*y`
"""
cudnnPoolingForward, cudnnPoolingForward!

# Public methods
cudnnPoolingForward(x; o...)     = cudnnPoolingForwardWithDefaults(x; o...)
cudnnPoolingForward!(y, x; o...) = cudnnPoolingForwardWithDefaults(x; y, o...)
cudnnPoolingForward(x, d::cudnnPoolingDescriptor; o...)     = cudnnPoolingForwardWithDefaults(x; poolingDesc=d, o...)
cudnnPoolingForward!(y, x, d::cudnnPoolingDescriptor; o...) = cudnnPoolingForwardWithDefaults(x; y, poolingDesc=d, o...)

# Private method
function cudnnPoolingForwardWithDefaults(
    x;                          # no type for x, could be AutoGrad.Value
    mode::cudnnPoolingMode_t = CUDNN_POOLING_MAX,
    nanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
    window::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = 2,
    padding::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = 0,
    stride::Union{Integer,Vector{<:Integer},Tuple{<:Integer,Vararg{Int}}} = window,
    poolingDesc::cudnnPoolingDescriptor = cudnnPoolingDescriptor(mode, nanOpt, Cint(max(2,ndims(x)-2)), pooldims(window,size(x)), pooldims(padding,size(x)), pooldims(stride,size(x))),
    format::cudnnTensorFormat_t = CUDNN_TENSOR_NCHW,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x; format),
    y = cudnnPoolingForwardOutput(x, xDesc, poolingDesc, format),
    yDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(y; format),
    alpha::Real = 1,
    beta::Real = 0,
)
    T = eltype(x)
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    cudnnPoolingForwardAD(x; poolingDesc, alpha, beta, xDesc, yDesc, y)
end

# Convert the integer, tuple or array to pooling dims compatible with array size
function pooldims(d, s::Dims{N}) where N
    if d isa Integer || length(d) == N-2
        Cint[reverse(min.(d,s[1:N-2]))...]
    else
        throw(DimensionMismatch("Cannot pool $(Base.dims2string(s)) array with $d pooldims."))
    end
end

pooldims(d, s::Dims{3}) = pooldims(d, (1,s...))
pooldims(d, s::Dims{2}) = pooldims(d, (1,1,s...))
pooldims(d, s::Dims{1}) = pooldims(d, (1,1,1,s...))
pooldims(d, s::Dims{0}) = pooldims(d, (1,1,1,1))

function cudnnPoolingForwardOutput(x, xDesc, poolingDesc, format)
    d = Array{Cint}(undef, max(4, ndims(x))) # d = [N,C,Yn,...,Y1] no matter what format
    cudnnGetPoolingNdForwardOutputDim(poolingDesc, xDesc, length(d), d)
    if length(d) > ndims(x) # This happens when x is (X,C,N), its TD is [N,C,X,1]
        @assert all(d[ndims(x)+1:end] .== 1)
        d = d[1:ndims(x)]
    end
    # ydims(NCHW)=(Y1,...,Yn,C,N) ydims(NHWC)=(C,Y1,...,Yn,N)
    ydims = (format === CUDNN_TENSOR_NCHW ? reverse(d) : (d[2],d[end:-1:3]...,d[1]))
    similar(x, ydims...)
end

# AD method
function cudnnPoolingForwardAD(x; poolingDesc, alpha, beta, xDesc, yDesc, y)
    cudnnPoolingForward(handle(), poolingDesc, alpha, xDesc, x, beta, yDesc, y)
    return y
end

# Deprecated methods
function cudnnPoolingForward(y::DenseCuArray{T,N}, x::DenseCuArray{T,N}, pdims::NNlib.PoolDims;
                             alpha=1, beta=0, mode=CUDNN_POOLING_MAX) where {T,N}
    @warn "`cudnnPoolingForward(y,x,d::PoolDims)` is deprecated, please use one of the methods in `@doc cudnnPoolingForward`." maxlog=1
    cudnnPoolingForward!(y, x; window=NNlib.kernel_size(pdims), padding=nnlibPadding(pdims), stride=NNlib.stride(pdims), mode, alpha, beta)
end


text:
using Test
using CSV
using DataFrames
using FlightMechanicsSimulator
using FlightMechanicsUtils

@testset "atmosphere_f16" begin
    df = DataFrame(CSV.File("data/adc.csv"))
    # df has duplicate altitudes with different vt
    df = unique(df, :alt)
    # Discontinuity alt (35000 ft) produces tests not passing because stevens atmosphere
    # is not continous at 35000 ft, so if conversion FT2M produces slightly different altitude,
    # a wrong atmosphere layer will be chosen
    df = df[(!isapprox).(df[:, :alt], 35000), :]
    for case in eachrow(df)
        T1, ρ1, a1, p1 = F16Stevens.atmosphere_f16(case.alt * FT2M)
        T2, ρ2, a2, p2 = F16Stevens.atmosphere(case.alt)
        @test isapprox(T1*KEL2RANK, T2)
        @test isapprox(ρ1*KGM32SLUGFT3, ρ2)
        @test isapprox(a1*M2FT, a2)
        @test isapprox(p1*PA2PSF, p2)
    end
end

# The atmosphere in Stevens does not provide exact results when tested against official tables
@testset "ISA1978-Morelli" begin
    # Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and
    # practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Appendix C, p. 574.
    # Test sea level
    T, ρ, a, p = F16Stevens.atmosphere(0.0)  # 0 ft
    @test_broken isapprox(T, 518.67, atol=0.005)  # Rankine
    @test_broken isapprox(p, 2116.2, atol=0.05)  # psf
    @test_broken isapprox(ρ, 0.0023769, atol=5e-8)  # slug/ft³
    @test_broken isapprox(a, 1116.44, atol=0.005)  # ft/s
    # Test in tropopause
    T, ρ, a, p = F16Stevens.atmosphere(50000.0)  # 0 ft
    @test_broken isapprox(T, 389.97, atol=0.005)  # Rankine
    @test_broken isapprox(p, 243.6, atol=0.05)  # psf
    @test_broken isapprox(ρ, 0.0003639, atol=5e-8)  # slug/ft³
    @test_broken isapprox(a, 968.08, atol=0.005)  # ft/s
end


text:
"""
# AccessControlSystem

## Description
> The state of an access control system is defined by a triple (S, O, A), where S is the set of subjects, O is the set of objects, and A is an access matrix, with rows corresponding to subjects and columns corresponding to objects. An entry A[s, o] is a set of rights.
(Ferraiolo at al., 2007)

## Bibliography
- Ferraiolo, D., Kuhn, D.R., Chandramouli, R., 2007. Role-based access control, 2nd ed. Artech House, Boston.

## TODO:
- Different data structures are possible to manage access control data, especially:
- BitArray{<:Bool, 3}. Pros: analytics. Cons: only small datasets, need to syncrhonize subject,object,right vectors.
- Sparse array: Pros & Cons: same except we can manage larger real-life datasets.
- Database table: Pros: large data. Cons: analytic perfs.
- We should consider both in-memory array storage and remote DB (JuliaDB?) in remote system for voluminous data.
- Possible approach: use a table-like data structure. And make it exportable to Array on-demand.
"""
mutable struct AccessControlSystem
    label::String;
    subject::Vector{<:Subject};
    object::Vector{<:Object};
    right::Vector{<:Right};
    access_matrix::BitArray{3};
end;
function AccessControlSystem(subject_length::Int = 0, object_length::Int = 0, right_length::Int = 0; label::String = "An ACS")
    subject = Subject(n=subject_length);
    right = Right(n=right_length);
    object = Object(n=object_length);
    access_matrix = falses(subject_length, object_length, right_length);
    acs = AccessControlSystem(label, subject, object, right, access_matrix);
    return acs;
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
mutable struct force_torque_t <: LCMType
    timestamp::Int64
    fx::Float64
    fy::Float64
    fz::Float64
    tx::Float64
    ty::Float64
    tz::Float64
end

@lcmtypesetup(force_torque_t)


text:
using Plots

# Model
include_model("bimanipulation_block_v2")

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[model.idx_u[1:4]] .= 10.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[model.idx_u[1:4]] .= -10.0

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
ϵ = 1.0e-8
q1 = [0.0, r, 0.0, -r - ϵ, r, r + ϵ, r]
x1 = [q1; q1]
ϕ_func(model, q1)
x_shift = 0.0
z_shift = 0.25
qT = [x_shift, r + z_shift, 0.0, -r + x_shift, r + z_shift, r + x_shift, r + z_shift]
ϕ_func(model, qT)
xT = [qT; qT]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

ϕ_func(model, zeros(7))
P_func(model, zeros(7))

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])
obj_velocity = velocity_objective(
    [t > T / 2 ? Diagonal(1.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3]))

obj_tracking = quadratic_tracking_objective(
    [Diagonal(0.0 * ones(model.n)) for t = 1:T],
    # [Diagonal(0.1 * ones(model.m)) for t = 1:T-1],
	[Diagonal([0.0 * ones(model.nu);
		zeros(model.nc);
		ones(model.nb);
		zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

function l_stage(x, u, t)
	J = 0.0
	p_box = view(x, 7 .+ (1:2))
	p_pusher1 = view(x, 7 .+ (4:5))
	p_pusher2 = view(x, 7 .+ (6:7))

	ϕ = ϕ_func(model, view(x, 7 .+ (1:7)))
	if true
		J += 1.0 * sum(ϕ[5:6].^2.0)
		# J += 100.0 * sum((p_box - p_pusher1).^2.0)
		# J += 100.0 * sum((p_box - p_pusher2).^2.0)
	end

	return J
end

obj_stage = nonlinear_stage_objective(l_stage, l_stage)

obj_penalty = PenaltyObjective(1.0e5, model.m)

obj_ctrl_velocity = control_velocity_objective(Diagonal([1.0 * ones(model.nu)..., 0.0e-3 * ones(model.nc + model.nb)..., zeros(model.m - model.nu - model.nc - model.nb)...]))

obj = MultiObjective([obj_tracking,
	obj_penalty, obj_velocity, obj_stage, obj_ctrl_velocity])

# Constraints
include_constraints(["contact", "stage"])
t_idx = vcat([t for t = 1:T-1])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])#, con_ctrl_comp, con_ctrl_lim])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.1 * randn(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-3, c_tol = 1.0e-3, max_iter = 2000)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
#open(vis)
visualize!(vis, model,
    q, u,
    Δt = h,
	r = model.block_dim[1])

plot(hcat(q̄...)', color = :red, width = 1.0, labels = "")
plot(hcat([ū..., ū[end]]...)[model.idx_u[1:2], :]', linetype = :steppost, width = 1.0, labels = ["x" "z"])
plot(hcat([ū..., ū[end]]...)[model.idx_u[3:4], :]', linetype = :steppost, width = 1.0, labels = ["x" "z"])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using NgSpice
using Test
using SafeTestsets

@safetestset "interface functions" begin
    include("ngspice_api.jl")
end

@testset "friend functions" begin
    include("wrapper_api.jl")
end

@safetestset "wrapper safetests" begin
    include("wrapper_safetests.jl")
end

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
# Copyright 2015, 2016, 2017, 2018 Martin Holters
# See accompanying license file.

using ACME

function diodeclipper(::Type{Circuit})
    @circuit begin
        j_in = voltagesource(), [-] ⟷ gnd
        r1 = resistor(1e3), [1] ⟷ j_in[+]
        c1 = capacitor(47e-9), [1] ⟷ r1[2], [2] ⟷ gnd
        d1 = diode(is=1e-15), [-] ⟷ gnd, [+] ⟷ r1[2]
        d2 = diode(is=1.8e-15), [-] ⟷ r1[2], [+] ⟷ gnd
        j_out = voltageprobe(), [-] ⟷ gnd, [+] ⟷ r1[2]
    end
end

diodeclipper(::Type{DiscreteModel}=DiscreteModel; fs=44100, solver=nothing) =
    solver === nothing ?
        DiscreteModel(diodeclipper(Circuit), 1//fs) :
        DiscreteModel(diodeclipper(Circuit), 1//fs, solver)


text:
##===----------------------------------------------------------------------===##
##                                   ACCELERATE.JL                            ##
##                                                                            ##
##===----------------------------------------------------------------------===##
##                                                                            ##
##  This file checks for the available frameworks and imports only the        ##
##  frameworks available; it also coordinates access to three frameworks      ##
##  through the unified Vectorize.f() syntax                                  ##
##                                                                            ##
##===----------------------------------------------------------------------===##
module Vectorize

using Libdl
using SpecialFunctions
using Statistics
export @vectorize, @replacebase

OS = Sys.KERNEL

## FUNCTION DICTIONARY
functions = Dict()

## Function for adding to function dictionary
function addfunction(d::Dict, k, v)
    if k in keys(d)
        append!(d[k], [v])
    else
        d[k] = []
        append!(d[k], [v])
    end
end

## START OSX
## On OS X, check if Accelerate can be found - future proofing
OS == :Darwin ? (ACCELERATE = true) : (ACCELERATE = false)
if ACCELERATE == true
    if isfile("/System/Library/Frameworks/Accelerate.framework/Accelerate")
        eval(:(include("Accelerate.jl")))
    else
        warn("Unable to locate the Apple Accelerate framework")
    end
end
## END OSX

## Find Yeppp library
libyeppp_ = Libdl.find_library(["libyeppp"])
if libyeppp_ != ""
    @eval const global libyeppp = libyeppp_
else
    currdir = @__FILE__
    @static if Sys.iswindows()
        bindir = currdir[1:end-16]*"deps\\src\\yeppp\\binaries\\"
    else
        bindir = currdir[1:end-16]*"deps/src/yeppp/binaries/"
    end

    if OS == :Darwin
        @eval const global libyeppp = bindir*"macosx/x86_64/libyeppp.dylib"
    elseif OS == :Linux
        @eval const global libyeppp = bindir*"linux/x86_64/libyeppp.so"
    elseif OS == :NT
        @eval const global libyeppp = bindir*"windows\\amd64\\yeppp.dll"
    end
end

if isfile(libyeppp)
    include("Yeppp.jl") # include Yeppp if present
end

## Find VML
const global libvml = Libdl.find_library(["libmkl_rt"], ["/opt/intel/mkl/lib"])
if libvml != ""
    include("VML.jl")
end

# vectorize macro
macro vectorize(ex)
    len = length(ex.args)
    if len == 2
        f = :(Vectorize.$(ex.args[1]))
        arg = ex.args[2]
        return esc(:($f($arg)))
    elseif len == 3
        f = :(Vectorize.$(ex.args[1]))
        arg1 = ex.args[2]
        arg2 = ex.args[3]
        return esc(:($f($arg1, $arg2)))
    else
        f = :(Vectorize.$(ex.args[1]))
        arg1 = ex.args[2]
        arg2 = ex.args[3]
        arg3 = ex.args[4]
        return esc(:($f($arg1, $arg2, $arg3)))
    end
end

# replacebase macro
"""
Replace all broadcasted base functions with benchmarked vectorized equivalents.
Alternatively, if the user may pass the functions to overload instead.
"""
macro replacebase(fs...)
    b = Expr(:block)
    for ((fvec, args), vectorized_f) in functions
        m, f = get_corresponding_f(fvec)
        if isempty(fs) || f in fs
            if m == :Base
                if ismutating(fvec)
                    Tdest = first(args)
                    Targs = makeargs(args[2:end])
                    e = quote
                        (Base.copyto!)(dest::$Tdest, bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($m.$f), $Targs}) where {Style, Axes, N} = (Vectorize.$fvec)(dest, bc.args...)
                    end
                else
                    Targs = makeargs(args)
                    e = quote
                        (Base.copy)(bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($m.$f), $Targs}) where {Style, Axes, N} = (Vectorize.$fvec)(bc.args...)
                    end
                end
                push!(b.args, e)
            end
        end
    end
    b
end

function get_corresponding_f(fvec)
    fstr = string(fvec)
    if ismutating(fvec)
        fstr = chop(fstr)
    end
    if fstr == "add"
        fstr = "+"
    elseif fstr == "sub"
        fstr = "-"
    elseif fstr == "mul"
        fstr = "*"
    elseif fstr == "div"
        fstr = "/"
    elseif fstr == "rec"
        fstr = "inv"
    elseif fstr == "pow"
        fstr = "^"
    end

    m = :Base
    if fstr in ("cdfnorminv", "erf", "erfc", "erfi", "erfinv", "erfcinv", "cdfnorm")
        m = :SpecialFunctions
    elseif fstr in ("mean", )
        m = :Statistics
    elseif fstr in ("dot", )
        m = :LinearAlgebra
    elseif fstr in ("pow3o2", "frac", "mulbyconj", "tanpi", "fdiv", "invsqrt",
                    "sqr", "pow2o3", "summag", "sumsqr", "invcbrt")
        m = :None
    end

    fsym = Meta.parse(fstr)
    return m, fsym
end
@inline ismutating(f) = last(string(f)) == '!'
function makeargs(Targs)
    arg_str = "Tuple{"
    for T in Targs
        arg_str *= "$T, "
    end
    arg_str = arg_str[1:end-2]
    arg_str *= "}"
    arg_str = replace(arg_str, " where N" => "")

    return Meta.parse(arg_str)
end

# Include optimized functions
include("Functions.jl")

end # module


text:
#!/usr/bin/env julia

using Pkg
Pkg.instantiate()
using SciMLBenchmarks

repo_dir = dirname(dirname(@__DIR__))

model = ARGS[1]
split_path = split(model, "/")
file = split_path[end]
folder = split_path[end-1]

project_abspath = joinpath(repo_dir, model)

@info("debugging!", model, file, folder, project_abspath, isfile(project_abspath), occursin(".jmd", file))

if !isfile(project_abspath)
    error("Invalid file $file")
end

if occursin(".jmd", file)
    @info("Rebuilding $folder/$file")
    SciMLBenchmarks.weave_file(folder, file)
elseif occursin(".toml", file)
    @info("Rebuilding $folder")
    SciMLBenchmarks.weave_folder(folder)
else
    error("Invalid file $file")
end


text:
module PkgDeps

using Pkg
using Pkg.Types: VersionNumber, VersionRange
using REPL
using TOML: parsefile
using UUIDs
using Compat

export PkgEntry, RegistryInstance
export NoUUIDMatch, PackageNotInRegistry
export users, reachable_registries
export direct_dependencies, dependencies

include("pkg_entry.jl")
include("registry_instance.jl")
include("exceptions.jl")
include("utilities.jl")

const GENERAL_REGISTRY = "General"

# borrowed from <https://github.com/JuliaRegistries/RegistryTools.jl/blob/77cae9ef6a075e1d6ec1592bc3e166234d3f01c8/src/builtin_pkgs.jl>
const stdlibs = isdefined(Pkg.Types, :stdlib) ? Pkg.Types.stdlib : Pkg.Types.stdlibs
const STDLIBS = stdlibs()

"""
    reachable_registries(registry_names::Array)
    reachable_registries(registry_name::String; depots::Union{String, Vector{String}}=Base.DEPOT_PATH)
    reachable_registries(; depots::Union{String, Vector{String}}=Base.DEPOT_PATH)

Get an array of found registries.

# Arguments
- `registry_names::Array`: List of registries to retrieve

# Keywords
- `depots::Union{String, Vector{String}}=Base.DEPOT_PATH`: Depots to look in for registries

# Returns
- `Array{RegistryInstance}`: List of all found registries
"""
function reachable_registries(
    registry_names::Array;
    depots::Union{String, Vector{String}}=Base.DEPOT_PATH,
    kwargs...
)
    registries = RegistryInstance[]

    if depots isa String
        depots = [depots]
    end

    for d in depots
        isdir(d) || continue
        reg_dir = joinpath(d, "registries")
        isdir(reg_dir) || continue

        for name in readdir(reg_dir)
            if isempty(registry_names) || name in registry_names
                file = joinpath(reg_dir, name, "Registry.toml")
                isfile(file) || continue
                push!(registries, RegistryInstance(joinpath(reg_dir, name)))
            end
        end
    end

    return registries
end
reachable_registries(registry_name::String; depots::Union{String, Vector{String}}=Base.DEPOT_PATH, kwargs...) = reachable_registries([registry_name]; depots=depots, kwargs...)
reachable_registries(; depots::Union{String, Vector{String}}=Base.DEPOT_PATH, kwargs...) = reachable_registries([]; depots=depots, kwargs...)

"""
    users(uuid::UUID; registries::Array{RegistryInstance}=reachable_registries())
    users(pkg_name::String, pkg_registry_name::String=GENERAL_REGISTRY; kwargs...)
    users(pkg_name::String, pkg_registry::RegistryInstance; kwargs...))

Find the users of a given package.

# Arguments
- `uuid::UUID`: UUID of the package.
- `pkg_name::String`: Find users of this package.
- `pkg_registry_name::String="$GENERAL_REGISTRY"`: Name of registry where `pkg_name` is
  registered. This is used to look up the UUID of `pkg_name`.

# Keywords
- `registries::Array{RegistryInstance}=reachable_registries()`: Registries to search for users.

# Returns
- `Array{String}`: All packages which are dependent on the given package.
"""
function users(
    uuid::UUID;
    registries::Array{RegistryInstance}=reachable_registries(),
)
    downstream_dependencies = String[]

    for rego in registries
        for (pkg, pkg_entry) in rego.pkgs
            deps = direct_dependencies(pkg_entry)
            if any(isequal(uuid), values(deps))
                push!(downstream_dependencies, pkg)
            end
        end
    end

    return downstream_dependencies
end

function users(pkg_name::String, pkg_registry_name::String=GENERAL_REGISTRY; kwargs...)
    uuid = _get_pkg_uuid(pkg_name, pkg_registry_name)
    return users(uuid; kwargs...)
end

function users(pkg_name::String, pkg_registry::RegistryInstance; kwargs...)
    uuid = _get_pkg_uuid(pkg_name, pkg_registry)
    return users(uuid; kwargs...)
end

"""
    direct_dependencies(pkg_name::String;
        registries::Array{RegistryInstance}=reachable_registries()) -> Dict{String, UUID}
    direct_dependencies(pkg_uuid::UUID; registries::Array{RegistryInstance}=reachable_registries()) -> Dict{String, UUID}

Returns the direct dependencies of the latest version of a given package
in the form of `Dict` with names as keys and UUIDs as values.

Standard libraries are assumed to not have any dependencies.
"""
direct_dependencies

function direct_dependencies(pkg_entry::PkgEntry)
    base_path = joinpath(pkg_entry.registry_path, pkg_entry.path)
    deps_file_path = joinpath(base_path, "Deps.toml")
    all_direct_dependencies = Dict{String, UUID}()

    if isfile(deps_file_path)
        deps_content = parsefile(deps_file_path)
        dependency_versions = collect(keys(deps_content))
        latest_version = _get_latest_version(base_path)

        # Use the latest_version of pkg, and check to see if pkg_name is in its dependents
        for version_range in dependency_versions
            if in(latest_version, VersionRange(version_range))
                for (k, v) in pairs(deps_content[version_range])
                    all_direct_dependencies[k] = UUID(v)
                end
            end
        end
    end

    return all_direct_dependencies
end

function direct_dependencies(pkg_name::String; registries::Array{RegistryInstance}=reachable_registries())
    if pkg_name in values(STDLIBS)
        return Dict{String, UUID}()
    end
    return direct_dependencies(_find_latest_pkg_entry(pkg_name, missing; registries))
end

function direct_dependencies(pkg_uuid::UUID; registries::Array{RegistryInstance}=reachable_registries())
    if pkg_uuid in keys(STDLIBS)
        return Dict{String, UUID}()
    end
    return direct_dependencies(_find_latest_pkg_entry(missing, pkg_uuid; registries))
end

"""
    dependencies(
        pkg_name::Union{AbstractString, Missing}, pkg_uuid::Union{Missing, UUID}=missing;
        registries::Array{RegistryInstance}=reachable_registries()
    )

Find all packages that the latest version of `pkg_name` depends on (directly or indirectly).

Note: each package in the dependency tree is assumed to be at the latest version; compat bounds
are ignored. Additionally, standard libraries are assumed to not have any dependencies.

# Arguments
- `pkg_name::AbstractString`: Name of the package
- `pkg_uuid::UUID`: UUID of the package, if available
# Keywords
- `registries::Array{RegistryInstance}=reachable_registries()`: Registries to look into
# Returns
- `Dict{String, UUID}()`: Packages which `pkg_name` depends on.
"""
function dependencies(
    pkg_name::Union{AbstractString, Missing}, pkg_uuid::Union{Missing, UUID}=missing;
    registries::Array{RegistryInstance}=reachable_registries()
)
    return _dependencies!(Dict{String, UUID}(), (pkg_name, pkg_uuid); registries)
end

function dependencies(
    uuid::UUID;
    registries::Array{RegistryInstance}=reachable_registries()
)
    return dependencies(missing, uuid; registries)
end

# recursive helper that accumulates into `deps_found`
function _dependencies!(
    deps_found::Dict{String, UUID},
    (name, uuid);
    registries::Array{RegistryInstance}=reachable_registries()
)
    if (uuid => name) in pairs(STDLIBS)
        return deps_found
    end

    direct_deps = direct_dependencies(_find_latest_pkg_entry(name, uuid; registries))

    for (dep_name, dep_uuid) in pairs(direct_deps)
        if dep_name ∈ keys(deps_found)
            if dep_uuid != deps_found[dep_name]
                error("Package (possibly transitively) depends on $(dep_name) twice, with different UUIDs: $(dep_uuid) and $(deps_found[dep_name])!")
            end
        else
            deps_found[dep_name] = dep_uuid
            _dependencies!(deps_found, (dep_name, dep_uuid); registries=registries)
        end
    end

    return deps_found
end

end  # module


text:
#

julia> using Checkers
julia> f(x) = x^2
julia> @test_formany -10<x<10, f(x) >= 0
#Test Passed
  Expression: (:((-10 < x < 10, f(x) >= 0)), :(mode = test_formany))

julia> @test_forall x in -10:10, f(x) >= 0
#Test Passed
  Expression: (x in -10:10, f(x) >= 0)
    @test_forall x in 0:10, f(x) >= 0

julia> @test_forall x in -10:10, f(x) > 0 #Should fail b/c f(0) = 0
#Test Failed
  Expression: (x in -10:10, f(x) > 0)

#tests the conditional property of Julia
 julia> @test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#100 test will be run
#Passes if 100 tests meet the crierea of holding true. Will fail becuase generating 100 pairs of x,y, some will bound to fail condition

#another test the conditional property of Julia
false.
#100,000 tests will performed until 100 passing test
julia> @test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#Same as top but will pass because 100000test will be performed.


text:
# Declare types and structs

"""
    Particle{T}

An abstract type with subtypes
- StaticParticle{T}
- DynamicParticle{T}

"""
abstract type Particle{T} end

"""
    StaticParticle{T} <: Particle{T}

- ```x``` Location of particle
- ```ℓw``` Log weight of particle (possibly unnormalised)
"""
struct StaticParticle{T} <: Particle{T}
    x::T # location
    ℓw::Real # log weight
end

"""
    DynamicParticle{T} <: Particle{T}
- ```x``` Location of particle
- ```ℓw``` Log weight of particle (possibly unnormalised)
- ```ℓdens``` Log density function
- ```τ``` Temperature (or index) for ℓdens
- ```cacheℓdens``` Cache of log density at current temperature (or index)
"""
struct DynamicParticle{T} <: Particle{T}
    x::T # location
    ℓw::Real # log weight
    ℓdens::Function # log density (functions by reference)
    τ::Real # temperature (or index) for ℓdens
    cacheℓdens::Real # cache of log density at current temperature (or index)
end

# Constructors - 'DynamicParticle'

"""
    Particle(x::T, w::Real, ℓdens::Function, τ::Real, logweight::Bool = false) where {T<:Any}

Explicit constructor for Particle with type ```DynamicParticle```.
Auto calculates cached log density ```cacheℓdens```

- ```x``` Location of particle
- ```w``` Weight of particle (possibly unnormalised)
- ```ℓdens``` Log density function
- ```τ``` Temperature (or index) for ℓdens
-- ```logweight``` Is ```w``` on log scale? Default is ```false```

"""
function Particle(x::T, w::Real, ℓdens::Function, τ::Real, logweight::Bool = false) where {T<:Any}
    if hasmethod(ℓdens, (T, Real))
        cacheℓdens = ℓdens(x, τ)
    elseif hasmethod(ℓdens, Tuple{T})
        cacheℓdens = ℓdens(x)
    else
        cacheℓdens = nothing
    end
    if logweight
        DynamicParticle(x, w, ℓdens, τ, real(cacheℓdens))
    else
        DynamicParticle(x, log(w), ℓdens, τ, real(cacheℓdens))
    end

end

"""
    Particle(initial::D, ℓdens::Function, w::Real = 0.0, τ::Real = 0.0, logweight::Bool = false) where {D<:Distribution}

Constructor Particle with type ```DynamicParticle``` from initial distribution ```init_distribution```.
Draws sample from ```init_distribution``` and
Auto calculates cached log density ```cacheℓdens```

- ```initial``` Initial distribution
- ```ℓdens``` Log density function
- ```w``` Weight of particle (possibly unnormalised)
- ```τ``` Temperature (or index) for ℓdens
- ```logweight``` Is ```w``` on log scale? Default is ```false```

"""
#Particle generated from intial distribution
function Particle(initial::D, ℓdens::Function, w::Real = 1.0, τ::Real = 0.0, logweight::Bool = false) where {D<:Distribution}
    x = rand(initial, 1)
    Particle(x, w, ℓdens, τ, logweight)
end

# Constructors - 'SimpleParticle'
"""
    Particle(x::T, w::Real, logweight::Bool = false) where {T<:Any}

Explicit constructor for Particle with type ```StaticParticle```.

- ```x``` Location of particle
- ```w``` Weight of particle (possibly unnormalised)
- ```logweight``` Is ```w``` on log scale? Default is ```false```
"""
function Particle(x::T, w::Real, logweight::Bool = false) where {T<:Any}
    if logweight
        StaticParticle(x, w)
    else
        StaticParticle(x, log(w))
    end
end

"""
    Particle(initial::D, w::Real = 1.0, logweight::Bool = false) where {D<:Distribution}

Constructor Particle with type ```StaticParticle``` from initial distribution ```init_distribution```.
Draws sample from ```init_distribution```

- ```initial``` Initial distribution
- ```w``` Weight of particle (possibly unnormalised)
- ```logweight``` Is ```w``` on log scale? Default is ```false```

"""
#Particle generated from intial distribution
function Particle(initial::D, w::Real = 1.0, logweight::Bool = false) where {D<:Distribution}
    x = rand(initial, 1)
    Particle(x, w, logweight)
end