text:
using .NestedSamplers

export NestedStatic

Base.@kwdef struct NestedStatic{B} <: AbstractNested
    nlive::Int = 500
    bounds::B=NestedSamplers.Bounds.MultiEllipsoid
    proposal::Symbol=:auto
end

"""
    sample(ns::NestedStatic, logj; nlive=400, kwargs...)
Takes a problem defined by the logj likelihood object and runs the
NestedSampler.jl nested sampler algorithm on it. This constructs
the approporate likliehood and prior transform for you if you stick to one
of the predefined models.

Returns a chain, state, names
"""
function StatsBase.sample(ns::NestedStatic, lj::Soss.ConditionalModel; kwargs...)
    lklhd, prt, tc, unflatten = _split_conditional(lj)

    sampler = Nested(dimension(tc),
                     ns.nlive,
                     bounds=ns.bounds,
                     proposal=ns.proposal)
    model = NestedModel(lklhd, prt)
    chain, state = sample(model, sampler; chain_type=Array, kwargs...)
    logz = state[:logz]
    logzerr = state[:logzerr]
    logl = state[:logl]
    stats = (logz=logz, logzerr=logzerr, logl=logl)
    return _create_tv(unflatten, @view(chain[:,1:end-1]), @view(chain[:,end]) ), stats
end


text:
output = joinpath(@__DIR__,"..","output")
mkpath(output)

@testset "import_dict: without free parameters" begin
  pmodel = Model.m
  D = @struct_2_dict(pmodel)

  sf = SingleFile(joinpath(output,"toy_model_all.jl"))
  export_dict(D, sf)

  contents = "pmodel.x = 500.0\npmodel.a.x = 500.0\npmodel.a.i = 500\npmodel.a.fb.x = 500.0\npmodel.a.fb.i = 500\n"
  open(sf.filename, "w") do io
    print(io, contents)
  end
  import_dict(D, sf, EntireStruct())
  @test D["pmodel.x"]      ≈ 500.0
  @test D["pmodel.x"]      ≈ 500.0
  @test D["pmodel.a.x"]    ≈ 500.0
  @test D["pmodel.a.i"]    ≈ 500
  @test D["pmodel.a.fb.x"] ≈ 500.0
  @test D["pmodel.a.fb.i"] ≈ 500
end

@testset "import_dict: with free parameters" begin
  pmodel = Model.m
  D = @struct_2_dict(pmodel)

  @free D["pmodel.x"] Distributions.Normal
  @free D["pmodel.a.x"]

  sf = SingleFile(joinpath(output,"toy_model_all.jl"))
  export_dict(D, sf)

  # contents = open(f->read(f, String), sf.filename)
  contents = "pmodel.a.fb.x = 500
pmodel.a.x = FreeParameter{Float64,Nothing,Nothing}(500, nothing, nothing)
pmodel.x = FreeParameter{Float64,UnionAll,Nothing}(500, Normal, nothing)
pmodel.a.i = 500
pmodel.a.fb.i = 500
"
  # Manually adjust file:
  open(sf.filename, "w") do io
    print(io, contents)
  end

  import_dict(D, sf, EntireStruct())
  @test FreeParameters.get_val(D["pmodel.x"]) ≈ 500
  @test FreeParameters.get_val(D["pmodel.a.x"]) ≈ 500
  @test D["pmodel.a.i"] ≈ 500
  @test D["pmodel.a.fb.x"] ≈ 500
  @test D["pmodel.a.fb.i"] ≈ 500
end

@testset "import_dict: with free parameters, subset" begin
  pmodel = Model.m
  D = @struct_2_dict(pmodel)

  @free D["pmodel.x"] Distributions.Normal
  @free D["pmodel.a.x"]

  sf = SingleFile(joinpath(output,"toy_model_all.jl"))
  export_dict(D, sf)

  # contents = open(f->read(f, String), sf.filename)
  contents = "pmodel.a.fb.x = 500
pmodel.a.x = FreeParameter{Float64,Nothing,Nothing}(500, nothing, nothing)
pmodel.x = FreeParameter{Float64,UnionAll,Nothing}(500, Normal, nothing)
pmodel.a.i = 500
pmodel.a.fb.i = 500
"
  # Manually adjust file:
  open(sf.filename, "w") do io
    print(io, contents)
  end

  import_dict(D, sf, FreeParametersOnly())
  @test FreeParameters.get_val(D["pmodel.x"]) ≈ 500
  @test FreeParameters.get_val(D["pmodel.a.x"]) ≈ 500
  @test !(D["pmodel.a.i"] ≈ 500)
  @test !(D["pmodel.a.fb.x"] ≈ 500)
  @test !(D["pmodel.a.fb.i"] ≈ 500)
end

#####
##### Complex model
#####

@testset "complex model: import_dict single file" begin
  pmodel = DycomsModel.model
  D = @struct_2_dict(pmodel)

  sf = SingleFile(joinpath(output,"real_model_fp.jl"))

  @free D["pmodel.turbulence.C_smag"] Distributions.Normal
  @free D["pmodel.boundarycondition.C_drag"]
  export_dict(D, sf, FreeParametersOnly())

  # contents = open(f->read(f, String), sf.filename)
  # Manually adjust file:
  contents = "pmodel.turbulence.C_smag = FreeParameter{Float32,UnionAll,Nothing}(500, Normal, nothing)\npmodel.boundarycondition.C_drag = FreeParameter{Float32,Nothing,Nothing}(500, nothing, nothing)\n"
  open(sf.filename, "w") do io
    print(io, contents)
  end

  import_dict(D, sf, FreeParametersOnly())
  @test FreeParameters.get_val(D["pmodel.turbulence.C_smag"]) ≈ 500
  @test FreeParameters.get_val(D["pmodel.boundarycondition.C_drag"]) ≈ 500
end

ext: jl


text:
using Random
using Logging, Plots
Random.seed!(7)

dt = 0.01
integration = :rk4

###################
## Parallel Park ##
###################

opts = SolverOptions()
opts.verbose = true
opts.cost_tolerance = 1e-7
opts.cost_tolerance_intermediate = 1e-6
opts.constraint_tolerance = 1e-7
opts.resolve_feasible = true
opts.outer_loop_update_type = :default
opts.R_infeasible = 10
opts.square_root = true

# Set up model, objective, and solver
model, obj = TrajectoryOptimization.Dynamics.dubinscar_parallelpark
solver = Solver(model, obj, integration=integration, dt=dt, opts=opts)
U0 = ones(solver.model.m,solver.N-1)
X0 = line_trajectory(solver)

res, stats = TrajectoryOptimization.solve(solver,U0)

solver2 = Solver(solver)
solver2.opts.al_type = :algencan
res2, stats2 = solve(solver2,U0)

plot()
plot_trajectory!(res2)

p = plot(stats["c_max"],yscale=:log10)
plot!(stats2["c_max"])
plot_vertical_lines!(p,stats2["outer_updates"])
stats2["major iterations"]
stats2


text:
module FefferyAntdComponents
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.1.0"

include("jl/'feffery'_antdaffix.jl")
include("jl/'feffery'_antdalert.jl")
include("jl/'feffery'_antdanchor.jl")
include("jl/'feffery'_antdavatar.jl")
include("jl/'feffery'_antdbacktop.jl")
include("jl/'feffery'_antdbadge.jl")
include("jl/'feffery'_antdbreadcrumb.jl")
include("jl/'feffery'_antdbutton.jl")
include("jl/'feffery'_antdcalendar.jl")
include("jl/'feffery'_antdcarousel.jl")
include("jl/'feffery'_antdcascader.jl")
include("jl/'feffery'_antdcheckbox.jl")
include("jl/'feffery'_antdcheckboxgroup.jl")
include("jl/'feffery'_antdcollapse.jl")
include("jl/'feffery'_antdcomment.jl")
include("jl/'feffery'_antdcountdown.jl")
include("jl/'feffery'_antddatepicker.jl")
include("jl/'feffery'_antddaterangepicker.jl")
include("jl/'feffery'_antddivider.jl")
include("jl/'feffery'_antddraggerupload.jl")
include("jl/'feffery'_antddrawer.jl")
include("jl/'feffery'_antddropdown.jl")
include("jl/'feffery'_antdempty.jl")
include("jl/'feffery'_antdicon.jl")
include("jl/'feffery'_antdimage.jl")
include("jl/'feffery'_antdinput.jl")
include("jl/'feffery'_antdinputnumber.jl")
include("jl/'feffery'_antdmentions.jl")
include("jl/'feffery'_antdmenu.jl")
include("jl/'feffery'_antdmessage.jl")
include("jl/'feffery'_antdmodal.jl")
include("jl/'feffery'_antdnotification.jl")
include("jl/'feffery'_antdpageheader.jl")
include("jl/'feffery'_antdpagination.jl")
include("jl/'feffery'_antdpopconfirm.jl")
include("jl/'feffery'_antdpopover.jl")
include("jl/'feffery'_antdprogress.jl")
include("jl/'feffery'_antdradiogroup.jl")
include("jl/'feffery'_antdrate.jl")
include("jl/'feffery'_antdresult.jl")
include("jl/'feffery'_antdribbon.jl")
include("jl/'feffery'_antdselect.jl")
include("jl/'feffery'_antdskeleton.jl")
include("jl/'feffery'_antdslider.jl")
include("jl/'feffery'_antdspace.jl")
include("jl/'feffery'_antdspin.jl")
include("jl/'feffery'_antdstatistic.jl")
include("jl/'feffery'_antdsteps.jl")
include("jl/'feffery'_antdswitch.jl")
include("jl/'feffery'_antdtable.jl")
include("jl/'feffery'_antdtag.jl")
include("jl/'feffery'_antdtimeline.jl")
include("jl/'feffery'_antdtimepicker.jl")
include("jl/'feffery'_antdtimerangepicker.jl")
include("jl/'feffery'_antdtooltip.jl")
include("jl/'feffery'_antdtransfer.jl")
include("jl/'feffery'_antdtree.jl")
include("jl/'feffery'_antdtreeselect.jl")
include("jl/'feffery'_antdupload.jl")
include("jl/'feffery'_antdcard.jl")
include("jl/'feffery'_antdcardgrid.jl")
include("jl/'feffery'_antddescriptionitem.jl")
include("jl/'feffery'_antddescriptions.jl")
include("jl/'feffery'_antdform.jl")
include("jl/'feffery'_antdformitem.jl")
include("jl/'feffery'_antdcol.jl")
include("jl/'feffery'_antdrow.jl")
include("jl/'feffery'_antdcontent.jl")
include("jl/'feffery'_antdfooter.jl")
include("jl/'feffery'_antdheader.jl")
include("jl/'feffery'_antdlayout.jl")
include("jl/'feffery'_antdsider.jl")
include("jl/'feffery'_antdtabpane.jl")
include("jl/'feffery'_antdtabs.jl")
include("jl/'feffery'_antdparagraph.jl")
include("jl/'feffery'_antdtext.jl")
include("jl/'feffery'_antdtitle.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "feffery_antd_components",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "feffery_antd_components.min.js",
    external_url = "https://unpkg.com/feffery_antd_components@0.1.0/feffery_antd_components/feffery_antd_components.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "feffery_antd_components.min.js.map",
    external_url = "https://unpkg.com/feffery_antd_components@0.1.0/feffery_antd_components/feffery_antd_components.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch2.dat")
include("../../Basin.jl")
pop!(ARGS)
#=
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch3.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch3.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch4.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch5.dat")
include("../../Basin.jl")
pop!(ARGS)
=#


text:
module CGAL
__precompile__(false)
include("load.jl")
end # module


text:
struct FBPFFTSquare{F<:AbstractCTFilter} <: AbstractFBP
    filter::F
end

FBPFFTSquare() = FBPFFTSquare(RamLak())
FBPFFTSquare(f::Function) = FBPFFTSquare(CTFilter(f))

_alg_method(::FBPFFTSquare) = IsIRadonSquare()
@_defiradonalgfn FBPFFTSquare fbp_fft_square

function fbp_fft_square end

@_defiradonfn fbp_fft_square begin
    κ::T = (nd - 1) / min(width(xs), width(ys))
    xxs = xs * (κ * ν)
    yys = ys * (κ * ν)
    temp_images = fill(deepcopy(tomog), 1)
    Threads.resize_nthreads!(temp_images)
    p = _iradon_progress(nϕ, progress)
    Threads.@threads for iϕ ∈ eachindex(scϕs)
        sϕ, cϕ = scϕs[iϕ]
        id = Threads.threadid()
        @inbounds img = temp_images[id]
        for ix ∈ eachindex(xs), iy ∈ eachindex(ys)
            x, y = xxs[ix], yys[iy]
            t::T = x * cϕ + y * sϕ + t₀
            if t ∈ 1..nd
                img[iy,ix] += interp(t, iϕ)
            end
        end
        next!(p)
    end
    foreach(temp_images) do x
        tomog .+= x
    end
    δt::T = π / length(scϕs) / nd * κ^2 * rows / cols
    tomog .*= δt
end

