text:
import Base.show
import DataStructures
import Random

struct Node{InputType, MetricReturnType}
    index::Int
    data::InputType
    radius::MetricReturnType
    min_dist::MetricReturnType
    max_dist::MetricReturnType
    left_child::Union{Node{InputType, MetricReturnType}, Nothing}
    right_child::Union{Node{InputType, MetricReturnType}, Nothing}
end

function Base.show(io::IO, n::Node)
    print("Node $(typeof(n).parameters[1]): $(n.data), index $(n.index) radius $(n.radius), depth $(_depth(n))")
end

function _depth(n)
    if n === nothing
        0
    else
        1 + max(_depth(n.left_child), _depth(n.right_child))
    end
end

"""
    VPTree(data::Vector{InputType}, metric::Function)

Construct Vantage Point Tree with a vector of `data` and given metric function `metric`.

## Example:

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
metric = (a, b) -> evaluate(Levenshtein(),a,b)
vptree = VPTree(data, metric)
query = "blau"
radius = 2
data[find(vptree, query, radius)]
# 2-element Array{String,1}:
#  "bla"
#  "blub"
n_neighbors = 3
data[find_nearest(vptree, query, n_neighbors)]
# 3-element Array{String,1}:
#  "baube"
#  "blub"
#  "bla"
```
"""
struct VPTree{InputType, MetricReturnType}
    data::Vector{InputType}
    metric::Function
    root::Node{InputType, MetricReturnType}
    MetricReturnType::DataType
    function VPTree(data::Vector{InputType}, metric::Function) where {InputType}
        @assert length(data) > 0 "Input data must contain at least one point."
        MetricReturnType = typeof(metric(data[1], data[1]))
        indexed_data = Random.shuffle!(collect(enumerate(data)))
        root = _construct_tree_rec!(indexed_data, metric, MetricReturnType)
        new{InputType, MetricReturnType}(data, metric, root, MetricReturnType)
    end
end

function VPTree(data::Vector{T}, metric::Function, MetricReturnType) where T
    VPTree(data, metric)
end

@deprecate VPTree(data::Vector, metric::Function, MetricReturnType::DataType) VPTree(data::Vector, metric::Function)

function _construct_tree_rec!(data::AbstractVector{Tuple{Int, InputType}}, metric, MetricReturnType) where InputType
    if isempty(data)
        return nothing
    end
    n_data = length(data)
    if n_data == 1
        return Node(data[1][1], data[1][2], zero(MetricReturnType), zero(MetricReturnType), zero(MetricReturnType), nothing, nothing)
    end
    vantage_point = data[1]
    rest = view(data, 2:length(data))
    distances = [metric(d[2], vantage_point[2]) for d in rest]
    i_middle = div(n_data - 1, 2) + 1
    distance_order = sortperm(distances, alg=PartialQuickSort(i_middle))
    permute!(rest, distance_order)

    left_rest = view(rest, 1:i_middle)
    left_node = _construct_tree_rec!(left_rest, metric, MetricReturnType)

    right_rest = view(rest, i_middle + 1:length(rest))
    right_node = _construct_tree_rec!(right_rest, metric, MetricReturnType)

    min_dist, max_dist = extrema(distances)
    radius = distances[distance_order[i_middle]]

    Node{InputType, MetricReturnType}(vantage_point[1], vantage_point[2], radius,  min_dist, max_dist, left_node, right_node)
end

"""
    find(vptree::VPTree{InputType, MetricReturnType}, query::InputType, radius::MetricReturnType)::Vector{Int}

Find all items in `vptree` within `radius` with respect to the metric defined in the VPTree.
Returns Indices into VPTree.data.

## Example

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
metric = (a, b) -> evaluate(Levenshtein(),a,b)
vptree = VPTree(data, metric)
query = "blau"
radius = 2
data[find(vptree, query, radius)]
# 2-element Array{String,1}:
#  "bla"
#  "blub"
```
"""
function find(vptree::VPTree{InputType, MetricReturnType}, query::InputType, radius::MetricReturnType)::Vector{Int} where {InputType, MetricReturnType}
    results = Vector{Int}()
    _find(vptree.root, query, radius, results, vptree.metric)
    results
end

function _find(vantage_point, query, radius, results, metric)
    distance = metric(vantage_point.data, query)
    if distance <= radius
        push!(results, vantage_point.index)
    end
    if distance - radius <= vantage_point.radius && vantage_point.left_child != nothing && distance + radius >= vantage_point.min_dist
        _find(vantage_point.left_child, query, radius, results, metric)
    end
    if distance + radius > vantage_point.radius && vantage_point.right_child != nothing && distance - radius <= vantage_point.max_dist
        _find(vantage_point.right_child, query, radius, results, metric)
    end
end

"""
    find_nearest(vptree::VPTree{InputType, MetricReturnType}, query::InputType, n_neighbors::Int)::Vector{Int}

Find `n_neighbors` items in `vptree` closest to `query` with respect to the metric defined in the VPTree.
Returns Indices into VPTree.data.

## Example:

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
metric = (a, b) -> evaluate(Levenshtein(),a,b)
vptree = VPTree(data, metric)
query = "blau"
n_neighbors = 3
data[find_nearest(vptree, query, n_neighbors)]
# 3-element Array{String,1}:
#  "baube"
#  "blub"
#  "bla"
```
"""
function find_nearest(vptree::VPTree{InputType, MetricReturnType}, query::InputType, n_neighbors::Int)::Vector{Int} where {InputType, MetricReturnType}
    @assert n_neighbors > 0 "Can't search for fewer than 1 neighbors"
    candidates = DataStructures.BinaryMaxHeap{Tuple{MetricReturnType, Int}}()
    _find_nearest(vptree.root, query, n_neighbors, candidates, vptree.metric)
    [t[2] for t in candidates.valtree]
end

function _find_nearest(vantage_point, query, n_neighbors, candidates, metric)
    distance = metric(vantage_point.data, query)
    radius = length(candidates) < n_neighbors ? typemax(typeof(vantage_point.radius)) : DataStructures.top(candidates)[1]
    if distance < radius
        push!(candidates, (distance, vantage_point.index))
        if length(candidates) > n_neighbors
            pop!(candidates)
        end
    end
    # Switch radius to one side to prevent overflow.
    if distance < vantage_point.radius
        # Switch radius to one side to prevent overflow.
        if distance - vantage_point.radius <= radius && vantage_point.left_child != nothing && distance - vantage_point.min_dist >= -radius
            _find_nearest(vantage_point.left_child, query, n_neighbors, candidates, metric)
        end
        if distance - vantage_point.radius > -radius && vantage_point.right_child != nothing && distance - vantage_point.max_dist <= radius
            _find_nearest(vantage_point.right_child, query, n_neighbors, candidates, metric)
        end
    else
        if distance - vantage_point.radius > -radius && vantage_point.right_child != nothing && distance - vantage_point.max_dist <= radius
            _find_nearest(vantage_point.right_child, query, n_neighbors, candidates, metric)
        end
        if distance - vantage_point.radius <= radius && vantage_point.left_child != nothing && distance - vantage_point.min_dist >= -radius
            _find_nearest(vantage_point.left_child, query, n_neighbors, candidates, metric)
        end
    end
end


text:
### Effective sample size (ESS)

actypes = (:bm, :imse, :ipse)

function ess(x::Vector{Float64}; vtype::Symbol=:imse, args...)
  @assert in(vtype, actypes) "Unknown ESS type $vtype"

  return length(x)*mcvar(x; vtype=:iid)/mcvar(x; vtype=vtype, args...)
end

ess(x::Matrix{Float64}, pars::Range=1:size(x, 2); vtype::Symbol=:imse, args...) =
  Float64[ess(x[:, pars[i]]; vtype=vtype, args...) for i = 1:length(pars)]

ess(x::Matrix{Float64}, par::Real; vtype::Symbol=:imse, args...) = ess(x, par:par; vtype=vtype, args...)

ess(c::MCChain, pars::Range=1:size(c.samples, 2); vtype::Symbol=:imse, args...) =
  ess(c.samples, pars; vtype=vtype, args...)

ess(c::MCChain, par::Real; vtype::Symbol=:imse, args...) = ess(c, par:par; vtype=vtype, args...)

ext: jl


text:
using JuMP
using Clp

m = Model(Clp.Optimizer)

@variable(m, x₁ ≥ 0)
@variable(m, x₂ ≥ 0)
@variable(m, x₃ ≥ 0)

@objective(m, Max, 2x₁  + 4x₂ + 3x₃)
@constraint(m, 3x₁  + 4x₂ + 2x₃ ≤ 60)
@constraint(m, 2x₁  + x₂ + 2x₃ ≤ 40)
@constraint(m, x₁  + 3x₂ + 2x₃ ≤ 80)

print(m)

optimize!(m)

status = termination_status(m)

println("Solution status: ", status)

println("Objective value: ", objective_value(m))
println("x₁ = ", value(x₁))
println("x₂ = ", value(x₂))
println("x₃ = ", value(x₃))


text:
module VAETools
    export get_data, Encoder, Decoder, reconstruct, model_loss

    using Flux
    using Flux.Data.MNIST
    using Flux.Data: DataLoader
    using Flux: onehotbatch
    using Flux.Losses: logitbinarycrossentropy

    function get_data(batch_size)
        imgs = MNIST.images()
        labels = MNIST.labels()

        xtrain = hcat(float.(reshape.(imgs, :))...)
        ytrain = onehotbatch(labels, 0:9)

        DataLoader((xtrain, ytrain), batchsize=batch_size, shuffle=true)
    end

    struct Encoder
        linear
        μ
        logσ
    end

    @functor Encoder

    Encoder(input_dim::Int, latent_dim::Int, hidden_dim::Int) =
    Encoder(
        Dense(input_dim, hidden_dim, tanh),     # linear
        Dense(hidden_dim, latent_dim),          # μ
        Dense(hidden_dim, latent_dim),          # logσ
    )

    function (encoder::Encoder)(x)
        h = encoder.linear(x)
        encoder.μ(h), encoder.logσ(h)
    end

    Decoder(input_dim::Int, latent_dim::Int, hidden_dim::Int) =
    Chain(
        Dense(latent_dim, hidden_dim, tanh),
        Dense(hidden_dim, input_dim)
    )

    function reconstruct(encoder, decoder, x, device)
        μ, logσ = encoder(x)
        z = μ + device(randn(Float32, size(logσ))) .* exp.(logσ)
        μ, logσ, decoder(z)
    end

    function model_loss(encoder, decoder, λ, x, device)
        μ, logσ, decoder_z = reconstuct(encoder, decoder, x, device)
        len = size(x)[end]
        # KL-divergence
        kl_q_p = 0.5f0 * sum(@. (exp(2f0 * logσ) + μ^2 -1f0 - 2f0 * logσ)) / len

        logp_x_z = -logitbinarycrossentropy(decoder_z, x, agg=sum) / len
        # regularization
        reg = λ * sum(x->sum(x.^2), Flux.params(decoder))

        -logp_x_z + kl_q_p + reg
    end

end

max_stars_repo_name: aniketjivani/variational-autoencoders

max_issues_repo_name: aniketjivani/variational-autoencoders

max_forks_repo_name: aniketjivani/variational-autoencoders


text:
@testset "Manopt.jl Error Measures" begin
    M = Sphere(2)
    N = PowerManifold(M, NestedPowerRepresentation(), 2)
    using Random: seed!
    seed!(42)
    d = Manifolds.uniform_distribution(M,[1.0,0.0,0.0])
    w = rand(d)
    x = rand(d)
    y = rand(d)
    z = rand(d)
    a = [w,x]
    b = [y,z]
    @test meanSquaredError(M,x,y) == distance(M,x,y)^2
    @test meanSquaredError(N,a,b) == 1/2*(distance(M,w,y)^2+distance(M,x,z)^2)
    @test meanAverageError(M,x,y) == distance(M,x,y)
    @test meanAverageError(N,a,b) == 1/2*sum( distance.(Ref(M), a, b) )
end


text:
using Documenter, VarianceCovarianceEstimators

makedocs(
    # format = :html,
    sitename = "VarianceCovarianceEstimators.jl",
    pages = [
        "index.md",
        "GettingStarted.md",
        "ModelAPI.md",
        "Examples.md",
        "References.md"
    ]
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo = "github.com/JuliaEconometrics/VarianceCovarianceEstimators.jl.git",
    julia  = "nightly"
)


text:
module LLVMBootstrap_jll
using Pkg, Pkg.BinaryPlatforms, Pkg.Artifacts, Libdl
import Base: UUID

# We put these inter-JLL-package API values here so that they are always defined, even if there
# is no underlying wrapper held within this JLL package.
const PATH_list = String[]
const LIBPATH_list = String[]

# Load Artifacts.toml file
artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")

# Extract all platforms
artifacts = Pkg.Artifacts.load_artifacts_toml(artifacts_toml; pkg_uuid=UUID("db172ae1-73ad-59d9-bc09-d7dba05a9256"))
platforms = [Pkg.Artifacts.unpack_platform(e, "LLVMBootstrap", artifacts_toml) for e in artifacts["LLVMBootstrap"]]

# Filter platforms based on what wrappers we've generated on-disk
platforms = filter(p -> isfile(joinpath(@__DIR__, "wrappers", triplet(p) * ".jl")), platforms)

# From the available options, choose the best platform
best_platform = select_platform(Dict(p => triplet(p) for p in platforms))

# Silently fail if there's no binaries for this platform
if best_platform === nothing
    @debug("Unable to load LLVMBootstrap; unsupported platform $(triplet(platform_key_abi()))")
else
    # Load the appropriate wrapper
    include(joinpath(@__DIR__, "wrappers", "$(best_platform).jl"))
end

end  # module LLVMBootstrap_jll


text:
using Libdl

const depsfile = joinpath(@__DIR__, "deps.jl")

