text:
struct Dynamics{T}
    evaluate::Any
    jacobian::Any
    hessian::Any
    num_next_state::Int
    num_state::Int
    num_action::Int
    num_parameter::Int
    num_jacobian::Int
    num_hessian::Int
    jacobian_sparsity::Vector{Vector{Int}}
    hessian_sparsity::Vector{Vector{Int}}
    evaluate_cache::Vector{T}
    jacobian_cache::Vector{T}
    hessian_cache::Vector{T}
end

function Dynamics(f::Function, num_next_state::Int, num_state::Int, num_action::Int;
    num_parameter::Int=0,
    evaluate_hessian=false)

    #TODO: option to load/save methods
    @variables y[1:num_next_state], x[1:num_state], u[1:num_action], w[1:num_parameter]
    evaluate = f(y, x, u, w)
    jac = Symbolics.sparsejacobian(evaluate, [x; u; y]);
    evaluate_func = eval(Symbolics.build_function(evaluate, y, x, u, w)[2]);
    jacobian_func = eval(Symbolics.build_function(jac.nzval, y, x, u, w)[2]);
    num_jacobian = length(jac.nzval)
    jacobian_sparsity = [findnz(jac)[1:2]...]
    if evaluate_hessian
        @variables λ[1:num_next_state]
        lag_con = dot(λ, evaluate)
        hessian = Symbolics.sparsehessian(lag_con, [x; u; y])
        hessian_func = eval(Symbolics.build_function(hessian.nzval, y, x, u, w, λ)[2])
        hessian_sparsity = [findnz(hessian)[1:2]...]
        num_hessian = length(hessian.nzval)
    else
        hessian_func = Expr(:null)
        hessian_sparsity = [Int[]]
        num_hessian = 0
    end

    return Dynamics(evaluate_func,
        jacobian_func,
        hessian_func,
        num_next_state,
        num_state,
        num_action,
        num_parameter,
        num_jacobian,
        num_hessian,
        jacobian_sparsity,
        hessian_sparsity,
        zeros(num_next_state),
        zeros(num_jacobian),
        zeros(num_hessian))
end

function Dynamics(constraint::Function, constraint_jacobian::Function, num_next_state::Int, num_state::Int, num_action::Int;
    num_parameter::Int=0)

    # jacobian function
    num_variables = num_state + num_action + num_next_state
    jacobian_func = (J, y, x, u, w) -> constraint_jacobian(reshape(view(J, :), num_next_state, num_variables), y, x, u, w)

    # number of Jacobian elements
    num_jacobian = num_next_state * num_variables

    # Jacobian sparsity
    row = Int[]
    col = Int[]
    for j = 1:num_variables
        for i = 1:num_next_state
            push!(row, i)
            push!(col, j)
        end
    end

    jacobian_sparsity = [row, col]

    # Hessian
    hessian_func = Expr(:null)
    hessian_sparsity = [Int[]]
    num_hessian = 0

    return Dynamics(
        constraint,
        jacobian_func,
        hessian_func,
        num_next_state,
        num_state,
        num_action,
        num_parameter,
        num_jacobian,
        num_hessian,
        jacobian_sparsity,
        hessian_sparsity,
        zeros(num_next_state),
        zeros(num_jacobian),
        zeros(num_hessian))
end

function constraints!(violations, indices, constraints::Vector{Dynamics{T}}, states, actions, parameters) where T
    for (t, con) in enumerate(constraints)
        con.evaluate(con.evaluate_cache, states[t+1], states[t], actions[t], parameters[t])
        @views violations[indices[t]] .= con.evaluate_cache
        fill!(con.evaluate_cache, 0.0) # TODO: confirm this is necessary
    end
end

function jacobian!(jacobians, indices, constraints::Vector{Dynamics{T}}, states, actions, parameters) where T
    for (t, con) in enumerate(constraints)
        con.jacobian(con.jacobian_cache, states[t+1], states[t], actions[t], parameters[t])
        @views jacobians[indices[t]] .= con.jacobian_cache
        fill!(con.jacobian_cache, 0.0) # TODO: confirm this is necessary
    end
end

function hessian_lagrangian!(hessians, indices, constraints::Vector{Dynamics{T}}, states, actions, parameters, duals) where T
    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_cache)
            con.hessian(con.hessian_cache, states[t+1], states[t], actions[t], parameters[t], duals[t])
            @views hessians[indices[t]] .+= con.hessian_cache
            fill!(con.hessian_cache, 0.0) # TODO: confirm this is necessary
        end
    end
end

function sparsity_jacobian(constraints::Vector{Dynamics{T}}, num_state::Vector{Int}, num_actions::Vector{Int};
    row_shift=0) where T

    row = Int[]
    col = Int[]
    for (t, con) in enumerate(constraints)
        col_shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_actions[1:t-1])) : 0)
        push!(row, (con.jacobian_sparsity[1] .+ row_shift)...)
        push!(col, (con.jacobian_sparsity[2] .+ col_shift)...)
        row_shift += con.num_next_state
    end

    return collect(zip(row, col))
end

function sparsity_hessian(constraints::Vector{Dynamics{T}}, num_state::Vector{Int}, num_actions::Vector{Int}) where T
    row = Int[]
    col = Int[]
    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_sparsity[1])
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_actions[1:t-1])) : 0)
            push!(row, (con.hessian_sparsity[1] .+ shift)...)
            push!(col, (con.hessian_sparsity[2] .+ shift)...)
        end
    end
    return collect(zip(row, col))
end

num_state_action_next_state(constraints::Vector{Dynamics{T}}) where T = sum([con.num_state + con.num_action for con in constraints]) + constraints[end].num_next_state
num_constraint(constraints::Vector{Dynamics{T}}) where T = sum([con.num_next_state for con in constraints])
num_jacobian(constraints::Vector{Dynamics{T}}) where T = sum([con.num_jacobian for con in constraints])
# num_hessian(constraints::Vector{Dynamics{T}}) where T = sum([con.num_hessian for con in constraints])

function constraint_indices(constraints::Vector{Dynamics{T}};
    shift=0) where T
    [collect(shift + (t > 1 ? sum([constraints[s].num_next_state for s = 1:(t-1)]) : 0) .+ (1:constraints[t].num_next_state)) for t = 1:length(constraints)]
end

function jacobian_indices(constraints::Vector{Dynamics{T}};
    shift=0) where T
    [collect(shift + (t > 1 ? sum([constraints[s].num_jacobian for s = 1:(t-1)]) : 0) .+ (1:constraints[t].num_jacobian)) for t = 1:length(constraints)]
end

function hessian_indices(constraints::Vector{Dynamics{T}}, key::Vector{Tuple{Int,Int}}, num_state::Vector{Int}, num_action::Vector{Int}) where T
    indices = Vector{Int}[]
    for (t, con) in enumerate(constraints)
        if !isempty(con.hessian_sparsity[1])
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
            row = collect(con.hessian_sparsity[1] .+ shift)
            col = collect(con.hessian_sparsity[2] .+ shift)
            rc = collect(zip(row, col))
            push!(indices, [findfirst(x -> x == i, key) for i in rc])
        else
            push!(indices, Int[])
        end
    end
    return indices
end

function state_indices(constraints::Vector{Dynamics{T}}) where T
    [[collect((t > 1 ? sum([constraints[s].num_state + constraints[s].num_action for s = 1:(t-1)]) : 0) .+ (1:constraints[t].num_state)) for t = 1:length(constraints)]...,
        collect(sum([constraints[s].num_state + constraints[s].num_action for s = 1:length(constraints)]) .+ (1:constraints[end].num_next_state))]
end

function action_indices(constraints::Vector{Dynamics{T}}) where T
    [collect((t > 1 ? sum([constraints[s].num_state + constraints[s].num_action for s = 1:(t-1)]) : 0) + constraints[t].num_state .+ (1:constraints[t].num_action)) for t = 1:length(constraints)]
end

function state_action_indices(constraints::Vector{Dynamics{T}}) where T
    [[collect((t > 1 ? sum([constraints[s].num_state + constraints[s].num_action for s = 1:(t-1)]) : 0) .+ (1:(+ constraints[t].num_state + constraints[t].num_action))) for t = 1:length(constraints)]...,
        collect(sum([constraints[s].num_state + constraints[s].num_action for s = 1:length(constraints)]) .+ (1:constraints[end].num_next_state))]
end

function state_action_next_state_indices(constraints::Vector{Dynamics{T}}) where T
    [collect((t > 1 ? sum([constraints[s].num_state + constraints[s].num_action for s = 1:(t-1)]) : 0) .+ (1:(+ constraints[t].num_state + constraints[t].num_action + constraints[t].num_next_state))) for t = 1:length(constraints)]
end

function dimensions(dynamics::Vector{Dynamics{T}};
    parameters=[0 for t = 1:(length(dynamics) + 1)]) where T
    states = [[d.num_state for d in dynamics]..., dynamics[end].num_next_state]
    actions = [[d.num_action for d in dynamics]..., 0]
    return states, actions, parameters
end


text:
@testset "Handles" begin
    @test vg.is_handle("VkInstance")
    @test !vg.is_handle("VkInstanceCreateInfo")
    @test vg.is_handle_with_create_info("VkInstance")
    @test !vg.is_handle_with_create_info("VkPhysicalDevice")
    @test vg.is_handle_destructible("VkDevice")
    @test !vg.is_handle_destructible("VkPhysicalDevice")
    @test vg.is_handle_with_multiple_create_info("VkPipeline")
    @test !vg.is_handle_with_multiple_create_info("VkPhysicalDevice")
    @test issubset(filter(!vg.isalias, unique(vg.vulkan_creation_info.name)), vg.vulkan_handles.name)
end

@testset "Type queries" begin
    @test vg.is_bitmask("VkQueryPoolCreateFlags")
    @test vg.is_enum("VkDebugUtilsMessageSeverityFlagBitsEXT")
    @test !vg.is_bitmask("VkInstance")
    @test vg.is_enum("VkResult")
    @test !vg.is_enum("VkInstance")
    @test vg.isalias("VkAccelerationStructureNV")
    @test !vg.isalias("VkAccelerationStructureKHR")
end

@testset "Array & count variables" begin
    @test vg.is_array_variable("ppEnabledLayerNames", "VkInstanceCreateInfo")
    @test !vg.is_array_variable("enabledLayerCount", "VkInstanceCreateInfo")
    @test !vg.is_count_variable("ppEnabledLayerNames", "VkInstanceCreateInfo")
    @test vg.is_count_variable("enabledLayerCount", "VkInstanceCreateInfo")
    @test vg.associated_array_variables("enabledLayerCount", "VkInstanceCreateInfo") == ["ppEnabledLayerNames"]
    @test vg.associated_count_variable("ppEnabledLayerNames", "VkInstanceCreateInfo") == "enabledLayerCount"
end

@testset "Defaults" begin
    @test vg.default("ppEnabledLayerNames", "Ptr{Cstring}") == "C_NULL"
    @test vg.default("pApplicationInfo", "Ptr{VkApplicationInfo}") == "C_NULL"
    @test vg.default("device", "VkDevice") == "C_NULL"
    @test vg.default("enabledLayerCount", "UInt32") == "0"
end

@testset "Count variables to be filled (enumerations)" begin
    @test vg.is_count_to_be_filled("pPhysicalDeviceCount", "vkEnumeratePhysicalDevices")
    @test vg.enumeration_command_count_variable("vkEnumeratePhysicalDevices").name == "pPhysicalDeviceCount"
    @test vg.enumeration_command_array_variable("vkEnumeratePhysicalDevices").name == "pPhysicalDevices"