# function fbp_fft_square(
#     sinog::AbstractMatrix{T},
#     xs::AbstractVector{U1},
#     ys::AbstractVector{U2},
#     ϕs::ClosedInterval = 0..2π;
#     background::Optional{U3} = nothing,
#     filter::Optional{F} = nothing,
#     interpolation::Optional{Interp} = nothing,
#     progress::Bool = true,
# ) where {
#     T <: Real,
#     U1 <: Real,
#     U2 <: Real,
#     U3 <: Real,
#     F <: AbstractCTFilter,
#     Interp <: AbstractInterp2DOrNone,
# }
#     nd, nϕ = size(sinog)
#     cols = length(xs)
#     rows = length(ys)
#     filtered = apply(maybe(Filters.RamLak(), filter)) do f
#         filter_freq = fft(sinog, 1) .* f(T, nd, nϕ)
#         ifft(filter_freq, 1) |> real
#     end
#     interpolation = maybe(interpolate, interpolation)
#     interp = interpolation(filtered)
#     t₀::T = T(nd + 1) / 2
#     l = min(rows, cols)
#     x₀, y₀ = (cols - l) ÷ 2, (rows - l) ÷ 2
#     scϕs = sincos.(linspace(T, ORI(ϕs), nϕ))
#     indices = Vector{NTuple{2,Int}}(undef, l^2)
#     xys = Vector{NTuple{2,T}}(undef, l^2)
#     @inbounds @simd for k ∈ eachindex(xys)
#         ix, iy = (k - 1) ÷ l + 1, (k - 1) % l + 1
#         indices[k] = x₀ + ix, y₀ + iy
#         xys[k] = xs[ix], ys[iy]
#     end
#     z::T = maybe(zero(T), background)
#     temp_images = fill(fill(z, rows, cols), 1)
#     Threads.resize_nthreads!(temp_images)
#     p = _iradon_progress(nϕ, progress)
#     Threads.@threads for iϕ ∈ eachindex(scϕs)
#         sϕ, cϕ = scϕs[iϕ]
#         id = Threads.threadid()
#         @inbounds img = temp_images[id]
#         @inbounds @simd for k ∈ eachindex(xys)
#             ix, iy = indices[k]
#             x, y = xys[k]
#             # To be consistent with our conventions should be '+'.
#             t::T = x * cϕ + y * sϕ + t₀
#             if t ∈ 1..nd
#                 img[iy,ix] += interp(t, iϕ)
#             end
#         end
#         next!(p)
#     end
#     tomog = similar(sinog, rows, cols)
#     CTTomogram(sum!(tomog, temp_images))
# end


text:
using Mimi
using MimiFUND
using Plots

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
module TrialsTests

using BenchmarkTools
using Test
using Statistics

#########
# Trial #
#########

trial1 = BenchmarkTools.Trial(BenchmarkTools.Parameters(evals = 2))
push!(trial1, 2, 1, 4, 5)
push!(trial1, 21, 0, 41, 51)

trial2 = BenchmarkTools.Trial(BenchmarkTools.Parameters(time_tolerance = 0.15))
push!(trial2, 21, 0, 41, 51)
push!(trial2, 2, 1, 4, 5)

push!(trial2, 21, 0, 41, 51)
@test length(trial2) == 3
deleteat!(trial2, 3)
@test length(trial1) == length(trial2) == 2
sort!(trial2)

@test trial1.params == BenchmarkTools.Parameters(evals = trial1.params.evals)
@test trial2.params == BenchmarkTools.Parameters(time_tolerance = trial2.params.time_tolerance)
@test trial1.times == trial2.times == [2.0, 21.0]
@test trial1.gctimes == trial2.gctimes == [1.0, 0.0]
@test trial1.memory == trial2.memory ==  4
@test trial1.allocs == trial2.allocs == 5

trial2.params = trial1.params

@test trial1 == trial2

@test trial1[2] == push!(BenchmarkTools.Trial(BenchmarkTools.Parameters(evals = 2)), 21, 0, 4, 5)
@test trial1[1:end] == trial1

@test time(trial1) == time(trial2) == 2.0
@test gctime(trial1) == gctime(trial2) == 1.0
@test memory(trial1) == memory(trial2) == trial1.memory
@test allocs(trial1) == allocs(trial2) == trial1.allocs
@test params(trial1) == params(trial2) == trial1.params

# outlier trimming
trial3 = BenchmarkTools.Trial(BenchmarkTools.Parameters(), [1, 2, 3, 10, 11],
                              [1, 1, 1, 1, 1], 1, 1)

trimtrial3 = rmskew(trial3)
rmskew!(trial3)

@test mean(trimtrial3) <= median(trimtrial3)
@test trimtrial3 == trial3

#################
# TrialEstimate #
#################

randtrial = BenchmarkTools.Trial(BenchmarkTools.Parameters())

for _ in 1:40
    push!(randtrial, rand(1:20), 1, 1, 1)
end

while mean(randtrial) <= median(randtrial)
    push!(randtrial, rand(10:20), 1, 1, 1)
end

rmskew!(randtrial)

tmin = minimum(randtrial)
tmed = median(randtrial)
tmean = mean(randtrial)
tmax = maximum(randtrial)

@test time(tmin) == time(randtrial)
@test gctime(tmin) == gctime(randtrial)
@test memory(tmin) == memory(tmed) == memory(tmean) == memory(tmax) == memory(randtrial)
@test allocs(tmin) == allocs(tmed) == allocs(tmean) == allocs(tmax) == allocs(randtrial)
@test params(tmin) == params(tmed) == params(tmean) == params(tmax) == params(randtrial)

@test tmin <= tmed
@test tmean <= tmed # this should be true since we called rmoutliers!(randtrial) earlier
@test tmed <= tmax

##############
# TrialRatio #
##############

randrange = 1.0:0.01:10.0
x, y = rand(randrange), rand(randrange)

@test (ratio(x, y) == x/y) && (ratio(y, x) == y/x)
@test (ratio(x, x) == 1.0) && (ratio(y, y) == 1.0)
@test ratio(0.0, 0.0) == 1.0

ta = BenchmarkTools.TrialEstimate(BenchmarkTools.Parameters(), rand(), rand(), rand(Int), rand(Int))
tb = BenchmarkTools.TrialEstimate(BenchmarkTools.Parameters(), rand(), rand(), rand(Int), rand(Int))
tr = ratio(ta, tb)

@test time(tr) == ratio(time(ta), time(tb))
@test gctime(tr) == ratio(gctime(ta), gctime(tb))
@test memory(tr) == ratio(memory(ta), memory(tb))
@test allocs(tr) == ratio(allocs(ta), allocs(tb))
@test params(tr) == params(ta) == params(tb)

@test BenchmarkTools.gcratio(ta) == ratio(gctime(ta), time(ta))
@test BenchmarkTools.gcratio(tb) == ratio(gctime(tb), time(tb))

##################
# TrialJudgement #
##################

ta = BenchmarkTools.TrialEstimate(BenchmarkTools.Parameters(time_tolerance = 0.50, memory_tolerance = 0.50), 0.49, 0.0, 2, 1)
tb = BenchmarkTools.TrialEstimate(BenchmarkTools.Parameters(time_tolerance = 0.05, memory_tolerance = 0.05), 1.00, 0.0, 1, 1)
tr = ratio(ta, tb)
tj_ab = judge(ta, tb)
tj_r = judge(tr)

@test ratio(tj_ab) == ratio(tj_r) == tr
@test time(tj_ab) == time(tj_r) == :improvement
@test memory(tj_ab) == memory(tj_r) == :regression
@test tj_ab == tj_r

tj_ab_2 = judge(ta, tb; time_tolerance = 2.0, memory_tolerance = 2.0)
tj_r_2 = judge(tr; time_tolerance = 2.0, memory_tolerance = 2.0)

@test tj_ab_2 == tj_r_2
@test ratio(tj_ab_2) == ratio(tj_r_2)
@test time(tj_ab_2) == time(tj_r_2) == :invariant
@test memory(tj_ab_2) == memory(tj_r_2) == :invariant

@test !(isinvariant(tj_ab))
@test !(isinvariant(tj_r))
@test isinvariant(tj_ab_2)
@test isinvariant(tj_r_2)

@test isregression(tj_ab)
@test isregression(tj_r)
@test !(isregression(tj_ab_2))
@test !(isregression(tj_r_2))

@test isimprovement(tj_ab)
@test isimprovement(tj_r)
@test !(isimprovement(tj_ab_2))
@test !(isimprovement(tj_r_2))

###################
# pretty printing #
###################

@test BenchmarkTools.prettypercent(.3120123) == "31.20%"

@test BenchmarkTools.prettydiff(0.0) == "-100.00%"
@test BenchmarkTools.prettydiff(1.0) == "+0.00%"
@test BenchmarkTools.prettydiff(2.0) == "+100.00%"

@test BenchmarkTools.prettytime(999) == "999.000 ns"
@test BenchmarkTools.prettytime(1000) == "1.000 μs"
@test BenchmarkTools.prettytime(999_999) == "999.999 μs"
@test BenchmarkTools.prettytime(1_000_000) == "1.000 ms"
@test BenchmarkTools.prettytime(999_999_999) == "1000.000 ms"
@test BenchmarkTools.prettytime(1_000_000_000) == "1.000 s"

@test BenchmarkTools.prettymemory(1023) == "1023 bytes"
@test BenchmarkTools.prettymemory(1024) == "1.00 KiB"
@test BenchmarkTools.prettymemory(1048575) == "1024.00 KiB"
@test BenchmarkTools.prettymemory(1048576) == "1.00 MiB"
@test BenchmarkTools.prettymemory(1073741823) == "1024.00 MiB"
@test BenchmarkTools.prettymemory(1073741824) == "1.00 GiB"

@test sprint(show, "text/plain", ta) == sprint(show, ta; context=:compact => false) == """
BenchmarkTools.TrialEstimate:
  time:             0.490 ns
  gctime:           0.000 ns (0.00%)
  memory:           2 bytes
  allocs:           1"""

@test sprint(show, ta) == "TrialEstimate(0.490 ns)"
@test sprint(
    show, ta;
    context = IOContext(
        devnull, :compact => true, :typeinfo => BenchmarkTools.TrialEstimate)
) == "0.490 ns"

@test sprint(show, [ta, tb]) == "BenchmarkTools.TrialEstimate[0.490 ns, 1.000 ns]"

@test sprint(show, "text/plain", [ta, tb]) == """
2-element Array{BenchmarkTools.TrialEstimate,1}:
 0.490 ns
 1.000 ns"""

end # module


text:
"""
    Base.:+(a::UncertainDataset, b::UncertainDataset; n = 10000)

Addition operator for two uncertain value datasets.

To obtain the sum for a pair of uncertain values `(xᵢ ∈ a, yᵢ ∈ b)`, both `u₁` and `u₂` are
resampled `n` times each, so `sample_xᵢ = resample(x, n)` and `sample_yᵢ = resample(y, n)`.

The i-th added uncertain value is a kernel density estimate to those sums, i.e.
`added_val = UncertainValue(UnivariateKDE, sample_xᵢ .+ sample_yᵢ)`.
"""
function Base.:+(a::UncertainDataset, b::UncertainDataset; n = 10000)
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] + b[i] for i = 1:N])
end

"""
    Base.:+(a::Real, b::UncertainDataset;
        n = 10000)

Addition operator for scalars and uncertain value datasets.

To obtain the sum for a pair of uncertain values `(xᵢ ∈ a, yᵢ ∈ b)`, both `u₁` and `u₂` are
resampled `n` times each, so `sample_xᵢ = resample(x, n)` and `sample_yᵢ = resample(y, n)`.

The i-th added uncertain value is a kernel density estimate to those sums, i.e.
`added_val = UncertainValue(UnivariateKDE, sample_xᵢ .+ sample_yᵢ)`.
"""
function Base.:+(a::T, b::UncertainDataset; n = 10000) where {T <: Number}
    N = length(a)
    UncertainDataset([a + b[i] for i = 1:N])
end

"""
    Base.:+(a::Vector{Real}, b::UncertainDataset; n = 10000)

Addition operator for scalars and uncertain value datasets.

To obtain the sum for a pair of uncertain values `(xᵢ ∈ a, yᵢ ∈ b)`, both `u₁` and `u₂` are
resampled `n` times each, so `sample_xᵢ = resample(x, n)` and `sample_yᵢ = resample(y, n)`.

The i-th added uncertain value is a kernel density estimate to those sums, i.e.
`added_val = UncertainValue(UnivariateKDE, sample_xᵢ .+ sample_yᵢ)`.
"""
function Base.:+(a::Vector{T}, b::UncertainDataset; n = 10000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] + b[i] for i = 1:N])
end

"""
    Base.:+(a::UncertainDataset, b::Real; n = 10000)

Addition scalars to uncertain value datasets.

To obtain the sum for a pair of uncertain values `(xᵢ ∈ a, yᵢ ∈ b)`, both `u₁` and `u₂` are
resampled `n` times each, so `sample_xᵢ = resample(x, n)` and `sample_yᵢ = resample(y, n)`.

The i-th added uncertain value is a kernel density estimate to those sums, i.e.
`added_val = UncertainValue(UnivariateKDE, sample_xᵢ .+ sample_yᵢ)`.
"""
function Base.:+(a::UncertainDataset, b::T; n = 10000) where {T <: Number}
    N = length(a)

    UncertainDataset([a[i] + b for i = 1:N])
end