const libpath = get(ENV, "JULIA_HDF5_LIBRARY_PATH", nothing)

if libpath === nothing
    # By default, use HDF5_jll
    open(depsfile, "w") do io
        print(io,
              """
              # This file is automatically generated
              # Do not edit
              using HDF5_jll
              check_deps() = nothing
              """
             )
    end
else
    libhdf5 = find_library("libhdf5", [libpath])
    libhdf5_hl = find_library("libhdf5_hl", [libpath])

    isempty(libhdf5) && error("libhdf5 not found in $libpath")
    isempty(libhdf5_hl) && error("libhdf5_hl not found in $libpath")

    libhdf5_size = filesize(dlpath(libhdf5))

    open(depsfile, "w") do io
        println(io,
                """
                # This file is automatically generated
                # Do not edit

                function check_deps()
                    if libhdf5_size != filesize(Libdl.dlpath(libhdf5))
                        error("HDF5 library has changed, re-run Pkg.build(\\\"HDF5\\\")")
                    end
                end
                """
               )
        println(io, :(const libhdf5 = $libhdf5))
        println(io, :(const libhdf5_hl = $libhdf5_hl))
        println(io, :(const libhdf5_size = $libhdf5_size))
    end
end


text:
using ArgParse
using JuliaFormatter

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--overwrite"
            action = :store_true
        "--indent"
            arg_type = Int
            default = 4
        "--margin"
            arg_type = Int
            default = 92
        "--for_in"
            action = :store_true
        "--whitespace_typedefs"
            action = :store_true
        "--whitespace_ops_in_indices"
            action = :store_true
        "--remove_extra_newlines"
            action = :store_true
        "--path"
            required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    # println("Parsed args:")
    # for (arg,val) in parsed_args
    #     println("  $arg  =>  $val\t", typeof(val))
    # end

    for (arg, val) in parsed_args
        a = Symbol(arg)
        @eval $a = $val
    end

    format(
        path,
        overwrite                   = overwrite,
        verbose                     = false,
        indent                      = indent,
        margin                      = margin,
        always_for_in               = for_in,
        whitespace_typedefs         = whitespace_typedefs,
        whitespace_ops_in_indices   = whitespace_ops_in_indices,
        remove_extra_newlines       = remove_extra_newlines,
        style                       = DefaultStyle()
    )
end

main()


text:
using Random
using Flux

struct Param
    beta
    lat
    nd
    volume
    tau
    nstep
    dt
    ntraj
    nrun
    nprint
    seed
    randinit
end

function Param(;
        beta = 6.0, lat = [64, 64],
        tau = 2.0, nstep = 50, ntraj = 256, nrun = 4, nprint = 256,
        seed = 11*13, randinit = false)
    nd = length(lat)
    volume = prod(lat)
    dt = tau / nstep
    Param(beta,lat,nd,volume,tau,nstep,dt,ntraj,nrun,nprint,seed,randinit)
end

summary(param::Param) = """latsize = $(param.lat)
    volume = $(param.volume)
    beta = $(param.beta)
    trajs = $(param.ntraj)
    tau = $(param.tau)
    steps = $(param.nstep)
    seed = $(param.seed)
    """
function uniquestr(param::Param)
    lat = join([string(x) for x in param.lat], ".")
    "out_l$(lat)_b$(param.beta)_n$(param.ntraj)_t$(param.tau)_s$(param.nstep).out"
end

action(param, f) = (-param.beta)*sum(cos.(plaqphase(f)))
force(param, f) = gradient(x -> action(param, x), f)[1]

plaqphase(f) = f[1,:,:] - f[2,:,:] - circshift(f[1,:,:], [0,-1]) + circshift(f[2,:,:], [-1,0])

topocharge(f) = round(sum(regularize(plaqphase(f))) / 2pi)
function regularize(f)
    p2 = 2pi
    # f_ = f .- pi
    p2 .* ((f .- pi) ./ p2 .- fld.((f .- pi), p2) .- 0.5)
end

function leapfrog(param, x, p)
    dt = param.dt
    x_ = x + 0.5dt .* p
    p_ = p + (-dt) .* force(param, x_)
    for i = 1:param.nstep
        x_ += dt .* p_
        p_ += (-dt) .* force(param, x_)
    end
    x_ += 0.5dt .* p_
    (x_, p_)
end
function hmc(param, x)
    p = randn(eltype(x), size(x))
    act0 = action(param, x) + 0.5*sum(p .* p)
    x_, p = leapfrog(param, x, p)
    x_ = regularize(x_)
    act = action(param, x_) + 0.5*sum(p .* p)
    prob = rand()
    dH = act-act0
    exp_mdH = exp(-dH)
    acc = prob < exp_mdH
    newx = acc ? x_ : x
    (dH, exp_mdH, acc, newx)
end

put(s) = write(stdout, s)

function main()
    param = Param(
        beta = 6.0, # 4.0
        lat = [64, 64], # [8, 8]
        tau = 2, # 0.3
        nstep = 50, # 3
        ntraj = 256, # 2**16 # 2**10 # 2**15
        nprint = 256,
        seed = 1331)

    Random.seed!(param.seed)

    field = param.randinit ? randn(Float64, (param.nd, param.lat...)) : zeros(Float64, (param.nd, param.lat...))  # mu, x, t

    open(uniquestr(param), "w") do O
        params = summary(param)
        write(O, params)
        put(params)
        plaq_, topo_ = (action(param, field) / (-param.beta*param.volume), topocharge(field))
        status = "Initial configuration:  plaq: $(plaq_)  topo: $(topo_)\n"
        write(O, status)
        put(status)
        ts = Float64[]
        for n = 1:param.nrun
            t = -time()
            for i = 1:param.ntraj
                dH, exp_mdH, acc, field = hmc(param, field)
                plaq = action(param, field) / (-param.beta*param.volume)
                topo = topocharge(field)
                ifacc = acc ? "ACCEPT" : "REJECT"
                status = "Traj: $((n-1)*param.ntraj+i) $ifacc:  dH: $dH  exp(-dH): $exp_mdH  plaq: $plaq  topo: $topo\n"
                write(O, status)
                if i % div(param.ntraj, param.nprint) == 0
                    put(status)
                end
            end
            t += time()
            append!(ts, t)
        end
        print("Run times: $ts\n")
        print("Per trajectory: $([t/param.ntraj for t in ts])\n")
    end
end

main()


text:
using luteos
using SymPy

# let # limit scope

dim = 3

# Material
mat = Material( ν=0.33, E=1.0 )