end

ext: jl


text:
include("Header.jl")
include("Frame.jl")

# 現在単位はÅ
struct Asd
    header   ::AbstractHeader
    frames   ::Array{Frame, 1}
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
# Runs the pipeline stage provided as argument

using Pkg

# setup environment
const projectdir = get(ENV, "PROJECTDIR", "/project")
const workdir = joinpath(projectdir, "data")

Pkg.activate(projectdir)

# include methods that implement the pipeline stages
include("utils.jl")

# look for the stage to run
const command = ARGS[1]

if command == "fetch_sources"
    fetch_sources()
    sleep(10)           # delay for demo purposes
elseif command == "extract_sources"
    extract_sources()
    sleep(10)           # delay for demo purposes
elseif command == "index_sources"
    index_sources()
    sleep(10)           # delay for demo purposes
elseif command == "run_search_server"
    run_search_server()
else
    error("unknown command $command")
end

max_stars_repo_name: tanmaykm/KubernetesIntegratedJulia

max_issues_repo_name: tanmaykm/KubernetesIntegratedJulia

max_forks_repo_name: tanmaykm/KubernetesIntegratedJulia


text:
"""
    z_up_transform(origin::Point3D, zaxis::FreeVector3D, from::CartesianFrame3D)

Return `Transform3D` from `from` to the frame in which both `origin` and `zaxis`
is expressed, such that z-axis of `from` is `zaxis`, and origin of `from` is
`origin`.
"""
function z_up_transform(origin::Point3D, zaxis::FreeVector3D, from::CartesianFrame3D)
    @framecheck origin.frame zaxis.frame
    to = origin.frame
    rotation = Rotations.rotation_between(SVector(0., 0., 1.), zaxis.v)
    translation = origin.v
    Transform3D(from, to, rotation, translation)
end

function forcebasis(μ::Float64, num_basis_vectors::Val{N}) where N
    Δθ = 2 * π / N
    basisvectors = ntuple(num_basis_vectors) do i
        θ = (i - 1) * Δθ
        normalize(SVector(μ * cos(θ), μ * sin(θ), 1.0))
    end
    hcat(basisvectors...)
end

nvars(model::Parametron.Model, ::Val{N}) where {N} = SVector(ntuple(_ -> Variable(model), Val(N)))

struct ContactPoint{N}
    normal_aligned_frame::CartesianFrame3D
    ρ::SVector{N, Variable} # basis vector multipliers
    force_local::FreeVector3D{SVector{3, Variable}} # contact force expressed in contact point's normal-aligned frame
    wrench_world::Wrench{Variable}
    position::Any # The position, normal, and μ fields are not concretely typed,
    normal::Any   # but these fields are not used in any of the contact handling algorithms.
    μ::Any        # Instead, they are just stored for debugging purposes.
    weight::typeof(Ref(1.0)) # In 0.7 and up, this is just Ref{Float64}
    maxnormalforce::typeof(Ref(1.0)) # In 0.7 and up, this is just Ref{Float64}

    function ContactPoint{N}(
            position::Union{Point3D, Parameter{<:Point3D}}, normal::Union{FreeVector3D, Parameter{<:FreeVector3D}}, μ::Union{Float64, Parameter{Float64}},
            state::MechanismState, model::Parametron.Model) where N
        # frames
        normal_aligned_frame = CartesianFrame3D()
        worldframe = root_frame(state.mechanism)

        # variables
        ρ = nvars(model, Val(N))
        force_local = FreeVector3D(normal_aligned_frame, nvars(model, Val(3)))
        wrench_world = Wrench(worldframe, nvars(model, Val(3)), nvars(model, Val(3)))

        ret = new{N}(normal_aligned_frame, ρ, force_local, wrench_world, position, normal, μ, Ref(0.0), Ref(0.0))

        # constraints
        basis = @expression(forcebasis(μ, Val(N)))
        maxnormalforce = let ret = ret
            Parameter(() -> ret.maxnormalforce[], model)
        end
        maxρ = @expression((maxnormalforce / (N * sqrt(μ*μ + 1))) * ones(N))  # μ*μ instead of μ^2 due to https://github.com/tkoolen/Parametron.jl/issues/66
        state_param = let state = state
            Parameter(() -> state, model)
        end
        toroot = @expression(transform_to_root(state_param, position.frame) * z_up_transform(position, normal, normal_aligned_frame))
        hat = RBD.Spatial.hat
        @constraint(model, force_local.v == basis * ρ)
        @constraint(model, ρ >= zeros(N))
        @constraint(model, ρ <= maxρ)
        @constraint(model, linear(wrench_world) == rotation(toroot) * force_local.v)
        @constraint(model, angular(wrench_world) == hat(translation(toroot)) * linear(wrench_world)) # TODO: ×
        ret
    end
end

disable!(point::ContactPoint) = point.maxnormalforce[] = 0
isenabled(point::ContactPoint) = point.maxnormalforce[] > 0

function objectiveterm(point::ContactPoint, model::Parametron.Model)
    weight = Parameter{Float64}(() -> point.weight[], model)
    f = point.force_local
    @expression weight * (f ⋅ f)
end

# TODO: type piracy:
function Parametron.value(m::Parametron.Model, wrench::Wrench{Variable})
    # Wrench(wrench.frame, value.(Ref(m), angular(wrench)), value.(Ref(m), linear(wrench))) # allocates on 0.6
    τ = angular(wrench)
    f = linear(wrench)
    @inbounds τval = SVector(value(m, τ[1]), value(m, τ[2]), value(m, τ[3]))
    @inbounds fval = SVector(value(m, f[1]), value(m, f[2]), value(m, f[3]))
    return Wrench(wrench.frame, τval, fval)
end


text:
const N_A = 6.02214076e23 # mol-1
const k_B = 1.380649e-23 # m2 kg s-2 K-1
const R̄   = N_A*k_B # m2 kg s-2 K-1 mol-1
const ħ = 1.054571817e-34 # J s
const h = 1.054571817e-34*2*π # J s


text:
"""
Estimate a linear model with high dimensional categorical variables / instrumental variables

### Arguments
* `df::AbstractDataFrame`
* `model::Model`: A model created using [`@model`](@ref)
* `save::Union{Bool, Symbol} = false`: Should residuals and eventual estimated fixed effects saved in a dataframe? Use `save = :residuals` to only save residuals. Use `save = :fe` to only save fixed effects.
* `method::Symbol = :lsmr`: Method to deman regressors. `:lsmr` is akin to conjugate gradient descent.  With parallel use `:lsmr_parallel`. To use multi threaded use `lsmr_threads`. Other choices are `:qr` and `:cholesky` (factorization methods)
* `contrasts::Dict = Dict()` An optional Dict of contrast codings for each categorical variable in the `formula`.  Any unspecified variables will have `DummyCoding`.
* `maxiter::Integer = 10000`: Maximum number of iterations
* `tol::Real =1e-8`: Tolerance

### Details
Models with instruments variables are estimated using 2SLS. `reg` tests for weak instruments by computing the Kleibergen-Paap rk Wald F statistic, a generalization of the Cragg-Donald Wald F statistic for non i.i.d. errors. The statistic is similar to the one returned by the Stata command `ivreg2`.

### Examples
```julia
using DataFrames, RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
df[:StateC] =  categorical(df[:State])
df[:YearC] =  categorical(df[:Year])
reg(df, @model(Sales ~ Price, fe = StateC + YearC))
reg(df, @model(Sales ~ NDI, fe = StateC + StateC&Year))
reg(df, @model(Sales ~ NDI, fe = StateC*Year))
reg(df, @model(Sales ~ (Price ~ Pimin)))
reg(df, @model(Sales ~ Price, weights = Pop))
reg(df, @model(Sales ~ NDI, subset = State .< 30))
reg(df, @model(Sales ~ NDI, vcov = robust))
reg(df, @model(Sales ~ NDI, vcov = cluster(StateC)))
reg(df, @model(Sales ~ NDI, vcov = cluster(StateC + YearC)))
reg(df, @model(Sales ~ YearC), contrasts = Dict(:YearC => DummyCoding(base = 80)))
```
"""
function reg(df::AbstractDataFrame, m::Model; kwargs...)
    reg(df, m.f; m.dict..., kwargs...)
end