"""
    Base.:+(a::UncertainDataset, b::Real; n = 10000)

Addition scalars to uncertain value datasets.

To obtain the sum for a pair of uncertain values `(xᵢ ∈ a, yᵢ ∈ b)`, both `u₁` and `u₂` are
resampled `n` times each, so `sample_xᵢ = resample(x, n)` and `sample_yᵢ = resample(y, n)`.

The i-th added uncertain value is a kernel density estimate to those sums, i.e.
`added_val = UncertainValue(UnivariateKDE, sample_xᵢ .+ sample_yᵢ)`.
"""
function Base.:+(a::UncertainDataset, b::Vector{T}; n = 10000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainDataset([a[i] + b[i] for i = 1:N])
end


text:
@doc raw"""
    SubGradientProblem <: Problem

A structure to store information about a subgradient based optimization problem

# Fields
* `manifold` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)
* `cost` – the function $F$ to be minimized
* `subgradient` – a function returning a subgradient $\partial F$ of $F$

# Constructor

    SubGradientProblem(M, f, ∂f)

Generate the [`Problem`] for a subgradient problem, i.e. a function `f` on the
manifold `M` and a function `∂f` that returns an element from the subdifferential
at a point.
"""
struct SubGradientProblem{T,mT<:AbstractManifold,C,S} <: Problem{T}
    M::mT
    cost::C
    subgradient!!::S
    function SubGradientProblem(
        M::mT,
        cost::C,
        subgrad::S;
        evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    ) where {mT<:AbstractManifold,C,S}
        return new{typeof(evaluation),mT,C,S}(M, cost, subgrad)
    end
end
"""
    get_subgradient(p, q)
    get_subgradient!(p, X, q)

Evaluate the (sub)gradient of a [`SubGradientProblem`](@ref)` p` at the point `q`.

The evaluation is done in place of `X` for the `!`-variant.
The `T=`[`AllocatingEvaluation`](@ref) problem might still allocate memory within.
When the non-mutating variant is called with a `T=`[`MutatingEvaluation`](@ref)
memory for the result is allocated.
"""
function get_subgradient(p::SubGradientProblem{AllocatingEvaluation}, q)
    return p.subgradient!!(p.M, q)
end
function get_subgradient(p::SubGradientProblem{MutatingEvaluation}, q)
    X = zero_vector(p.M, q)
    return p.subgradient!!(p.M, X, q)
end
function get_subgradient!(p::SubGradientProblem{AllocatingEvaluation}, X, q)
    return copyto!(p.M, X, p.subgradient!!(p.M, q))
end
function get_subgradient!(p::SubGradientProblem{MutatingEvaluation}, X, q)
    return p.subgradient!!(p.M, X, q)
end

"""
    SubGradientMethodOptions <: Options
stories option values for a [`subgradient_method`](@ref) solver

# Fields
* `retraction_method` – the retration to use within
* `stepsize` – a [`Stepsize`](@ref)
* `stop` – a [`StoppingCriterion`](@ref)
* `x` – (initial or current) value the algorithm is at
* `x_optimal` – optimal value
* `∂` the current element from the possible subgradients at `x` that is used
"""
mutable struct SubGradientMethodOptions{
    TR<:AbstractRetractionMethod,TS<:Stepsize,TSC<:StoppingCriterion,P,T
} <: Options where {P,T}
    retraction_method::TR
    stepsize::TS
    stop::TSC
    x::P
    x_optimal::P
    ∂::T
    function SubGradientMethodOptions(
        M::TM,
        x::P;
        stopping_criterion::SC=StopAfterIteration(5000),
        stepsize::S=ConstantStepsize(M),
        subgrad::T=zero_vector(M, x),
        retraction_method::TR=default_retraction_method(M),
    ) where {
        TM<:AbstractManifold,
        P,
        T,
        SC<:StoppingCriterion,
        S<:Stepsize,
        TR<:AbstractRetractionMethod,
    }
        return new{TR,S,SC,P,T}(
            retraction_method, stepsize, stopping_criterion, x, deepcopy(x), subgrad
        )
    end
    @deprecate SubGradientMethodOptions(M, x, S, s, r=ExponentialRetraction()) SubGradientMethodOptions(
        M, x; stopping_criterion=S, stepsize=s, retraction_method=r
    )
end


text:
@testset "src/graph/graph.jl|output.jl" begin
    @testset "write_visualization" begin
        tmp_directory = mktempdir()
        write_visualization(network, joinpath(tmp_directory, "van_zyl"), del_files = false)
        @test isfile(joinpath(tmp_directory, "van_zyl_graph.pdf"))
        @test isfile(joinpath(tmp_directory, "van_zyl_cbar.pdf"))
        @test isfile(joinpath(tmp_directory, "van_zyl_graph_w_cb.pdf"))
        write_multi_time_viz(network_mn, solution, joinpath(tmp_directory, "van_zyl_wm"))
        @test isfile(joinpath(tmp_directory, "van_zyl_wm_graph_w_cb.pdf"))
        rm(tmp_directory, recursive = true)
    end
end


text:
using Flux: OneHotVector, onehot

import Base

struct OneHotArray{N, A<: AbstractArray{OneHotVector}} <: AbstractArray{Bool, N}
    dims::Int
    data::A
    OneHotArray(dims, data::A) where A = new{ndims(A)+1, A}(dims, data)
end

Base.size(xs::OneHotArray) = (Int64(xs.dims), size(xs.data)...)

Base.getindex(o::OneHotArray, i::Integer) = o[CartesianIndices(o)[i]]
Base.getindex(o::OneHotArray, i::Integer, j::Integer) = o.data[j][i]
Base.getindex(o::OneHotArray, i::Integer, j::Integer, I::Vararg{Int, N}) where N = o.data[j, I...][i]
Base.getindex(o::OneHotArray, ::Colon, i::Integer) = o.data[i]
Base.getindex(o::OneHotArray, ::Colon, i::AbstractArray) = OneHotArray(o.dims, o.data[i])

OneHotArray(xs::A) where A = OneHotArray{ndims(A)+1, A}(Int(xs[1].of), xs)
OneHotArray(nums::Int, xs::AbstractArray{Int}) = OneHotArray(nums, map(i->onehot(i, 1:nums), xs))

onehotarray(num::Int, xs::AbstractArray{Int}) = OneHotArray(num, indices2onehot(num, xs))

import Adapt: adapt, adapt_structure

adapt_structure(T, xs::OneHotArray) = OneHotArray(xs.dims, adapt(T, xs.data))

tofloat(::Type{F}, o::OneHotArray{N, <:AbstractArray}) where {F<:AbstractFloat,N} = Array{F}(o)

Base.convert(::Type{OneHotVector}, x::Int) = OneHotVector(x & 0xffffffff, x >> 32)
Base.convert(::Type{Int}, x::OneHotVector) = Int(x.ix)

Base.hcat(x::OneHotArray, xs::OneHotArray...) = begin
    !all(isequal(x.dims), map(o->o.dims, xs)) && throw(DimensionMismatch("OneHot dimension are not all the same"))
    OneHotArray(x.dims, vcat(x.data, map(o->o.data, xs)...))
end

Base.cat(x::OneHotArray, xs::OneHotArray...; dims::Int) = begin
    !all(isequal(x.dims), map(o->o.dims, xs)) && throw(DimensionMismatch("OneHot dimension are not all the same"))
    dims < 2 && error("OneHotArray concatenation at dimension $dims is not allow")
    OneHotArray(x.dims, cat(x.data, map(o->o.data, xs)...; dims=dims-1))
end

"turn one hot encoding to indices"
onehot2indices(xs::OneHotArray) = onehot2indices(xs.data)

#cpu onehot to indices
onehot2indices(x::AbstractArray{OneHotVector}) = map(i->Int(i.ix), x)

#cpu indices to onehot
indices2onehot(nums::Int, xs::AbstractArray{Int}) = map(i->onehot(i, 1:nums), xs)

_labelindex(num::Int, x::AbstractArray) = x .+ (num << 32)

using ZygoteRules: @adjoint, AContext, pullback
import ZygoteRules: _pullback

import Base: *
*(A::AbstractMatrix{T}, b::OneHotArray{N}) where {T,N} = gather(A, b)

_pullback(cx::AContext, ::typeof(*), A::AbstractMatrix{T}, b::OneHotArray{N}) where {T, N} = _pullback(cx, gather, A, b)


text:
using LittleScienceTools.Random
using LittleScienceTools.Ising
# using Erdos
using Base.Test

g = Graph(10,20)
Js = random_couplings(g)
@test length(Js) == nv(g)
for i=1:nv(g)
    @test length(Js[i]) == degree(g, i)
end
jdict = Dict{Tuple{Int,Int}, Int}()
for i=1:nv(g)
    for (k,j) in enumerate(neighbors(g,i))
        if i < j
            jdict[(i,j)] = Js[i][k]
        end
    end
end
for i=1:nv(g)
    for (k,j) in enumerate(neighbors(g,i))
        if i > j
            @test Js[i][k] == jdict[(j,i)]
        end
    end
end

include("mincut.jl")
include("tau-eo.jl")

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
module Beams

## Package imports
#==========================================================================================#

using SparseArrays
using StaticArrays
using LinearAlgebra
using SplitApplyCombine

abstract type AbstractBeam end

struct Beam{T <: AbstractBeam}
    section :: Vector{T}
end

## Material definition
#==========================================================================================#

abstract type AbstractMaterial end

struct Material{T} <: AbstractMaterial
    elastic_modulus :: T
    shear_modulus   :: T
    yield_stress    :: T
    density         :: T
end

"""
    Material(E, J, σ_max, ρ)
    Material(; elastic_modulus, shear_modulus, yield_stress, density)

Define a `Material` with positive, real elastic modulus ``E``, shear modulus ``G``, yield stress ``σ_{\\max}``, and density ``ρ``.

The default assignments for the named variables are set to the properties of aluminium?

# Arguments
- `elastic_modulus :: Real = 85e9`: Elastic modulus.
- `shear_modulus :: Real = 25e9`: Shear modulus.
- `yield_stress :: Real = 350e6`: Yield stress.
- `density :: Real = 1.6e3`: Density.
"""
function Material(E, G, σ_max, ρ)
    T = promote_type(eltype(E), eltype(G), eltype(σ_max), eltype(ρ))
    @assert E > 0. && G > 0. && σ_max > 0. && ρ > 0. "All quantities must be positive."
    Material{T}(E, G, σ_max, ρ)
end

Material(; elastic_modulus = 85e9, shear_modulus = 25e9, yield_stress = 350e6, density = 1.6e3) = Material(elastic_modulus, shear_modulus, yield_stress, density)

elastic_modulus(mat :: Material) = mat.elastic_modulus
shear_modulus(mat :: Material)   = mat.shear_modulus
yield_stress(mat :: Material)    = mat.yield_stress
density(mat :: Material)         = mat.density

## Tube definition
#==========================================================================================#

"""
    Tube(material :: Material, length, radius, thickness)

Define a hollow tube of fixed radius with a given material, length, and thickness.
"""
struct Tube{T <: Real} <: AbstractBeam
    material  :: Material{T}
    length    :: T
    radius    :: T
    thickness :: T
end

function Tube(material :: Material, length, radius, thickness)
    T = promote_type(eltype(length), eltype(radius), eltype(thickness))
    @assert length > 0. && radius > 0. && thickness > 0. "Length, outer radius and thickness must be positive."
    Material{T}(material, length, radius, thickness)
end

Base.length(tube :: Tube) = tube.length
material(tube :: Tube)    = tube.material
radius(tube :: Tube)      = tube.radius
thickness(tube :: Tube)   = tube.thickness

radii(tube :: Tube) = radius(tube) - thickness(tube), radius(tube)

function area(tube :: Tube)
    r1, r2 = radii(tube)
    π * (r2^2 - r1^2)
end

function moment_of_inertia(tube :: Tube)
    r1, r2 = radii(tube)
    π/4 * (r2^4 - r1^4)
end

function polar_moment_of_inertia(tube :: Tube)
    r1, r2 = radii(tube)
    π/2 * (r2^4 - r1^4)
end

## Stress calculations (NEEDS CHECKING)
#==========================================================================================#

principal_stress(E, L, R, dx, dθ_yz) = E * (dx / L + R * dθ_yz / L)
torsional_stress(G, L, R, dθ_x)    = G * R * dθ_x / L
von_mises_stress(σ_xx, σ_xt)      = √(σ_xx^2 + 3σ_xt^2)

function von_mises_stress(tube :: Tube, ds, θs)
    E = (elastic_modulus ∘ material)(tube)
    G = (shear_modulus ∘ material)(tube)
    R = radius(tube)
    L = length(tube)

    dx     = ds[1]
    dθ_yz  = norm(θs[2:end])
    σ_xx_1 = principal_stress(E, L, R,  dx, dθ_yz)
    σ_xx_2 = principal_stress(E, L, R, -dx, dθ_yz)
    σ_xt   = torsional_stress(G, L, R, θs[1])

    von_mises_stress.(SVector(σ_xx_1, σ_xx_2), σ_xt)
end

## Stiffness matrices
#==========================================================================================#

# E * I / L^3 .* [ 12 .* [ 1. -1. ; -1.  1.]  6 .* [-1. -1.; 1. 1.] ;
#                   6 .* [-1.  1. ; -1.  1.]  2 .* [ 2.  1.; 1. 2.] ]

## Coefficient matrices
J_coeffs(A, B, L) = A * B / L .* [ 1 -1; -1  1 ]

k_coeffs(k0, k1, k2) = @SMatrix [  k0 -k1 -k0 -k1 ;
                                  -k1 2k2  k1  k2 ;
                                  -k0  k1  k0  k1 ;
                                  -k1  k2  k1 2k2 ]

Iyy_coeffs(E, I, L) = let k0 = 12, k1 =  6L, k2 = 2L^2; E * I / L^3 * k_coeffs(k0, k1, k2) end
Izz_coeffs(E, I, L) = let k0 = 12, k1 = -6L, k2 = 2L^2; E * I / L^3 * k_coeffs(k0, k1, k2) end

## Composite deflection stiffness matrix of adjacent beam elements
function bending_stiffness_matrix(Es, Is, Ls, direction = :z)
    n = length(Es)
    mat = @MMatrix zeros(2(n+1), 2(n+1))
    block_size = CartesianIndices((4,4))

    for (c_ind, E, I, L) in zip(CartesianIndex.(0:2:2(n-1), 0:2:2(n-1)), Es, Is, Ls)
        coeffs = ifelse(direction == :z, Izz_coeffs(E, I, L), Iyy_coeffs(E, I, L))
        @. mat[c_ind + block_size] += coeffs
    end

    sparse(mat)
end

## Composite torsional stiffness matrix for adjacent beam elements
function axial_stiffness_matrix(Es, Is, Ls)
    n = length(Es)
    mat = @MMatrix zeros(n+1, n+1)
    block_size = CartesianIndices((2,2))

    for (c_ind, E, I, L) in zip(CartesianIndex.(0:n-1, 0:n-1), Es, Is, Ls)
        coeffs = J_coeffs(E, I, L)
        @. mat[c_ind + block_size] += coeffs
    end

    sparse(mat)
end

"""
    tube_stiffness_matrix(E, G, A, Iyy, Izz, J, L)
    tube_stiffness_matrix(E, G, A, Iyy, Izz, J, L, num)
    tube_stiffness_matrix(tube :: Tube)
    tube_stiffness_matrix(tube :: Vector{Tube})
    tube_stiffness_matrix(x :: Matrix{Real})

Generate the stiffness matrix using the properties of a tube.

The required properties are the elastic modulus ``E``, shear modulus ``G``, area ``A``, moments of inertia about the ``y-`` and ``z-`` axes ``I_{yy}, I_{zz}``, polar moment of inertia ``J``, length ``L``. A composite stiffness matrix is generated with a specified number of elements.
"""
function tube_stiffness_matrix end

## Full stiffness matrix for one element
tube_stiffness_matrix(E, G, A, Iy, Iz, J, L) = blockdiag((sparse ∘ Iyy_coeffs)(E, Iy, L), (sparse ∘ Izz_coeffs)(E, Iz, L), (sparse ∘ J_coeffs)(E, A, L), (sparse ∘ J_coeffs)(G, J, L))

## Composite stiffness matrix for adjacent tube finite elements
function tube_stiffness_matrix(Es :: Vector{T}, Gs :: Vector{T}, As :: Vector{T}, Iys :: Vector{T}, Izs :: Vector{T}, Js :: Vector{T}, Ls :: Vector{T}) where T <: Real
    # Check if sizes match
    @assert length(Es) == (length ∘ zip)(Gs, Iys, Izs, Js, Ls) "Lengths of coefficients must be the same as the dimension."

    Izs = bending_stiffness_matrix(Es, Izs, Ls, :z)
    Iys = bending_stiffness_matrix(Es, Iys, Ls, :y)
    As  = axial_stiffness_matrix(Es, As, Ls)
    Js  = axial_stiffness_matrix(Gs, Js, Ls)

    blockdiag(Iys, Iys, As, Js)
end

tube_stiffness_matrix(E, G, A, Iy, Iz, J, L, num :: Integer) = tube_stiffness_matrix(fill(E, num), fill(G, num), fill(A, num), fill(Iy, num), fill(Iz, num), fill(J, num), fill(L / num, num))

tube_stiffness_matrix(tube :: Tube) = tube_stiffness_matrix(elastic_modulus(material(tube)), shear_modulus(material(tube)), area(tube), moment_of_inertia(tube), moment_of_inertia(tube), polar_moment_of_inertia(tube), length(tube))

function tube_stiffness_matrix(x :: Matrix{<: Real})
    @assert size(x)[2] == 7 "Input must have 7 columns."
    @views tube_stiffness_matrix(x[:,1], x[:,2], x[:,3], x[:,4], x[:,5], x[:,6], x[:,7])
end

tube_stiffness_matrix(tubes :: Vector{<: Tube}) = tube_stiffness_matrix((elastic_modulus ∘ material).(tubes), (shear_modulus ∘ material).(tubes), area.(tubes), moment_of_inertia.(tubes), moment_of_inertia.(tubes), polar_moment_of_inertia.(tubes), length.(tubes))

function sparse_stiffness_matrix(num_Ks, D_start, D_end, D_mid, D_12, D_21)
    # Build sparse matrix
    stiffness = spzeros(6 * (num_Ks + 2), 6 * (num_Ks + 2))
    @views stiffness[7:12,7:12]           = D_start
    @views stiffness[end-5:end,end-5:end] = D_end

    for m in 1:num_Ks
        mid_inds = 6(m+1)+1:6(m+1)+6
        off_inds = 6(m)+1:6(m)+6
        if m < num_Ks
            @views stiffness[mid_inds,mid_inds] = D_mid[:,:,m]
        end
        @views stiffness[off_inds,mid_inds] = D_12[:,:,m]
        @views stiffness[mid_inds,off_inds] = D_21[:,:,m]
    end

    stiffness
end

function build_stiffness_matrix(D, constraint_indices)
    # Number of elements
    num_Ks = @views length(D[1,1,:])

    # First element
    D_start = @views D[1:6,1:6,1]

    # Last element
    D_end   = @views D[7:end,7:end,end]

    # Summing adjacent elements corresponding to diagonals
    D_mid   = @views D[7:end,7:end,1:end-1] + D[1:6,1:6,2:end]

    # Upper and lower diagonals
    D_12    = @views D[1:6,7:end,:]
    D_21    = @views D[7:end,1:6,:]

    stiffness = sparse_stiffness_matrix(num_Ks, D_start, D_end, D_mid, D_12, D_21)

    # Fixed boundary condition by constraining the locations.
    cons = 6 * (constraint_indices .+ 1)
    arr = 1:6

    col = arr
    row = combinedimsview(map(con -> con .+ arr, cons))

    @views stiffness[CartesianIndex.(col, row)] .= 1e9
    @views stiffness[CartesianIndex.(row, col)] .= 1e9

    stiffness
end

## Cantilever setup
function solve_cantilever_beam(Ks, loads, constraint_indices)
    # Create the stiffness matrix from the array of individual stiffnesses
    # Also specifies the constraint location
    K = build_stiffness_matrix(Ks, constraint_indices)

    # Build force vector with constraint
    f = [ zeros(6); vec(loads) ]

    # Solve FEM system
    x = K \ f

    # Throw away the junk values for the constraint
    @views reshape(x[7:end], 6, length(Ks[1,1,:]) + 1)
end

end


text:
# A lot of this module is adapted from Tracker.jl and ReverseDiff.jl
# ReverseDiff.jl is not actively developed but it would be nice to move the code in this
# module to ReverseDiff at some point.

###############
## logsumexp ##
###############

logsumexp(x::TrackedArray; dims=:) = track(logsumexp, x, dims = dims)
@grad function logsumexp(x::AbstractArray; dims)
    x_value = value(x)
    lse = logsumexp(x_value; dims=dims)
    return lse, Δ -> (Δ .* exp.(x_value .- lse),)
end

############
## linalg ##
############

function LinearAlgebra.cholesky(A::Symmetric{<:Any, <:TrackedMatrix}; check=true)
    uplo = A.uplo == 'U' ? (:U) : (:L)
    factors, info = symm_turing_chol(parent(A), check, uplo)
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end
function LinearAlgebra.cholesky(A::TrackedMatrix; check=true)
    factors, info = turing_chol(A, check)
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end

function symm_turing_chol(x::TrackedArray{V,D}, check, uplo) where {V,D}
    tp = tape(x)
    x_value = value(x)
    (factors,info), back = DistributionsAD.symm_turing_chol_back(x_value, check, uplo)
    C = Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
    out = track(C.factors, D, tp)
    record!(tp, SpecialInstruction, symm_turing_chol, (x, check, uplo), out, (back, issuccess(C)))
    return out, C.info
end
function turing_chol(x::TrackedArray{V,D}, check) where {V,D}
    tp = tape(x)
    x_value = value(x)
    (factors,info), back = DistributionsAD.turing_chol_back(x_value, check)
    C = Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
    out = track(C.factors, D, tp)
    record!(tp, SpecialInstruction, turing_chol, (x, check), out, (back, issuccess(C)))
    return out, C.info
end

for f in (:turing_chol, :symm_turing_chol)
    @eval begin
        @noinline function ReverseDiff.special_reverse_exec!(
            instruction::SpecialInstruction{typeof($f)},
        )
            output = instruction.output
            instruction.cache[2] || throw(PosDefException(C.info))
            input = instruction.input
            input_deriv = deriv(input[1])
            P = instruction.cache[1]
            input_deriv .+= P((factors = deriv(output),))[1]
            unseed!(output)
            return nothing
        end
    end
end

@noinline function ReverseDiff.special_forward_exec!(
    instruction::SpecialInstruction{typeof(turing_chol)},
)
    output, input = instruction.output, instruction.input
    factors = turing_chol(value.(input)...)[1]
    value!(output, factors)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(
    instruction::SpecialInstruction{typeof(symm_turing_chol)},
)
    output, input = instruction.output, instruction.input
    factors = symm_turing_chol(value.(input)...)[1]
    value!(output, factors)
    return nothing
end


text:
using PyCall, BenchmarkTools, DataStructures
using PyCall: PyArray_Info

results = OrderedDict{String,Any}()

let
    np = pyimport("numpy")
    nprand = np["random"]["rand"]
    # nparray_pyo(x) = pycall(np["array"], PyObject, x)
    # pytestarray(sz::Int...) = pycall(np["reshape"], PyObject, nparray_pyo(1:prod(sz)), sz)

    # no convert baseline
    nprand_pyo(sz...)   = pycall(nprand, PyObject, sz...)

    for arr_size in [(2,2), (100,100)]
        pyo_arr = nprand_pyo(arr_size...)
        results["nprand_pyo$arr_size"] = @benchmark $nprand_pyo($arr_size...)
        println("nprand_pyo $arr_size:\n"); display(results["nprand_pyo$arr_size"])
        println("--------------------------------------------------")

        results["convert_pyarr$arr_size"] = @benchmark $convert(PyArray, $pyo_arr)
        println("convert_pyarr $arr_size:\n"); display(results["convert_pyarr$arr_size"])
        println("--------------------------------------------------")

        results["PyArray_Info$arr_size"] = @benchmark $PyArray_Info($pyo_arr)
        println("PyArray_Info $arr_size:\n"); display(results["PyArray_Info$arr_size"])
        println("--------------------------------------------------")

        results["convert_pyarrbuf$arr_size"] = @benchmark $PyArray($pyo_arr)
        println("convert_pyarrbuf $arr_size:\n"); display(results["convert_pyarrbuf$arr_size"])
        println("--------------------------------------------------")

        results["convert_arr$arr_size"] = @benchmark convert(Array, $pyo_arr)
        println("convert_arr $arr_size:\n"); display(results["convert_arr$arr_size"])
        println("--------------------------------------------------")

        results["convert_arrbuf$arr_size"] = @benchmark $NoCopyArray($pyo_arr)
        println("convert_arrbuf $arr_size:\n"); display(results["convert_arrbuf$arr_size"])
        println("--------------------------------------------------")

        pyarr = convert(PyArray, pyo_arr)
        results["setdata!$arr_size"] = @benchmark $setdata!($pyarr, $pyo_arr)
        println("setdata!:\n"); display(results["setdata!$arr_size"])
        println("--------------------------------------------------")

        pyarr = convert(PyArray, pyo_arr)
        pybuf=PyBuffer()
        results["setdata! bufprealloc$arr_size"] =
            @benchmark $setdata!($pyarr, $pyo_arr, $pybuf)
        println("setdata! bufprealloc:\n"); display(results["setdata! bufprealloc$arr_size"])
        println("--------------------------------------------------")
    end
end
println()
println("Mean times")
println("----------")
foreach((r)->println(rpad(r[1],27), ": ", mean(r[2])), results)


text:
abstract type KLIEPSolver end

struct CD_KLIEP <: KLIEPSolver end

### interface

KLIEP(Ψx, Ψy, ::CD_KLIEP) = KLIEP!(SparseIterate(size(Ψx, 1)), Ψx, Ψy)

function KLIEP!(x::SparseIterate, Ψx, Ψy)
    f = CDKLIEPLoss(Ψx, Ψy)
    g = ProxZero()
    coordinateDescent!(x, f, g)
end

function spKLIEP!(x::SparseIterate, Ψx, Ψy, λ)
    f = CDKLIEPLoss(Ψx, Ψy)
    g = ProxL1(λ)
    coordinateDescent!(x, f, g)
end

spKLIEP!(x::SparseIterate, f::CDKLIEPLoss, g::ProxL1) = coordinateDescent!(x, f, g)

function _compute_loadings(x, Ψx, Ψy)
    p, nx = size(Ψx)
    ny = size(Ψy, 2)
    Ψyr = transpose( rhat(x, Ψy) .* transpose(Ψy) )
    s = zeros(p)
    for k = 1:p
        s[k] = sqrt( (var(Ψx[k, :]) / nx) + (var(Ψyr[k, :]) / ny) )
    end
    s
end

function spKLIEP_scaled!(x::SparseIterate, Ψx, Ψy, λ)
    f = CDKLIEPLoss(Ψx, Ψy)

    γ = _compute_loadings(x, Ψx, Ψy)
    g = ProxL1(λ, γ)
    coordinateDescent!(x, f, g)

    γ = _compute_loadings(x, Ψx, Ψy)
    g = ProxL1(λ, γ)
    coordinateDescent!(x, f, g, CDOptions(; warmStart=true))
end

spKLIEP(Ψx, Ψy, λ, ::CD_KLIEP; loadings=true) =
    (loadings ? spKLIEP_scaled!(SparseIterate(size(Ψx, 1)), Ψx, Ψy, λ) : spKLIEP!(SparseIterate(size(Ψx, 1)), Ψx, Ψy, λ))

# Invert row of a Hessian using the scaled lasso
function Hinv_row(H, row, λ0)
    p = size(H, 1)
    e = zeros(p)
    e[row] = -1.

    x = SparseIterate(p)
    x[row] = 1.
    f = CDQuadraticLoss(H, e)

    # init sigma
    σ = sqrt( dot(x, H * x) )
    λ = e .+ 1.

    for iter = 1:10
        coordinateDescent!(x, f, ProxL1(λ0 * σ, λ))
        σnew = sqrt( dot(x, H * x) )
        if abs(σnew - σ) / σ < 1e-3
            break
        end
        σ = σnew
    end

    x
end


text:
@testset "34.find-first-and-last-position-of-element-in-sorted-array.jl" begin
    @test search_range([5, 7, 7, 8, 8, 10], 8) == [3, 4]
    @test search_range([5, 7, 7, 8, 8, 10], 6) == [-1, -1]
    @test search_range(Int[], 0) == [-1, -1]
    @test search_range([1, 1, 1, 2, 3, 3, 4, 4], 1) == [0, 2]
end


text:
using Distributed
using Printf

mutable struct Node{T<:AbstractFloat}
    split::T
    low::T
    up::T
    logit::Union{Nothing, T}
    left::Union{Node, Nothing}
    right::Union{Node, Nothing}
    function Node(
        split::T, low::T, up::T, logit::Union{T, Nothing} = nothing
    ) where T<:AbstractFloat
        new{T}(split, low, up, logit, nothing, nothing)
    end
end

##
isleaf(v::Node) = isnothing(v.left) && isnothing(v.right)
descend(v::Node{T}, x::T) where T<:AbstractFloat = (x ≤ v.split) ? (v.left, 0) : (v.right, 1)

# tostr(x::Union{Node, Float64, Nothing}; add_children=false) = begin
#     if isnothing(x)
#         return "(*)"
#     elseif isa(x, Float64)
#         return @sprintf("%.2f", x)
#     else
#         s = "Node($(tostr(x.split)), $(tostr(x.logit)))\n"
#         if add_children
#             s *= "|--$(tostr(x.left))\n"
#             s *= "|--$(tostr(x.right))"
#         end
#         return s
#     end
# end
# Base.show(io::IO, v::Node) = show(io, "Node(id=$(objectid(v)))")
# Base.print(io::IO, v::Node) = print(io, tostr(v, add_children=true))

function getval(v::Node{T}, x::T) where {T<:AbstractFloat}
    @assert !isnothing(v.logit) "can only evaluate nodes with values, use bfs_fill_tree"
    (x ≤ v.split) ? v.logit : -v. logit
end

function uniform_binary_splits(xmin::T, xmax::T, depth::Int) where T <: AbstractFloat
    lows = T[]
    mids = T[]
    ups = T[]
    levs = Int[]
    splitvals = collect(LinRange(xmin, xmax, 2^depth + 1))
    for d in 1:depth
        width = 2 ^ (depth - d)
        num_intervals = 2 ^ (d - 1)
        for j in 1:num_intervals
            start = 2 * (j - 1) * width + 1
            push!(lows, splitvals[start])
            push!(mids, splitvals[start + width])
            push!(ups, splitvals[start + 2width])
            push!(levs, d - 1)
        end
    end
    (lows=lows, mids=mids, ups=ups, levs=levs, splitvals=splitvals)
end

function eval_logprob(root::Node{T}, x::T) where T<:AbstractFloat
    @assert !isleaf(root) "must start from non-leaf"
    logprob = zero(T)
    v = root
    while !isnothing(v)
        β = getval(v, x)
        if β ≥ zero(T)
            logprob += -log(one(T) + exp(-β))
        else
            logprob += β - log(one(T)+ exp(β))
        end
        v, _ = descend(v, x)
    end
    return logprob
end

function eval_logdens(root::Node{T}, x::T) where T<:AbstractFloat
    @assert !isleaf(root) "must start from non-leaf"
    logprob = zero(T)
    v = root
    prev = v
    descdir = 0
    while !isnothing(v)
        β = getval(v, x)
        if β ≥ zero(T)
            logprob += -log(one(T) + exp(-β))
        else
            logprob += β - log(one(T)+ exp(β))
        end
        prev = v
        v, descdir = descend(v, x)
    end
    δ = (descdir == 0) ? prev.split - prev.low : prev.up - prev.split
    return logprob - log(δ)
end

##
function print_bfs(v::Node)
    visited = Queue{Node}()
    enqueue!(visited, v)
    while !isempty(visited)
        curr = dequeue!(visited)
        println(curr)
        !isnothing(curr.left) && enqueue!(visited, curr.left)
        !isnothing(curr.right) && enqueue!(visited, curr.right)
    end
end

##

function make_tree_from_bfs(
    lows::Vector{T},
    mid::Vector{T},
    ups::Vector{T},
    logits::Vector{S} = fill(nothing, length(mid))
) where {T<:AbstractFloat, S<:Union{Nothing, AbstractFloat}}
    num_nodes = length(mid)
    next_parent = Queue{Node}()
    root = Node(mid[1], lows[1], ups[1], logits[1])
    curr = root  # current parent
    for i in 2:num_nodes
        v = Node(mid[i], lows[i], ups[i], logits[i])
        while !(ups[i] ≈ curr.split || lows[i] ≈ curr.split)
            @assert !isempty(next_parent)  "disconnected tree!"
            curr = dequeue!(next_parent)
        end
        if ups[i] ≈ curr.split
            curr.left = v
        elseif lows[i] ≈ curr.split
            curr.right = v
        end
        enqueue!(next_parent, v)
    end
    return root
end

function make_tree(
    lows::Vector{T},
    mid::Vector{T},
    ups::Vector{T},
    logits::Vector{S} = fill(nothing, length(mid))
) where {T<:AbstractFloat, S<:Union{Nothing, AbstractFloat}}
    num_nodes = length(mid)
    nodes = [Node(mid[i], lows[i], ups[i], logits[i]) for i in 1:num_nodes]
    is_root = Set{Int}(1:num_nodes)
    for i in 1:num_nodes
        # serch for parent upwards
        child = nodes[i]
        child_lower, child_upper = lows[i], ups[i]
        for j in 1:num_nodes
            parent = nodes[j]
            parent_mid = mid[j]
            parent_lower, parent_upper = lows[j], ups[j]
            if child_lower ≈ parent_mid && child_upper ≈ parent_upper
                parent.right = child
                pop!(is_root, i)
                break
            elseif child_upper ≈ parent_mid && child_lower ≈ parent_lower
                parent.left = child
                pop!(is_root, i)
                break
            end
        end
    end
    if isempty(is_root)
        error("circular graph")
    elseif length(is_root) > 1
        error("multiple roots")
    end
    nodes[pop!(is_root)]
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

include("./models2D.jl")
include("./plot2D.jl")

using CairoMakie
CairoMakie.activate!(type = "svg")
filetype = ".svg"
# CairoMakie.activate!(type = "png")
# filetype = ".png"

function fig_corridor_therm2D_eps(eps_lab)
    # test blocking with shadow plot
    m = model_labyrinth([0.4, -0.4, 0.4, 0.2, 0.5, -0.3, -0.2], 0.1, 0.02)
    # blocking with tilewalk
    vfmat = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat)
    n = 30
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    @time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
    @time calculating_vf!(m, vfmat, normit = true)
    # solve Qp
    temp = zeros(m.nelem,1)
    set_bc_part!(m, temp, 1, 600)
    set_bc_part!(m, temp, 2:m.npar, 400)
    # eps = 1.0
    epsilon = zeros(m.nelem,1)
    set_bc_part!(m, epsilon, 1, 1.0)
    set_bc_part!(m, epsilon, 2:(m.npar-1), eps_lab)
    set_bc_part!(m, epsilon, m.npar, 1.0)
    @time Qp, G = tempsolver(m, vfmat, temp, epsilon)
    Qp_parts = [sum(Qp[m.elem2par[i].first:m.elem2par[i].last,1]) for i = 1:m.npar]
    println("heat flux start of labyrith: ", Qp_parts[1])
    println("heat flux end of labyrith: ", Qp_parts[end])
    # lab_min = minimum(Qp[m.elem2par[2].first:m.elem2par[end].last,1])
    # lab_max = maximum(Qp[m.elem2par[2].first:m.elem2par[end].last,1])
    # println("Qp of labyrith between: ", lab_min, " and ", lab_max)
    # heat flux density
    area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
    qp_area = Qp[:] ./ area[:]
    qp_area_kw = qp_area ./ 1000
    lab_min = minimum(qp_area_kw[m.elem2par[2].first:m.elem2par[end].last,1])
    lab_max = maximum(qp_area_kw[m.elem2par[2].first:m.elem2par[end].last,1])
    println("qp_area_kw of labyrith between: ", lab_min, " and ", lab_max)
    # modifying qw for displaying only incoming heat flux
    qp_area_kw_mod = qp_area_kw
    qp_area_kw_mod[qp_area_kw.>=0.0] .= 0.0
    qp_area_kw_mod .*= (-1)
    # 2D plot
    fig = Figure(resolution = (270, 400), font = "Arial", fontsize = 10)
    ax = fig[1, 1] = Axis(fig)
    linewidth = 1.0
    setup_axis(ax, linewidth)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    # changing legend in therm2D plotting to appear underneath
    cbar = plot_model_with_value(fig, ax, m, qp_area_kw_mod, "Incoming heat flux density in kW/m", cmap = :viridis, linewidth = linewidth*1.5, showcbar = true)
    # cbar.ticks = -50:20:50
    # cbar.limits = (lab_min,lab_max)
    cbar.limits = (0,0.4)
    cbar.ticks = 0.0:0.1:0.4
    cbar.vertical = false
    cbar.width = 120
    cbar.height = 12
    cbar.labelpadding = 1.0 # default 5.0
    cbar.ticklabelpad = 1.0 # default 3.0
    cbar.ticksize = 2.0 # default 6.0
    cbar.tellwidth = false # erlaubt es das colorbar und plot andere groeßen zueinander haben
    # plot_model_elements(fig, ax, m, m.elem2par[1].first:m.elem2par[1].last, color = :red, linewidth = linewidth*1.5)
    # display(fig)
    eps_save = replace(string(eps_lab),"." => "n")
    save("fig_c3_therm2D_eps"*eps_save*"_cb"*filetype, fig)
    # for i = 1:m.nnodes
    #     println(m.nodes[i].x, ";", m.nodes[i].y)
    # end
    # for i = 1:m.nelem
    #     println(m.elem[i].node1, ";", m.elem[i].node2, ";", m.elem[i].com, ";", m.elem[i].nvec, ";", m.elem[i].area)
    # end
    # for i = 1:m.nelem
    #     println(temp[i,:], ";", epsilon[i,:], ";", qp_area_kw_mod[i,:])
    # end