## Set up problem
#   Get functions for exact solution
function ExactSol3D( Cstiff )

  dim = 3

  @syms x1 x2 x3

  u1 = 10 * (x2 - x2.^2) .* (x3 - x3.^2) .* sin.(pi*x1) .* (1 - x1) .* (1 - x2/2) .* (1 - x3/4)
  u2 =  2 * (x1 - x1.^2) .* (x3 - x3.^2) .* sin.(pi*x2) .* (1 - x2) .* (1 - x1/2) .* (1 - x3/4)
  u3 =  5 * (x1 - x1.^2) .* (x2 - x2.^2) .* sin.(pi*x3) .* (1 - x3) .* (1 - x1/2) .* (1 - x2/4)

  ∂u1∂x1 = diff( u1, x1 )
  ∂u1∂x2 = diff( u1, x2 )
  ∂u1∂x3 = diff( u1, x3 )
  ∂u2∂x1 = diff( u2, x1 )
  ∂u2∂x2 = diff( u2, x2 )
  ∂u2∂x3 = diff( u2, x3 )
  ∂u3∂x1 = diff( u3, x1 )
  ∂u3∂x2 = diff( u3, x2 )
  ∂u3∂x3 = diff( u3, x3 )

  ∇u = [ ∂u1∂x1 ∂u2∂x1 ∂u3∂x1; ∂u1∂x2 ∂u2∂x2 ∂u3∂x2; ∂u1∂x3 ∂u2∂x3 ∂u3∂x3 ]

  ϵ = 1/2 * ( ∇u + ∇u.' )

  σ = fill( 0*x1, dim, dim )

  for ii in 1:dim, jj in 1:dim, kk in 1:dim, ll in 1:dim
    σ[ii,jj] += (Cstiff[ii,jj,kk,ll] * ϵ[kk,ll])[1,1,1,1]
  end

  F  = fill( 0*x1, dim )

  ∇σ = fill( 0*x1, dim, dim, dim )

  for ii in 1:dim, jj in 1:dim
    ∇σ[ii,jj,1] = diff( σ[ii,jj], x1 )
    ∇σ[ii,jj,2] = diff( σ[ii,jj], x2 )
    ∇σ[ii,jj,3] = diff( σ[ii,jj], x3 )
  end

  for ii in 1:dim, jj in 1:dim
    F[ii] -= ∇σ[ii,jj,jj]
  end

  u1func_org = lambdify( u1 )
  u2func_org = lambdify( u2 )
  u3func_org = lambdify( u3 )

  F1func_org  = lambdify( F[1] )
  F2func_org  = lambdify( F[2] )
  F3func_org  = lambdify( F[3] )

  σ11func_org = lambdify( σ[1,1] )
  σ12func_org = lambdify( σ[1,2] )
  σ13func_org = lambdify( σ[1,3] )
  σ22func_org = lambdify( σ[2,2] )
  σ23func_org = lambdify( σ[2,3] )
  σ33func_org = lambdify( σ[3,3] )

  ϵ11func_org = lambdify( ϵ[1,1] )
  ϵ12func_org = lambdify( ϵ[1,2] )
  ϵ13func_org = lambdify( ϵ[1,3] )
  ϵ22func_org = lambdify( ϵ[2,2] )
  ϵ23func_org = lambdify( ϵ[2,3] )
  ϵ33func_org = lambdify( ϵ[3,3] )

  # make sure you can call it with matrices
  u1func = x -> u1func_org(x[:,1], x[:,2], x[:,3])
  u2func = x -> u2func_org(x[:,1], x[:,2], x[:,3])
  u3func = x -> u3func_org(x[:,1], x[:,2], x[:,3])

  F1func = x -> F1func_org(x[:,1], x[:,2], x[:,3])
  F2func = x -> F2func_org(x[:,1], x[:,2], x[:,3])
  F3func = x -> F3func_org(x[:,1], x[:,2], x[:,3])

  σ11func = x -> σ11func_org(x[:,1], x[:,2], x[:,3])
  σ12func = x -> σ12func_org(x[:,1], x[:,2], x[:,3])
  σ13func = x -> σ13func_org(x[:,1], x[:,2], x[:,3])
  σ22func = x -> σ22func_org(x[:,1], x[:,2], x[:,3])
  σ23func = x -> σ23func_org(x[:,1], x[:,2], x[:,3])
  σ33func = x -> σ33func_org(x[:,1], x[:,2], x[:,3])

  ϵ11func = x -> ϵ11func_org(x[:,1], x[:,2], x[:,3])
  ϵ12func = x -> ϵ12func_org(x[:,1], x[:,2], x[:,3])
  ϵ13func = x -> ϵ13func_org(x[:,1], x[:,2], x[:,3])
  ϵ22func = x -> ϵ22func_org(x[:,1], x[:,2], x[:,3])
  ϵ23func = x -> ϵ23func_org(x[:,1], x[:,2], x[:,3])
  ϵ33func = x -> ϵ33func_org(x[:,1], x[:,2], x[:,3])

  return( u1func, u2func, u3func,
          F1func, F2func, F3func,
          σ11func, σ12func, σ13func, σ22func, σ23func, σ33func,
          ϵ11func, ϵ12func, ϵ13func, ϵ22func, ϵ23func, ϵ33func )

end

( u1func, u2func, u3func,
  F1func, F2func, F3func,
  σ11func, σ12func, σ13func, σ22func, σ23func, σ33func,
  ϵ11func, ϵ12func, ϵ13func, ϵ22func, ϵ23func, ϵ33func ) = ExactSol3D( mat.Cstiff[dim] )

#   Setup boundary conditions
function funcB( p::Array{Float64} )
  return fill(0.0, size(p,1), 3)
end
bctype = [1,1,1,1,1,1] # All Dirichlet
#   Setup source function
source = (p) -> [F1func(p) F2func(p) F3func(p)]

P = 1 # Polynomial order
N = 3 # Grid size

@time mesh   = Mesh3D( "cube", P, N = N)
@time master = Master3D( P )

@time compJacob!( mesh, master )

prob = Problem( @sprintf("Reg %i %i", P, N), source, bctype, 0,
        [funcB, funcB, funcB, funcB, funcB, funcB] )

(uhathTri, uh, σh, ϵh) = hdgSolveElas( master, mesh, mat, prob )

# end # limit scope


text:
using AlmostBandedMatrices, BandedMatrices, FillArrays, Compat.Test, Compat.LinearAlgebra
using Compat
import BandedMatrices: _BandedMatrix

# Constructions:
@testset "construction" begin
    let A
end


text:
import Base
import AbstractAlgebra, Nemo
#import Oscar

include("utils.jl")
include("AbstractTypes.jl")

#------------------------------------------------------------------------------

function form_derivative(varname::String, order::Integer)
    return "$(varname)^($order)"
end

function decompose_derivative(dername::String)
    regex = r"^(.*)\^\((\d+)\)$"
    m = match(regex, dername)
    if m === nothing
        throw(DomainError("$dername is not a name of a derivative"))
    end
    return (String(m.captures[1]), parse(Int, m.captures[2]))
end

#------------------------------------------------------------------------------

mutable struct DiffPolyRing <: DiffRing
    base_ring::AbstractAlgebra.Ring
    poly_ring::MPolyRing
    max_ord::Integer
    varnames::Array{String, 1}
    derivation::Dict{Any, Any}

    function DiffPolyRing(R::AbstractAlgebra.Ring, varnames::Array{String, 1}, max_ord=20)
        all_varnames = [form_derivative(v, ord) for v in varnames for ord in 0:max_ord]
        poly_ring, _ = AbstractAlgebra.PolynomialRing(R, all_varnames)
        derivation = Dict()
        for v in varnames
            for ord in 0:(max_ord - 1)
                derivation[str_to_var(form_derivative(v, ord), poly_ring)] =
                    str_to_var(form_derivative(v, ord + 1), poly_ring)
            end
        end
        return new(R, poly_ring, max_ord, varnames, derivation)
    end
end

function DifferentialPolynomialRing(R::AbstractAlgebra.Ring, varnames::Array{String, 1}, max_ord=20)
    R = DiffPolyRing(R, varnames, max_ord)
    return R, Tuple(gens(R))
end

#------------------------------------------------------------------------------

mutable struct DiffPoly <: DiffRingElem
    parent::DiffPolyRing
    algdata::AbstractAlgebra.MPolyElem

    function DiffPoly(R::DiffPolyRing, alg_poly::AbstractAlgebra.MPolyElem)
        return new(R, parent_ring_change(alg_poly, R.poly_ring))
    end
end

function Base.parent(a::DiffPoly)
    return a.parent
end

function algdata(a::DiffPoly)
    return a.algdata
end

AbstractAlgebra.elem_type(::Type{DiffPolyRing}) = DiffRingElem

AbstractAlgebra.parent_type(::Type{DiffRingElem}) = DiffPolyRing

AbstractAlgebra.parent_type(::Type{DiffPoly}) = DiffPolyRing

function get_der(R::DiffPolyRing, vname::String, ord::Integer)
    if ord > R.max_ord
        throw(DomainError("Requested order $ord exceed the maximal order of the ring $(R.max_ord)"))
    end
    return DiffPoly(R, str_to_var(form_derivative(vname, ord), R.poly_ring))
end

#------------------------------------------------------------------------------

function AbstractAlgebra.gens(R::DiffPolyRing)
    return [DiffPoly(R, str_to_var(form_derivative(v, 0), R.poly_ring)) for v in R.varnames]
end

#------------------------------------------------------------------------------

mutable struct DiffRational <: DiffFieldElem
    parent::DiffPolyRing
    algdata::Generic.Frac{<: AbstractAlgebra.MPolyElem}

    function DiffRational(R::DiffPolyRing, alg_rational::Generic.Frac{<: AbstractAlgebra.MPolyElem})
        p, q = map(x -> parent_ring_change(x, R.poly_ring), [numerator(alg_rational), denominator(alg_rational)])
        return new(R, p // q)
    end
end

function Base.parent(a::DiffRational)
    return a.parent
end

function algdata(a::DiffRational)
    return a.algdata
end

function Base.numerator(a::DiffRational)
    return DiffPoly(parent(a), numerator(algdata(a)))
end

function Base.denominator(a::DiffRational)
    return DiffPoly(parent(a), denominator(algdata(a)))
end

AbstractAlgebra.parent_type(::Type{DiffRational}) = DiffPolyRing

#------------------------------------------------------------------------------

function (R::DiffPolyRing)(b)
    return DiffPoly(R, R.poly_ring(b))
end

function (R::DiffPolyRing)(b::Generic.Frac)
    return DiffRational(R, b)
end

function (R::DiffPolyRing)(b::DiffPoly)
    return DiffPoly(R, algdata(b))
end

function (R::DiffPolyRing)(b::DiffRational)
    return DiffRational(R, algdata(b))
end

function (R::DiffPolyRing)()
    return zero(R)
end

function Base.zero(R::DiffPolyRing)
    return R(0)
end

function Base.one(R::DiffPolyRing)
    return R(1)
end

#------------------------------------------------------------------------------

function Base.show(io::IO, R::DiffPolyRing)
    print(io, "Differential polynomial ring over $(R.base_ring) in " * join(R.varnames, ", "))
end

function Base.show(io::IO, p::DiffRingElem)
    show(io, algdata(p))
end

#------------------------------------------------------------------------------

function check_parent(a::DiffRingElem, b::DiffRingElem)
    if parent(a) != parent(b)
        throw(DomainError("$a and $b are from different rings: $(parent(a)) and $(parent(b)), respectively"))
    end
end

#------------------------------------------------------------------------------

function Base.iszero(a::DiffRingElem)
    return algdata(a) == 0
end

#------------------------------------------------------------------------------

function Base.:+(a::DiffRingElem, b::DiffRingElem)
    check_parent(a, b)
    return parent(a)(algdata(a) + algdata(b))
end

function Base.:+(a::DiffRingElem, b)
    return parent(a)(algdata(a) + b)
end

function Base.:+(a, b::DiffRingElem)
    return parent(b)(algdata(b) + a)
end

function Base.:+(a::DiffRingElem, b::Union{RingElem, AbstractFloat, Integer, Rational, Nemo.fmpq})
    return parent(a)(algdata(a) + b)
end

function AbstractAlgebra.addeq!(a::DiffRingElem, b::DiffRingElem)
    a = a + b
end

#------------------------------------------------------------------------------

function Base.:*(a::DiffRingElem, b::DiffRingElem)
    check_parent(a, b)
    return parent(a)(algdata(a) * algdata(b))
end

function Base.:*(a::RingElem, b::DiffRingElem)
    if typeof(a) <: DiffRingElem
        return parent(a)(algdata(a) * algdata(b))
    end
    return parent(b)(a * algdata(b))
end

#function Base.:*(a, b::DiffRingElem)
#    return parent(b)(algdata(b) * a)
#end

function Base.:*(a::Union{AbstractFloat, Integer, Rational}, b::DiffRingElem)
    return parent(b)(a * algdata(b))
end

function AbstractAlgebra.mul!(a::DiffRingElem, b::DiffRingElem, c::DiffRingElem)
    a = b * c
end

#------------------------------------------------------------------------------

function Base.:-(a::DiffRingElem, b::DiffRingElem)
    check_parent(a, b)
    return parent(a)(algdata(a) - algdata(b))
end

function Base.:-(a::DiffRingElem, b)
    return parent(a)(algdata(a) - b)
end

function Base.:-(a, b::DiffRingElem)
    return parent(b)(-algdata(b) + a)
end

function Base.:-(a::DiffRingElem, b::Union{RingElem, AbstractFloat, Integer, Rational, Nemo.fmpq})
    return parent(a)(algdata(a) - b)
end

#------------------------------------------------------------------------------

function Base.:-(a::DiffRingElem)
    return parent(a)(-algdata(a))
end

#------------------------------------------------------------------------------

function Base.:^(a::DiffRingElem, i::Integer)
    return parent(a)(algdata(a)^i)
end

#------------------------------------------------------------------------------

function Base.:(==)(a::DiffRingElem, b::DiffRingElem)
    return algdata(a) == algdata(b)
end

#------------------------------------------------------------------------------

function Base.:(//)(a::DiffRingElem, b::DiffRingElem)
    check_parent(a, b)
    return parent(a)(algdata(a) // algdata(b))
end

function Base.:(//)(a::DiffRingElem, b)
    return parent(a)(algdata(a) // b)
end

function Base.:(//)(a, b::DiffRingElem)
    return parent(b)(a // algdata(b))
end

# Fixing ambiguity
function Base.:(//)(a::Union{Integer, Rational}, b::DiffRingElem)
    return parent(b)(a // algdata(b))
end

#------------------------------------------------------------------------------

function Base.hash(a::DiffRingElem)
    return hash(algdata(a))
end

#------------------------------------------------------------------------------

function d_aux(p::MPolyElem, der::Dict{Any, Any})
    result = zero(parent(p))
    for v in vars(p)
        if !(v in keys(der))
            throw(DomainError("No derivative defined for $v. Most likely you have exceeded the maximal order."))
        end
        result += der[v] * derivative(p, v)
    end
    return result
end

function d(a::DiffPoly)
    return DiffPoly(parent(a), d_aux(algdata(a), parent(a).derivation))
end

function d(a::DiffRational)
    p, q = numerator(algdata(a)), denominator(algdata(a))
    dp, dq = map(x -> d_aux(x, parent(a).derivation), [p, q])
    return DiffRational(parent(a), (dp * q - p * dq) // q^2)
end

function d(a::DiffRingElem, ord::Integer)
    if ord == 0
        return a
    end
    return d(d(a), ord - 1)
end

#------------------------------------------------------------------------------

function AbstractAlgebra.degree(p::DiffPoly, v::DiffPoly)
    return degree(algdata(p), algdata(v))
end

function AbstractAlgebra.total_degree(p::DiffPoly)
    return total_degree(algdata(p))
end

function AbstractAlgebra.gcd(a::DiffPoly, b::DiffPoly)
    return DiffPoly(parent(a), gcd(algdata(a), algdata(b)))
end

function AbstractAlgebra.coeffs(a::DiffPoly)
    return coeffs(algdata(a))
end

function AbstractAlgebra.divexact(p::DiffPoly, q::DiffPoly)
    return DiffPoly(parent(p), divexact(algdata(p), algdata(q)))
end

#------------------------------------------------------------------------------

function order(p::DiffPoly, v::DiffPoly)
    R = parent(p)
    varname = decompose_derivative(var_to_str(algdata(v)))[1]
    orders = []
    for v in vars(algdata(p))
        vname, ord = decompose_derivative(var_to_str(v))
        if vname == varname
            push!(orders, ord)
        end
    end
    if length(orders) == 0
        return -1
    end
    return max(orders...)
end

function order(r::DiffRational, v::DiffPoly)
    return max(order(numerator(r), v), order(denominator(r), v))
end

#------------------------------------------------------------------------------

function leader(p::DiffPoly, v::DiffPoly)
    """
    Finds the highest derivative of variable v appearing in p
    """
    R = parent(p)
    h = order(p, v)
    if h == -1
        throw(DomainError("Variable $v does not apperar in $p"))
    end
    return get_der(R, decompose_derivative(var_to_str(algdata(v)))[1], h)
end

function leader(p::DiffPoly)
    if length(parent(p).varnames) != 1
        throw(DomainError("Expected a unvariate differential polynomial but got $p"))
    end
    return leader(p, gens(parent(p))[1])
end

#------------------------------------------------------------------------------

function initial(p::DiffPoly, v::DiffPoly)
    """
    Finds the initail of p with respect to variable v
    """
    ld = leader(p, v)
    d = degree(p, ld)
    initial = extract_coefficients(algdata(p), [algdata(ld)])[ [d] ]
    return DiffPoly(parent(p), initial)
end

function initial(p::DiffPoly)
    if length(parent(p).varnames) != 1
        throw(DomainError("Expected a unvariate differential polynomial but got $p"))
    end
    return initial(p, gens(parent(p))[1])
end

#------------------------------------------------------------------------------

function separant(p::DiffPoly, v::DiffPoly)
    """
    Finds the separant of p with respect to variable v
    """
    ld = leader(p, v)
    return DiffPoly(parent(p), derivative(algdata(p), algdata(ld)))
end

function separant(p::DiffPoly)
    if length(parent(p).varnames) != 1
        throw(DomainError("Expected a unvariate differential polynomial but got $p"))
    end
    return separant(p, gens(parent(p))[1])
end

#------------------------------------------------------------------------------

function diffred(g::DiffPoly, f::DiffPoly, v::DiffPoly)
    """
    Performs a differential reduction of g with respect to f considered as polynomials v
    Returns a triple (q, a, b), where
      - q is the remainder
      - a is the power of the separant
      - b is the power of the initial
    """
    check_parent(f, g)
    result = g
    a, b = 0, 0

    h = order(f, v)
    ld = leader(f, v)
    deg = degree(f, ld)
    sep = separant(f, v)
    init = initial(f, v)

    while order(result, v) > h
        H = order(result, v)
        ld_res = leader(result, v)
        D = degree(result, ld_res)
        result = result * sep - initial(result, v) * ld_res^(D - 1) * d(f, H - h)
        a += 1
    end

    while degree(result, ld) > deg
        D = degree(result, ld)
        result = result * init - initial(result, v) * ld^(D - deg) * f
        b += 1
    end

    return (result, a, b)
end

function diffred(g::DiffPoly, f::DiffPoly)
    check_parent(f, g)
    if length(parent(g).varnames) != 1
        throw(DomainError("Expected a unvariate differential polynomial but got $p"))
    end
    return diffred(g, f, gens(parent(g))[1])
end

#------------------------------------------------------------------------------

function evaluate_alg(a::DiffPoly, eval_point::Array{<: Any, 1})
    return evaluate(algdata(a), eval_point)
end

function evaluate_alg(a::DiffRational, eval_point::Array{<: Any, 1})
    return evaluate(numerator(algdata(a)), eval_point) // evaluate(denominator(algdata(a)), eval_point)
end

function Generic.evaluate(a::DiffRingElem, eval_point::Array{<: Any, 1})
    R = parent(a)
    eval_dict = Dict()
    for (i, v) in enumerate(gens(R))
        eval_dict[v] = eval_point[i]
        for ord in 1:order(a, v)
            prev = eval_dict[v]
            v = d(v)
            eval_dict[v] = d(prev)
        end
    end
    eval_point_alg = [get(eval_dict, R(v), zero(R)) for v in gens(R.poly_ring)]
    return evaluate_alg(a, eval_point_alg)
end

function Generic.evaluate(a::DiffRingElem, to_eval::Array{<: Any, 1}, values::Array{<: Any, 1})
    eval_dict = Dict()
    R = parent(a)
    for g in gens(R.poly_ring)
        dg = R(g)
        eval_dict[dg] = dg
    end
    for (e, v) in zip(to_eval, values)
        while true
            eval_dict[e] = v
            try
                e = d(e)
                v = d(v)
            catch e
                break
            end
        end
    end
    eval_point_alg = [get(eval_dict, R(v), zero(R)) for v in gens(R.poly_ring)]
    return evaluate_alg(a, eval_point_alg)
end

#------------------------------------------------------------------------------

# orders must be sorted in the nondecreasing order
function wronskian(polys::Array{DiffPoly, 1}, orders::Array{<: Integer, 1})
    R = parent(polys[1])
    n = length(polys)
    S = AbstractAlgebra.MatrixSpace(R, n, n)
    W = zero(S)
    for (i, p) in enumerate(polys)
        W[i, 1] = d(p, orders[1])
        for j in 2:n
            W[i, j] = d(W[i, j - 1], orders[j] - orders[j - 1])
        end
    end
    return W
end

function wronskian(polys::Array{DiffPoly, 1})
    return wronskian(polys, [i - 1 for i in 1:length(polys)])
end

function wronskian(polys::Array{DiffPoly, 1}, ind_to_omit::Integer)
    return wronskian(polys, [i for i in 0:length(polys) if i != ind_to_omit])
end

max_stars_repo_name: pogudingleb/DifferentialAlgebra

max_issues_repo_name: pogudingleb/DifferentialAlgebra

max_forks_repo_name: pogudingleb/DifferentialAlgebra


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

using Markdown
using InteractiveUtils

# ╔═╡ d757a3a8-1ee3-4afd-b0bd-9d8429250f96
begin
    import Pkg
    Pkg.develop("PoreMatMod")
end

# ╔═╡ 322d1bd1-86ae-47e7-b49d-38687eb3885a
using PoreMatMod

# ╔═╡ 494d2664-b5cf-43fe-850b-538a090dd0e8
using PoreMatMod.ExampleHelpers

# ╔═╡ 095d18b8-b8ab-4bd8-89dc-4382f9375b42
begin
    parent = Crystal("UiO-66.cif")
    for x in parent.atoms.coords.xf
        if x == 1.
            x -= eps(Float64)
        elseif x == 0.
            x += eps(Float64)
        end
    end
end

# ╔═╡ ac94ddb2-314c-49a7-b0bb-4371f93c0429
infer_bonds!(parent, true)

# ╔═╡ ecb86319-1c5c-48f8-98b9-8bccca5e9953
SBU = moiety("SBU.xyz")

# ╔═╡ d7b5abf5-13fb-4644-9502-6c54106be51a
bimetallic_SBU = moiety("bimetallic_SBU.xyz")

# ╔═╡ 44b0cc75-036b-468f-9da7-04c2d97536a9
begin
    multi_metal = replace(parent, SBU => bimetallic_SBU)
    translate_by!(multi_metal.atoms.coords, Frac([-0.25, -0.25, 0.]))
end

# ╔═╡ 780437f4-3a42-4aa2-b2e8-599d8d232dc2
write_cif(multi_metal, "multi-metal.cif")

# ╔═╡ ad888be3-2123-4d2f-865a-ff7a63cebf9e
begin
    shifted_uio66 = deepcopy(parent)
    translate_by!(shifted_uio66.atoms.coords, Frac([-0.25, -0.25, 0.]))
    write_cif(shifted_uio66, "shifted_parent.cif")
end

# ╔═╡ Cell order:
# ╠═d757a3a8-1ee3-4afd-b0bd-9d8429250f96
# ╠═322d1bd1-86ae-47e7-b49d-38687eb3885a
# ╠═494d2664-b5cf-43fe-850b-538a090dd0e8
# ╠═095d18b8-b8ab-4bd8-89dc-4382f9375b42
# ╠═ac94ddb2-314c-49a7-b0bb-4371f93c0429
# ╠═ecb86319-1c5c-48f8-98b9-8bccca5e9953
# ╠═d7b5abf5-13fb-4644-9502-6c54106be51a
# ╠═44b0cc75-036b-468f-9da7-04c2d97536a9
# ╠═780437f4-3a42-4aa2-b2e8-599d8d232dc2
# ╠═ad888be3-2123-4d2f-865a-ff7a63cebf9e


text:
function add_triangle!(point::Vector{SVector{3,Float64}}, i::Vector{SVector{3,Int64}}, a::SVector{3,Float64},
        b::SVector{3,Float64}, c::SVector{3,Float64})

    push!(i, SVector{3,Int64}(1, 2, 3) .+ length(point))
    push!(point, a, b, c)
    return nothing
end

function add_rectangle!(point::Vector{SVector{3,Float64}}, i::Vector{SVector{3,Int64}}, a⁻::SVector{3,Float64},
        a⁺::SVector{3,Float64}, b⁻::SVector{3,Float64}, b⁺::SVector{3,Float64})

    add_triangle!(point, i, a⁻, a⁺, b⁺)
    add_triangle!(point, i, a⁻, b⁺, b⁻)
end

function add_slice(point::Vector{SVector{3,Float64}}, i::Vector{SVector{3,Int64}}, v::NTuple{N,SVector{3,Float64}},
        θ⁰::Float64, θ¹::Float64) where {N}
    # adds a rotationally symmetric slice of point sequence v to the vector of points (point) and vector of indices i

    function add_rot!(point::Vector{SVector{3,Float64}}, i::Vector{SVector{3,Int64}}, a::SVector{3,Float64},
            b::SVector{3,Float64}, RZ_0::RotZ{Float64}, RZ_1::RotZ{Float64})

        a⁻ = RZ_0 * a
        a⁺ = RZ_1 * a
        b⁻ = RZ_0 * b
        b⁺ = RZ_1 * b

        # TODO: is this necessary given mesh repair_mesh?
        if (a[1] == 0.0)  # && (a[2] == 0.0)  # point has no xy extent
            add_triangle!(point, i, a, b⁺, b⁻)
        elseif (b[1] == 0.0)  # && (b[2] == 0.0)
            add_triangle!(point, i, a⁻, a⁺, b)
        else
            add_rectangle!(point, i, a⁻, a⁺, b⁻, b⁺)
        end
    end

    RZ_0 = RotZ(θ⁰)
    RZ_1 = RotZ(θ¹)
    for k = 1:(N - 1)
        add_rot!(point, i, v[k], v[k + 1], RZ_0, RZ_1)
    end
end

function obj_from_point_sequence(point_vec_2D::Vector{SVector{2,Float64}}, n_theta::Int64=10)
    # creates a rotationally symmetric object by rotating a vector of 2D points (x,z) about the z axis to create a
    # vector of 3D points (x,y,z)

    ### FIXES degenerate points near the axis of rotation
    for k = 1:length(point_vec_2D)
        item_1, item_2 = point_vec_2D[k]
        tol = 1.0e-12
        if item_1 <= -tol
            error("negative extent")
        # elseif item_1 == 0.0
        #     #
        # elseif item_1 <= tol
        elseif 0.0 < item_1 <= tol
            @warn("point near rotation axis set to 0.0")
            point_vec_2D[k] = SVector{2,Float64}(0.0, item_2)
        end
    end

    point = Vector{SVector{3,Float64}}()
    i = Vector{SVector{3,Int64}}()
    point_vec = [SVector{3,Float64}(k[1], 0.0, k[2]) for k = point_vec_2D]

    theta = collect(LinRange{Float64}(0, 2*pi, n_theta + 1)) .+ (pi/2)
    for k = 1:n_theta
        add_slice(point, i, Tuple(point_vec), theta[k], theta[mod1(k + 1, n_theta)])
    end

    hex_mesh = eMesh(point, i, nothing, nothing)
    mesh_repair!(hex_mesh)
    return hex_mesh
end


text:
demand = nhhy * (1-unemp)*(1-tau)* ((prodlb+produb)/2)/nf

space = GridSpace((k1,k2))

covidmodel = ABM(CovidAgent.GenericAgent,space; scheduler = fastest)

#create emtpy list of employees for all firms
workerlist =  Array{Int}[]
for n=1:nf
    push!(workerlist, Int[])
end
#unemployd list
unemplist = [Dict{Int,CovidAgent.GenericAgent}() for tt=1:nsec, i=1:k1, j=1:k2] # only dimension nsec because there are no unemployd public workers
shorttimelist = [Dict{Int,CovidAgent.GenericAgent}() for tt=1:nsec, i=1:k1, j=1:k2]

#generate list of firms in each location
locf = [Int[] for tt = 1:nsec+1, i=1:k1, j=1:k2] # separated by type of firm
locallf = [Int[] for i=1:k1, j=1:k2] # all types

# first setor
for n =1: nfsec[1]
    k1help = 1+fld(mod(n-1,k1*k2),k2)
    k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
    push!(locf[1,k1help,k2help],n)
    push!(locallf[k1help,k2help],n)
end

for tt = 2:nsec
    for n = sum(nfsec[1:tt-1])+1: sum(nfsec[1:tt])
        k1help = 1+fld(mod(n-1,k1*k2),k2)
        k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
        push!(locf[tt,k1help,k2help],n)
        push!(locallf[k1help,k2help],n)
    end
end
for n = sum(nfsec[1:nsec])+1: nf
    k1help = 1+fld(mod(n-1,k1*k2),k2)
    k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
    push!(locf[nsec+1,k1help,k2help],n)
    push!(locallf[k1help,k2help],n)
end

#generate list of hh in each location
lochh = [Dict{Int,CovidAgent.GenericAgent}() for tt = 1:2, i=1:k1, j=1:k2]

#list of households
hh  = [collect(1:nhhy), collect(nhhy+1:nhhy+nhho)]

#list of firms
firms = collect(1:nf)

global empcount = zeros(nsec+1)
global unempcount = zeros(nsec)
global shorttimecount = zeros(nsec)
global oldcount = 0

# add young households
for n =1:nhhy
    k1help = 1+fld(mod(n-1,k1*k2),k2)
    k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
    if rand() < unemp
        emp=0
        tt = conselect(fracemp) # no unemployed for public sector
        unempcount[tt] += 1
    else
        emp = rand(locallf[k1help,k2help])
        push!(workerlist[emp],n)
        if emp <= nfsec[1]
            tt = 1
        end
        for tth=2:nsec+1
            if sum(nfsec[1:tth-1]) < emp && emp <= sum(nfsec[1:tth])
                tt = tth
            end
        end

        empcount[tt] +=1
    end
    if rand() < homeofficefrac[tt] # can this worker do home-office ?
        ho = true
    else
        ho = false
    end
    shopday = rand(0:6,nsec)
    agent = CovidAgent.GenericAgent(n,(k1help,k2help),tt,1,1,false,0,0,emp,false,0,ho,Dict{Int,CovidAgent.GenericAgent}(),Dict{Int,CovidAgent.GenericAgent}(),incomehh[tt],savinghh[tt],consbud,prevconsbud,shopday,shoprepini,0,0,0,0,0,pastrevini, 0,0)

    add_agent_pos!(agent,covidmodel)
    lochh[1,k1help,k2help][n] = agent
    if emp == 0
        unemplist[tt,k1help,k2help][n] = agent
    end
end

# add old households
for n = nhhy+1:nhhy+nhho
    global oldcount += 1
    k1help = 1+fld(mod(n-1,k1*k2),k2)
    k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
    shopday = rand(0:6,nsec)
    agent = CovidAgent.GenericAgent(n,(k1help,k2help),1,2,1,false,0,0,0,false,0,false,Dict{Int,CovidAgent.GenericAgent}(),Dict{Int,CovidAgent.GenericAgent}(),incomehho,savinghho,consbud,prevconsbud,shopday,shoprepini,0,0,0,0,0,pastrevini,0,0)

    add_agent_pos!(agent,covidmodel)
    lochh[2,k1help,k2help][n] = agent
end

# add firms
#first sector
for n = 1: nfsec[1]
    k1help = 1+fld(mod(n-1,k1*k2),k2)
    k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
    product =  (prodlb[1]+ (produb[1]-prodlb[1])*rand())*baseprod[1] # determine productivity
    shopday = [0,0,0]

    worker = Dict{Int,CovidAgent.GenericAgent}()
    for id in workerlist[n]
        worker[id] = covidmodel[id]
    end

    agent = CovidAgent.GenericAgent(nhhy+nhho+n,(k1help,k2help),nsec+2,1,1,false,0,0,0,false,0,false,worker,Dict{Int,CovidAgent.GenericAgent}(),0,savingf,consbud,prevconsbud,shopday,shoprepini,product,0,demand[1],demandexpini[1],0,pastrevini,etmin[1],0)
    add_agent_pos!(agent,covidmodel)
end
#other sectors
for tt = 2:nsec+1
    for n = sum(nfsec[1:tt-1])+1: sum(nfsec[1:tt])
        k1help = 1+fld(mod(n-1,k1*k2),k2)
        k2help = 1+mod(n-1,k1*k2)-(k1help-1)*k2 # distribute evenly on grid
        product =  (prodlb[tt]+ (produb[tt]-prodlb[tt])*rand())*baseprod[tt] # determine productivity
        shopday = [0,0,0]

        worker = Dict{Int,CovidAgent.GenericAgent}()
        for id in workerlist[n]
            worker[id] = covidmodel[id]
        end

        agent = CovidAgent.GenericAgent(nhhy+nhho+n,(k1help,k2help),nsec+1+tt,1,1,false,0,0,0,false,0,false,worker,Dict{Int,CovidAgent.GenericAgent}(),0,savingf,consbud,prevconsbud,shopday,shoprepini,product,0,demand[tt],demandexpini[tt],0,pastrevini, etmin[tt], 0)
        add_agent_pos!(agent,covidmodel)
    end
end

divperhh = 0
weeklyconsumption = zeros(nsec, k1, k2)

covidmodel_copy = copy(covidmodel.agents)


text:
s = readline()
p = readline()

s = s^3
if occursin(p, s)
    println("Yes")
else
    println("No")
end


text:
#### Turbulence closures
using DocStringExtensions
using CLIMA.PlanetParameters
using CLIMA.SubgridScaleParameters
export ConstantViscosityWithDivergence, SmagorinskyLilly, Vreman, AnisoMinDiss

abstract type TurbulenceClosure end

vars_state(::TurbulenceClosure, FT) = @vars()
vars_gradient(::TurbulenceClosure, FT) = @vars()
vars_diffusive(::TurbulenceClosure, FT) = @vars()
vars_aux(::TurbulenceClosure, FT) = @vars()

function atmos_init_aux!(::TurbulenceClosure, ::AtmosModel, aux::Vars, geom::LocalGeometry)
end
function atmos_nodal_update_aux!(::TurbulenceClosure, ::AtmosModel, state::Vars, aux::Vars, t::Real)
end
function diffusive!(::TurbulenceClosure, diffusive, ∇transform, state, aux, t, ν)
end
function gradvariables!(::TurbulenceClosure, transform::Vars, state::Vars, aux::Vars, t::Real)
end

"""
  PrincipalInvariants{FT}

Calculates principal invariants of a tensor. Returns struct with fields first,second,third
referring to each of the invariants.
"""
struct PrincipalInvariants{FT}
  first::FT
  second::FT
  third::FT
end
function compute_principal_invariants(X::StaticArray{Tuple{3,3}})
  first = tr(X)
  second = 1/2 *((tr(X))^2 - tr(X .^ 2))
  third = det(X)
  return PrincipalInvariants{eltype(X)}(first,second,third)
end

"""
    ConstantViscosityWithDivergence <: TurbulenceClosure

Turbulence with constant dynamic viscosity (`ρν`). Divergence terms are included in the momentum flux tensor.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantViscosityWithDivergence{FT} <: TurbulenceClosure
  "Dynamic Viscosity [kg/m/s]"
  ρν::FT
end
function dynamic_viscosity_tensor(m::ConstantViscosityWithDivergence, S,
  state::Vars, diffusive::Vars, ∇transform::Grad, aux::Vars, t::Real)
  return m.ρν
end
function scaled_momentum_flux_tensor(m::ConstantViscosityWithDivergence, ρν, S)
  @inbounds trS = tr(S)
  return (-2*ρν) * S + (2*ρν/3)*trS * I
end

"""
    SmagorinskyLilly <: TurbulenceClosure

  § 1.3.2 in CliMA documentation

  article{doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2,
  author = {Smagorinksy, J.},
  title = {General circulation experiments with the primitive equations},
  journal = {Monthly Weather Review},
  volume = {91},
  number = {3},
  pages = {99-164},
  year = {1963},
  doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  URL = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  eprint = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
  }

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SmagorinskyLilly{FT} <: TurbulenceClosure
  "Smagorinsky Coefficient [dimensionless]"
  C_smag::FT
end

vars_aux(::SmagorinskyLilly,T) = @vars(Δ::T)
vars_gradient(::SmagorinskyLilly,T) = @vars(θ_v::T)

function atmos_init_aux!(::SmagorinskyLilly, ::AtmosModel, aux::Vars, geom::LocalGeometry)
  aux.turbulence.Δ = lengthscale(geom)
end

function gradvariables!(m::SmagorinskyLilly, transform::Vars, state::Vars, aux::Vars, t::Real)
  transform.turbulence.θ_v = aux.moisture.θ_v
end

"""
  buoyancy_correction(normSij, θᵥ, dθᵥdz)
  return buoyancy_factor, scaling coefficient for Standard Smagorinsky Model
  in stratified flows

Compute the buoyancy adjustment coefficient for stratified flows
given the strain rate tensor inner product |S| ≡ SijSij ≡ normSij,
local virtual potential temperature θᵥ and the vertical potential
temperature gradient dθvdz.

Brunt-Vaisala frequency N² defined as in equation (1b) in
  Durran, D.R. and J.B. Klemp, 1982:
  On the Effects of Moisture on the Brunt-Väisälä Frequency.
  J. Atmos. Sci., 39, 2152–2158,
  https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2

Ri = N² / (2*normSij)
Ri = gravity / θᵥ * ∂θᵥ∂z / 2 |S_{ij}|

§1.3.2 in CliMA documentation.

article{doi:10.1111/j.2153-3490.1962.tb00128.x,
author = {LILLY, D. K.},
title = {On the numerical simulation of buoyant convection},
journal = {Tellus},
volume = {14},
number = {2},
pages = {148-172},
doi = {10.1111/j.2153-3490.1962.tb00128.x},
url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1962.tb00128.x},
eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1962.tb00128.x},
year = {1962}
}
"""
function squared_buoyancy_correction(normS, ∇transform::Grad, aux::Vars)
  ∂θ∂Φ = dot(∇transform.turbulence.θ_v, aux.orientation.∇Φ)
  N² = ∂θ∂Φ / aux.moisture.θ_v
  Richardson = N² / (normS^2 + eps(normS))
  sqrt(clamp(1 - Richardson*inv_Pr_turb, 0, 1))
end

function strain_rate_magnitude(S::SHermitianCompact{3,FT,6}) where {FT}
  sqrt(2*S[1,1]^2 + 4*S[2,1]^2 + 4*S[3,1]^2 + 2*S[2,2]^2 + 4*S[3,2]^2 + 2*S[3,3]^2)
end

function dynamic_viscosity_tensor(m::SmagorinskyLilly, S, state::Vars, diffusive::Vars, ∇transform::Grad, aux::Vars, t::Real)
  # strain rate tensor norm
  # Notation: normS ≡ norm2S = √(2S:S)
  # ρν = (Cₛ * Δ * f_b)² * √(2S:S)
  FT = eltype(state)
  @inbounds normS = strain_rate_magnitude(S)
  f_b² = squared_buoyancy_correction(normS, ∇transform, aux)
  # Return Buoyancy-adjusted Smagorinsky Coefficient (ρ scaled)
  return state.ρ * normS * f_b² * FT(m.C_smag * aux.turbulence.Δ)^2
end
function scaled_momentum_flux_tensor(m::SmagorinskyLilly, ρν, S)
  (-2*ρν) * S
end

"""
  Vreman{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient
is specified and used to compute the equivalent Vreman coefficient.

1) ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
2) βij = Δ²αₘᵢαₘⱼ
3) Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²
βᵢⱼ is symmetric, positive-definite.
If Δᵢ = Δ, then β = Δ²αᵀα

@article{Vreman2004,
  title={An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications},
  author={Vreman, AW},
  journal={Physics of fluids},
  volume={16},
  number={10},
  pages={3670--3681},
  year={2004},
  publisher={AIP}
}

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Vreman{FT} <: TurbulenceClosure
  "Smagorinsky Coefficient [dimensionless]"
  C_smag::FT
end
vars_aux(::Vreman,FT) = @vars(Δ::FT)
vars_gradient(::Vreman,FT) = @vars(θ_v::FT)
function atmos_init_aux!(::Vreman, ::AtmosModel, aux::Vars, geom::LocalGeometry)
  aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(m::Vreman, transform::Vars, state::Vars, aux::Vars, t::Real)
  transform.turbulence.θ_v = aux.moisture.θ_v
end
function dynamic_viscosity_tensor(m::Vreman, S, state::Vars, diffusive::Vars, ∇transform::Grad, aux::Vars, t::Real)
  FT = eltype(state)
  ∇u = ∇transform.u
  αijαij = sum(∇u .^ 2)
  @inbounds normS = strain_rate_magnitude(S)
  f_b² = squared_buoyancy_correction(normS, ∇transform, aux)
  βij = f_b² * (aux.turbulence.Δ)^2 * (∇u' * ∇u)
  Bβinvariants = compute_principal_invariants(βij)
  @inbounds Bβ = Bβinvariants.second
  return state.ρ * max(0,m.C_smag^2 * 2.5 * sqrt(abs(Bβ/(αijαij+eps(FT)))))
end
function scaled_momentum_flux_tensor(m::Vreman, ρν, S)
  (-2*ρν) * S
end

"""
  AnisoMinDiss{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient
is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the
eigenvalue problem for the Laplacian operator).

@article{doi:10.1063/1.4928700,
author = {Rozema,Wybe  and Bae,Hyun J.  and Moin,Parviz  and Verstappen,Roel },
title = {Minimum-dissipation models for large-eddy simulation},
journal = {Physics of Fluids},
volume = {27},
number = {8},
pages = {085107},
year = {2015},
doi = {10.1063/1.4928700},
URL = {https://aip.scitation.org/doi/abs/10.1063/1.4928700},
eprint = {https://aip.scitation.org/doi/pdf/10.1063/1.4928700}
}
-------------------------------------------------------------------------------------
# TODO: Future versions will include modifications of Abkar(2016), Verstappen(2018)
@article{PhysRevFluids.1.041701,
title = {Minimum-dissipation scalar transport model for large-eddy simulation of turbulent flows},
author = {Abkar, Mahdi and Bae, Hyun J. and Moin, Parviz},
journal = {Phys. Rev. Fluids},
volume = {1},
issue = {4},
pages = {041701},
numpages = {10},
year = {2016},
month = {Aug},
publisher = {American Physical Society},
doi = {10.1103/PhysRevFluids.1.041701},
url = {https://link.aps.org/doi/10.1103/PhysRevFluids.1.041701}
}

"""
struct AnisoMinDiss{FT} <: TurbulenceClosure
  C_poincare::FT
end
vars_aux(::AnisoMinDiss,T) = @vars(Δ::T)
vars_gradient(::AnisoMinDiss,T) = @vars(θ_v::T)
function atmos_init_aux!(::AnisoMinDiss, ::AtmosModel, aux::Vars, geom::LocalGeometry)
  aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(m::AnisoMinDiss, transform::Vars, state::Vars, aux::Vars, t::Real)
  transform.turbulence.θ_v = aux.moisture.θ_v
end
function dynamic_viscosity_tensor(m::AnisoMinDiss, S, state::Vars, diffusive::Vars, ∇transform::Grad, aux::Vars, t::Real)
  FT = eltype(state)
  ∇u = ∇transform.u
  αijαij = dot(∇u,∇u)
  coeff = (aux.turbulence.Δ * m.C_poincare)^2
  βij = -(∇u' * ∇u)
  ν_e = max(0,coeff * (dot(βij, S) / (αijαij + eps(FT))))
  return state.ρ * ν_e
end
function scaled_momentum_flux_tensor(m::AnisoMinDiss, ρν, S)
  (-2*ρν) * S
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
@testset "deprecated" begin
    a = allowmissing(1.0:1.0:20.0)
    a[[2, 3, 7]] .= missing
    mask = map(!ismissing, a)

    @testset "Drop" begin
        result = impute(a, :drop; limit=0.2)
        expected = copy(a)
        deleteat!(expected, [2, 3, 7])

        @test result == expected

        # Mutating method
        a2 = copy(a)
        Impute.drop!(a2; limit=0.2)
        @test a2 == expected
    end

    @testset "Interpolate" begin
        result = impute(a, :interp; limit=0.2)
        @test result == collect(1.0:1.0:20)
        @test result == interp(a)

        # Test in-place method
        a2 = copy(a)
        Impute.interp!(a2; limit=0.2)
        @test a2 == result

        # Test interpolation between identical points
        b = ones(Union{Float64, Missing}, 20)
        b[[2, 3, 7]] .= missing
        @test interp(b) == ones(Union{Float64, Missing}, 20)

        # Test interpolation at endpoints
        b = ones(Union{Float64, Missing}, 20)
        b[[1, 3, 20]] .= missing
        result = interp(b)
        @test ismissing(result[1])
        @test ismissing(result[20])
    end

    @testset "Fill" begin
        @testset "Value" begin
            fill_val = -1.0
            result = impute(a, :fill, fill_val; limit=0.2)
            expected = copy(a)
            expected[[2, 3, 7]] .= fill_val

            @test result == expected
        end

        @testset "Mean" begin
            result = impute(a, :fill; limit=0.2)
            expected = copy(a)
            expected[[2, 3, 7]] .= mean(a[mask])

            @test result == expected

            a2 = copy(a)
            Impute.fill!(a2; limit=0.2)
            @test a2 == result
        end
    end

    @testset "LOCF" begin
        result = impute(a, :locf; limit=0.2)
        expected = copy(a)
        expected[2] = 1.0
        expected[3] = 1.0
        expected[7] = 6.0

        @test result == expected
        a2 = copy(a)
        impute!(a2, :locf; limit=0.2)
        @test a2 == result
    end

    @testset "NOCB" begin
        result = impute(a, :nocb; limit=0.2)
        expected = copy(a)
        expected[2] = 4.0
        expected[3] = 4.0
        expected[7] = 8.0

        @test result == expected
        a2 = copy(a)
        Impute.nocb!(a2; limit=0.2)
        @test a2 == result
    end

    @testset "DataFrame" begin
        data = dataset("boot", "neuro")
        df = impute(data, :interp; limit=1.0)
    end

    @testset "Matrix" begin
        data = Matrix(dataset("boot", "neuro"))

        @testset "Drop" begin
            result = Iterators.drop(data)
            @test size(result, 1) == 4
        end

        @testset "Fill" begin
            result = impute(data, :fill, 0.0; limit=1.0)
            @test size(result) == size(data)
        end
    end

    @testset "Not enough data" begin
        @test_throws ImputeError impute(a, :drop)
    end

    @testset "Chain" begin
        orig = dataset("boot", "neuro")

        @testset "DataFrame" begin
            result = chain(
                orig,
                Impute.Interpolate(),
                Impute.LOCF(),
                Impute.NOCB();
                limit=1.0
            )

            @test size(result) == size(orig)
            # Confirm that we don't have any more missing values
            @test all(!ismissing, Matrix(result))
        end

        @testset "Column Table" begin
            data = Tables.columntable(orig)
            result = chain(
                data,
                Impute.Interpolate(),
                Impute.LOCF(),
                Impute.NOCB();
                limit=1.0
            ) |> Tables.matrix

            @test size(result) == size(orig)
            # Confirm that we don't have any more missing values
            @test all(!ismissing, result)
        end

        @testset "Matrix" begin
            data = Matrix(orig)
            result = chain(
                data,
                Impute.Interpolate(),
                Impute.LOCF(),
                Impute.NOCB();
                limit=1.0
            )

            @test size(result) == size(data)
            # Confirm that we don't have any more missing values
            @test all(!ismissing, result)
        end
    end

    @testset "Alternate missing functions" begin
        data1 = dataset("boot", "neuro")                # Missing values with `missing`
        data2 = impute(data1, :fill, NaN; limit=1.0)     # Missing values with `NaN`

        @test impute(data1, :drop; limit=1.0) == dropmissing(data1)

        result1 = chain(data1, Impute.Interpolate(), Impute.Drop(); limit=1.0)
        result2 = chain(data2, isnan, Impute.Interpolate(), Impute.Drop(); limit=1.0)
        @test result1 == result2

        @test Impute.drop(data1; limit=1.0) == impute(data2, isnan, :drop; limit=1.0)
    end
end


text:
module SyntheticDataGeneration

export
    MarkerPositionGenerationOptions,
    PoseDataGenerationOptions,
    generate_marker_positions,
    generate_joint_offset,
    generate_pose_data,
    generate_calibration_problem

using MotionCaptureJointCalibration
using RigidBodyDynamics
using StaticArrays
using Parameters
using Base.Iterators

using Random: randperm, rand!
using MotionCaptureJointCalibration: Point3DS

@with_kw struct MarkerPositionGenerationOptions
    marker_measurement_stddev::Float64 = 1e-5
    marker_offset_max::Float64 = 0.1
    num_markers::Int = 4
    num_measured_markers::Int = 3
end

@with_kw struct PoseDataGenerationOptions
    num_poses::Int = 25
    motion_capture_noise_stddev::Float64 = 1e-6
    joint_configuration_noise_stddev::Float64 = 1e-4
    occlusion_probability::Float64 = 0.025
end

function generate_marker_positions(bodies::AbstractVector{<:RigidBody},
        options::MarkerPositionGenerationOptions = MarkerPositionGenerationOptions())
    # Marker positions in body frame
    B = eltype(bodies)
    T = Float64
    ground_truth_marker_positions = Dict(b => Vector{Point3DS{T}}() for b in bodies)
    measured_marker_positions = Dict(b => Vector{Tuple{Point3DS{T}, Point3DS{T}}}() for b in bodies)
    for body in bodies
        frame = default_frame(body) # TODO: use some other frame to improve test coverage
        measured_marker_inds = randperm(options.num_markers)[1 : options.num_measured_markers]
        for i = 1 : options.num_markers
            ground_truth = Point3D(frame, (rand(SVector{3}) - 0.5) * 2 * options.marker_offset_max)
            push!(ground_truth_marker_positions[body], ground_truth)
            bounds = if i ∈ measured_marker_inds
                measurement_error = FreeVector3D(frame, options.marker_measurement_stddev * randn(SVector{3}))
                measured = ground_truth + measurement_error
                measured, measured
            else
                lower = Point3D(frame, fill(-0.2, SVector{3}))
                upper = Point3D(frame, fill(0.2, SVector{3}))
                lower, upper
            end
            push!(measured_marker_positions[body], bounds)
        end
    end
    ground_truth_marker_positions, measured_marker_positions
end

function generate_joint_offset(joint::Joint, max_offset::Number)
    (rand(num_positions(joint)) .- 0.5) .* max_offset .* 2
end

function generate_joint_offsets(joints::AbstractVector{<:Joint}, max_offset::Number)
    Dict(j => generate_joint_offset(j, max_offset) for j in joints)
end

function generate_pose_data(
        state::MechanismState{X, M, C},
        ground_truth_marker_positions::AbstractDict{<:RigidBody{M}, <:AbstractVector{Point3DS{T}}},
        ground_truth_offsets::AbstractDict{<:Joint{M}, <:AbstractVector{T}},
        free_joints::AbstractVector{<:Joint{M}},
        options::PoseDataGenerationOptions = PoseDataGenerationOptions()) where {X, M, C, T}
    ground_truth_pose_data = Vector{PoseData{C}}()
    measured_pose_data = Vector{PoseData{C}}()
    mechanism = state.mechanism
    for i = 1 : options.num_poses
        rand!(state)

        # Joint configurations
        q_ground_truth = copy(configuration(state))
        q_measured = copy(q_ground_truth)
        for (joint, offset) in ground_truth_offsets
            qjoint = q_measured[joint]
            qjoint .+= offset
        end
        for joint in free_joints
            qjoint = q_measured[joint]
            zero_configuration!(qjoint, joint)
        end
        for joint in setdiff(tree_joints(mechanism), free_joints)
            qjoint = q_measured[joint]
            qjoint .+= options.joint_configuration_noise_stddev * randn(num_positions(joint))
        end

        # Markers
        S = promote_type(C, T)
        markerbodies = keys(ground_truth_marker_positions)
        ground_truth_marker_data = Dict(b => Vector{Point3DS{S}}() for b in markerbodies)
        measured_marker_data = Dict(b => Vector{Point3DS{S}}() for b in markerbodies)
        for body in markerbodies
            toworld = transform_to_root(state, body)
            for (j, marker_body) in enumerate(ground_truth_marker_positions[body])
                marker_world = transform(marker_body, toworld)
                push!(ground_truth_marker_data[body], marker_world)
                occluded = rand() < options.occlusion_probability
                measured = if occluded
                    Point3D(marker_world.frame, fill(S(NaN), SVector{3}))
                else
                    noise = FreeVector3D(marker_world.frame, options.motion_capture_noise_stddev * randn(SVector{3}))
                    marker_world + noise
                end
                push!(measured_marker_data[body], measured)
            end
        end
        push!(ground_truth_pose_data, PoseData(q_ground_truth, ground_truth_marker_data))
        push!(measured_pose_data, PoseData(q_measured, measured_marker_data))
    end
    ground_truth_pose_data, measured_pose_data
end

function generate_calibration_problem(state::MechanismState{T}, body_weights::Dict{RigidBody{T}, T};
        marker_options::MarkerPositionGenerationOptions = MarkerPositionGenerationOptions(),
        pose_options::PoseDataGenerationOptions = PoseDataGenerationOptions()) where {T}
    bodies = collect(keys(body_weights))
    mechanism = state.mechanism
    correction_joints = unique(flatten([collect(RigidBodyDynamics.path(mechanism, body1, body2)) for (body1, body2) in product(bodies, bodies)]))
    calibration_param_bounds = Dict{Joint{T}, Vector{Tuple{T, T}}}(j => fill((-0.05, 0.05), num_positions(j)) for j in correction_joints)
    free_joint_configuration_bounds = Dict{Joint{T}, Vector{Tuple{T, T}}}(
        j => fill((-1., 1.), num_positions(j)) for j in tree_joints(mechanism) if isfloating(j))
    free_joints = collect(keys(free_joint_configuration_bounds))
    ground_truth_marker_positions, measured_marker_positions = generate_marker_positions(bodies, marker_options)
    ground_truth_offsets = Dict{Joint{T}, Vector{T}}(j => generate_joint_offset(j, 1e-2) for j in correction_joints)
    ground_truth_pose_data, measured_pose_data =
        generate_pose_data(state, ground_truth_marker_positions, ground_truth_offsets, free_joints, pose_options)
    problem = CalibrationProblem(
        mechanism,
        calibration_param_bounds,
        free_joint_configuration_bounds,
        measured_marker_positions,
        measured_pose_data, body_weights)
    configurations = [data.configuration for data in ground_truth_pose_data]
    ground_truth = CalibrationResult(:Optimal, 0., ground_truth_offsets, configurations, ground_truth_marker_positions)
    problem, ground_truth
end

end # module


text:
export
    default_grid_size, plot_function_around, plot_local_slices

"""
    $SIGNATURES

Find integers `w` and `h` such that `w * h ≥ n` and `abs2(w/h-aspect_ratio)` is
minimized. Return `(w, h)`.

Useful for determining grid size for `n` subplots.

The algorithm uses a heuristic, making an initial guess as the closest integer
to `√(aspect_ratio/r)`, and minimizing the objective within the nearest `± d`
integers.
"""
function default_grid_size(n; aspect_ratio = 4//3, d = 2)
    hs = filter(h -> h ∈ 1:n, @. round(√(n/aspect_ratio) + (-d:d)))
    ws = cld.(n, hs)
    _, ix = findmin(@. abs2((ws / hs) - aspect_ratio))
    ws[ix], hs[ix]
end

"""
    $SIGNATURES

Return `PGFPlotsX.Plot` elements visualizing the function ``f`` on ``x ± Δ``,
using `n` points. For plotting directly, use

```julia
Axis(plot_function_around(f, x, Δ)...)
```
"""
function plot_function_around(f, x::Real, Δ::Real; n = 100)
    xs = linspace(x - Δ, x + Δ, n)
    @pgf (Plot(Table(xs, f.(xs))),
          Plot({ only_marks }, Table([x], [f(x)])))
end

"""
    $SIGNATURES

Plot slices of the ``ℝⁿ→ℝ`` around `x ± Δs`. Return a `GroupPlot`, with subplots
arranged using `default_grid_size` using `aspect_ratio`. `xlabelfun` generates
labels for coordinates from the coordinate index, while `ylabel` provides the
`y` label.
"""
function plot_local_slices(f, x::AbstractVector, Δs::AbstractVector;
                           aspect_ratio = 4//3,
                           grid_size = default_grid_size(length(x)),
                           xlabelfun = i -> "\$x_{$(i)}\$",
                           ylabel = raw"$f$")
    @argcheck length(x) == length(Δs)
    w, h = grid_size
    plots = []
    for (i, Δ) in enumerate(Δs)
        push!(plots, @pgf { xlabel = xlabelfun(i) })
        append!(plots, plot_function_around(function_slice(f, x, i),
                                            zero(eltype(x)), Δ))
    end
    @pgf GroupPlot({ group_style = { group_size = "$w by $h"},
                     ylabel = ylabel},
                   plots...)
end


text:
"""
**Laplace likelihood**

Laplace likelihood for regression: ``\\frac{1}{2\\beta}\\exp\\left(-\\frac{|y-f|}{\\beta}\\right)``
see [wiki page](https://en.wikipedia.org/wiki/Laplace_distribution)

```julia
LaplaceLikelihood(β::T=1.0)  #  Laplace likelihood with scale β
```

---

For the analytical solution, it is augmented via:
```math
p(y|f,\\omega) = \\mathcal{N}(y|f,\\omega^{-1})
```
where ``\\omega \\sim \\text{Exp}\\left(\\omega \\mid \\frac{1}{2 \\beta^2}\\right)``, and Exp is the [Exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution)
We approximate ``q(\\omega) = \\mathcal{GIG}\\left(\\omega \\mid a,b,p\\right)
"""
struct LaplaceLikelihood{T<:Real} <: RegressionLikelihood{T}
    β::LatentArray{T}
    a::LatentArray{T}
    p::LatentArray{T}
    b::LatentArray{Vector{T}}
    θ::LatentArray{Vector{T}} #Expected value of ω
    function LaplaceLikelihood{T}(β::AbstractVector{T}) where {T<:Real}
        new{T}(β,β.^-2,0.5*ones(size(β)))
    end
    function LaplaceLikelihood{T}(β::AbstractVector{T},b::AbstractVector{<:AbstractVector{T}},θ::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
        new{T}(β,β.^(-2),0.5*ones(size(β)),b,θ)
    end
end

function LaplaceLikelihood(β::T=1.0) where {T<:Real}
    LaplaceLikelihood{T}([β])
end

function init_likelihood(likelihood::LaplaceLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        LaplaceLikelihood{T}(
        likelihood.β,
        [abs2.(T.(rand(T,nSamplesUsed))) for _ in 1:nLatent],
        [zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        LaplaceLikelihood{T}(likelihood.β)
    end
end

function pdf(l::LaplaceLikelihood,y::Real,f::Real)
    Distributions.pdf(Laplace(f,l.β[1]),y) #WARNING multioutput invalid
end

function Base.show(io::IO,model::LaplaceLikelihood{T}) where T
    print(io,"Laplace likelihood")
end

function compute_proba(l::LaplaceLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    return μ,max.(σ²,0.0).+ 2*l.β[1]^2
end

## Local Updates ##

function local_updates!(model::Union{VGP{T,<:LaplaceLikelihood,<:AnalyticVI},VStP{T,<:LaplaceLikelihood,<:AnalyticVI}}) where {T}
    model.likelihood.b .= broadcast((Σ,μ,y)->(Σ+abs2.(μ-y)),diag.(model.Σ),model.μ,model.y)
    model.likelihood.θ .= broadcast((a,b)->sqrt(a)./sqrt.(b),model.likelihood.a,model.likelihood.b)
end

function local_updates!(model::SVGP{T,<:LaplaceLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.b .= broadcast((K̃,κ,Σ,μ,y)->(K̃ + opt_diag(κ*Σ,κ) + abs2.(κ*μ-y)),model.K̃,model.κ,model.Σ,model.μ,model.inference.y)
    model.likelihood.θ .= broadcast((a,b)->sqrt(a)./sqrt.(b),model.likelihood.a,model.likelihood.b)
end

function sample_local!(model::VGP{T,<:LaplaceLikelihood,<:GibbsSampling}) where {T}
    model.likelihood.b .= broadcast((μ,y,β)->rand.(GeneralizedInverseGaussian.(1/β^2,abs2.(μ-y),0.5)),model.μ,model.inference.y,model.likelihood.β)
    model.likelihood.θ .= broadcast(b->1.0./b,model.likelihood.b)
    return nothing
end

@inline ∇E_μ(model::AbstractGP{T,<:LaplaceLikelihood,<:GibbsorVI}) where {T} =  hadamard.(model.likelihood.θ,model.inference.y)
@inline ∇E_μ(model::AbstractGP{T,<:LaplaceLikelihood,<:GibbsorVI},i::Int) where {T} =  model.likelihood.θ[i].*model.inference.y[i]

@inline ∇E_Σ(model::AbstractGP{T,<:LaplaceLikelihood,<:GibbsorVI}) where {T} = 0.5*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:LaplaceLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5*model.likelihood.θ[i]

## ELBO ##
function ELBO(model::AbstractGP{T,<:LaplaceLikelihood,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - GIGExpKL(model) - GaussianKL(model)
end

function expecLogLikelihood(model::VGP{T,LaplaceLikelihood{T},AnalyticVI{T}}) where {T}
    tot = -0.5*model.nLatent*model.nSample*log(twoπ)
    tot += 0.5.*sum(broadcast(θ->sum(log.(θ)),model.likelihood.θ))
    tot += -0.5.*sum(broadcast((θ,Σ,μ,y)->dot(θ,(Σ+abs2.(μ)-2.0*μ.*y-abs2.(y))),model.likelihood.θ,diag.(model.Σ),model.μ,model.y))
    return tot
end

function expecLogLikelihood(model::SVGP{T,LaplaceLikelihood{T},AnalyticVI{T}}) where {T}
    tot = -0.5*model.nLatent*model.inference.nSamplesUsed*log(twoπ)
    tot += 0.5.*sum(broadcast(θ->sum(log.(θ)),model.likelihood.θ))
    tot += -0.5.*sum(broadcast((θ,K̃,κ,Σ,μ,y)->dot(θ,(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)-2.0*(κ*μ).*y-abs2.(y))),model.likelihood.θ,model.K̃,model.κ,model.Σ,model.μ,model.inference.y))
    return model.inference.ρ*tot
end

function GIGExpKL(model::AbstractGP{T,<:LaplaceLikelihood}) where {T}
    GIGEntropy(model)-expecExponentialGIG(model)
end

function GIGEntropy(model::AbstractGP{T,<:LaplaceLikelihood}) where {T}
    model.inference.ρ*sum(broadcast(GIGEntropy,model.likelihood.a,model.likelihood.b,model.likelihood.p))
end

function expecExponentialGIG(model::AbstractGP{T,<:LaplaceLikelihood}) where {T}
    sum(broadcast((β,a,b)->sum(-log(2*β^2).-0.5*(a.*sqrt.(b)+b.*sqrt(a))./(a.*b*β^2)),model.likelihood.β,model.likelihood.a,model.likelihood.b))
end

## PDF and Log PDF Gradients ##

function grad_quad(likelihood::LaplaceLikelihood{T},y::Real,μ::Real,σ²::Real,inference::Inference) where {T<:Real}
    nodes = inference.nodes*sqrt2*sqrt(σ²) .+ μ
    Edlogpdf = dot(inference.weights,grad_log_pdf.(likelihood,y,nodes))
    Ed²logpdf =  (1/sqrt(twoπ*σ²))/(likelihood.β[1]^2)
    return -Edlogpdf::T, Ed²logpdf::T
end

@inline grad_log_pdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real} = sign(y-f)./l.β[1]

function gradpdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real}
    grad_log_pdf_μ(l,y,f)*pdf(l,y,f)
end

function hessiandiagpdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real}
    pdf(l,y,f)/(l.β[1]^2)
end

@inline hessian_log_pdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real} = zero(T)


text:
module UnregisteredWithProject
export f
f(x) = x

end

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
# Pre-image Computation:
# ======================

abstract InferenceAlgorithm

abstract PartitionAlgorithm <: InferenceAlgorithm
abstract SamplingAlgorithm <: InferenceAlgorithm
abstract MCMCAlgorithm <: SamplingAlgorithm

"A stop function used as dummy stopping criteria in while loops"
neverstop(_...) = true

"Is this box small (enough to be accepted)"
function issmall(box::Boxes, precision::Float64)
  for dim in dims(box)
    (measure(box[dim]) > precision)&& return false
  end
  return true
end

include("refinement/partition.jl")
include("refinement/chain.jl")
include("refinement/mcmc.jl")

# Specific algorithms
include("refinement/ams.jl")
include("refinement/bfs.jl")


text:
# functions for HDF5
"""
    save_qmsim_data_hdf5(pop,hdf5file)

Save a genome structure to a file (HDF5), which will replace the old file.
"""
function save_qmsim_data_hdf5(pop,hdf5file)
   map = pop.map
   h5open(hdf5file,"w") do fid
      # map
      create_group(fid, "map")
      write(fid, "map/nchr", map.nchr)
      write(fid, "map/totalSNP", map.totalSNP)
      write(fid, "map/totalQTL", map.totalQTL)
      create_group(fid, "map/chr")
      for i=1:map.nchr
         n = @sprintf("%d",i)
         create_group(fid, "map/chr/"*n)
         write(fid, "map/chr/"*n*"/nLoci", map.chr[i].nLoci)
         write(fid, "map/chr/"*n*"/nSNP", map.chr[i].nSNP)
         write(fid, "map/chr/"*n*"/nQTL", map.chr[i].nQTL)
         write(fid, "map/chr/"*n*"/seqQTL", map.chr[i].seqQTL)
         write(fid, "map/chr/"*n*"/pos", map.chr[i].pos)
         write(fid, "map/chr/"*n*"/maxAllele", map.chr[i].maxAllele)
         write(fid, "map/chr/"*n*"/naQTL", map.chr[i].naQTL)
         write(fid, "map/chr/"*n*"/effQTL", map.chr[i].effQTL)
      end

      # genotypes
      create_group(fid, "genome")
      gr = fid["genome"]
      packedsz = get_packed_genome_size(pop.map)
      write(gr, "packedsz", packedsz)
      packed = zeros(Int8,packedsz)
      nind = length(pop.individual)
      write(gr, "nind", nind)
      dset = create_dataset(gr, "genotypes", datatype(Int8), ((packedsz,nind),(packedsz,-1)), chunk=(packedsz,1))
      for k=1:nind
         pack_genome!(pop.map, pop.individual[k],packed)
         dset[1:packedsz,k] = packed
      end
   end
   return nothing
end

"""
    create_qmsim_map_hdf5(map,hdf5file)

Create an HDF5 file with `map`.
A dataset for genotypes will be empty
"""
function create_qmsim_map_hdf5(map,hdf5file)
   pop = QMSimPopulationGenome(map,Vector{QMSimIndividualGenome}(undef,0))
   save_qmsim_data_hdf5(pop,hdf5file)
end

"""
    add_qmsim_individual_hdf5(map,hdf5file)

Adds a genotyped individual to pre-created HDF5 file.
"""
function add_qmsim_individual_hdf5(map,hdf5file,individual)
   h5open(hdf5file,"r+") do fid
      gr = fid["genome"]
      packedsz = get_packed_genome_size(map)
      nind = read(gr, "nind")
      dset = gr["genotypes"]
      (dims,max_dims) = HDF5.get_extent_dims(dset)
      if nind>=dims[2]
         # doubled size
         HDF5.set_extent_dims(dset, (packedsz,max(10,dims[2])*2))
      end
      packed = zeros(Int8,packedsz)
      pack_genome!(map, individual, packed)
      dset[1:packedsz,nind+1] = packed
      q = gr["nind"]
      write(gr["nind"], nind+1)
   end
end

"""
    pop = read_qmsim_map_hdf5(hdf5file)

Read a map from `hdf5file`, which has been created by `save_qmsim_data`.
"""
function read_qmsim_map_hdf5(hdf5file)
   nchr = 0
   totalSNP = 0
   totalQTL = 0
   h5open(hdf5file,"r") do fid
      gr = fid["map"]
      nchr = read(gr,"nchr")
      totalSNP = read(gr,"totalSNP")
      totalQTL = read(gr,"totalQTL")
   end
   chr = Vector{QMSimChromosomeMap}(undef,nchr)
   h5open(hdf5file,"r") do fid
      gr = fid["map"]
      for i=1:nchr
         n = @sprintf("%d",i)
         cgr = fid["map/chr/"*n]
         nLoci = read(cgr,"nLoci")
         nSNP = read(cgr,"nSNP")
         nQTL = read(cgr,"nQTL")
         seqQTL = read(cgr,"seqQTL")
         pos = read(cgr,"pos")
         maxAllele = read(cgr,"maxAllele")
         naQTL = read(cgr,"naQTL")
         effQTL = read(cgr,"effQTL")
         chr[i] = QMSimChromosomeMap(nLoci,nSNP,nQTL,seqQTL,pos,maxAllele,naQTL,effQTL)
      end
   end
   return QMSimMap(nchr,totalSNP,totalQTL,chr)
end

"""
    genotype = read_qmsim_individual_hdf5(map,hdf5file,id)

Read a genotype of particular individual (`id`) from `hdf5file`, which has been created by `save_qmsim_data`.
"""
function read_qmsim_individual_hdf5(gmap,hdf5file,idx)
   chromosome_set = read_qmsim_chromosome_set_hdf5(gmap,hdf5file,idx)
   tbv = get_true_breeding_value(gmap,chromosome_set)
   return QMSimIndividualGenome(tbv,chromosome_set)
end

"""
    genotype = read_qmsim_chromosome_set_hdf5(map,hdf5file,id)

Read a set of chromosomes of particular individual (`id`) from `hdf5file`, which has been created by `save_qmsim_data`.
"""
function read_qmsim_chromosome_set_hdf5(gmap,hdf5file,idx)
   chromosome_set = Vector{QMSimChromosomeGenome}()
   h5open(hdf5file,"r") do fid
      gr = fid["genome"]
      #packedsz = read(gr, "packedsz")
      #nind = read(gr, "nind")
      dset = gr["genotypes"]
      packed = dset[1:end,idx]
      chromosome_set = unpack_genome(gmap,packed)
   end
   return chromosome_set
end

"""
    genome = read_qmsim_all_genotypes_hdf5(map,hdf5file)

Read all genotypes from `hdf5file`, which has been created by `save_qmsim_data`.
"""
function read_qmsim_all_genotypes_hdf5(gmap,hdf5file)
   # number of genotypes in the file
   nind = 0
   h5open(hdf5file,"r") do fid
      gr = fid["genome"]
      nind = read(gr, "nind")
   end

   # read the genotypes
   individual = Vector{QMSimIndividualGenome}(undef, nind)
   h5open(hdf5file,"r") do fid
      gr = fid["genome"]
      packedsz = read(gr, "packedsz")
      nind = read(gr, "nind")
      dset = gr["genotypes"]
      packed = zeros(Int8,packedsz)
      for i=1:nind
         packed .= dset[1:end,i]
         chromosome_set = unpack_genome(gmap,packed)
         tbv = get_true_breeding_value(gmap,chromosome_set)
         individual[i] = QMSimIndividualGenome(tbv, chromosome_set)
      end
   end
   return individual
end

"""
    data = read_qmsim_data_hdf5(hdf5file)

Read the enture data from `hdf5file`, which has been created by `save_qmsim_data`.
"""
function read_qmsim_data_hdf5(hdf5file)
   gmap = read_qmsim_map_hdf5(hdf5file)
   individual = QMSimFiles.read_qmsim_all_genotypes_hdf5(gmap,hdf5file)
   return QMSimPopulationGenome(gmap,individual)
end

"""
    (nind,nloci) = get_qmsim_genotype_size(gmap,hdf5file)

Get the number of individuals and loci.
"""
function get_qmsim_genotype_size(hdf5file)
   nind = 0
   ngen = 0
   nloci = 0
   h5open(hdf5file,"r") do fid
      gr = fid["genome"]
      packedsz = read(gr, "packedsz")
      nind = read(gr, "nind")
      (ngen,nloci) = size(gr["genotypes"])
   end
   return (nind,nloci)
end

function get_packed_genome_size(gmap)
   return 2*(gmap.totalSNP + gmap.totalQTL)
end

function pack_genome(gmap::QMSimMap,g::QMSimIndividualGenome)
   sz = get_packed_genome_size(gmap)
   packed = zeros(Int8,sz)
   pack_genome!(gmap,g,packed)
   return packed
end

function pack_genome(gmap::QMSimMap,chromosome_set::Vector{QMSimChromosomeGenome})
   sz = get_packed_genome_size(gmap)
   packed = zeros(Int8,sz)
   pack_genome!(gmap,chromosome_set,packed)
   return packed
end

function pack_genome!(gmap::QMSimMap,g::QMSimIndividualGenome,packed::Vector{Int8})
   pack_genome!(gmap,g.chr,packed)
   return nothing
end

function pack_genome!(gmap::QMSimMap,chr::Vector{QMSimChromosomeGenome},packed::Vector{Int8})
   fst = 1
   lst = 0
   for i=1:gmap.nchr
      fst = lst + 1
      lst = fst + gmap.chr[i].nLoci - 1
      packed[fst:lst] = chr[i].gp[1:end]
      fst = lst + 1
      lst = fst + gmap.chr[i].nLoci - 1
      packed[fst:lst] = chr[i].gm[1:end]
   end
   return nothing
end

function unpack_genome(gmap::QMSimMap,packed::Vector{Int8})
   chromosome_set = generate_chromosome_set(gmap)
   fst = 1
   lst = 0
   for i=1:gmap.nchr
      fst = lst + 1
      lst = fst + gmap.chr[i].nLoci - 1
      chromosome_set[i].gp[1:end] = packed[fst:lst]
      fst = lst + 1
      lst = fst + gmap.chr[i].nLoci - 1
      chromosome_set[i].gm[1:end] = packed[fst:lst]
   end
   return chromosome_set
end


text:
using Base: String, Int
import Base.+

########## Command ##########

@enum Direction forward down up

struct Command{Direction}

    coeff::Int

end

function Command(entry::String)

    direction_string, coeff_string = split(entry)

    direction_string == "forward" ? direction = forward :
    direction_string == "down" ? direction = down :
    direction_string == "up" ? direction = up :
    error("wrong command direction")

    Command{direction}(parse(Int, coeff_string))

end

coeff(c::Command) = c.coeff

########## Coordinates ##########

struct Coordinates

    x::Int
    y::Int

end

Coordinates() = Coordinates(0, 0)

x(c::Coordinates) = c.x
y(c::Coordinates) = c.y

+(a::Coordinates, b::Coordinates) = Coordinates(x(a) + x(b), y(a) + y(b))

########## Submarines_type 1 ##########

mutable struct Submarine_type1

    position::Coordinates

end

Submarine_type1() = Submarine_type1(Coordinates())

position(s::Submarine_type1) = s.position

instruct!(submarine::Submarine_type1, command::Command{forward}) = (submarine.position += Coordinates(coeff(command), 0))
instruct!(submarine::Submarine_type1, command::Command{down}) = (submarine.position += Coordinates(0, coeff(command)))
instruct!(submarine::Submarine_type1, command::Command{up}) = (submarine.position += Coordinates(0, -coeff(command)))

########## Submarines_type 2 ##########

mutable struct Submarine_type2

    position::Coordinates
    aim::Int

end

Submarine_type2() = Submarine_type2(Coordinates(), 0)

position(s::Submarine_type2) = s.position
aim(s::Submarine_type2) = s.aim

instruct!(submarine::Submarine_type2, command::Command{forward}) = (submarine.position += Coordinates(coeff(command), aim(submarine) * coeff(command)))
instruct!(submarine::Submarine_type2, command::Command{down}) = (submarine.aim += coeff(command))
instruct!(submarine::Submarine_type2, command::Command{up}) = (submarine.aim -= coeff(command))

########## Script ##########

function main()

    # open file and store content into a string
    input = open("inputs//day2.txt", "r") do io
        read(io, String)
    end

    # split the string into a vector and parse it as Commands
    data = Command.(String.(split(input, "\r\n")))

    # instantiate submarine
    submarine1 = Submarine_type1()
    submarine2 = Submarine_type2()

    # send all instructions
    for command in data
        instruct!(submarine1, command)
        instruct!(submarine2, command)
    end

    println("part1: ", x(position(submarine1)) * y(position(submarine1)))
    println("part2: ", x(position(submarine2)) * y(position(submarine2)))

end

main()


text:
module C

using Compat
export PetscInt
using MPI

# These three constants where originally in defs.jl.
# I had to changed them here because I swapped the
# order of the next two includes
const MPI_COMM_SELF=MPI.COMM_SELF
const MPI_Comm=MPI.Comm
const comm_type=MPI.Comm

include("libPETSc_commonRealDouble.jl")
include("defs.jl")
if have_petsc[1]
  include("PETScRealDouble.jl")
end
if have_petsc[2]
  include("PETScRealSingle.jl")
end
if have_petsc[3]
  include("PETScComplexDouble.jl")
end
include("error.jl")
include("defs2.jl")
include("c_funcs.jl")
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Libtask_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Libtask")
JLLWrappers.@generate_main_file("Libtask", UUID("3ae2931a-708c-5973-9c38-ccf7496fb450"))
end  # module Libtask_jll


text:
"""
Retrieve the output of a screen as a data frame. The group is GENERAL by default.
The private flag is only required if the name designates as well a public screen.
See the Equity Screening page in the Bloomberg terminal (EQS <GO>).

N.B. If the calculation takes too long the function will timeout (after 2 min by default).
In that case the result has to be retrieved (to avoid issues with later requests) doing
DataFrames.DataFrame(Bloomie.process_tree(Bloomie.retrieve_response()))

#Examples
```julia
beqs("Global Automotive Parts Retailers")
beqs("Fisher Inspired Screen",group="Guru Screens")
```
"""
function beqs(name::String;group::String="General",private::Bool=false,timeout::Int=120)
    global session
    global services
    service = services["refdata"]
    request = @with_pointer request blpapi_Service_createRequest(service,request,"BeqsRequest")
    elements = blpapi_Request_elements(request)
    try
        blpapi_Element_setElementString(elements,"screenName",C_NULL,name)
        if private
            blpapi_Element_setElementString(elements,"screenType",C_NULL,"PRIVATE")
        else
            blpapi_Element_setElementString(elements,"screenType",C_NULL,"GLOBAL")
        end
        blpapi_Element_setElementString(elements,"Group",C_NULL,group)
        correlation_id = blpapi_CorrelationId(0,0)
        blpapi_Session_sendRequest(session,request,Ref(correlation_id),C_NULL,C_NULL,C_NULL,0)
    finally
        blpapi_Request_destroy(request)
    end
    DataFrames.DataFrame(process_tree(retrieve_response(timeout)))
end


text:
# logarithmic spirals

spiraldata = [
    (10,  0.05),
    (4,   0.10),
    (0.5, 0.17)]

grid = GridRect(O - (200, 0), 175, 50)
@png begin
    background("grey15")
    for aspiral in spiraldata
        @layer begin
            sethue("antiquewhite")
            setcolor(sethue("antiquewhite"))
            Luxor.translate(nextgridpoint(grid))
            setline(.5)
            spiral(first(aspiral), last(aspiral), log=true, period=10π, :stroke)
            fontsize(16)
            label(string(aspiral), :S, offset=100)
        end
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
Recipes = Dict{AbstractString, Recipe}()

macro IN(iname, iqty)
	:(r.ins[$iname] = $iqty / r.time)
end

macro OUT(iname, iqty)
	:(r.outs[$iname] = $iqty / r.time)
end

macro FOUTS(K) :(sum([@FOUT(F, $K) for F in @Recipes])) end
macro FINS(K) :(sum([@FIN(F, $K) for F in @Recipes])) end

macro FCON(R, I) :(getvalue(Factories[$R]) * get(Recipes[$R].ins, $I, 0)) end
macro FCONS(K) :(sum([@FCON(F, $K) for F in @Recipes])) end

macro Recipes() :(collect(keys(Recipes))) end
macro FIN(F, R) :(Factories[$F] * get(Recipes[$F].ins, $R, 0)) end
macro FOUT(F, R) :(Factories[$F] * get(Recipes[$F].outs, $R, 0))end
macro FOUTS(K) :(sum([@FOUT(F, $K) for F in @Recipes])) end
macro FINS(K) :(sum([@FIN(F, $K) for F in @Recipes])) end

macro minsec(n) :((60($n), ($n), 60 * $n * mins)) end
macro PerSec(R)
	return quote
		r = get(Recipes, $R, nothing)
		if r == nothing
			1
		else
			1 / Recipes[$R].time
		end
	end
end

r = Recipes["high-tech-science-pack"] = Recipe(14, 5.5)
@IN "battery" 1
@IN "processing-unit" 3
@IN "speed-module" 1
@IN "copper-cable" 30
@OUT "high-tech-science-pack" 2.8

r = Recipes["battery"] = Recipe(5, 2.5)
@IN "sulfuric-acid" 20
@IN "iron-plate" 1
@IN "copper-plate" 1
@OUT "battery" 1.2

r = Recipes["processing-unit"] = Recipe(10, 5.5)
@IN "electronic-circuit" 20
@IN "advanced-circuit" 2
@IN "sulfuric-acid" 5
@OUT "processing-unit" 1.4

r = Recipes["advanced-circuit"] = Recipe(6, 5.5)
@IN "electronic-circuit" 2
@IN "plastic-bar" 2
@IN "copper-cable" 4
@OUT "advanced-circuit" 1.4

r = Recipes["electronic-circuit"] = Recipe(0.5, 5.5)
@IN "iron-plate" 1
@IN "copper-cable" 3
@OUT "electronic-circuit" 1.4

r = Recipes["speed-module"] = Recipe(15, 8.75)
@IN "advanced-circuit" 5
@IN "electronic-circuit" 5
@OUT "speed-module" 1

r = Recipes["copper-cable"] = Recipe(0.5, 5.5)
@IN "copper-plate" 1
@OUT "copper-cable" 2.8

r = Recipes["science-pack-3"] = Recipe(12, 5.5)
@IN "advanced-circuit" 1
@IN "engine-unit" 1
@IN "electric-mining-drill" 1
@OUT "science-pack-3" 1

r = Recipes["engine-unit"] = Recipe(10, 5.5)
@IN "steel-plate" 1
@IN "iron-gear-wheel" 1
@IN "pipe" 2
@OUT "engine-unit" 1

r = Recipes["steel-plate"] = Recipe(17.5, 9.4)
@IN "iron-plate" 5
@OUT "steel-plate" 1.2

r = Recipes["iron-gear-wheel"] = Recipe(0.5, 5.5)
@IN "iron-plate" 1
@OUT "iron-gear-wheel" 1.4

r = Recipes["pipe"] = Recipe(0.5, 5.5)
@IN "iron-plate" 1
@OUT "pipe" 1.4

r = Recipes["electric-mining-drill"] = Recipe(2, 8.75)
@IN "electronic-circuit" 3
@IN "iron-gear-wheel" 5
@IN "iron-plate" 10
@OUT "electric-mining-drill" 1

r = Recipes["plastic-bar"] = Recipe(1, 2.5)
@IN "coal" 1
@IN "petroleum-gas" 20
@OUT "plastic-bar" 1.3

r = Recipes["science-pack-2"] = Recipe(6, 5.5)
@IN "inserter" 1
@IN "transport-belt" 1
@OUT "science-pack-2" 1.4

r = Recipes["inserter"] = Recipe(0.5, 8.75)
@IN "electronic-circuit" 1
@IN "iron-gear-wheel" 1
@IN "iron-plate" 1
@OUT "inserter" 1

r = Recipes["transport-belt"] = Recipe(0.5, 8.75)
@IN "iron-gear-wheel" 1
@IN "iron-plate" 1
@OUT "transport-belt" 2

r = Recipes["science-pack-1"] = Recipe(5, 5.5)
@IN "copper-plate" 1
@IN "iron-gear-wheel" 1
@OUT "science-pack-1" 1.4

r = Recipes["production-science-pack"] = Recipe(14, 5.5)
@IN "electric-engine-unit" 1
@IN "assembling-machine-1" 1
@IN "electric-furnace" 1
@OUT "production-science-pack" 2

r = Recipes["electric-engine-unit"] = Recipe(10, 5.5)
@IN "engine-unit" 1
@IN "electronic-circuit" 2
@IN "lubricant" 15
@OUT "electric-engine-unit" 1.4

r = Recipes["assembling-machine-1"] = Recipe(0.5, 8.75)
@IN "electronic-circuit" 3
@IN "iron-gear-wheel" 5
@IN "iron-plate" 9
@OUT "assembling-machine-1" 1

r = Recipes["electric-furnace"] = Recipe(5, 8.75)
@IN "steel-plate" 10
@IN "advanced-circuit" 5
@IN "stone-brick" 10
@OUT "electric-furnace" 1

r = Recipes["stone-brick"] = Recipe(3.5, 9.4)
@IN "stone" 2
@OUT "stone-brick" 1

r = Recipes["military-science-pack"] = Recipe(10, 5.5)
@IN "piercing-rounds-magazine" 1
@IN "grenade" 1
@IN "gun-turret" 1
@OUT "military-science-pack" 2

r = Recipes["piercing-rounds-magazine"] = Recipe(3, 8.75)
@IN "firearm-magazine" 1
@IN "steel-plate" 1
@IN "copper-plate" 5
@OUT "piercing-rounds-magazine" 1

r = Recipes["firearm-magazine"] = Recipe(1, 8.75)
@IN "iron-plate" 4
@OUT "firearm-magazine" 1

r = Recipes["grenade"] = Recipe(8, 8.75)
@IN "coal" 10
@IN "iron-plate" 5
@OUT "grenade" 1

r = Recipes["gun-turret"] = Recipe(8, 8.75)
@IN "copper-plate" 10
@IN "iron-gear-wheel" 10
@IN "iron-plate" 20
@OUT "gun-turret" 1

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
using DifferentialEquations
function f(du, u, p, t)
    du[1] = u[2]
    du[2] = -p
end

function condition(u, t, integrator)
    u[1]
end

function affect!(integrator)
    integrator.u[2] = -integrator.u[2]
    integrator.u[2] *= 0.86
end
cb = ContinuousCallback(condition, affect!)
tspan = (0.0, 40.0)
p = 9.8
u0 = [100.0, -6]
prob = ODEProblem(f, u0, tspan, p)
ball_sol = solve(prob, Tsit5(), callback = cb)

using Luxor
demomovie = Movie(400, 400, "bouncing ball")
function frame(scene::Scene, framenumber)
    ballradius = scene.movie.width / 8
    ground = scene.movie.height / 2 * 0.95
    height = scene.movie.height / 2
    i = rescale(framenumber, 1, scene.framerange.stop, tspan[1], tspan[end])
    yval = clamp(ball_sol(i)[1], 0.0, height)
    ypos = rescale(yval, 0, height / 2, ground - ballradius, -height / 2)
    xpos = rescale(
        framenumber,
        1,
        scene.framerange.stop,
        -scene.movie.width / 2 + 100,
        scene.movie.width / 2 - 100,
    )
    sethue("orange")
    circle(Point(xpos, ypos), ballradius, :fill)
    sethue("black")
    setline(5)
    rule(Point(0, ground))
end

destination = "simple-bouncing-ball.gif"

animate(
    demomovie,
    [
     Scene(demomovie, (s, f) -> background("white"), 1:100),
     Scene(demomovie, frame, 1:100),
    ],
    creategif = true,
    pathname = destination,
)

include("viewer.jl"); Viewer(destination)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using PSMV
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here

n = 2*10^5

A = sprand(n, n, 0.0005)
x = rand(n)
C = PSMV.MultithreadedTransMatVec(A, A')
@show @test norm(A*x - C*x) < 1e-10

C = PSMV.MultithreadedMatVec(A, 1)
@show @test norm(A*x - C*x) < 1e-10

max_stars_repo_name: camruggles/PSVM

max_issues_repo_name: camruggles/PSVM

max_forks_repo_name: camruggles/PSVM