function reg(df::AbstractDataFrame, f::Formula;
    fe::Union{Symbol, Expr, Nothing} = nothing,
    vcov::Union{Symbol, Expr, Nothing} = :(simple()),
    weights::Union{Symbol, Expr, Nothing} = nothing,
    subset::Union{Symbol, Expr, Nothing} = nothing,
    maxiter::Integer = 10000, contrasts::Dict = Dict(),
    tol::Real= 1e-8, df_add::Integer = 0,
    save::Union{Bool, Symbol} = false,  method::Symbol = :lsmr, drop_singletons = true
   )
    feformula = fe
    if isa(vcov, Symbol)
        vcovformula = VcovFormula(Val{vcov})
    else
        vcovformula = VcovFormula(Val{vcov.args[1]}, (vcov.args[i] for i in 2:length(vcov.args))...)
    end

    ##############################################################################
    ##
    ## Parse formula
    ##
    ##############################################################################
    rf = deepcopy(f)
    (has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose_iv!(rf)
    rt = Terms(rf)
    has_absorb = feformula != nothing
    if has_absorb
        # check depth 1 symbols in original formula are all CategoricalVector
        if isa(feformula, Symbol)
            x = feformula
            !isa(df[x], CategoricalVector) && error("$x should be CategoricalVector")
        elseif feformula.args[1] == :+
            x = feformula.args
            for i in 2:length(x)
                isa(x[i], Symbol) && !isa(df[x[i]], CategoricalVector) && error("$(x[i]) should be CategoricalVector")
            end
        end
    end
    has_weights = (weights != nothing)

    ##############################################################################
    ##
    ## Save keyword argument
    ##
    ##############################################################################
    if !isa(save, Bool)
        if save ∉ (:residuals, :fe)
            error("the save keyword argument must be a Bool or a Symbol equal to :residuals or :fe")
        end
    end
    save_residuals = (save == :residuals) | (save == true)
    save_fe = (save == :fe) | ((save == true) & has_absorb)

    ##############################################################################
    ##
    ## Construct new dataframe after removing missing values
    ##
    ##############################################################################

    # create a dataframe without missing values & negative weights
    vars = allvars(rf)
    iv_vars = allvars(iv_formula)
    endo_vars = allvars(endo_formula)
    absorb_vars = allvars(feformula)
    vcov_vars = allvars(vcovformula)

    # create a dataframe without missing values & negative weights
    all_vars = vcat(vars, vcov_vars, absorb_vars, endo_vars, iv_vars)
    all_vars = unique(Symbol.(all_vars))

    esample = completecases(df[all_vars])

    if has_weights
        esample .&= isnaorneg(df[weights])
    end
    if subset != nothing
        subset = eval(evaluate_subset(df, subset))
        if length(subset) != size(df, 1)
            error("df has $(size(df, 1)) rows but the subset vector has $(length(subset)) elements")
        end
        esample .&= convert(BitArray, subset)
    end

    if has_absorb
        fes, ids = parse_fixedeffect(df, Terms(@eval(@formula(nothing ~ $(feformula)))))
        if drop_singletons
            for fe in fes
                @show sum(esample)
                remove_singletons!(esample, fe)
                @show sum(esample)
            end
        end
    end

    nobs = sum(esample)
    (nobs > 0) || error("sample is empty")

    # Compute weights
    sqrtw = get_weights(df, esample, weights)

    # Compute pfe, a FixedEffectMatrix
    has_intercept = rt.intercept
    if has_absorb
        # in case some FixedEffect does not have interaction, remove the intercept
        if any([isa(fe.interaction, Ones) for fe in fes])
            rt.intercept = false
            has_intercept = true
        end
        fes = FixedEffect[_subset(fe, esample) for fe in fes]
        pfe = FixedEffectMatrix(fes, sqrtw, Val{method})
    end

    # Compute data for std errors
    vcov_method_data = VcovMethod(df[esample, unique(Symbol.(vcov_vars))], vcovformula)

    ##############################################################################
    ##
    ## Dataframe --> Matrix
    ##
    ##############################################################################

    mf = ModelFrame2(rt, df, esample; contrasts = contrasts)

    # Obtain y
    # for a Vector{Float64}, conver(Vector{Float64}, y) aliases y
    y = convert(Vector{Float64}, model_response(mf)[:])
    yname = rt.eterms[1]
    y .= y .* sqrtw

    # Obtain X
    coef_names = coefnames(mf)
    has_exo = !isempty(mf.terms.terms) | mf.terms.intercept
    if has_exo
        Xexo = ModelMatrix(mf).m
        Xexo .= Xexo .* sqrtw
    else
        Xexo = Matrix{Float64}(undef, nobs, 0)
    end

    if has_iv
        mf = ModelFrame2(endo_terms, df, esample)
        coef_names = vcat(coef_names, coefnames(mf))
        Xendo = ModelMatrix(mf).m
        Xendo .= Xendo .* sqrtw

        mf = ModelFrame2(iv_terms, df, esample)
        Z = ModelMatrix(mf).m
        Z .= Z .* sqrtw
    else
        Xendo = Matrix{Float64}(undef, nobs, 0)
        Z = Matrix{Float64}(undef, nobs, 0)
    end

    # compute tss now before potentially demeaning y
    tss = compute_tss(y, has_intercept, sqrtw)

    if has_absorb
        # used to compute tss even without save_fe
        if save_fe
            oldy = deepcopy(y)
            oldX = hcat(Xexo, Xendo)
        end

        # initialize iterations and converged
        iterations = Int[]
        convergeds = Bool[]

        y, b, c = solve_residuals!(y, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        Xexo, b, c = solve_residuals!(Xexo, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        Xendo, b, c = solve_residuals!(Xendo, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        Z, b, c = solve_residuals!(Z, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        iterations = maximum(iterations)
        converged = all(convergeds)
        if converged == false
            @warn "convergence not achieved in $(iterations) iterations; try increasing maxiter or decreasing tol."
        end
    end

    ##############################################################################
    ##
    ## Get Linearly Independent Components of Matrix
    ##
    ##############################################################################

    # Compute linearly independent columns + create the Xhat matrix
    if has_iv
        if size(Z, 2) < size(Xendo, 2)
            error("Model not identified. There must be at least as many ivs as endogeneneous variables")
        end
        # get linearly independent columns
        # note that I do it after residualizing
        baseall = basecol(Z, Xexo, Xendo)
        basecolXexo = baseall[(size(Z, 2)+1):(size(Z, 2) + size(Xexo, 2))]
        basecolXendo = baseall[(size(Z, 2) + size(Xexo, 2) + 1):end]
        Z = getcols(Z, baseall[1:size(Z, 2)])
        Xexo = getcols(Xexo, basecolXexo)
        Xendo = getcols(Xendo, basecolXendo)
        basecoef = vcat(basecolXexo, basecolXendo)

        # Build
        X = hcat(Xexo, Xendo)
        newZ = hcat(Xexo, Z)
        crossz = cholesky!(Symmetric(newZ' * newZ))
        Pi = crossz \ (newZ' * Xendo)
        Xhat = hcat(Xexo, newZ * Pi)

        # prepare residuals used for first stage F statistic
        ## partial out Xendo in place wrt (Xexo, Z)
        Xendo_res = gemm!('N', 'N', -1.0, newZ, Pi, 1.0, Xendo)
        ## partial out Z in place wrt Xexo
        Pi2 = cholesky!(Symmetric(Xexo' * Xexo)) \ (Xexo' * Z)
        Z_res = gemm!('N', 'N', -1.0, Xexo, Pi2, 1.0, Z)

    else
        # get linearly independent columns
        basecolXexo = basecol(Xexo)
        Xexo = getcols(Xexo, basecolXexo)
        Xhat = Xexo
        X = Xexo
        basecoef = basecolXexo
    end

    ##############################################################################
    ##
    ## Do the regression
    ##
    ##############################################################################

    crossx =  cholesky!(Symmetric(Xhat' * Xhat))
    coef = crossx \ (Xhat' * y)
    residuals = y - X * coef

    ##############################################################################
    ##
    ## Optionally save objects in a new dataframe
    ##
    ##############################################################################

    augmentdf = DataFrame()
    if save_residuals
        augmentdf[:residuals] =  Vector{Union{Missing, Float64}}(missing, length(esample))
        augmentdf[esample, :residuals] = residuals ./ sqrtw
    end
    if save_fe
        oldX = getcols(oldX, basecoef)
        newfes, b, c = solve_coefficients!(oldy - oldX * coef, pfe; tol = tol, maxiter = maxiter)
        for j in 1:length(fes)
            augmentdf[ids[j]] = Vector{Union{Float64, Missing}}(missing, length(esample))
            augmentdf[esample, ids[j]] = newfes[j]
        end
    end

    ##############################################################################
    ##
    ## Test Statistics
    ##
    ##############################################################################

    # Compute degrees of freedom
    df_intercept = 0
    if has_absorb || rt.intercept
        df_intercept = 1
    end
    df_absorb = 0
    if has_absorb
        for fe in fes
            # adjust degree of freedom only if fe is not fully nested in a cluster variable:
            if isa(vcovformula, VcovClusterFormula)
                if any(isnested(fe, vcov_method_data.clusters[v]) for v in names(vcov_method_data.clusters))
                    df_absorb = 1 # if fe is nested you still lose 1 degree of freedom
                    break
                end
            end
            #only count groups that exists
            df_absorb += length(Set(fe.refs))
        end
    end
    dof_residual = max(1, nobs - size(X, 2) - df_absorb - df_add)

    # Compute rss, tss, r2, r2 adjusted
    rss = sum(abs2, residuals)
    mss = tss - rss
    r2 = 1 - rss / tss
    adjr2 = 1 - rss / tss * (nobs - has_intercept) / dof_residual
    if has_absorb
        r2_within = 1 - rss / compute_tss(y, rt.intercept, sqrtw)
    end

    # Compute standard error
    vcov_data = VcovData(Xhat, crossx, residuals, dof_residual)
    matrix_vcov = vcov!(vcov_method_data, vcov_data)

    # Compute Fstat
    (F, p) = compute_Fstat(coef, matrix_vcov, nobs, rt.intercept, vcov_method_data, vcov_data)

    # Compute Fstat of First Stage
    if has_iv
        Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
        (F_kp, p_kp) = ranktest!(Xendo_res, Z_res, Pip,
                                  vcov_method_data, size(X, 2), df_absorb)
    end

    ##############################################################################
    ##
    ## Return regression result
    ##
    ##############################################################################

    # add omitted variables
    if !all(basecoef)
        newcoef = fill(zero(Float64), length(basecoef))
        newmatrix_vcov = fill(NaN, (length(basecoef), length(basecoef)))
        newindex = [searchsortedfirst(cumsum(basecoef), i) for i in 1:length(coef)]
        for i in 1:length(coef)
            newcoef[newindex[i]] = coef[i]
            for j in 1:length(coef)
                newmatrix_vcov[newindex[i], newindex[j]] = matrix_vcov[i, j]
            end
        end
        coef = newcoef
        matrix_vcov = newmatrix_vcov
    end

    # return
    if !has_iv && !has_absorb
        return RegressionResult(coef, matrix_vcov, esample, augmentdf,
                                coef_names, yname, f, nobs, dof_residual,
                                rss, tss, r2, adjr2, F, p)
    elseif has_iv && !has_absorb
        return RegressionResultIV(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, nobs, dof_residual,
                                  rss, tss,  r2, adjr2, F, p, F_kp, p_kp)
    elseif !has_iv && has_absorb
        return RegressionResultFE(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, feformula, nobs, dof_residual,
                                  rss, tss, r2, adjr2, r2_within, F, p, iterations, converged)
    elseif has_iv && has_absorb
        return RegressionResultFEIV(coef, matrix_vcov, esample, augmentdf,
                                   coef_names, yname, f, feformula, nobs, dof_residual,
                                   rss, tss, r2, adjr2, r2_within, F, p, F_kp, p_kp,
                                   iterations, converged)
    end
end

##############################################################################
##
## Fstat
##
##############################################################################

function compute_Fstat(coef::Vector{Float64}, matrix_vcov::Matrix{Float64},
    nobs::Int, hasintercept::Bool,
    vcov_method_data::AbstractVcovMethod, vcov_data::VcovData)
    coefF = copy(coef)
    # TODO: check I can't do better
    length(coef) == hasintercept && return NaN, NaN
    if hasintercept
        coefF = coefF[2:end]
        matrix_vcov = matrix_vcov[2:end, 2:end]
    end
    F = (Diagonal(coefF) * (matrix_vcov \ Diagonal(coefF)))[1]
    df_ans = df_FStat(vcov_method_data, vcov_data, hasintercept)
    dist = FDist(nobs - hasintercept, max(df_ans, 1))
    return F, ccdf(dist, F)
end

function compute_tss(y::Vector{Float64}, hasintercept::Bool, sqrtw::AbstractVector)
    if hasintercept
        tss = zero(Float64)
        m = (mean(y) / sum(sqrtw) * length(y))::Float64
        @inbounds @simd for i in 1:length(y)
            tss += abs2(y[i] - sqrtw[i] * m)
        end
    else
        tss = sum(abs2, y)
    end
    return tss
end

##############################################################################
##
## Syntax without keywords
##
##############################################################################
function evaluate_subset(df, ex::Expr)
    if ex.head == :call
        return Expr(ex.head, ex.args[1], (evaluate_subset(df, ex.args[i]) for i in 2:length(ex.args))...)
    else
        return Expr(ex.head, (evaluate_subset(df, ex.args[i]) for i in 1:length(ex.args))...)
    end
end
evaluate_subset(df, ex::Symbol) = df[ex]
evaluate_subset(df, ex)  = ex


text:
using BinDeps

@BinDeps.setup

mpg123 = library_dependency("libmpg123")
lame = library_dependency("libmp3lame")

provides(AptGet, "libmpg123-dev", mpg123)
provides(AptGet, "libmp3lame-dev", lame)

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides(Homebrew.HB, "mpg123", mpg123)
    provides(Homebrew.HB, "lame", lame)
end

@BinDeps.install Dict("libmpg123" => :libmpg123, "libmp3lame" => :libmp3lame)


text:
__precompile__(true)
"""
    UnitfulBrewing

Module extending Unitful.jl with beer brewing units.
"""
module UnitfulBrewing

using Unitful

using UnitfulEquivalences: Equivalence, dimtype, @eqrelation
import UnitfulEquivalences: edconvert

export DensityConcentration, SugarGravity, SugarGravityQuad

# New dimensions
@dimension 𝐂𝐋   "𝐂𝐋"     Color
@dimension 𝐃𝐏   "𝐃𝐏"     DiastaticPower
@dimension 𝐁𝐓   "𝐁𝐓"     Bitterness
@dimension 𝐒𝐂   "𝐒𝐂"     SugarContents

# Time units: adding beerjson symbols as alias to Unitful symbols
const sec = Unitful.s
const min = Unitful.minute
const day = Unitful.d
const week = Unitful.wk

# US volumes
@unit gal       "gal"       Gallon      231*(Unitful.inch)^3    false
@unit qt        "qt"        Quart       gal//4                  false
@unit pt        "pt"        Pint        qt//2                   false
@unit floz      "floz"      FluidOunce  pt//16                  false
@unit bbl       "bbl"       Barrel      42gal                   false

# Imperial volumes
@unit ifloz     "ifloz"     ImperialFluidOunce  28.4130625*(Unitful.mm) false
@unit gi        "gi"        Gill                5ifloz                  false
@unit ipt       "ipt"       ImperialPint        20ifloz                 false
@unit iqt       "iqt"       ImperialQuart       2ipt                    false
@unit igal      "igal"      ImperialGallon      8ipt                    false
@unit ibbl      "ibbl"      ImperialBarrel      36igal                  false

# Culinary volumes
@unit cup       "cup"       Cup         pt//2                   false
@unit tbsp      "tbsp"      Tablespoon  cup//16                 false
@unit tsp       "tsp"       Teaspoon    tbsp//3                 false

# Sugar content
@refunit °P     "°P"        Plato               𝐒𝐂          false
@unit Brix      "Brix"      Brix                1°P         false
@unit Balling   "Balling"   Balling             1°P         false
const P = °P

# Specific gravity
# uconvert(Unitful.NoUnits, 1.010u"sg") == 1.01
# uconvert(u"permille", 1.040u"sg") == 1040.0 ‰
# uconvert(u"gu", 1.040u"sg") == 40.0 gu
# uconvert(u"sg", 40u"gu") == 1.04 sg
@unit       sg  "sg"        SpecificGravity     1.0         false
@affineunit gu  "gu"        1000.0 * Unitful.permille
const gp = gu # gravity points

# Diastatic Power
@refunit        °Lintner    "°Lintner"  Lintner     𝐃𝐏                  false
@unit           °WK_aux     "°WK_aux"   WK_aux      (10//35)°Lintner    false
@affineunit     °WK         "°WK"       16°WK_aux
const Lintner = °Lintner
const WK = °WK

# Color units
@refunit    SRM     "SRM"       SRM                 𝐂𝐋                  false
@unit       EBC     "EBC"       EBC                 (100//197)SRM       false
@unit       L_aux   "L_aux"     L_aux               (13546//10000)SRM    false
@affineunit °L      "°L"        -(7600//13546)L_aux
const srm = SRM
const Lovi = °L
const Lovibond = °L
const ebc = EBC

# Carbonation units

# International Bitterness Unit

@refunit    IBU     "IBU"       InternationalBitternessUnit 𝐁𝐓      false

# Concentration units
#=
Unitful already has percent, permille and pertenthousand,
so we only add ppm, ppb, and ppt
@unit percent         "%"    Percent         1//100             false
@unit permille        "‰"    Permille        1//1000            false
@unit pertenthousand  "‱"    Pertenthousand  1//10000           false
=#
@unit       ppm     "ppm"       PartsPerMillion     1//10^6         false
@unit       ppb     "ppb"       PartsPerBillion     1//10^9         false
@unit       ppt     "ppt"       PartsPerTrillion    1//10^9         false

# pH logarithmic scale
@logscale pH⁺    "pH⁺"       powerofHydrogen      10      10      false
const pwrH = pH⁺ # either pH\^+<ESC> or, with dead keys, pH\^<SPACE>+<ESC>
#=
Using symbols pwrH and pH⁺ since there is a dimensional symbol pH already
defined in Unitful:
```julia-repl
julia> typeof(Unitful.pH)
Unitful.FreeUnits{(pH,),𝐋² 𝐌 𝐈⁻² 𝐓⁻²,nothing}
```
=#

## Define the conversion functions between Plato and gravity units
"""
    gu_to_plato(gu::Number)

Convert a value in Gravity Units to degrees Plato according
to the quadratic formula

    Plato = 0.25802gu - 0.00020535gu^2,

which is equivalent to the formula for specific gravity

    Plato = 668.72 * sg - 463.37 - 205.35 * sg^2,

with
    gu = 1000 ( sg - 1.000 ).
"""
gu_to_plato(gu::Number) = 0.25802gu - 0.00020535gu^2

"""
    plato_to_gu(p::Number)

Convert a value in degrees Plato to Gravity Units by inverting the quadratic
formula for degrees Plato, so that

    gu = e - sqrt(e^2 - g * Plato)) if ^2 - g * p >= 0 else gu = e

where
    e = 0.25802 / 0.00020535 / 2 = 628.2444606768931
and
    g = 1 / 0.00020535 = 4869.734599464329

The value gu = d when p > d/2 is just for definiteness of the function since
the in this range the conversion is meaningless.
"""
function plato_to_gu(p::Number)
    e = 628.2444606768931
    g = 4869.734599464329
    d = e^2 - g*p
    if d >= 0
        e - sqrt(d)
    else
        e
    end
end

# Define the equivalences

"""
    DensityConcentration()

Equivalence to convert between Density and Concentration for water-based
solutions, with relatively small quantities of solutes, as in water treatment,
so that 1 mg/L is equivalent to 1 ppm (parts per million).

# Examples

```jldoctest
julia> uconvert(u"mg/l", 10u"ppm", DensityConcentration())
10 mg L⁻¹
julia> uconvert(u"ppm", 1u"g/l", DensityConcentration())
1000 ppm
```
"""
struct DensityConcentration <: Equivalence end

@eqrelation DensityConcentration Unitful.Density / Unitful.DimensionlessQuantity = 1u"kg/l"

"""
    SugarGravity()

Equivalence to convert between Sugar Contents and Specific Gravity quantities.

Convert between degrees Plato and specific gravity according to the equation

    Plato = 259 * (1 - 1/sg),

# Examples

```jldoctest
julia> uconvert(u"°P", 1.040u"sg", SugarGravity())
9.961538461538483 °P
julia> uconvert(u"sg", 15u"°P", SugarGravity())
1.0614754098360655 sg
julia> uconvert(u"gu", 12u"°P", SugarGravity())
48.58299595141693 gu
julia> uconvert(u"°P", 40u"gu", SugarGravity())
9.961538461538483 °P
```
"""
struct SugarGravity <: Equivalence end

edconvert(::dimtype(SugarContents), x::Unitful.DimensionlessQuantity,
    ::SugarGravity) = 259 * (1 - 1 / uconvert(sg, x).val) * °P

edconvert(::dimtype(Unitful.DimensionlessQuantity), x::SugarContents,
    ::SugarGravity) = 259 / (259 - uconvert(°P, x).val) * sg

"""
    SugarGravityQuad()

Equivalence to convert between Sugar Contents and Specific Gravity quantities.

Convert between degrees Plato and gravity units according to the quadratic
equation

    Plato = 0.25802gu - 0.00020535gu^2,

which is equivalent to the formula for specific gravity

    Plato = 668.72 * sg - 463.37 - 205.35 * sg^2,

with

    gu = 1000 ( sg - 1.000 ).

# Examples

```jldoctest
julia> uconvert(u"°P", 1.040u"sg", SugarGravityQuad())
9.992240000000002 °P
julia> uconvert(u"sg", 15u"°P", SugarGravityQuad())
1.0611068377146748 sg
julia> uconvert(u"gu", 12u"°P", SugarGravityQuad())
48.370088784473296 gu
julia> uconvert(u"°P", 40u"gu", SugarGravityQuad())
9.992240000000002 °P
```
"""
struct SugarGravityQuad <: Equivalence end

edconvert(::dimtype(SugarContents), x::Unitful.DimensionlessQuantity,
    ::SugarGravityQuad) = gu_to_plato(uconvert(gu, x).val) * °P

edconvert(::dimtype(Unitful.DimensionlessQuantity), x::SugarContents,
    ::SugarGravityQuad) = plato_to_gu(uconvert(°P, x).val) * gu

# The function below is just so I get things straight
function show_quantity_info(x::Quantity{T,D,U}) where {T,D,U}
    println("Here is the result of (T, D, U, U()) for $x:")
    return T, D, U, U()
end

# Register the above units and dimensions in Unitful
const localpromotion = Unitful.promotion # only needed with new dimensions
function __init__()
    Unitful.register(UnitfulBrewing) # needed for new Units
    merge!(Unitful.promotion, localpromotion) # only needed with new dimensions
end

end # module


text:
#Setup pre-evalutated utility matrix and the utility-maximum vector
function setuputilityarrays(a::Vector, w::Vector, utility::Function)
    cardinality_a = length(a)
    cardinality_w = length(w)

    #pre-compute utilities, find maxima
    U_pre = zeros(cardinality_a, cardinality_w)
    Umax = zeros(cardinality_w)
    for i in 1:cardinality_w
        for j in 1:cardinality_a
            U_pre[j,i]=utility(a[j], w[i])
        end
        Umax[i] = maximum(U_pre[:,i])
    end

    return U_pre, Umax
end

#This function computes p_boltz = 1/Z * p0 * exp(β*ΔU),
#where Z is the normalization constant (partition function)
### arguments:
#p0 ... prior distribution (vector of length N)
#β ... inverse temperature (scalar)
#ΔU ... utility or potential-difference (vector of length N)
### returns:
#p_boltz ... 1/Z * p0 * exp(β*ΔU)
function boltzmanndist(p0::Vector, β, ΔU::Vector)
    p_boltz = p0.*exp(β.*ΔU)
    p_boltz /= sum(p_boltz)
    return p_boltz

    #=
    #do the computation with more precision (but slower)
    deltaU = convert(Array{BigFloat}, ΔU)
    p_boltz = p0.*exp(β.*ΔU)
    p_boltz = p_boltz/sum(p_boltz)
    return float64(p_boltz)
    =#

end

#This function performs Blahut-Arimoto iterations
function BAiterations(pa_init::Vector, β, U_pre::Matrix, pw::Vector, ε_conv::Real, maxiter::Integer;
                       compute_performance::Bool=false, performance_per_iteration::Bool=false,
                       performance_as_dataframe::Bool=false)
    pa_new = pa_init
    card_a = size(U_pre,1)
    card_w = size(U_pre,2)
    pagw_new = zeros(card_a,card_w)

    #if performance measures don't need to be returned, don't compute them per iteration
    if compute_performance==false
        performance_per_iteration = false
    end

    #preallocate if necessary
    if performance_per_iteration
        I_i = zeros(maxiter)
        Ha_i = zeros(maxiter)
        Hagw_i = zeros(maxiter)
        EU_i = zeros(maxiter)
        RDobj_i = zeros(maxiter)
    end

    #main iteration
    iter = 0 #initialize counter, so it persists beyond the loop
    for iter in 1:maxiter
        pa = deepcopy(pa_new)  #make sure not to just copy the reference
        #pa_new = zeros(card_a)
        pagw = deepcopy(pagw_new)
        pagw_new = zeros(card_a,card_w)

        for k in 1:card_w
            #update p(a|o)
            pagw_new[:,k] = boltzmanndist(pa,β,vec(U_pre[:,k]))
            #update p(a)
            pa_new = pa_new + vec((pagw_new[:,k]')*pw[k])
        end

        #add small value to pa_new to make sure there are no zero-entries
        #due to limited numerical precision, then re-normalize
        pa_new += eps()
        pa_new /= sum(pa_new)

        #compute entropic quantities (if requested with additional parameter)
        if performance_per_iteration
            I_i[iter], Ha_i[iter], Hagw_i[iter], EU_i[iter], RDobj_i[iter] = analyzeBAsolution(pw, pa_new, pagw_new, U_pre, β)
        end

        #check for convergence
        if norm(pagw-pagw_new) < ε_conv
            break
        end
    end

    #check if iteration limit has been reached (before convergence)
    if iter == maxiter
        warn("[BAiterations] maximum iteration reached - returning... (results might be inaccurate)")
    end

    #return results
    if compute_performance == false
        return pagw_new, vec(pa_new)  #the squeeze will turn pa into a vector again
    else
        if performance_per_iteration == false
            #compute performance measures for final solution
            I, Ha, Hagw, EU, RDobj = analyzeBAsolution(pw, pa_new, pagw_new, U_pre, β)
        else
            #"cut" valid results from preallocated vector
            I = I_i[1:iter]
            Ha = Ha_i[1:iter]
            Hagw = Hagw_i[1:iter]
            EU = EU_i[1:iter]
            RDobj = RDobj_i[1:iter]
        end

        #if needed, transform to data frame
        if performance_as_dataframe == false
            return pagw_new, vec(pa_new), I, Ha, Hagw, EU, RDobj
        else
            performance_df = performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj)
            return pagw_new, vec(pa_new), performance_df
        end
    end

end


text:
####################
#
# Polyline
#
####################

# Union of Shapefile Polyline types
ShpPolyline = Union{Shapefile.Polyline, Shapefile.PolylineZ, Shapefile.PolylineM}

#=
For an array of ESRIShape points
1. convert to a different coordinate system (optional)
2. convert to a Compose line

Note:
Units in Compose can be either:
- Context - cx, cy - using the axis coordinates
- Relative - w (width), h (height) - using the size of the image
- Absolute - pt, mm, cm, inch
If you just provide numbers e.g. if pts below was an array of tuples of Float64
the units would default to context units.
=#
function ESRItoComposeLine(points::AbstractArray{T, 1},
                           convertcoords) where {T<:ShpPoint}
    pts = Array{Tuple{Measures.Length{:cx, Float64}, Measures.Length{:cy, Float64}}}(undef, 0)
    for pt in points
        x, y = convertcoords(pt.x, pt.y)
        push!(pts, (x * cx, y * cy))
    end
    Compose.line(pts)
end

# Plot a single part polyline
function compose_single(polyline::T, canvas, convertcoords, line_width,
                        line_color) where {T<:ShpPolyline}
    compose(canvas, (context(), ESRItoComposeLine(polyline.points, convertcoords),
                     linewidth(line_width), stroke(line_color)))
end

# Plot a multi part polyline
function compose_multi(polyline::T, canvas, convertcoords, line_width,
                       line_color) where {T<:ShpPolyline}
    start = 1
    finish = length(polyline.points)
    for idx in polyline.parts[2:end]
        canvas = compose(canvas, (context(),
                                  ESRItoComposeLine(polyline.points[start:idx], convertcoords),
                                  linewidth(line_width), stroke(line_color)))
        start = idx + 1
    end
    canvas = compose(canvas, (context(),
                              ESRItoComposeLine(polyline.points[start:finish], convertcoords),
                              linewidth(line_width), stroke(line_color)))
end

# Plot an array of polylines
function draw_shp(shapes::AbstractArray{T, 1}, canvas, convertcoords,
                  line_width, line_color, fill_color,
                  radius) where {T<:ShpPolyline}
    for polyline in shapes
        if length(polyline.parts) == 1
            canvas = compose_single(polyline, canvas, convertcoords, line_width, line_color)
        else
            canvas = compose_multi(polyline, canvas, convertcoords, line_width, line_color)
        end
    end
    canvas
end


text:
#!/usr/bin/env julia

using Luxor
using Test
using Random
Random.seed!(42)

function bezierarithmetic()
    bps = BezierPathSegment(ngon(O, 100, 4, vertices=true)...)
    @test length(bps) == 4
    @test Luxor.get_bezier_length(bps) ≈ 281.4703644162024
    @test Luxor.get_bezier_length(bps, steps=100) ≈ 282.8313796714682
    l, r = splitbezier(bps, 0.5)
    @test last(l) ≈ first(r)
    @test first(l) ≈ first(bps)
    @test last(r) ≈ last(bps)
end

function bezierpathtest(fname)
    Random.seed!(3)
    currentwidth = 1200
    currentheight = 1200
    Drawing(currentwidth, currentheight, fname)
    origin()
    background(setgray(0.8)...)
    tiles = Tiler(currentwidth, currentheight, 4, 4, margin=20)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            pts = randompointarray(Point(-tiles.tilewidth/2, -tiles.tilewidth/2), Point(tiles.tilewidth/2, tiles.tilewidth/2), rand(3:5))
            setopacity(0.7)
            sethue("black")
            prettypoly(pts, :stroke, close=true)
            randomhue()
            drawbezierpath(makebezierpath(pts, smoothing=rand(0:0.1:1)), :fill)
        end
    end
    @test finish() == true
end

function pathtobezierpath(fname)
    Random.seed!(3)
    currentwidth = 800
    currentheight = 800
    Drawing(currentwidth, currentheight, fname)
    origin()
    background("white")
    st = "julia"
    setopacity(0.7)
    thefontsize = 200
    sethue("red")
    fontsize(thefontsize)
    move(Point(-350, 0))
    textpath(st)
    nbps = pathtobezierpaths()
    for nbp in nbps
        setline(.25)
        sethue("grey50")
        drawbezierpath(nbp, :stroke)
        # exaggerate position of Bezier control points
        for p in nbp
            sethue("magenta")
            exthandle1 = between(p[2], p[1], -3)
            exthandle2 = between(p[3], p[4], -3)
            circle(exthandle1, 1, :fill)
            circle(exthandle2, 1, :fill)
            line(p[1], exthandle1, :stroke)
            line(exthandle2, p[4], :stroke)
            if p[1] != p[4]
                sethue("black")
                circle(p[1], 1, :fill)
                circle(p[4], 1, :fill)
            end
        end
    end
    @test finish() == true
end

bezierarithmetic()

fname = "bezierpath.png"
bezierpathtest(fname)

println("...bezier test: output in $(fname)")

fname = "pathtobezierpath.svg"
pathtobezierpath(fname)

println("...finished bezier tests: output in $(fname)")


text:
function build_channel(ch::IonChannel; name = get_name(ch))
    @parameters t
    D = Differential(t)
    @variables V(t) Ca(t)
    eqs, states, parameters, current, defaultmap = channel_dynamics(ch, V, Ca, D, t)
    return ODESystem(eqs, t, [V, Ca, states...], [parameters...];
                                            observed=current,
                                            name = name,
                                            defaults = defaultmap)
end

function build_channel(syn::Synapse; name = get_name(syn))
    @parameters t
    D = Differential(t)
    @variables Vpre(t) Vpost(t)
    eqs, states, parameters, current, defaultmap = channel_dynamics(syn, Vpre, Vpost, D, t)
    return ODESystem(eqs, t, [Vpre, Vpost, states...], [parameters...];    observed = current,
                        name = name,
                        defaults = defaultmap)
end

# 10 nF/mm² for default specific membrane capacitance
function build_neuron(channels;  capacitance = 10, V_init = -60, Ca_init = 0.05, name = :unidentified_neuron, num_inputs = 0)

    @parameters t cm
    states = @variables V(t) Ca(t)
    syns = [Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(:Isyn, i)))(t) for i in 1:num_inputs]

    D = Differential(t)
    channel_systems = build_channel.(channels)

    summed_membrane_currents = sum(ionic_current(ch, ch_sys) for (ch, ch_sys) in zip(channels, channel_systems))

    summed_calcium_flux = sum(calcium_current(ch, ch_sys) for (ch, ch_sys) in zip(channels, channel_systems))

    connections = cat(
                    [voltage_hook(V, ch, chs) for (ch, chs) in zip(channels, channel_systems)],
                    [calcium_hook(Ca, ch, chs) for (ch, chs) in zip(channels, channel_systems)],
                    dims=1
                    )

    diffeqs =   cat(
                    [D(V) ~ (summed_membrane_currents + sum(syns))/cm],
                [D(Ca) ~ (summed_calcium_flux*.0939 - Ca + .05)/200],
                dims=1)

    eqs = cat(connections, diffeqs; dims=1)
    eqs = filter(x -> !(x == (Num(0) ~ Num(0))), eqs)

    all_states = [states..., syns...]

    neur = ODESystem(eqs, t, all_states, [cm];
                    systems = channel_systems,
                    defaults = [V => V_init, Ca => Ca_init, cm => capacitance],
                    name = name
                    )

    return neur
end

function build_group(arr_of_neurons; name = :group)
    @parameters t
    eqs = Array{Equation,1}(undef,0)
    ODESystem(eqs, t, [], []; name=name, systems = arr_of_neurons)
end

"""
Issue: need to do all connections in one ḡChol

"""
function add_connection(group, pre, post, syn::Synapse; name = ModelingToolkit.getname(group), i=1)

    @parameters t
    prename, postname = ModelingToolkit.getname.([pre, post])
    synapse = build_channel(syn; name = Symbol(prename, :to, postname, get_name(syn)))

    oldeqs = ModelingToolkit.get_eqs(group)
    neweqs = [
            get_pre(syn, pre) ~ my_pre(syn, synapse),
            get_post(syn, post) ~ my_post(syn, synapse),
            ionic_current(syn, synapse) ~ getproperty(post, Symbol(post_connector(syn), i))
    ]

    eqs = cat(oldeqs, neweqs, dims=1)
    systems = cat(ModelingToolkit.get_systems(group), synapse, dims=1)

    return connected = ODESystem(eqs, t, [], [];
                                name=name,
                                systems = systems)
end

max_stars_repo_name: sg-s/NeuronBuilder

max_issues_repo_name: sg-s/NeuronBuilder

max_forks_repo_name: sg-s/NeuronBuilder


text:
using Test, MLToolkit.PlotsX, Plots
pyplot()
theme(:ggplot2; size=(400, 300))

@testset "Plots" begin
    @testset "autoget_lims" begin
        x = [0 1; 0 1]
        xlims, ylims, zlims = autoget_lims(x)
        @test xlims == [-0.1, 1.1]
        @test ylims == [-0.1, 1.1]
        @test zlims == nothing

        x = [0 1; 0 1; 0 1]
        xlims, ylims, zlims = autoget_lims(x)
        @test xlims == [-0.1, 1.1]
        @test ylims == [-0.1, 1.1]
        @test zlims == [-0.1, 1.1]
    end

    @testset "make_imggrid" begin
        d = 5
        rand_img = rand(d, d, 100)
        n_rows = n_cols = 3
        gap = 1
        imggrid = make_imggrid(rand_img, n_rows, n_cols)
        @test all(imggrid[1,:] .== 0.5)
        @test all(imggrid[end,:] .== 0.5)
        @test all(imggrid[:,1] .== 0.5)
        @test all(imggrid[:,end] .== 0.5)
        @test size(imggrid) == (
            n_rows * d + gap * (n_rows + 1),
            n_cols * d + gap * (n_cols + 1)
        )
    end

    @testset "count_leadingzeros" begin
        @test count_leadingzeros([0, 0, 1]) == 2
    end

    function test_save(p, fname; do_savetex=true)
        save(joinpath(@__DIR__, "$fname.png"), p)
        do_savetex && save(joinpath(@__DIR__, "$fname.tex"), p)
    end

    @testset "OneDimFunction" begin
        using Distributions
        dist = Normal(0, 1)

        p = plot(dist)

        test_save(p, "OneDimFunction")
    end

    @testset "ErrorBarLines" begin
        x = 1:10
        ys = [sin.(x) + rand(10) for i in 1:3]
        lines = ErrorBarLines(x, ys)

        p = plot(lines; nstd=2, label="Noisy sin(x) with 2 std", legend=:topright)

        test_save(p, "ErrorBarLines")
    end

    @testset "TwoYAxesLines" begin
        x = collect(1:0.5:10)
        y1 = sin.(x)
        y2 = x .^ 2

        lines = TwoYAxesLines(x, y1, y2)
        p = plot(lines; xlabel="x", ylabel1="sin(x)", ylabel2="x^2")

        test_save(p, "TwoYAxesLines")
    end

    @testset "TwoDimFunction" begin
        using Distributions
        dist = MvNormal(zeros(2), 1)

        p = plot(dist)

        test_save(p, "TwoDimFunction")
    end

    @testset "TwoDimPath" begin
        x = collect(1:10) + rand(10)
        y = collect(1:10) + rand(10)
        path = TwoDimPath(x, y)

        p = plot(path)

        test_save(p, "TwoDimPath")
    end

    @testset "LabelledScatter" begin
        ls = LabelledScatter((x=rand(2, 100), y=rand(2, 100)))

        p = plot(ls)

        test_save(p, "LabelledScatter")
    end

    @testset "ImageGrid" begin
        ig = ImageGrid(rand(32, 32, 3, 100))

        p = plot(ig)

        test_save(p, "ImageGrid")
    end

    @testset "FeatureActivations" begin
        fa = FeatureActivations(rand(Bool, 20, 50))

        p = plot(fa)

        test_save(p, "FeatureActivations")
    end
end


text:
# # Bouncing ball

# ## Model

# We model the bouncing ball as a hybrid automaton with one location and a self-loop.
# See for example [^LG09] pp. 79-83.

module bouncing_ball #jl

using ReachabilityAnalysis, ModelingToolkit

function _bouncing_ball()

    var = @variables x v

    ## "falling" mode with invariant x >= 0
    invariant = HalfSpace(x ≥ 0, var)
    flow = @system(z' = [0.0 1.0; 0.0 0.0] * z + [0.0, -9.81], z ∈ invariant)

    ## guard x == 0 && v ≤ 0
    guard = HPolyhedron([0 ≤ x, x ≤ 0, v ≤ 0], var)

    ## reset map v⁺ := -cv
    assignment = ConstrainedLinearMap([1.0 0.0; 0.0 -0.75], guard)

    ## initial-value problem
    return HybridSystem(flow, assignment)
end

function model(X0)
    H = _bouncing_ball()
    return @ivp(H, z(0) ∈ X0)
end

end #jl

# ## References

# [^LG09]: Le Guernic, Colas. Reachability analysis of hybrid systems with linear continuous dynamics. Diss. 2009.

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# ---
# title: 462. Minimum Moves to Equal Array Elements II
# id: problem462
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/minimum-moves-to-equal-array-elements-ii/description/>
# hidden: true
# ---
#
# Given a **non-empty** integer array, find the minimum number of moves required
# to make all array elements equal, where a move is incrementing a selected
# element by 1 or decrementing a selected element by 1.
#
# You may assume the array's length is at most 10,000.
#
# **Example:**
#
#
#
#     Input:
#     [1,2,3]
#
#     Output:
#     2
#
#     Explanation:
#     Only two moves are needed (remember each move increments or decrements one element):
#
#     [1,2,3]  =>  [2,2,3]  =>  [2,2,2]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
A = [  1.0  -2.0   3.0  -4.0  5.0;
       2.0  -1.0   0.0  -1.0  3.0;
      -1.0   0.0   4.0  -3.0  2.0;
      -1.0  -1.0  -1.0   1.0  3.0]
b = [1.0, 2.0, 3.0, 4.0]

m, n = size(A)

f = Translate(SqrNormL2(), -b)
f2 = LeastSquares(A, b)
lam = 0.1*norm(A'*b, Inf)
g = NormL1(lam)

x_star = [-3.877278911564627e-01, 0, 0, 2.174149659863943e-02, 6.168435374149660e-01]

## Nonfast/Nonadaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 150
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# Nonfast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 300
println(sol)

# Fast/Nonadaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 100
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# Fast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, adaptive=true, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 200
println(sol)

# ZeroFPR/Nonadaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15
println(sol)

#testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# ZeroFPR/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15

# PANOC/Nonadaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

## PANOC/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)

# Douglas-Rachford

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.DRS(x0; f=f2, g=g, gamma=10.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 30
println(sol)

#testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

#####################
# Complex Variables #
#####################

Random.seed!(123)
m, n = 10,5

A = randn(m,n)+im*randn(m,n)
b = randn(m)+im*randn(m)

f = Translate(SqrNormL2(), -b)
f2 = LeastSquares(A, b)
lam = 0.01*norm(A'*b, Inf)
g = NormL1(lam)

x0 = zeros(n)+im*zeros(n)
it, x_star, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2, verbose = 0)

## Nonfast/Nonadaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 200
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# Nonfast/Adaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 250
println(sol)

# Fast/Nonadaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 120
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# Fast/Adaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, adaptive=true, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 120
println(sol)

# ZeroFPR/Nonadaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15
println(sol)

#testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

# ZeroFPR/Adaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15

# PANOC/Nonadaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)

# testing solver already at solution
@time it, x = ProximalAlgorithms.run!(sol)
@test it == 1

## PANOC/Adaptive

x0 = zeros(n)+im*zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)

#############################################
# Real Variables  mapped to Complex Numbers #
#############################################

Random.seed!(123)
n = 2^6
A = AbstractOperators.DFT(n)[1:div(n,2)]      # overcomplete dictionary

x = sprandn(n,0.5)
b = fft(x)[1:div(n,2)]

#f = Translate(LogisticLoss(ones(n)), -b)
f = Translate(SqrNormL2(), -b)
lam = 0.01*norm(A'*b,Inf)
g = NormL1(lam)

x0 = zeros(n)
it, x_star, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, verbose = 0, tol = 1e-8)

# Nonfast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 50
println(sol)

# Fast/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=A, g=g, tol=1e-6, adaptive=true, fast=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 50
println(sol)

# ZeroFPR/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15
println(sol)

# PANOC/Adaptive

x0 = zeros(n)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=A, g=g, tol=1e-6, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)


text:
# a naive implementation of a terminal color code printer with a Pluto/HTML integration
# there is probably a nice Julia package for that, TODO: use it or at least refactor this piece of code

const PlainOrHtmlMIME = Union{MIME"text/plain",MIME"text/html"}

# It should be more tied to the type of diff to support background color and so-on
@enum Colors Green Red Reset

####
# Terminal colors
# For the diff printer interface implement with_color(f, io, m, col)
####

const color_codes =
    Dict{Colors,String}(Green => "\e[32m", Red => "\e[31m", Reset => "\e[39m")

function with_color(f::Function, io::IO, ::MIME"text/plain", col::Colors)
    write(io, color_codes[col])
    res = f()
    write(io, color_codes[Reset])
    res
end

####
# Pluto colors
# TODO: escaping of the user inputs using HyperTextLiteral?
####

# Colors are shade 400 of red/green from tailwindcss's default color theme
const hex_color_codes = Dict{Colors,String}(
    Green => "#34D399",
    Red => "#F87171",
    # no reset
)

function with_color(f::Function, io::IO, ::MIME"text/html", col::Colors)
    write(io, "")
    res = f()
    write(io, "")
    res
end

####

function color(io::IO, m::PlainOrHtmlMIME, col::Colors, text)
    with_color(() -> write(io, text), io, m, col)
end

green(io::IO, m::PlainOrHtmlMIME, text) = color(io, m, Green, text)
red(io::IO, m::PlainOrHtmlMIME, text) = color(io, m, Red, text)

with_green(f::Function, io::IO, m::PlainOrHtmlMIME) = with_color(f, io, m, Green)
with_red(f::Function, io::IO, m::PlainOrHtmlMIME) = with_color(f, io, m, Red)

Base.show(io::IO, diff::T) where {T<:AbstractDiff} = Base.show(io, MIME"text/plain"(), diff)
Base.show(io::IO, m::MIME"text/plain", diff::T) where {T<:AbstractDiff} = disp(io, m, diff)
Base.show(io::IO, m::MIME"text/html", diff::T) where {T<:AbstractDiff} = disp(io, m, diff)

disp(io::IO, m::PlainOrHtmlMIME, ins::Insertion) = green(io, m, ins.text)
disp(io::IO, m::PlainOrHtmlMIME, del::Deletion) = red(io, m, del.text)
disp(io::IO, ::PlainOrHtmlMIME, eq::Equality) = write(io, eq.text)

abstract type AbstractDiffViewer end

# how am i supposed to do these things ?
Base.show(io::IO, adv::T) where {T<:AbstractDiffViewer} =
    Base.show(io, MIME"text/plain"(), adv)
Base.show(io::IO, m::MIME"text/plain", adv::T) where {T<:AbstractDiffViewer} =
    disp(io, m, adv)
function Base.show(io::IO, m::MIME"text/html", adv::T) where {T<:AbstractDiffViewer}
    write(io, "
")
    disp(io, m, adv)
    write(io, "
")
end

struct SimpleDiffViewer <: AbstractDiffViewer
    diffs::Vector{<:AbstractDiff}
end

function disp(io::IO, m::PlainOrHtmlMIME, dv::SimpleDiffViewer)
    for diff in dv.diffs
        show(io, m, diff)
    end
end

newline(io::IO, ::MIME"text/plain") = write(io, '\n')
newline(io::IO, ::MIME"text/html") = write(io, "
")
tabulate(io::IO, ::MIME"text/plain") = write(io, '\t')
tabulate(io::IO, ::MIME"text/html") = write(io, "")

"""
Show both initial inputs on top of each other

```julia-repl
julia> PlotuDiff.stringdiff(\"\"\"
       This is the text before
       \"\"\", \"\"\"
       This is
       the text
       after
       \"\"\", viewer = PlotuDiff.SimpleLineDiffViewer)
        - This is the text before
        -
        + This is
        + the text
        + after
        +
```
"""
struct SimpleLineDiffViewer <: AbstractDiffViewer
    diffs::Vector{<:AbstractDiff}
end

function _split_on_newline(io::IO, m::MIME, s::AbstractString, startline)
    parts = split(s, '\n'; keepempty = true)
    for part in parts
        write(io, part)

        newline(io, m)
        startline()
    end
end
function _split_on_newline(io::IO, m::MIME, c::Char, startline)
    if c == '\n'
        newline(io, m)
        startline()
    else
        error("_split_on_newline called with $c")
    end
end

function disp(io::IO, m::PlainOrHtmlMIME, ldv::SimpleLineDiffViewer)
    with_red(io, m) do
        function startline()
            tabulate(io, m)
            write(io, "- ")
        end
        startline()

        for diff in ldv.diffs
            diff isa Insertion && continue
            if '\n' ∈ diff.text
                _split_on_newline(io, m, diff.text, startline)
                write(io, diff.text)
            end
        end
    end
    newline(io, m)

    with_green(io, m) do
        function startline()
            tabulate(io, m)
            write(io, "+ ")
        end
        startline()

        for diff in ldv.diffs
            diff isa Deletion && continue
            if '\n' ∈ diff.text
                _split_on_newline(io, m, diff.text, startline)
            else
                write(io, diff.text)
            end
        end
    end
end

function split_diff(diff::T, sep = '\n') where {T<:AbstractDiff}
    out = AbstractDiff[]
    for (i, part) in enumerate(split(diff.text, sep, keepempty = true))
        if i != 1
            push!(out, T(string(sep)))
        end

        if !isempty(part)
            push!(out, T(part))
        end
    end
    out
end

"""
    split_newlines(diffs::Vector{<:AbstractDiff})

Splits the diff on new lines. Expects the diffs to have already been merged.

```julia-repl
julia> split_newlines([Insertion("new\\nlines")])
[Insertion("new"), Insertion('\\n'), Insertion("\\lines")]
```
"""
function split_newlines(diffs::Vector{<:AbstractDiff})
    isempty(diffs) && return diffs

    out = AbstractDiff[]

    for diff in diffs
        for subdiff in split_diff(diff, '\n')
            push!(out, subdiff)
        end
    end

    out
end

"""
    LineViewer(diffs::Vector{<:AbstractDiff})

This viewer show the differences between texts line by line starting with
deletions.
"""
struct LineViewer <: AbstractDiffViewer
    diffs::Vector{<:AbstractDiff}
end

function disp(io::IO, m::PlainOrHtmlMIME, lv::LineViewer)
    diggs = merge_diffs(lv.diffs) .|> to_string |> split_newlines

    start_of_line = i = 1

    while i < length(diggs)
        start_of_line = i

        is_insertion = false
        is_deletion = false

        while diggs[i].text != "\n" && i < length(diggs)
            is_insertion |= diggs[i] isa Insertion
            is_deletion |= diggs[i] isa Deletion

            i += 1
        end

        if is_deletion
            diffs = filter(d -> !(d isa Insertion), diggs[start_of_line:(i - 1)])
            for diff in diffs
                show(io, m, diff)
            end
            newline(io, m)
        end

        if is_insertion
            diffs = filter(d -> !(d isa Deletion), diggs[start_of_line:(i - 1)])
            for diff in diffs
                show(io, m, diff)
            end
            newline(io, m)
        end

        if !is_insertion && !is_deletion
            value =
                filter(d -> !(d isa Deletion), diggs[start_of_line:(i - 1)]) .|>
                text |>
                join
            show(io, m, Equality(value))
            newline(io, m)
        end

        i += 1
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using FEMBase

"""
    craig_bampton(K, M, r, l, n)

Reduce the stiffness and mass matrices with the Craig-Bampton method.
K = original stiffness matrix, M = original mass matrix,
r = retained DOF:s, l = internal DOF:s, n = the number of modes to keep.

"""
function craig_bampton(K, M, r, l, n)
    Krr = K[r,r]; Krl = K[r,l]
    Klr = K[l,r]; Kll = K[l,l]
    Mrr = M[r,r]; Mrl = M[r,l]
    Mlr = M[l,r]; Mll = M[l,l]
    X = zeros(size(Kll,1), n)
    if issparse(K) && issparse(M)
        SparseArrays.dropzeros!(Kll)
        nz = get_nonzero_rows(Kll)
        w, X[nz,:] = eigs(Kll[nz,nz], Mll[nz,nz]; nev=n, which=:SM)
        kll = (Kll + Kll')/2
    else
        w2 = eigvals(Kll,Mll)
        w = w2[1:n]
        X1 = eigvecs(Kll,Mll)
        X = X1[:,1:n]
    end
    Kmm = X'*Kll*X
    P = Krl*X
    Q = Mrl*X
    invw = diagm(inv.(w))
    Kbb = Krr - P*invw*P'
    Kmb = zeros(n, length(r))
    Kbm = zeros(length(r), n)
    MU = X'*Mll*X
    Mbb = Mrr - Q*invw*P' - P*invw*Q' + P*invw*MU*invw*P'
    Mmb = Q' - MU*invw*P'
    Mbm = Q - P*invw*MU
    Mmm = MU
    # Final system:
    Mred = [Mbb Mbm; Mmb Mmm]
    Kred = [Kbb Kbm; Kmb Kmm]
    Z = 10e-10
    Mred[abs.(Mred) .< Z] = 0.0
    Kred[abs.(Kred) .< Z] = 0.0
    return Mred, Kred
end

"""
    craig_bampton(K, M, c, r, l, n)

Reduce the stiffness, damping and mass matrices with the Craig-Bampton method.
K = original stiffness matrix, M = original mass matrix, c = damping ratio,
r = retained DOF:s, l = internal DOF:s, n = the number of modes to keep.

"""
function craig_bampton(K, M, c, r, l, n)
    Krr = K[r,r]; Krl = K[r,l]
    Klr = K[l,r]; Kll = K[l,l]
    Mrr = M[r,r]; Mrl = M[r,l]
    Mlr = M[l,r]; Mll = M[l,l]
    if issparse(K) && issparse(M)
        w2, X1 = eigs(Kll, Mll)
    else
        w2 = eigvals(Kll,Mll)
        X1 = eigvecs(Kll,Mll)
    end
    X = X1[:,1:n]
    V = X1'*Kll*X1; Z = 10e-6
    Kmm = X'*Kll*X
    b = -X*inv(Kmm)*X'*Klr
    B = -X1*inv(V)*X1'*Klr
    Kbb = Krr + Krl*B
    Kbm = (Krl + b'*Kll)*X
    Kmb = X'*(Klr + Kll*b)
    Mbb = Mrr + Mrl*B + B'*Mlr + B'*Mll*B
    Mmm = X'*Mll*X
    Mbm = (Mrl + b'*Mll)*X
    Mmb = X'*(Mlr + Mll*b)
    Cbb = zeros(Kbb)
    Cbm = zeros(Kbm)
    Cmb = zeros(Kmb)
    Kmm[abs.(Kmm) .< Z] = 0.0
    Cmm = 2*c*sqrt.(Kmm)
    Mred = [Mbb Mbm; Mmb Mmm]
    Mred[abs.(Mred) .< Z] = 0.0
    Kred = [Kbb Kbm; Kmb Kmm]
    Kred[abs.(Kred) .< Z] = 0.0
    Cred = [Cbb Cbm; Cmb Cmm]
    Cred[abs.(Cred) .< Z] = 0.0
    return Mred, Cred, Kred
end


text:
using VisClaw
using Test
using Printf

# make a list
exdir = joinpath(dirname(pathof(VisClaw)), "../example")
filelist = readdir(exdir)
filter!(x->occursin(".jl", x), filelist)
#map(s->filter!(x->!occursin(s, x), filelist), ["check", "fgmax"])
map(s->filter!(x->!occursin(s, x), filelist), ["check"])
#filter!(x->occursin("current_Plots_ike", x), filelist)

# number
nf = length(filelist)
println(@sprintf("%d", nf)*" files are going to be tested...")

using GR: GR
GR.inline("png")

@testset "VisClaw.jl" begin
    # Write your own tests here.
    # for loop
    for f in filelist
        println(f)
        @test !isa(try include(joinpath(exdir,f)) catch ex ex end, Exception)
        #=
        try
            @test_nowarn include(joinpath(exdir,f))
        catch e
            try
                include(joinpath(exdir,f))
            catch e
                println("Failed "*f)

                bt = backtrace()
                msg = sprint(showerror, e, bt)
                println(msg)
                println()
                continue
            end
            continue
        end
        =#
    end
end


text:
#
# Use this in an active julia kernel:
#
# include("prepareresults.jl")
# df = prepareresults()
#
# This way, a nicely aggregated DataFrame is returned.
#
using DataFrames

function prepareresults(folder::ASCIIString="./seven_vs_all")

    # concatenate DataFrames
    df = DataFrame()
    for seed in readdir(folder)
        if isdir(folder*"/"*seed)
            accfile = folder*"/"*seed*"/acc.csv"
            if isfile(accfile)
                println("Reading $accfile...")
                df = vcat(df, readtable(accfile))
            end
        end
    end

    # only consider LogReg and Boosting (ensemble vote)
    df = df[convert(Array{Bool,1}, [ c == "LogReg" || c == "Boosting B1-3" for c in df[:classifier]]),:]

    # rename columns so it fits on screen
    rename!(df, [:batchsize, :classifier, :trainsize, :trainacc, :testacc, :iterrate],
                [:batch,     :learner,    :numex,     :train,    :test,    :rate])

    # shorten classifier name
    df[:learner] = [ (learner == "LogReg") ? "LogReg" : "Boosting" for learner in df[:learner] ]

    # print some info
    meanex = mean(df[df[:learner] .== "Boosting",:][:numex])
    stdex  = stdm(df[df[:learner] .== "Boosting",:][:numex], meanex)
    println("\nBoosting used $meanex examples on average (std: $stdex)\n")

    # compute std and mean
    df = aggregate(df, [:gdopt, :ls, :batch, :data, :learner], [mean, std])

    # select and reorder columns
    df = df[:, [:learner, :gdopt, :batch,
                :test_mean, :test_std, :time_mean, :time_std,
                :rate_mean, :rate_std, :train_mean, :train_std ]]

    # sort
    sort!(df, rev=[true, false, false, false], cols=[:test_mean, :test_std, :time_mean, :time_std])

    return df

end

function latex(df::DataFrame)

    _format(v::Real)           = @sprintf "%.3f" v
    _format(v::Int)            = @sprintf "%d"   v
    _format(v::AbstractString) = replace(v, "_", "")
    _format(v::Symbol)         = replace(string(v), "_", "")

    println("\\documentclass[margin={10pt 10pt 10pt 10pt}]{standalone}")
    println("\\usepackage[utf8]{inputenc}")
    println("\\usepackage{color}")
    println("\\usepackage{amsmath}")
    println("\\usepackage{booktabs}\n")
    println("\\begin{document}\n")

    println("\\begin{tabular}{",
            apply(string, map(n->"l", names(df))), "}")
    println("    ", join(map(n->_format(n), names(df)), " & "), " \\\\")
    println("    \\midrule")
    for r = eachrow(df)
        println("    ", join(map(n->_format(r[n]), names(r)), " & "), " \\\\")
    end
    println("    \\bottomrule")
    println("\\end{tabular}")

    println("\n\\end{document}")

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import FMIImport: fmi3StatusError

myFMU = fmi3Load("BouncingBall", "ModelicaReferenceFMUs", "0.0.14")

### CASE A: Print log ###
inst = fmi3InstantiateCoSimulation!(myFMU; loggingOn=true)
@test inst != 0

open(joinpath(pwd(), "stdout"), "w") do out
    open(joinpath(pwd(), "stderr"), "w") do err
        redirect_stdout(out) do
            redirect_stderr(err) do
                @test fmi3ExitInitializationMode(inst) == fmi3StatusError
            end
        end
    end
end
# ToDo: this test is wrong / not working (capture doesn't work for color output)
#output = read(joinpath(pwd(), "stdout"), String)
#@test output == ""
#output = read(joinpath(pwd(), "stderr"), String)
#@test output == ""

### CASE B: Print log, but catch infos ###

inst = fmi3InstantiateCoSimulation!(myFMU; loggingOn=true, logStatusError=false)
@test inst != 0

open(joinpath(pwd(), "stdout"), "w") do out
    open(joinpath(pwd(), "stderr"), "w") do err
        redirect_stdout(out) do
            redirect_stderr(err) do
                @test fmi3ExitInitializationMode(inst) == fmi3StatusError
            end
        end
    end
end
output = read(joinpath(pwd(), "stdout"), String)
@test output == ""

if VERSION >= v"1.7.0"
    output = read(joinpath(pwd(), "stderr"), String)
    @test startswith(output, "┌ Warning: fmi3ExitInitializationMode(...): Needs to be called in state `fmi3InstanceStateInitializationMode`.\n")
end

### CASE C: Disable Log ###

inst = fmi3InstantiateCoSimulation!(myFMU; loggingOn=false)
@test inst != 0

open(joinpath(pwd(), "stdout"), "w") do out
    open(joinpath(pwd(), "stderr"), "w") do err
        redirect_stdout(out) do
            redirect_stderr(err) do
                @test fmi3ExitInitializationMode(inst) == fmi3StatusError
            end
        end
    end
end
# ToDo: this test is wrong / not working (capture doesn't work for color output)
#output = read(joinpath(pwd(), "stdout"), String)
#@test output == ""
#output = read(joinpath(pwd(), "stderr"), String)
#@test output == ""

# cleanup
fmi3Unload(myFMU)


text:
abstract type GraphicsArea end
abstract type Canvas <: GraphicsArea end

origin(c::Canvas) = (origin_x(c), origin_y(c))
Base.size(c::Canvas) = (width(c), height(c))
pixel_size(c::Canvas) = (pixel_width(c), pixel_height(c))

function Base.print(io::IO, c::GraphicsArea)
    for row in 1:nrows(c)
        printrow(io, c, row)
        print(io, "\n")
    end
end

function Base.show(io::IO, c::GraphicsArea)
    b = border_dashed
    border_length = ncols(c)
    print_border_top(io, "", border_length, :solid)
    print(io, "\n")
    for row in 1:nrows(c)
        printstyled(io, b[:l]; color = :white)
        printrow(io, c, row)
        printstyled(io, b[:r], "\n";color = :white)
    end
    print_border_bottom(io, "", border_length, :solid)
    print(io, "\n")
end

function pixel!(c::Canvas, pixel_x::Int, pixel_y::Int; color::Symbol = :white)
    pixel!(c, pixel_x, pixel_y, color)
end

function points!(c::Canvas, plot_x::Number, plot_y::Number, color::Symbol)
    origin_x(c) <= plot_x < origin_x(c) + width(c) || return c
    origin_y(c) <= plot_y < origin_y(c) + height(c) || return c
    plot_offset_x = plot_x - origin_x(c)
    pixel_x = plot_offset_x / width(c) * pixel_width(c)
    plot_offset_y = plot_y - origin_y(c)
    pixel_y = pixel_height(c) - plot_offset_y / height(c) * pixel_height(c)
    pixel!(c, floor(Int, pixel_x), floor(Int, pixel_y), color)
end

function points!(c::Canvas, plot_x::Number, plot_y::Number; color::Symbol = :white)
    points!(c, plot_x, plot_y, color)
end

function points!(c::Canvas, X::AbstractVector, Y::AbstractVector, color::Symbol)
    length(X) == length(Y) || throw(DimensionMismatch("X and Y must be the same length"))
    for i in 1:length(X)
        points!(c, X[i], Y[i], color)
    end
    c
end

function points!(c::Canvas, X::AbstractVector, Y::AbstractVector; color::Symbol = :white)
    points!(c, X, Y, color)
end

# Implementation of the digital differential analyser (DDA)
function lines!(c::Canvas, x1::Number, y1::Number, x2::Number, y2::Number, color::Symbol)
    if (x1 < origin_x(c) && x2 < origin_x(c)) ||
        (x1 > origin_x(c) + width(c) && x2 > origin_x(c) + width(c))
        return c
    end
    if (y1 < origin_y(c) && y2 < origin_y(c)) ||
        (y1 > origin_y(c) + height(c) && y2 > origin_y(c) + height(c))
        return c
    end
    toff = x1 - origin_x(c)
    px1 = toff / width(c) * pixel_width(c)
    toff = x2 - origin_x(c)
    px2 = toff / width(c) * pixel_width(c)
    toff = y1 - origin_y(c)
    py1 = pixel_height(c) - toff / height(c) * pixel_height(c)
    toff = y2 - origin_y(c)
    py2 = pixel_height(c) - toff / height(c) * pixel_height(c)
    dx = px2 - px1
    dy = py2 - py1
    nsteps = abs(dx) > abs(dy) ? abs(dx) : abs(dy)
    inc_x = dx / nsteps
    inc_y = dy / nsteps
    cur_x = px1
    cur_y = py1
    fpw = Float64(pixel_width(c))
    fph = Float64(pixel_height(c))
    pixel!(c, floor(Int, cur_x), floor(Int, cur_y), color)
    for i = 1:nsteps
        cur_x += inc_x
        cur_y += inc_y
        pixel!(c, floor(Int, cur_x), floor(Int, cur_y), color)
    end
    c
end

function lines!(c::Canvas, x1::Number, y1::Number, x2::Number, y2::Number; color::Symbol = :white)
    lines!(c, x1, y1, x2, y2, color)
end

function lines!(c::Canvas, X::AbstractVector, Y::AbstractVector, color::Symbol)
    length(X) == length(Y) || throw(DimensionMismatch("X and Y must be the same length"))
    for i in 1:(length(X)-1)
        lines!(c, X[i], Y[i], X[i+1], Y[i+1], color)
    end
    c
end

function lines!(c::Canvas, X::AbstractVector, Y::AbstractVector; color::Symbol = :white)
    lines!(c, X, Y, color)
end


text:
#############################################################################
# ProxAL/ExaTron Example File
# This example runs ProxAL with ExaTron as a backend and outputs a profile
# file. PProf has to be installed in the global environment.
#############################################################################

using ProxAL
using DelimitedFiles, Printf
using LinearAlgebra, JuMP, Ipopt
using CatViews
using CUDA
using MPI
using Logging
using LazyArtifacts

MPI.Init()

case = "case9"
demandfiles = "$(case)_oneweek_168"
# Large 2000 bus case
# case = "case_ACTIVSg2000"
# demandfiles = "$(case)"

# choose one of the following (K*T subproblems in each case)
if length(ARGS) == 0
    (T, K) = (2, 1)
elseif length(ARGS) == 4
    case = ARGS[1]
    demandfiles = ARGS[2]
    T = parse(Int, ARGS[3])
    K = parse(Int, ARGS[4])
else
    println("Usage: [mpiexec -n nprocs] julia --project examples/exatron.jl [case demandfiles T K]")
    println("")
    println("       (case,demandfiles,T,K) defaults to (case9,case9_oneweek_168,2,1)")
    exit()
end

# choose backend
# backend = ProxAL.JuMPBackend()
# With ExaTronBackend(), CUDADevice will used
backend = ProxAL.ExaTronBackend()

# Load case
const DATA_DIR = joinpath(artifact"ExaData", "ExaData")
case_file = joinpath(DATA_DIR, "$(case).m")
load_file = joinpath(DATA_DIR, "mp_demand", demandfiles)

# Model/formulation settings
modelinfo = ModelInfo()
modelinfo.num_time_periods = T
modelinfo.load_scale = 1.0
modelinfo.ramp_scale = 0.2
modelinfo.corr_scale = 0.5
modelinfo.allow_obj_gencost = true
modelinfo.allow_constr_infeas = true
modelinfo.weight_constr_infeas = 1e8
modelinfo.time_link_constr_type = :penalty
modelinfo.ctgs_link_constr_type = :corrective_penalty
modelinfo.allow_line_limits = false
modelinfo.case_name = case
modelinfo.num_ctgs = K

# Algorithm settings
algparams = AlgParams()
algparams.verbose = 1
algparams.tol = 1e-3
algparams.decompCtgs = (K > 0)
algparams.iterlim = 100
if isa(backend, ProxAL.ExaTronBackend)
    algparams.device = ProxAL.CUDADevice
end
algparams.optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0) #,  "tol" => 1e-1*algparams.tol)
algparams.tron_rho_pq=5*1e4
algparams.tron_rho_pa=5*1e5
algparams.tron_outer_iterlim=30
algparams.tron_inner_iterlim=2000
algparams.tron_scale=1e-5
algparams.mode = :coldstart
algparams.init_opf = false
tron_outer_eps = 1e-6

ranks = MPI.Comm_size(MPI.COMM_WORLD)
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
   println("ProxAL/ExaTron $ranks ranks, $T periods, $K contingencies")
end
cur_logger = global_logger(NullLogger())
elapsed_t = @elapsed begin
  redirect_stdout(devnull) do
    global nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, backend)
  end
end
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    global_logger(cur_logger)
    println("Creating problem: $elapsed_t")
    println("Benchmark Start")
    np = MPI.Comm_size(MPI.COMM_WORLD)
    elapsed_t = @elapsed begin
      info = ProxAL.optimize!(nlp)
    end
    println("AugLag iterations: $(info.iter) with $np ranks in $elapsed_t seconds")
else
  info = ProxAL.optimize!(nlp)
end

MPI.Finalize()

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
module TestTableTraits

using Test, DataFrames, IteratorInterfaceExtensions, TableTraits, DataValues

struct ColumnSource
end

TableTraits.supports_get_columns_copy_using_missing(::ColumnSource) = true
TableTraits.isiterabletable(::ColumnSource) = true

function TableTraits.get_columns_copy_using_missing(x::ColumnSource)
    return (a=[1 ,2, 3], b=[4.0, 5.0, 6.0], c=["A", "B", "C"])
end

@testset "TableTraits" begin
    df = DataFrame(a=[1, 2 ,3], b=[1.0, missing, 3.0])
    @test IteratorInterfaceExtensions.isiterable(df)
    @test TableTraits.isiterabletable(df)
    @test collect(IteratorInterfaceExtensions.getiterator(df)) ==
        [(a=1, b=DataValue(1.0)), (a=2, b=DataValue{Float64}()), (a=3, b=DataValue(3.0))]
    sdf = view(df, 1:2, :)
    @test IteratorInterfaceExtensions.isiterable(sdf)
    @test TableTraits.isiterabletable(sdf)
    @test collect(IteratorInterfaceExtensions.getiterator(sdf)) ==
        [(a=1, b=DataValue(1.0)), (a=2, b=DataValue{Float64}())]

    df = DataFrame(ColumnSource())

    @test size(df)==(3, 3)
    @test df.a==[1, 2, 3]
    @test df.b==[4.0, 5.0, 6.0]
    @test df.c==["A", "B", "C"]
end

end # module


text:
# Autogenerated wrapper script for Giflib_jll for aarch64-apple-darwin
export libgif

JLLWrappers.@generate_wrapper_header("Giflib")
JLLWrappers.@declare_library_product(libgif, "@rpath/libgif.7.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgif,
        "lib/libgif.7.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()