end

# fig_corridor_therm2D_eps(0.0001)
# fig_corridor_therm2D_eps(0.05)
# fig_corridor_therm2D_eps(0.75)

# include("./test/run_corridor.jl")


text:
# methods for estimating dtw #

abstract type DTWDistance{D <: Union{Function, Distances.PreMetric}} end

"""
    struct DTW{D} <: DTWDistance{D}

# Arguments:
- `radius`: The maximum allowed deviation of the matching path from the diagonal
- `dist`: Inner distance
- `transportcost` If >1, an additional penalty factor for non-diagonal moves is added.
"""
Base.@kwdef struct DTW{D} <: DTWDistance{D}
    "The maximum allowed deviation of the matching path from the diagonal"
    radius::Int
    dist::D = SqEuclidean()
    "If >1, an additional penalty factor for non-diagonal moves is added."
    transportcost::Float64 = 1.0
    DTW(r,dist=SqEuclidean(),transportcost=1) = new{typeof(dist)}(r,dist,transportcost)
end

"""
    struct SoftDTW{D, T} <: DTWDistance{D}

# Arguments:
- `γ`: smoothing parameter
- `dist`
- `transportcost`
"""
Base.@kwdef struct SoftDTW{D,T} <: DTWDistance{D}
    γ::T
    "The maximum allowed deviation of the matching path from the diagonal"
    dist::D = SqEuclidean()
    "If >1, an additional penalty factor for non-diagonal moves is added."
    transportcost::Float64 = 1.0
    SoftDTW(γ, dist=SqEuclidean(),transportcost=1) = new{typeof(dist), typeof(γ)}(γ,dist,transportcost)
end

"""
    struct FastDTW{D} <: DTWDistance{D}

- `radius`
- `dist` inner distance
"""
Base.@kwdef struct FastDTW{D} <: DTWDistance{D}
    radius::Int
    dist::D = SqEuclidean()
    FastDTW(r, dist=SqEuclidean()) = new{typeof(dist)}(r, dist)
end

Distances.evaluate(d::DTW, x, y) = dtw_cost(x, y, d.dist, d.radius)
Distances.evaluate(d::SoftDTW, x, y) = soft_dtw_cost(x, y, d.dist, γ=d.γ)
Distances.evaluate(d::FastDTW, x, y) =
    fastdtw(x, y, d.dist, d.radius)[1]

distpath(d::DTW, x, y) = dtw(x, y, d.dist)
distpath(d::DTW, x, y, i2min::AbstractVector, i2max::AbstractVector) =
    dtw(x, y, i2min, i2max, d.dist)
distpath(d::FastDTW, x, y) = fastdtw(x, y, d.dist, d.radius)

(d::DTWDistance)(x,y) = Distances.evaluate(d,x,y)

"""
    distance_profile(d::DTWDistance, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where S

Optimized method for computing the distance profile using DTW distances. kwargs are sent to [`dtwnn`](@ref).
"""
function SlidingDistancesBase.distance_profile(d::DTWDistance, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where S
    m = lastlength(Q)
    n = lastlength(T)
    n >= m || throw(ArgumentError("Q cannot be longer than T"))
    l = n-m+1
    res = dtwnn(Q, T, d.dist, d.radius; saveall=true, kwargs...)
    res.dists
end


text:
module Coords

import Base: ==

using AstroTime: Epoch, TimeScale

import ..position,
    ..state,
    ..velocity

using ..Bodies:
    CelestialBody,
    earth,
    grav_param

using ..Ephemerides: AbstractEphemeris

using ..Frames:
    AbstractFrame,
    icrf

import ..TwoBody:
    cartesian,
    keplerian,
    period

using StaticArrays: SVector

include("states.jl")
include("conversions.jl")
include("bodyfixed.jl")
include("events.jl")
include("trajectories.jl")

end


text:
using OpenCL

function testRelease()
   const cl = OpenCL
   function test()
       device, ctx, queue = cl.create_compute_context()
       local data = Array(Float32, 3)
       data[1] = 0.1
       data[2] = 0.22
       data[3] = 0.333

       local X = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=data)
       local result = Array(Float32, 3)
       cl.enqueue_read_buffer(queue, X, result, unsigned(0), nothing, true)
       println(data, "\n-----------\n", result)
       @assert(data == result)
   end

   try
      test()
      gc()
   catch err
      gc()
      throw(err)
   end
end


text:
# # Cameras

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/cameras.ipynb)

# In this section, we will showcase the cameras that are supported by
# `AutomotiveVisualization.jl`.
#
# All cameras make use of the `CameraState` object to represent the internal camera state.
#
#md # ```@docs
#md # CameraState
#md # ```
#
# The following camera types are supported by `AutomotiveVisualization`:
# `StaticCamera`, `TargetFollowCamera`, `SceneFollowCamera`, `ZoomingCamera`,
# `ComposedCamera`.
# The keyword arguments accepted by the `CameraState` object can also be passed to
# the constructors of these cameras, and they will be forwarded to the internal
# camera state. Now, we will show how the different cameras behave based on the
# simple stadium example.
#
# First, we prepare a sequence of scenes visualize

using AutomotiveSimulator
using AutomotiveVisualization
AutomotiveVisualization.colortheme["background"] = colorant"white"; # hide

nlanes = 4
nveh = 5
nticks = 100
timestep = 0.1

roadway = gen_stadium_roadway(nlanes)

scene = Scene([Entity(VehicleState(
    Frenet(roadway[LaneTag(1,rand(1:nlanes))], 10.0*i), roadway, 4.0 + 2.0rand()
), VehicleDef(), i) for i in 1:nveh])

models = Dict((
    i => LatLonSeparableDriver(ProportionalLaneTracker(), IntelligentDriverModel())
    for i in 1:nveh
))
set_desired_speed!.(values(models), 8.0 .+ 4.0rand(nveh))

scenes = simulate(scene, roadway, models, nticks, timestep)
nothing # hide

# The following helper function takes care of generating an animation from a
# sequence of scenes, using a given `camera`.
# Note that for `camera=nothing`, the render function attempts to fit all the
# renderable objects on the canvas.

using Reel

function animate(roadway, scenes, camera=nothing)
    animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
        i = Int(floor(t/dt)) + 1
        update_camera!(camera, scenes[i])
        renderables = [
            roadway, scenes[i],
            IDOverlay(scene=scenes[i], x_off=-2, y_off=1),
        ]
        render(renderables, camera=camera)
    end
    return animation
end
nothing # hide

# The call to `update_camera!(camera, scenes[i])` is important as it allows the
# internal camera state to be updated based on the positions of the vehicles.

# ## Default camera: fit to content
#
# The default camera uses the render instructions that are stored inside the
# `RenderModel` in order to compute a bounding box which holds all renderable
# objects.

animation = animate(roadway, scenes);
write("stadium_cam_default.gif", animation) # hide

#md # ![default camera](stadium_cam_default.gif)

# Note that some renderable objects (such as custom objects) may not be
# considered by the default fit to content camera when the bounding box is
# computed. Therefore, they may be cut off.

# ## Static Camera
#
# A static camera maintains its position, zoom and rotation. It is the only camera
# that doesn't change its state when calling `update_camera!()`.

animation = animate(roadway, scenes, StaticCamera(zoom=5.))
write("stadium_cam_static_default.gif", animation) # hide

#md # ![static default camera](stadium_cam_static_default.gif)

# We can pass in zoom, position and rotation parameters in order to change the camera state.

animation = animate(
    roadway, scenes,
    StaticCamera(position=VecE2(20.,10.), rotation=.2, zoom=8.)
)
write("stadium_cam_static_custom.gif", animation) # hide

#md # ![static custom camera](stadium_cam_static_custom.gif)

# ## Target Follow Camera
#
# In many applications, we are interested in tracking one particular vehicle.
# To do this, we can use a `TargetFollowCamera` which follows a vehicle with a
# given ID. Let's follow the vehicle with ID `3`:

animation = animate(roadway, scenes, TargetFollowCamera(3, zoom=5.))
write("stadium_cam_target.gif", animation) # hide

#md # ![target follow camera](stadium_cam_target.gif)

# Sometimes, it is desireable fix the x or y position of the camera.
# For example, we can use the `TargetFollowCamera` to track a vehicle
# horizontally (in x direction), while keeping the vertical position
# (y direction) fixed.
# To do this, we simply provide a fixed value for the corresponding keyword
# argument `y`.

animation = animate(roadway, scenes, TargetFollowCamera(3, zoom=10., y=40.))
write("stadium_cam_target_y.gif", animation) # hide

#md # ![target follow camera, fixed y](stadium_cam_target_y.gif)

# ## Scene Follow Camera
#
# In contrast to the `TargetFollowCamera`, the `SceneFollowCamera` tracks the
# center of mass of all the vehicles in a scene. The `padding` parameter specifies
# how man metres of padding there should be around the vehicles in the scene

animation = animate(roadway, scenes, SceneFollowCamera(padding=4.))
write("stadium_cam_scene.gif", animation) # hide

#md # ![scene follow camera](stadium_cam_scene.gif)

# If we desire to follow the scene, but fix the zoom level,
# we can do so by passing in the corresponding keyword argument `zoom`.

animation = animate(roadway, scenes, SceneFollowCamera(zoom=20.))
write("stadium_cam_scene_zoom.gif", animation) # hide

#md # ![scene follow camera](stadium_cam_scene_zoom.gif)

# ## Zooming Camera
#
# The zooming camera changes its zoom level by an increment `dz` anytime the
# `update_camera!()` function is called, until a zoom level of `zoom_target` is reached.

animation = animate(roadway, scenes, ZoomingCamera(zoom_target=20., dz=0.4, zoom=2.))
write("stadium_cam_zooming.gif", animation) # hide

#md # ![zooming camera](stadium_cam_zooming.gif)

# ## Composed Camera
#
# Finally, the `ComposedCamera` can be used to combine different camera behaviors.
# For example, one can use the `TargetFollowCamera` together with the
# `ZoomingCamera` in order to follow a vehicle and zoom in on it.
# The constructor of the `CompostedCamera` takes an array of cameras which are
# applied in the same sequence as they are provided in the array.

animation = animate(roadway, scenes, ComposedCamera([
    TargetFollowCamera(3),
    ZoomingCamera(zoom_target=25., dz=0.4, zoom=1.),
]))
write("stadium_cam_composed.gif", animation) # hide

#md # ![composed camera](stadium_cam_composed.gif)

# ## Custom Camera
#
# Finally, we will show how to build a custom camera. We will build a camera whose
# angle is always aligned with the axis of the target vehicle with ID `target_id`

mutable struct TargetRotateCamera{I} <: Camera
    state::CameraState
    target_id::I
end
TargetRotateCamera(target_id; kwargs...) = TargetRotateCamera(CameraState(;kwargs...), target_id)
nothing # hide

# In the `update_camera!()` function, we need to set the camera rotation to be the
# same as the vehicle yaw angle

function AutomotiveVisualization.update_camera!(camera::TargetRotateCamera{I}, scene::Scene{E}) where {I,E<:Entity}
    target = get_by_id(scene, camera.target_id)
    yaw = posg(target.state)[3]
    set_camera!(camera.state, rotation=-yaw)
    return camera.state
end
nothing # hide

# Finally, we can use our custom camera just like the other cameras

animation = animate(roadway, scenes, TargetRotateCamera(3, zoom=5.))
write("stadium_cam_custom.gif", animation) # hide

#md # ![custom camera](stadium_cam_custom.gif)

# Also, we can combine the camera with other cameras to our liking

animation = animate(roadway, scenes, ComposedCamera([
    TargetFollowCamera(3),
    TargetRotateCamera(3),
    ZoomingCamera(zoom_target=25., dz=0.4, zoom=1.),
]))
write("stadium_cam_composed_custom.gif", animation) # hide

#md # ![composed camera](stadium_cam_composed_custom.gif)


text:
export
    TransformProperties, calc_transform!

using ..Math

mutable struct TransformProperties{T <: AbstractFloat}
    position::Point{T}
    rotation::T       # radians
    scale::Point{T}

    aft::Math.AffineTransform{T}
    inverse::Math.AffineTransform{T}

    function TransformProperties{T}() where {T <: AbstractFloat}
        new(Point{T}(),
            0.0,
            Point{T}(1.0, 1.0),
            Math.AffineTransform{T}(),
            Math.AffineTransform{T}())
    end
end

function set_position!(prop::TransformProperties{T}, x::T, y::T) where {T <: AbstractFloat}
    Geometry.set!(prop.position, x, y);
end

function rotation_in_degrees!(prop::TransformProperties{T}) where {T <: AbstractFloat}
    rad2deg(prop.rotation)
end

function set_rotation_in_degrees!(prop::TransformProperties{T}, rotation::T) where {T <: AbstractFloat}
    prop.rotation = deg2rad(rotation);
end

function uniform_scale(prop::TransformProperties{T}) where {T <: AbstractFloat}
    @assert(prop.scale.x == prop.scale.y)

    prop.scale.x
end

function set_scale!(prop::TransformProperties{T}, s::T) where {T <: AbstractFloat}
    Geometry.set!(prop.scale, s, s);
end

function set_nonuniform_scale!(prop::TransformProperties{T}, sx::T, sy::T) where {T <: AbstractFloat}
    Geometry.set!(prop.scale, sx, sy);
end

# Mostly used by Filter type nodes
function calc_filtered_transform!(prop::TransformProperties{T},
        exclude_translation::Bool, exclude_rotation::Bool, exclude_scale::Bool,
        aft::Math.AffineTransform{T},
    ) where {T <: AbstractFloat}
    to_identity!(aft)

    if !exclude_translation
        make_translate!(aft, prop.position.x, prop.position.y)
    end

    if !exclude_rotation && prop.rotation ≠ 0.0
        rotate!(aft, prop.rotation)
    end

    if !exclude_scale && (prop.scale.x ≠ 1.0 || prop.scale.y ≠ 1.0)
        scale!(aft, prop.scale.x, prop.scale.y)
    end

    nothing
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function mhash(A::Array{Int64, 2})
  return reduce(*, map(string, A))
end

function nonempty(A::Array{Int64, 2})
  return prod(vec(sum(A, 1)) .+ vec(sum(A, 2))) > 0
end

function hashmat(x)
  V = round(Int64, sqrt(length(x)))
  A = reshape(x, (V, V))
  if nonempty(A)
    return mhash(A)
  end
end

function unhashmat(x)
  c = Int64[]
  for i in 1:length(x)
    push!(c, parse(Int64, x[i]))
  end
  V = round(Int64, sqrt(length(x)))
  return reshape(vec(c), (V, V))
end

function from_binarytree(h, s)
   tree = zeros(Int64, 2^(h+1)-2)
   for i in 1:length(tree)
      if iseven(i)
         tree[i] = 1
      end
   end
   terminal_nodes = (length(tree)-2^h+1):length(tree)
   paths = []
   for n in terminal_nodes
      i = h
      position = n
      sequence = zeros(Int64, h)
      while(i > 0)
         sequence[i] = tree[position]
         position = round(Int,floor((position - 1)/2))
         i = i-1
      end
      if sum(sequence) == s
         push!(paths, sequence)
      end
   end
   return paths
end


text:
@testset "Problem" begin
    # object creation
    lp = Problem{Int}()
    A = sprand(Int, 3, 5, 0.5)
    b = rand(Int, 3)
    c = rand(Int, 5)
    lo = rand(Int, 5)
    hi = rand(Int, 6) # too long
    @test_throws DimensionMismatch load_problem!(lp, A, b, c, lo, hi)
end


text:
# TODO
#=
Requirements:
- incomplete KSAM, i.e. (KEY_S, ignore scancode, key_down, ignore modifier)
- info message for help screen
- function call on success

=#

struct Key
    skip::NTuple{4, Bool}
    ksam::NTuple{4, Int64}
end

"""
    Key(key[; scancode, action, modifier])

Creates a Key for key bindings. Every argument is optional and won't be checked
if it isn't defined here.
"""
function Key(
        key::Int64 = -10;
        scancode::Int64 = -10,
        action::Int64 = -10,
        modifier::Int64 = -10
    )
    Key(
        (key == -10, scancode == -10, action == -10, modifier == -10),
        (key, scancode, action, modifier)
    )
end

function ==(key::Key, ksam::NTuple{4, Int64})
    (key.skip[1] || (key.ksam[1] == ksam[1])) &&
    (key.skip[2] || (key.ksam[2] == ksam[2])) &&
    (key.skip[3] || (key.ksam[3] == ksam[3])) &&
    (key.skip[4] || (key.ksam[4] == ksam[4]))
end

# This shouldn't be global
struct Hotkeys
    window::GLWindow.Screen
    hotkeys::Dict{Key, Tuple{String, Function}}
    key_listener_signal::Signal
    signals::Vector{Signal}
    destructor_signal::Signal
end

function Hotkeys(window::GLWindow.Screen)
    hotkeys = Dict{Key, Tuple{String, Function}}()
    key_listener_signal = map(window.inputs[:keyboard_buttons]) do ksam
        #!value(window.inputs[:window_open]) && return nothing
        for key in keys(hotkeys)
            if key == ksam
                hotkeys[key][2](ksam)
                break
            end
        end
        nothing
    end
    signals = Signal[]
    destructor_signal = map(window.inputs[:window_open]) do is_open
        if !is_open
            for s in signals[end:-1:1]
                close(s)
            end
            close(key_listener_signal)
            empty!(signals)
            empty!(hotkeys)
            println("Signals closed.")
        end
        nothing
    end
    Hotkeys(window, hotkeys, key_listener_signal, signals, destructor_signal)
end

function push!(h::Hotkeys, key::Key, info::String, action::Function)
    push!(h.hotkeys, Pair(key, (info, action)))
end

"""
    init_screenshots!(Hotkeys[, key = Key(GLFW.KEY_F2)])

Sets *key* as a hotkey for screenshots.
"""
function init_screenshots!(
        h::Hotkeys,
        key::Key=Key(GLFW.KEY_F2, action = GLFW.RELEASE)
    )
    push!(h, key, "Screenshot",
        _ -> begin
            name = string(time_ns())*".png"
            println("Screenshot saved as " * name)
            screenshot(h.window, path=name)
        end
    )
    nothing
end

"""
    init_QWEASD_controls!(Hotkeys, cam[, step::Float32=1°])

Sets up Q, W, E, A, S, D as camera controls keys. A and D perform left and right
rotations, W and S do up and down rotation, Q and E do rolls. *step* can be used
to change the rotational speed.
"""
function init_QWEASD_controls!(
        h::Hotkeys,
        cam::Camera,
        step::Signal{Float32} = Signal(Float32(pi/180))
    )

    pre() = begin
        eyepos = value(cam.eyeposition)
        lookat = value(cam.lookat)
        up = value(cam.up)
        dir = eyepos - lookat
        distance = norm(dir)
        dir_norm = normalize(dir)
        right = normalize(cross(dir, up))
        # lookat, up, dir, dir_norm, right
        return dir_norm, right, up, lookat, dir
    end

    post(rot, lookat, up, dir) = begin
        push!(cam.eyeposition, lookat + rot * dir)
        push!(cam.up, normalize(rot * up))
        return nothing
    end

    push!(h, Key(GLFW.KEY_Q), "Anticlockwise roll", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(dir_norm, -value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)

    push!(h, Key(GLFW.KEY_E), "Clockwise roll", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(dir_norm, value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)

    push!(h, Key(GLFW.KEY_W), "Upwards tilt", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(right, -value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)

    push!(h, Key(GLFW.KEY_S), "Downwards tilt", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(right, value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)

    push!(h, Key(GLFW.KEY_A), "Left rotation", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(up, value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)

    push!(h, Key(GLFW.KEY_D), "Right rotation", ksam -> begin
        if (ksam[3] == GLFW.PRESS) || (ksam[3] == GLFW.REPEAT)
            dir_norm, right, up, lookat, dir = pre()
            rot = Mat3f0(Quaternions.qrotation(up, -value(step)))
            post(rot, lookat, up, dir)
            nothing
        end
    end)
    nothing
end

function init_angles(h::Hotkeys)
    step_signal = Signal(Float32(pi/180))
    push!(h.signals, step_signal)

    push!(
        h,
        Key(GLFW.KEY_1, action = GLFW.PRESS),
        "Set angle step to 1°.",
        _ -> push!(step_signal, Float32(pi/180))
    )

    for i in 2:9
        push!(
            h,
            Key(48+i, action = GLFW.PRESS),
            "Set angle step to 1°.",
            _ -> push!(step_signal, Float32(2pi/i))
        )
    end

    step_signal
end


text:
"""
_compute_mutliplier: Additional function to estimate Lagrange multiplier of the problems
    (guarantee if LICQ holds)

`_compute_mutliplier(pb :: AbstractNLPModel, x :: T, gx :: T, cx :: T, Jx :: MT; active_prec_c :: Real = 1e-6, active_prec_b :: Real = 1e-6)`
"""
function _compute_mutliplier(
  pb::AbstractNLPModel,
  x::T,
  gx::T,
  cx::T,
  Jx::MT;
  active_prec_c::Real = 1e-6,
  active_prec_b::Real = 1e-6,
) where {MT, T}
  n = length(x)
  nc = cx == _init_field(typeof(cx)) ? 0 : length(cx)

  #active res_bounds
  Ib = findall(x -> (norm(x) <= active_prec_b), min(abs.(x - pb.meta.lvar), abs.(x - pb.meta.uvar)))
  if nc != 0
    #active constraints
    Ic = findall(
      x -> (norm(x) <= active_prec_c),
      min(abs.(cx - pb.meta.ucon), abs.(cx - pb.meta.lcon)),
    )

    Jc = hcat(Matrix(1.0I, n, n)[:, Ib], Jx'[:, Ic])
  else
    Ic = []
    Jc = hcat(Matrix(1.0I, n, n)[:, Ib])
  end

  mu, lambda = zeros(n), zeros(nc)
  if (Ib != []) || (Ic != [])
    l = Jc \ (-gx)
    mu[Ib], lambda[Ic] = l[1:length(Ib)], l[(length(Ib) + 1):length(l)]
  end

  return mu, lambda
end


text:
# this file contains the code to wrap hidapi using Clang.jl

using hidapi_jll
using Clang

mkpath(joinpath(@__DIR__,"..","gen"))

wc = Clang.init(; headers = [joinpath(hidapi_jll.artifact_dir,"include","hidapi","hidapi.h")],
            output_file = joinpath(@__DIR__,"..","gen","hidapi_api.jl"),
            common_file = joinpath(@__DIR__,"..","gen","hidapi_common.jl"),
            clang_includes = vcat(CLANG_INCLUDE),
            clang_args = map(x->"-I"*x, find_std_headers()),
            header_wrapped = (root, current)->root == current,
            header_library = x->"hidapi",
            clang_diagnostics = true,
            )

run(wc)

# certain adjustments need to be done afterwards:
# -remove ctypes.jl (not needed in this package)
# -remove the LibTemplate (we are relying on the _jll to provide the library)
# -remove "const HID_API_EXPORT_CALL = HID_API_EXPORT" from hidapi_common.jl


text:
module ReactiveMPTestingHelpers

using Test
using ReactiveMP

import ReactiveMP: deep_eltype

@testset "Helpers" begin

    @testset "deep_eltype" begin

        for type in [ Float32, Float64, Complex{Float64}, BigFloat ]

            @test deep_eltype(type) === type
            @test deep_eltype(zero(type)) === type

            vector = zeros(type, 10)
            matrix = zeros(type, 10, 10)
            vector_of_vectors  = [ vector, vector ]
            vector_of_matrices = [ matrix, matrix ]
            matrix_of_vector   = [ vector vector; vector vector ]
            matrix_of_matrices = [ matrix matrix; matrix matrix ]

            @test deep_eltype(vector) === type
            @test deep_eltype(matrix) === type
            @test deep_eltype(vector_of_vectors) === type
            @test deep_eltype(vector_of_matrices) === type
            @test deep_eltype(matrix_of_vector) === type
            @test deep_eltype(matrix_of_matrices) === type

        end

    end

end

end


text:
using Revise
module BayesMNIST

using Flux, Flux.Data.MNIST, Statistics, LinearAlgebra, StatsBase
using Flux: onehotbatch, onecold, crossentropy, throttle, @epochs
using Flux: kldivergence
using CUDAapi, CuArrays, Random, Distributions, DataFrames#, Plots
using Revise

#load the data
#the purpose is to explore flux in a Bayesian context

Random.seed!(11)

  include("MNISTModels.jl")
  #include("TuringBNN.jl")

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using ADerrors
using ADjson
using Test

# Check if file produced by pyerrors can be read.
tmp = load_json("./data/pyerrors_out")
display(tmp)
println()

# Check if file with a constant produced by pyerrors can be read.
tmp =load_json("./data/covobs_out")
display(tmp)
println()

# Check if file containing and array produced by pyerrors can be read.
tmp = load_json("./data/array_out")
display(tmp)
println()

# Check if file containing a correlator produced by pyerrors can be read.
tmp = load_json("./data/corr_out")
display(tmp)
println()

# Check if file containing a gapped observable produced by pyerrors can be read.
tmp = load_json("./data/gapped_out")
display(tmp)
println()

# Write file to disk, read it and check equality
test = uwreal(rand(4000), "Test ensemble")
dump_to_json(test, "test_file", "File containing an observable defined on 4000 random samples")
test2= load_json("test_file")
println()

iamzero = test - test2
uwerr(iamzero)
@test isapprox(iamzero.mean, 0.0; atol = 1e-14)
@test isapprox(iamzero.err, 0.0; atol = 1e-14)

# Write file containing constant to disk, read it and check equality
test *= uwreal([1.1, 0.2], "Constant")
dump_to_json(test, "constant_file", "File containing Monte carlo data and a constant with error")
test2= load_json("constant_file")
println()

iamzero = test - test2
uwerr(iamzero)
@test isapprox(iamzero.mean, 0.0; atol = 1e-14)
@test isapprox(iamzero.err, 0.0; atol = 1e-14)

# Test exceptions
gapped_observable = uwreal(Vector{Float64}(collect(1:2:20)), "go", collect(1:2:20), 20)
@test_throws ErrorException dump_to_json(gapped_observable, "this_should_not_work")

@test_throws ErrorException load_json("this_file_does_not_exist")


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

using Documenter, FEMDynamics

makedocs(modules=[FEMDynamics],
         format = :html,
         checkdocs = :all,
         sitename = "FEMDynamics.jl",
         pages = ["index.md"]
        )


text:
using SHIPs, JuLIP, BenchmarkTools, LinearAlgebra

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

trans = PolyTransform(2, 1.0)
fcut = PolyCutoff2s(2, 0.5, 3.0)
BB = [ SHIPBasis(SparseSHIP(2, 20, wL=1.0), trans, fcut),
       SHIPBasis(SparseSHIP(3, 16, wL=1.5), trans, fcut),
       SHIPBasis(SparseSHIP(4, 12, wL=1.5), trans, fcut),
       SHIPBasis(SparseSHIP(5, 10, wL=1.5), trans, fcut) ]

Nat = 30
Rs, Zs, z0 = randR(Nat)
btmp = SHIPs.alloc_temp(BB[1])
@info("profile precomputation of A")
@btime SHIPs.precompute_A!($btmp, $(BB[1]), $Rs, $Zs, 1)

@info("profile ship-basis and fast-ship site energies")
for n = 2:5
   @info("  body-order $(n+1):")
   Rs, Zs, z0 = randR(Nat)
   B = BB[n-1]
   coeffs = randcoeffs(B)
   🚢 = SHIP(B, coeffs)
   b = SHIPs.alloc_B(B)
   btmp = SHIPs.alloc_temp(B, length(Rs))
   tmp = SHIPs.alloc_temp(🚢, length(Rs))
   @info("     evaluate a site energy:")
   print("         SHIPBasis: "); @btime SHIPs.eval_basis!($b, $btmp, $B, $Rs, $Zs, $z0)
   print("         SHIP     : "); @btime SHIPs.evaluate!($tmp, $🚢, $Rs, $Zs, $z0)

   dEs = zeros(JVecF, length(Rs))
   db = SHIPs.alloc_dB(B, length(Rs))
   dbtmp = SHIPs.alloc_temp_d(B, length(Rs))
   tmp = SHIPs.alloc_temp_d(🚢, length(Rs))

   @info("     site energy gradient:")
   print("         SHIPBasis: "); @btime SHIPs.eval_basis_d!($db, $dbtmp, $B, $Rs, $Zs, $z0)
   print("         SHIP     : "); @btime SHIPs.evaluate_d!($dEs, $tmp, $🚢, $Rs, $Zs, $z0)
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
function read_atoms(fn::IO)
    sane_readline() = strip(strip_split(readline(fn), ":")[1], '\'')
    nspecies = parse(Int, sane_readline())
    species_counter = 0

    reading_atoms = false
    atname = :nothing
    natoms = 0
    atcounter = 0
    positions = Point3{Float64}[]
    bfcmts = Vec3{Float64}[]
    atoms = Tuple{Symbol,Point3{Float64},Vec3{Float64}}[]
    while species_counter < nspecies
        l = sane_readline()
        if !reading_atoms
            atname = Symbol(split(l, ".")[1])
            natoms = parse(Int, sane_readline())
            reading_atoms = true
        elseif atcounter < natoms
            sline = split(l)
            push!(positions, Point3{Float64}(parse.(Float64, sline[1:3])))
            push!(bfcmts, Vec3{Float64}(parse.(Float64, sline[4:end])))
            atcounter += 1
        elseif isempty(l)
            for (p, mag) in zip(positions, bfcmts)
                push!(atoms, (atname, p, mag))
            end
            atname          = :nothing
            reading_atoms   = false
            natoms          = 0
            atcounter       = 0
            positions       = Point3{Float64}[]
            bfcmts          = Vec3{Float64}[]
            species_counter += 1
        end
    end
    return atoms
end

function elk_parse_DFTU(atoms::Vector{Atom}, blocknames_flaglines::Dict{Symbol,Any})
    lines = pop!(blocknames_flaglines, Symbol("dft+u"))
    t_l = split(lines[1])
    species_start = 3
    if length(t_l) == 1
        dftu_type = parse(Int, t_l[1])
        inpdftu   = parse(Int, lines[2])
    else
        dftu_type, inpdftu = parse.(Int, t_l)
        species_start = 2
    end
    species_names = name.(unique(atoms))
    for line in lines[species_start:end]
        sline         = split(line)
        species_id, l = parse.(Int, sline[1:2])
        U, J0         = parse.(Float64, sline[3:4])
        for at in filter(x -> x.name == species_names[species_id], atoms)
            at.dftu = DFTU{Float64}(; U = U, J0 = J0, l = l)
        end
    end
    return dftu_type, inpdftu
end

"""
    elk_read_calculation(filename; execs=[Exec(exec="elk")], run=true, structure_name="noname")

Reads an Elk calculation file. The `ELK_EXEC` inside execs gets used to find which flags are allowed in this calculation file, and convert the read values to the correct Types.
Returns a `Calculation{Elk}` and the `Structure` that is found in the calculation.
"""
function elk_read_calculation(fn::String; exec = Exec(; exec = "elk"), run = true,
                              structure_name = "noname")
    blocknames_flaglines = Dict{Symbol,Any}()
    atoms = Atom[]
    calculations = Calculation[]
    dir, file = splitdir(fn)

    #TODO it's a bit messy that the atoms block is not handled like this, not sure if we want to clean this up
    #blocks that require specific parsing
    special_blocks = (:avec, :tasks, Symbol("dft+u"))

    open(fn, "r") do f
        sane_readline() = strip_split(readline(f), ':')[1]
        while !eof(f)
            line = sane_readline()
            isempty(line) || line[1] == '!' && continue
            if !isempty(line)
                blockname = Symbol(line)
                info = elk_block_info(blockname)
                if blockname == :atoms
                    atoms = read_atoms(f)
                elseif blockname == :wannierExtra
                    @warn "Please supply the .win file instead of wannierExtra block for extracting the Wannier Calculation."

                    wflags, wdata, ab, cb, proj_block = wan_read_calculation(f)
                    continue
                    # push!(calculations, Calculation{Wannier90}("wannier", dir, wflags, wdata, execs, run))
                elseif info !== nothing
                    blocklines = String[]
                    line = sane_readline()
                    while !isempty(line)
                        push!(blocklines,
                              strip(replace(strip(split(line, ':')[1]), "'" => ""), '.'))
                        line = sane_readline()
                    end
                    if !in(blockname, special_blocks)
                        blocknames_flaglines[blockname] = parse_block_flags(blockname,
                                                                            collect(Iterators.flatten(split.(filter(x -> !isempty(x),
                                                                                blocklines)))))
                    else
                        blocknames_flaglines[blockname] = filter(x -> !isempty(x),
                                                                 blocklines)
                    end
                end
            end
        end
    end
    # for i in
    #cell
    scale = haskey(blocknames_flaglines, :scale) ? last(pop!(blocknames_flaglines[:scale])) : 1.0
    cell  = uconvert.(Ang, Mat3(reshape(scale .* parse.(Float64, collect(Iterators.flatten(split.(pop!(blocknames_flaglines, :avec))))), 3, 3) .* 1a₀))'

    newats = [Atom(; name = x[1], element = x[1].element, position_cryst = x[2],
                   position_cart = cell * x[2], magnetization = x[3]) for x in atoms]
    #structure #TODO make positions be in lattices coordinates
    structure = Structure(cell, newats)

    #different tasks
    tasks_     = pop!(blocknames_flaglines, :tasks)
    should_run = map(x -> !occursin("!", x), tasks_)
    tasks      = map(x -> strip(strip(x, '!')), tasks_)

    #The wan tasks can prbably be skipped completely and added depending on the wannier calculation
    wan_tasks = String[]
    for (x, run) in zip(tasks, should_run)
        if x ∈ ["601", "602", "603", "604", "605"]
            push!(wan_tasks, x)
        else
            push!(calculations,
                  calculation_from_task(x, blocknames_flaglines, dir, execs, run))
        end
    end
    if haskey(blocknames_flaglines, :wannier) && !isempty(wan_tasks)
        flags = pop!(blocknames_flaglines, :wannier)
        # flags[:elk2wan_tasks] = wan_tasks
        push!(calculations,
              Calculation{Elk}(; name = "elk2wannier", flags = flags,
                               exec = exec, run = true))
    end
    for f in (:ngrid, :vkloff, :plot1d, :plot2d, :plot3d)
        haskey(blocknames_flaglines, f) && pop!(blocknames_flaglines, f)
    end

    #elk2wannier

    flags = Dict{Symbol,Any}()

    if haskey(blocknames_flaglines, Symbol("dft+u"))
        flags[:dftu], flags[:inpdftu] = elk_parse_DFTU(structure.atoms,
                                                       blocknames_flaglines)
    end

    for (k, v) in blocknames_flaglines
        merge!(flags, v)
    end

    for i in calculations
        if isempty(i.flags)
            i.flags = deepcopy(flags)
        end
    end
    return calculations, structure
end

function find_data(flags, blocknames_flaglines)
    data = InputData[]
    for f in flags
        if haskey(blocknames_flaglines, f)
            push!(data, InputData(f, :nothing, blocknames_flaglines[f]))
        end
    end
    return data
end

function calculation_from_task(task, blocknames_flaglines, dir, exec, run)
    if task ∈ ["0", "1", "2", "3", "5"]
        data = find_data((:ngridk, :vkloff), blocknames_flaglines)
    elseif task ∈ ["20", "21"]
        data = find_data((:plot1d, :plot2d, :plot3d), blocknames_flaglines)
    end
    return Calculation{Elk}(; name = task, data = data, exec = exec, run = run)
end

function parse(::Type{UnitRange{Int}}, l::AbstractString)
    for delim in [":", ",", "-"]
        if occursin(delim, l)
            return UnitRange(parse.(Int, split(l, delim))...)
        end
    end
end

function parse_block_flags(blockname::Symbol, lines::Vector{<:AbstractString})
    flags = Dict{Symbol,Any}()
    i = 1
    totlen = length(lines)
    for flaginfo in elk_block_info(blockname).flags
        flagname = flaginfo.name
        i > length(lines) && break
        typ = eltype(flaginfo)
        if typ == String || typ == Vector{String}
            flags[flagname] = lines[i]
            i += 1
        elseif typ == Symbol
            flags[flagname] = Symbol(lines[i])
            i += 1
        else
            eltyp = eltype(typ)
            # The assumption is that the infinite repeating values are at the end of a block
            if flaginfo == last(elk_block_info(blockname).flags)
                typelen = length(typ)
                parsed = parse.(eltyp, lines[i:end])
                flags[flagname] = typ[]
                @assert length(parsed) % typelen == 0
                for j in 1:typelen:length(parsed)
                    push!(flags[flagname], typ(parsed[j:j+typelen-1]...))
                end
            elseif typ <: UnitRange
                flags[flagname] = parse(typ, lines[i])
                i += 1
            else
                typelen = length(typ)
                parsed = parse.(eltyp, lines[i:i+typelen-1])
                flags[flagname] = length(parsed) == 1 ? parsed[1] : typ(parsed)
                i += typelen
            end
        end
    end
    return flags
end

function elk_write_structure(f, structure)
    write(f, "atoms\n\t")
    species = unique(map(x -> x.name, structure.atoms))
    nspecies = length(species)
    write(f, "$nspecies\n")
    for s in species
        write(f, "\t'$s.in'\n")
        ats = filter(x -> x.name == s, structure.atoms)
        write(f, "\t$(length(ats))\n")
        for a in ats
            write(f, "\t")
            for i in a.position_cryst
                write(f, "$(round(i, digits=8)) ")
            end
            for i in magnetization(a)
                write(f, "$i ")
            end
            write(f, "\n")
        end
        write(f, "\n")
    end
    write(f, "\n")
    write(f, "avec\n\t")

    for i in eachindex(structure.cell')
        c = ustrip(uconvert(a₀, structure.cell[i]))
        if i % 3 == 0
            write(f, "$c\n\t")
        else
            write(f, "$c ")
        end
    end
    return write(f, "\n")
end

function construct_flag_blocks(calculations::Vector{Calculation{Elk}})
    all_flags = Dict{Symbol,Any}()
    for i in calculations
        merge!(all_flags, i.flags)
        for d in data(i)
            merge!(all_flags, d.data)
        end
    end
    block_flags = Dict{Symbol,Any}()

    for b in ELK_CONTROLBLOCKS
        bvals = Any[]
        for f in b.flags
            if haskey(all_flags, f.name)
                push!(bvals, pop!(all_flags, f.name))
            end
        end
        if !isempty(bvals)
            block_flags[b.name] = bvals
        end
    end
    return block_flags
end

function elk_write_flag_val(f, val)
    write(f, "\t")
    if isa(val, UnitRange)
        write(f, "$(first(val))-$(last(val))")
    elseif isa(val, AbstractVector)
        for i in val
            write(f, "$i ")
        end
    else
        if isa(val, AbstractString)
            write(f, "'$val'")
        else
            write(f, "$val")
        end
    end
    return write(f, "\n")
end

function elk_write_DFTU(f::IO, structure::DFC.Structure, dftu_vals::Vector)
    write(f, "dft+u\n")
    for v in dftu_vals
        elk_write_flag_val(f, v)
    end
    species = unique(structure.atoms)
    for (i, s) in enumerate(species)
        write(f, "\t$i $(string(Elk, dftu(s)))\n")
    end
    return write(f, "\n")
end

function save(calculations::Vector{Calculation{Elk}}, structure::Structure, dir)
    tasks = map(x -> x.run ? "$(x.name)" : "!$(x.name)",
                filter(x -> x.name != "elk2wannier", calculations))
    elk2wan_calculation = getfirst(x -> x.name == "elk2wannier", calculations)
    if elk2wan_calculation != nothing && haskey(elk2wan_calculation.flags, :elk2wan_tasks)
        append!(tasks, pop!(elk2wan_calculation.flags, :elk2wan_tasks))
    end
    flags = construct_flag_blocks(calculations)
    open(joinpath(dir, "elk.in"), "w") do f
        write(f, "tasks\n")
        for t in tasks
            write(f, "\t$t\n")
        end
        write(f, "\n")

        elk_write_structure(f, structure)

        haskey(flags, Symbol("dft+u")) &&
            elk_write_DFTU(f, structure, pop!(flags, Symbol("dft+u")))

        #remaining generic blocks
        for (k, v) in flags
            write(f, "$k\n")
            for v_ in values(v)
                if isa(v_, Vector{<:AbstractVector})
                    for v__ in v_
                        elk_write_flag_val(f, v__)
                    end
                else
                    elk_write_flag_val(f, v_)
                end
            end
            write(f, "\n")
        end
    end
end