text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.4
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using Plots
using Optim

traceon = Optim.Options(
    store_trace = true,
    extended_trace = true
)

# model
Distributions.TDist(μ, ρ, ν) = LocationScale(μ, ρ, TDist(ν))

# test data
X = [-0.01, 0.01, 1.0]

# estimations with different initial values
r1 = optimize(x -> -loglikelihood(TDist(x[1], eps() + 10^x[2], eps() + 10^x[3]), X), [0.33, 0.0, 0.16], LBFGS(), traceon)
r2 = optimize(x -> -loglikelihood(TDist(x[1], eps() + 10^x[2], eps() + 10^x[3]), X), [0.33, 0.0, 0.18], LBFGS(), traceon)
r3 = optimize(x -> -loglikelihood(TDist(x[1], eps() + 10^x[2], eps() + 10^x[3]), X), [0.33, 0.0, 0.27], LBFGS(), traceon);

# %%
function plot_result(r, X)
    x = r.initial_x
    m = r.minimizer
    l = -minimum(r)
    t = "initial = $x\nfinal = $(round.(m; digits=3))\nfinal loglikelihood = $(round(l; digits=3))"
    xtick = axes(X, 1)
    normal = fit_mle(Normal, X)

    plot(; legend=:topleft)
    scatter!(xtick, X; label="data", xtick, xlim=(first(xtick)-0.5, last(xtick)+0.5))
    hline!([m[1]]; label="result of t-dist. model", lw=1.5)
    hline!([mean(normal)]; label="result of normal dist. model", ls=:dash, lw=1.5)
    title!(t; titlefontsize=10)
end

R1 = plot_result(r1, X)
R2 = plot_result(r2, X)
R3 = plot_result(r3, X)
plot(R1, R2, R3; size=(800, 700), xtick=1:3)

# %%
function plot_trace(r)
    x = r.initial_x
    c = hcat(x, (t.metadata["x"] for t in r.trace)...)
    l = -minimum(r)
    m = r.minimizer
    t = "initial = $x\nfinal = $(round.(m; digits=3))\nfinal loglikelihood = $(round(l; digits=3))"

    plot(; legend=:topleft)
    scatter3d!([x[1]], [x[2]], [x[3]]; label="initial", c=:blue, ms=3)
    plot!(c[1,:], c[2,:], c[3,:]; label="path", c=:darkcyan)
    scatter3d!([m[1]], [m[2]], [m[3]]; label="final", c=:yellow, ms=3)
    plot!(; xlabel="μ", ylabel="log₁₀ρ", zlabel="log₁₀ν")
    title!(t; titlefontsize=10)
end

P1 = plot_trace(r1)
P2 = plot_trace(r2)
P3 = plot_trace(r3)
plot(P1, P2, P3; size=(800, 700), tickfontsize=6)

# %%
function f(μ, log10ν, X)
    log10rhos = range(-15, 10; length=200)
    -maximum(loglikelihood(TDist(μ, 10^log10ρ, 10^log10ν), X) for log10ρ in log10rhos)
end

function plot_trace2d(r)
    x = r.initial_x
    c = hcat(x, (t.metadata["x"] for t in r.trace)...)
    l = -minimum(r)
    m = r.minimizer
    t = "initial = $x\nfinal = $(round.(m; digits=3))\nfinal loglikelihood = $(round(l; digits=3))"

    mus = range(-0.1, 0.4; length=200)
    log10nus = range(-2.5, 8.0; length=200)
    z = f.(mus', log10nus, Ref(X))

    plot(; legend=:topleft, colorbar=false)
    plot!(; xlim=extrema(mus), ylim=extrema(log10nus))
    heatmap!(mus, log10nus, z; clim=(-1, 5), c=reverse(cgrad(:CMRmap)))
    scatter!([x[1]], [x[3]]; label="initial", ms=4, c=:blue)
    plot!(c[1,:], c[3,:]; label="path", c=:cyan, lw=1.5)
    scatter!([m[1]], [m[3]]; label="final", ms=4, c=:yellow)
    plot!(; xlabel="μ", ylabel="log₁₀ν")
    title!(t; titlefontsize=10)
end

Q1 = plot_trace2d(r1)
Q2 = plot_trace2d(r2)
Q3 = plot_trace2d(r3)
plot(Q1, Q2, Q3; size=(800, 700))

# %%

max_stars_repo_path: 0024/TDist model is not log-convex - LBFGS.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
#=
There exists a staircase with N steps, and you can climb up either 1 or 2 steps at a time. Given N, write a function that returns the number of unique ways you can climb the staircase. The order of the steps matters.

For example, if N is 4, then there are 5 unique ways:

    1, 1, 1, 1
    2, 1, 1
    1, 2, 1
    1, 1, 2
    2, 2

What if, instead of being able to climb 1 or 2 steps at a time, you could climb any number from a set of positive integers X? For example, if X = {1, 3, 5}, you could climb 1, 3, or 5 steps at a time.
=#

using Test
include("Solutions/problem12_count_climbing_ways.jl")

@test climbingways(4,2) == 5

@test climbingways(4,3) == 7

#=
if N = 4, and X = 3, then there is 7 ways
    1, 1, 1, 1
    2, 1, 1
    1, 2, 1
    1, 1, 2
    2, 2
    3, 1
    1, 3
=#

ext: jl

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using ClimateBase, Test, Dates
using Statistics

Time = ClimateBase.Time

# TODO: Further test spatial averaging by making one hemisphere 1 and other 0

# Create the artificial dimensional array A that will be used in tests
function monthly_insolation(t::TimeType, args...)
    d = ClimateBase.monthspan(t)
    mean(insolation(τ, args...) for τ in d)
end

lats = -86:4:86
lons = collect(0.5:10:360)
t = Date(2000, 3, 15):Month(1):Date(2020, 3, 15)
tdense = Date(2000, 3, 15):Day(1):Date(2020, 3, 15)

d = (Lon(lons), Lat(lats), Time(t))
A = zeros([length(x) for x in (lons, lats, t)]...)
B = copy(A)

# generate solar rad
for i in 1:length(lats)
    θ = lats[i]
    s = [monthly_insolation(ti, θ) for ti in t]
    for j in 1:length(lons)
        c = cosd(lons[j])
        A[j, i, :] .= c.*s
        B[j, i, :] .= s
    end
end

A = ClimArray(A, d; name = "lon-variation")
B = ClimArray(B, d; name = "lon-constant")

# %%

@testset "Dropping dimensions" begin
    dt = timemean(A)
    @test dt isa ClimArray
    @test !hasdim(dt, Time)
    @test hasdim(dt, Lon)
    @test hasdim(dt, Lat)
    dl = zonalmean(A)
    @test dl isa ClimArray
    @test !hasdim(dl, Lon)
    @test hasdim(dl, Time)
    @test hasdim(dl, Lat)
    da = latmean(A)
    @test da isa ClimArray
    @test !hasdim(da, Lat)
    @test hasdim(da, Lon)
    @test hasdim(da, Time)
    ds = spacemean(A)
    @test ds isa ClimArray
    @test !hasdim(ds, Lon)
    @test !hasdim(ds, Lat)
    @test hasdim(ds, Time)
end

@testset "Default physical weights" begin
    x, y = hemispheric_means(B)
    @test timemean(x) != dropagg(mean, x, Time)
    d1 = timemean(x) - timemean(y)
    @test abs(d1) < 0.01
    d2 = dropagg(mean, x, Time) - dropagg(mean, y, Time)
    @test d2 < -0.1
    x = B[:, :, 1]
    @test zonalmean(latmean(x)) ≈ latmean(zonalmean(x)) ≈ spacemean(x)
    @test spacemean(x) - mean(x) > 50
    @test all(y -> abs(y) < 1e-8, spacemean(A))
end

@testset "Temporal weighting" begin
    # spatiotemporal weights:
    W = zero(A)
    W[Time(5)] .= 1.0
    res = timeagg(mean, A, W)
    @test all(res .≈ A[Time(5)])
    @test dims(A, Lon) == dims(res, Lon)
    res = timeagg(std, A, W)
    @test all(x -> isapprox(x, 0; atol = 1e-8), res)

    # just time weight
    w = zeros(length(t))
    w[5] = 1.0
    res = timeagg(mean, A, w)
    @test all(res .≈ A[Time(5)])
    @test dims(A, Lon) == dims(res, Lon)

    # TODO: more tests needed here, e.g. for timeagg(mean, t, a, w)
end

@testset "Advanced temporal manipulation" begin
    tdense = Date(2000, 3, 1):Day(1):Date(2020, 3, 31)
    tyearly = Date(2000, 3, 1):Year(1):Date(2020, 3, 31)
    mdates = unique!([(year(d), month(d)) for d in tdense])
    ydates = unique!([year(d) for d in tdense])
    tranges = temporalrange(tdense, Dates.month)
    yranges = temporalrange(tdense, Dates.year)
    @test length(tranges) == length(mdates)
    @test length(yranges) == length(ydates)

    @test temporal_sampling(t) == :monthly
    @test temporal_sampling(tdense) == :daily
    @test temporal_sampling(tyearly) == :yearly

    # test yearly temporal weights (i.e., no weighting)
    X = ClimArray(rand(3,3), (Lon(1:3), Time(tyearly[1:3])))
    W = [0, 1, 0]
    @test vec(timemean(X).data) == vec(mean(X.data; dims = 2))
    @test timemean(X, W).data == X.data[:, 2]

    C = ClimArray(zeros(length(lats), length(tdense)), (Lat(lats), Time(tdense)))

    # First version: just count number of days
    for j in 1:length(tdense)
        for i in 1:length(lats)
            C[i, j] = daysinmonth(tdense[j])
        end
    end
    Cm = monthlyagg(C)

    @test length(unique(Array(Cm))) == 4 # there are four unique number of days
    # test that each value when rounded to an integer is an integer (for first slice only
    # all remaining slices are the same)
    for e in Cm[:, 1]
        @test round(Int, e) == e
    end
    @test step(dims(Cm, Time).val) == Month(1)
    @test temporal_sampling(dims(Cm, Time).val) == :monthly

    for j in 1:length(tdense)
        for i in 1:length(lats)
            C[i, j] = daysinyear(tdense[j])
        end
    end
    Cy = yearlyagg(C)
    @test length(unique(Array(Cy))) == 2
    for e in Cy
        @test round(Int, e) == e
    end
    @test step(dims(Cy, Time).val) == Year(1)
    @test temporal_sampling(dims(Cy, Time).val) == :yearly

    # Second version: test actual physics
    for j in 1:length(tdense)
        for i in 1:length(lats)
            C[i, j] = insolation(tdense[j], lats[i])
        end
    end
    Bz = zonalmean(B)
    Cm = monthlyagg(C)
    @test all(Cm .≈ Bz)
end

@testset "Spatial weighting" begin
    W = zero(A)
    W[Lon(5)] .= 1.0
    res = spaceagg(mean, A, W)
    @test all(res .≈ latmean(A)[Lon(5)])
    @test dims(A, Time) == dims(res, Time)
    res = spaceagg(std, A, W)
    @test all(x -> x>0, res)

    W = zeros(length(lons), length(lats))
    W[5, :] .= 1.0
    W = ClimArray(W, dims(A, (Lon, Lat)))
    res = spaceagg(mean, A, W)
    @test all(res .≈ latmean(A)[Lon(5)])
    @test dims(A, Time) == dims(res, Time)
end

@testset "Insolation/Hemispheres" begin
    @test !any(x -> x < 0, B)
    x = timemean(B)
    x, y = hemispheric_functions(zonalmean(x))
    # @test x isa ClimArray
    x[Lat(1)] == y[Lat(1)]
    for j in 2:size(x, Lat)
        @test x[j] < x[j-1]
        @test y[j] < y[j-1]
    end
end


text:
using Dates

geolist = [Date(2019, 1, 1), Date(2019, 1, 3), Date(2019, 1, 6), Date(2019, 1, 10)]
intlist = [(geolist[i], geolist[j]) for i = 1:length(geolist)-1 for j = i+1:length(geolist)]


text:
function branchTree(quadTree)
  x = quadTree.boundary.x
  y = quadTree.boundary.y
  width = quadTree.boundary.width/2
  height = quadTree.boundary.height/2
  branch = quadTree.branch + 1

  trBounds = Rectangle(x + width, y - height, width, height)
  topRight = QuadTree(trBounds, quadTree.capacity, branch)

  brBounds = Rectangle(x + width, y + height, width, height)
  bottomRight = QuadTree(brBounds, quadTree.capacity, branch)

  blBounds = Rectangle(x - width, y + height, width, height)
  bottomLeft = QuadTree(blBounds, quadTree.capacity, branch)

  tlBounds = Rectangle(x - width, y - height, width, height)
  topLeft = QuadTree(tlBounds, quadTree.capacity, branch)

  return [topRight, bottomRight, bottomLeft, topLeft]
end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
using BcdV3
using StaticArrays
using MimuConstants
using LinearAlgebra
using Printf
using Dates

# Noise less simulations
TT = Float64
active_imus = [3,4,5,6, 27,28, 29, 30]
r_all = get_InitialPosMIMU4444BT_corrected()
r_m = r_all[:, active_imus]
r_ref = r_m[:,1]
r_m .-= r_ref
r_all .-= r_ref

sig_a = get_rms_noise_acc()
sig_g = get_rms_noise_gyro()
Na = size(r_m, 2)
Ng = 8

N = 5000 # Dyn
Ns = []
show_every = 100
N_bcd  = 1000
tol_dlogp = 1.0e-9

Nt = Na + Ng
g_mag = 9.81

b = [@SVector randn(3) for k = 1:Nt]
T = [I + 1.0e-3*@SMatrix randn(3,3) for k = 1:Nt]
T[1] = I + 1.0e-3*SMatrix{3,3, Float64}(triu(randn(3,3)))
r = [SVector{3,TT}(r_m[:,k]) for k = 1:Na]

w = [deg2rad(1000.0) * @SVector randn(3) for _ = 1:N]
w_dot = [deg2rad(40e3) * @SVector randn(3) for _ = 1:N]
s = [20.0 * @SVector randn(3) for _ = 1:N]
w_m = hcat(Array.(w)...)

y_d_true = Array{SVector{3,TT}}(undef, Nt,N)
y_d_true_2 = Array{SVector{3,TT}}(undef, Nt,N)
for n = 1:N
    for k = 1:Na
        y_d_true[k,n] = T[k]*(Ω(w[n])^2*r[k] + Ω(w_dot[n])*r[k] + s[n]) + b[k]
        y_d_true_2[k,n] = T[k]*(Ω(2*w[n])^2*r[k] + Ω(2*w_dot[n])*r[k] + 2*s[n]) + b[k]
    end
    for k in range(Na+1, length=Ng)
        y_d_true[k,n] = T[k]*w[n] + b[k]
        y_d_true_2[k,n] = T[k]*w[n]*2 + b[k]
    end
end

Q_inv = [SMatrix{3,3,Float64}(I)/( k <= Na ? sig_a^2 :  sig_g^2) for k = 1:Nt]

b0 = [b_n + 1.0e-2*randn(3) for  b_n in b]
T0 = [T[k] + 1.0e-3*(k == 1 ? triu(randn(3,3)) : randn(3,3))  for k = 1:Nt]
r0_m = copy(r_m)
r0_m[:,2:end] += 1.0e-3*randn(3,Na-1)
r0 = [SVector{3}(r0_m[:,k]) for k = 1:Na]

y = Array{SVector{3,TT}}(undef, Nt,N)
y2 = Array{SVector{3,TT}}(undef, Nt,N)
for n = 1:N
    for k = 1:Na
        e = sig_a * randn(3)
        y[k,n] = y_d_true[k,n] + e
        y2[k,n] = y_d_true_2[k,n] + e
    end
    for k in range(Na+1, length=Ng)
        e = sig_g * randn(3)
        y[k,n] = y_d_true[k,n] + e
        y2[k,n] = y_d_true_2[k,n] + e
    end
end
dw = zeros(3, N)
dw2 = zeros(3, N)

w0 = deepcopy(w)
r_hat, T_hat, b_hat, eta = bcd!(
    copy(r), copy(T), copy(b), y, Q_inv, Ns, g_mag, Val{3*(Na+Ng)}(),Val{Na}(),
    1.0e-12, 100, 1.0e-8, 200, tol_dlogp, 1.0-14, 1; warn = false, w0 = w0
)
w_hat_m = hcat(Array.(eta[:w])...)
dw .+= (w_m  - w_hat_m).^2
mse1 = sum(dw)/length(dw)
rmse1 = sqrt(mse1)

r_hat, T_hat, b_hat, eta = bcd!(
    copy(r), copy(T), copy(b), y2, Q_inv, Ns, g_mag, Val{3*(Na+Ng)}(),Val{Na}(),
    1.0e-12, 100, 1.0e-8, 200, tol_dlogp, 1.0-14, 1; warn = false, w0 = w0.*2
)
w_hat_m = hcat(Array.(eta[:w])...)
dw2 .+= (w_m.*2  - w_hat_m).^2
mse2 = sum(dw2)/length(dw2)
rmse2 = sqrt(mse2)

# rmse is lower for second


text:
function helper_test()
	@testset "Helpers" begin
        @test isequal(regressor_name(:gpa, 2), :afqt2)

	end
end

@testset "Helpers" begin
    helper_test()
end

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


text:
# TODO: column generation for all these functions.

function compute_routing(rd::RoutingData, ::Load, ::MinimumMaximum, ::FormulationType, ::Val{false}, ::Automatic, ::NoUncertaintyHandling, ::NoUncertainty)
    start = time_ns()

    # Create the problem.
    m = _create_model(rd)
    rm = basic_routing_model_unitary(m, rd)

    # Objective function.
    @variable(m, mu >= 0)
    rm.mu = mu
    @objective(m, Min, mu)

    # Constraints.
    mu_capacity_constraints(rm, rd.traffic_matrix)

    time_create_master_model_ms = (time_ns() - start) / 1_000_000.

    # Solve the problem.
    start = time_ns()
    optimize!(m)
    stop = time_ns()

    return RoutingSolution(rd,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_create_master_model_ms=time_create_master_model_ms,
                           time_solve_ms=(stop - start) / 1_000_000.,
                           objectives=objective_value(m),
                           matrices=rd.traffic_matrix,
                           routings=Routing(rd, value.(rm.routing)),
                           master_model=rm)
end

_warn(edge_obj::EdgeWiseObjectiveFunction, sense::String) =
    @warn "The objective function $edge_obj does not support $(sense)imisation. Proceed with caution."
_warn_min(edge_obj::EdgeWiseObjectiveFunction) =
    supports_min(edge_obj) || _warn(edge_obj, "min")
_warn_max(edge_obj::EdgeWiseObjectiveFunction) =
    supports_max(edge_obj) || _warn(edge_obj, "max")

_ensure_objective_compatibility(edge_obj::EdgeWiseObjectiveFunction, ::Union{MinimumTotal, MinimumMaximum}) = _warn_min(edge_obj)
_ensure_objective_compatibility(edge_obj::EdgeWiseObjectiveFunction, ::Union{MaximumTotal, MaximumMinimum}) = _warn_max(edge_obj)

function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::Union{MinimumTotal, MaximumTotal}, ::FormulationType, ::Val{false}, ::Automatic, ::NoUncertaintyHandling, ::NoUncertainty)
    _ensure_objective_compatibility(edge_obj, agg_obj)

    start = time_ns()

    # Create the problem.
    m = Model(rd.solver)
    set_silent(m)
    rm = basic_routing_model_unitary(m, rd)

    # Constraints.
    capacity_constraints(rm, rd.traffic_matrix)

    # Objective function.
    obj = sum(objective_edge_expression(rm, edge_obj, e) for e in edges(rd))
    if agg_obj == MinimumTotal()
        @objective(m, Min, obj)
    else
        @assert agg_obj == MaximumTotal()
        @objective(m, Max, obj)
    end

    # Done!
    optimize!(m)
    stop = time_ns()

    return RoutingSolution(rd,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_solve_ms=(stop - start) / 1_000_000.,
                           objectives=objective_value(m),
                           matrices=rd.traffic_matrix,
                           routings=Routing(rd, value.(rm.routing)),
                           master_model=rm)
end

function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::Union{MinimumMaximum, MaximumMinimum}, ::FormulationType, ::Val{false}, ::Automatic, ::NoUncertaintyHandling, ::NoUncertainty)
    _ensure_objective_compatibility(edge_obj, agg_obj)

    start = time_ns()

    # Create the problem.
    m = _create_model(rd)
    rm = basic_routing_model_unitary(m, rd)

    # Constraints.
    capacity_constraints(rm, rd.traffic_matrix)

    # Objective function. Don't define a lower or an upper bound, in case the
    # individual terms do not have such a bound (like α-fairness).
    @variable(m, obj)
    for e in edges(rd)
        expr = objective_edge_expression(rm, edge_obj, e)
        if agg_obj == MinimumMaximum()
            @constraint(m, obj >= expr)
        else
            @assert agg_obj == MaximumMinimum()
            @constraint(m, obj <= expr)
        end
    end

    if agg_obj == MinimumMaximum()
        @objective(m, Min, obj)
    else
        @assert agg_obj == MaximumMinimum()
        @objective(m, Max, obj)
    end

    # Done!
    optimize!(m)
    stop = time_ns()

    return RoutingSolution(rd,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_solve_ms=(stop - start) / 1_000_000.,
                           objectives=objective_value(m),
                           matrices=rd.traffic_matrix,
                           routings=Routing(rd, value.(rm.routing)),
                           master_model=rm)
end


text:
import .StaticArrays: StaticArrayStyle, StaticArray
mtb_copy(bc::Broadcasted{<:StaticArrayStyle}) = copy(bc) # I believe no one need it
mtb_copyto!(dest::StaticArray, bc::Broadcasted{Nothing}) = copyto!(dest, bc) # I believe no one need it
_similar(bc::Broadcasted{<:StaticArrayStyle}, T) = throw("StaticArrays is not supported for non-inplace @tab")


text:
module TestExamples1DEulerMulti

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "Compressible Euler Multicomponent" begin

  @trixi_testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [1.54891912e-01, 4.45430525e-01, 1.70222013e-02, 3.40444026e-02, 6.80888053e-02],
      linf = [2.71276498e-01, 9.95140742e-01, 3.93069410e-02, 7.86138820e-02, 1.57227764e-01])
  end

  @trixi_testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [1.53387916e-01, 4.41585576e-01, 3.93605635e-02, 7.87211270e-02],
      linf = [2.49632117e-01, 7.21088064e-01, 6.38328770e-02, 1.27665754e-01])
  end

  @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"),
      l2   = [8.57523604e-05, 1.63878043e-04, 1.94126993e-05, 3.88253986e-05],
      linf = [3.05932773e-04, 6.24480393e-04, 7.25312144e-05, 1.45062429e-04])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, 6.186482304747219e-6, 1.2372964609494437e-5],
      linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, 2.603347164065184e-5, 5.206694328130368e-5])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [1.88845048e-05, 5.49106005e-05, 9.42673716e-07, 1.88534743e-06, 3.77069486e-06, 7.54138973e-06],
      linf = [1.16223512e-04, 3.07922197e-04, 3.21774233e-06, 6.43548465e-06, 1.28709693e-05, 2.57419386e-05],
      volume_flux = flux_chandrashekar)
  end

  @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_two_interacting_blast_waves.jl"),
      l2   = [1.28886761e+00, 8.27133526e+01, 3.50680272e-03, 1.36987844e-02, 1.91795185e-02],
      linf = [2.96413045e+01, 1.32258448e+03, 9.19191937e-02, 3.10929710e-01, 4.41798976e-01],
      tspan = (0.0, 0.0001))
  end

end

end # module


text:
"Context `Filter` method. The default returns NEUTRAL."
filter{E <: Event}(flt::Filter, evnt::E) = FilterResult.NEUTRAL

"Appender `Filter` method. The default returns NEUTRAL."
filter(flt::Filter, level::Level.EventLevel, marker::MARKER, msg) = FilterResult.NEUTRAL

start(flt::FILTER) = !isnull(flt) && start(get(flt))
stop(flt::FILTER) = !isnull(flt) && stop(get(flt))

include("filters/composite.jl")
include("filters/marker.jl")
include("filters/threshold.jl")


text:
# File test/runtests.jl
module ScruffTest

using Test

@testset "ScruffTests" begin
    include("test_core.jl")
    include("test_sfuncs.jl")
    include("test_score.jl")
    include("test_utils.jl")
    include("test_ve.jl")
    include("test_lsfi.jl")
    include("test_bp.jl")
    include("test_importance.jl")
    include("test_net.jl")
    include("test_filter.jl")
    # These don't test anything, but we want to make sure any changes haven't broken the examples
    redirect_stdout(devnull) do
        include("../docs/examples/novelty_example.jl")
        include("../docs/examples/novelty_lazy.jl")
        include("../docs/examples/novelty_filtering.jl")
        include("../docs/examples/rembrandt_example.jl")
        include("../docs/examples/soccer_example.jl")
    end
end

end


text:
using Kelpie
using Documenter

DocMeta.setdocmeta!(Kelpie, :DocTestSetup, :(using Kelpie); recursive=true)

makedocs(;
    modules=[Kelpie],
    authors="Thomas A. Christensen II <25492070+MillironX@users.noreply.github.com> and contributors",
    repo="https://github.com/MillironX/Kelpie.jl/blob/{commit}{path}#{line}",
    sitename="Kelpie.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MillironX.github.io/Kelpie.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/MillironX/Kelpie.jl", devbranch="master")


text:
export Dim
export check_dim, multiplice_matrix_chain

struct Dim
    name::String
    row::Int64
    col::Int64
end

function check_dim(chain::Vector{Dim})
    for i = 1:(length(chain)-1)
        (chain[i].col != chain[i+1].row) && (return false)
    end

    return true
end

function to_dims(chain::Vector{Dim})
    dims = []
    for dim in chain
        push!(dims, dim.row)
    end
    push!(dims, chain[end].col)

    return dims
end

function print_chain(
    chain::Vector{Dim},
    references::Array,
    i::Int64, j::Int64;
    order=""
)
    if i == j
        order *= "$(chain[i].name)"
    else
        order *= "("
        order = print_chain(chain, references, i, references[i, j], order=order)
        order = print_chain(chain, references, references[i, j]+1, j, order=order)
        order *= ")"
    end

    return order
end

function multiplice_matrix_chain(chain::Vector{Dim})
    if !check_dim(chain) throw(DimensionMismatch) end

    dims = to_dims(chain)
    n = length(chain)

    m = zeros(Int, n, n)
    references = zeros(Int, n, n)

    for d = 2:n
        for i = 1:(n-d+1)
            j = i+d-1
            m[i, j] = typemax(Int64)
            for k = i:(j-1)
                temp = m[i, k] + m[k+1, j] + dims[i]*dims[k+1]*dims[j+1]
                (temp < m[i, j]) && (m[i, j] = temp; references[i, j] = k)
            end
        end
    end

    order = print_chain(chain, references, 1, n)
    complexity = m[1, n]

    return order, complexity
end


text:
using Images
using Dates

function init_CCobj()
  CCobj = CCstats(zeros(Int, 1,1,1),  # Bools
  zeros(Int, 1,1,1), # Connection Cube
  trues(3,3,3),
  0, # numExtremes
  0, # minAffectedVoxels
  [0], # labels
  [0], # perm
  [0], # affV
  [0], # affA
  [0], # affT
  [0.0], #affV
  [0.0], #affA [km^2]
  [0.0], # affT [days]
  [0.0], # Tmin
  [0.0], # Tmax
  [0.0], # Tcenters
  [0.0], # LATmin
  [0.0], # LATmax
  [0.0], # LATcenters
  [0.0], # LONmin
  [0.0], # LONmax
  [0.0], # LONcenters
  500.0:0.00001:500.0,  # Latitudes
  500.0:0.00001:500.0, # Longitudes
  [0.0], #Times
  [0.0],
  [0.0],
  Array{Dict{Tuple{Int64,Int64},Int64},1}(undef,0),
  Array{Dict{Tuple{Int64,Int64,Int64},Float64},1}(undef,0),
  Dict{String, VariableStats}() # Variable Statistics
  );
  return(CCobj)
end

export init_CCobj

function get_affectedTimes(CCobj::CCstats, TimeStepLength)
  CCobj.affectedTime = zeros(Float64, length(CCobj.Tmins))
  for i = 1:length(CCobj.Tmins)
   CCobj.affectedTime[i] = Float64(CCobj.Tmaxs[i] - CCobj.Tmins[i]) / (3600 * 24 * 1000) + TimeStepLength
  end
  CCobj.affectedTime
end

function getLonLatTimebounds(CCobj::CCstats)
  (minbounds, maxbounds) = getBounds(CCobj.LabelCube, CCobj.numExtremes)

for j = 1:3
  replace!(i -> i > size(CCobj.LabelCube, j) ? size(CCobj.LabelCube, j) : i, view(minbounds,j,:))
  replace!(i -> i < 1 ? 1 : i, view(minbounds,j,:))
  replace!(i -> i > size(CCobj.LabelCube, j) ? size(CCobj.LabelCube, j) : i, view(maxbounds,j,:))
  replace!(i -> i < 1 ? 1 : i, view(maxbounds,j,:))
end

  #println(extrema(minbounds[1,:]))
  CCobj.LONmins = CCobj.Longitudes[minbounds[1,:]]
  CCobj.LATmins = CCobj.Latitudes[minbounds[2,:]]
  CCobj.Tmins = CCobj.Times[minbounds[3,:]]

  CCobj.LONmaxs = CCobj.Longitudes[maxbounds[1,:]]
  CCobj.LATmaxs = CCobj.Latitudes[maxbounds[2,:]]
  CCobj.Tmaxs = CCobj.Times[maxbounds[3,:]]

  CCobj.affectedTimeSteps = maxbounds[3,:] .- minbounds[3,:] .+ 1

  return(minbounds, maxbounds)
end

function countVoxels(LabelCube::AbstractArray{<:Integer, N}, numExtremes::Int = 0) where {N}
  if numExtremes == 0
    numExtremes = maximum(LabelCube)
  end
  lAr=zeros(Int,numExtremes)
  for i=1:length(LabelCube)
    j=LabelCube[i]
    if j>0
      lAr[j]=lAr[j]+1
    end
  end
  return lAr
end

# get latitude index of a floatrange of latitudes
function getLatidx(Latitudes::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, lat::Float64)
    latidx = 0
    counter = 0
    for i in Latitudes
    counter += 1
        if i == lat
            latidx = counter
        end
    end
    if latidx == 0 error("lat = $lat not in Latitudes = $Latitudes") end
    return latidx
end

# get pixel Area of a given latitude
function pixelArea(CCobj::CCstats, lat::Float64)
    latidx = getLatidx(CCobj.Latitudes, lat)
    km_const = 111.324 #2pi / 360 * 6378.388 (Earth radius in km)
    latgridsize = abs(convert(Float64, CCobj.Latitudes.step))
    longridsize = abs(convert(Float64, CCobj.Longitudes.step))
    dlon = km_const * longridsize
    dlat = km_const * CCobj.CosLatitudes[latidx] * latgridsize
    Area = dlon * dlat
    return Area
end

function pixelArea(CCobj::CCstats, latidx::Int)
    #latidx = getLatidx(CCobj.Latitudes, lat)
    km_const = 111.324 #2pi / 360 * 6378.388 (Earth radius in km)
    latgridsize = abs(convert(Float64, CCobj.Latitudes.step))
    longridsize = abs(convert(Float64, CCobj.Longitudes.step))
    dlon = km_const * longridsize
    dlat = km_const * CCobj.CosLatitudes[latidx] * latgridsize
    Area = dlon * dlat
    return Area
end

function get_affectedLonLatCounter(CCobj::CCstats)
  CCobj.affectedLonLatIdxCounts = Array{Dict{Tuple{Int,Int}, Int},1}(undef, CCobj.numExtremes)
  get_affectedLonLatCounter!(CCobj.affectedLonLatIdxCounts, CCobj)
  return(CCobj.affectedLonLatIdxCounts)
end

export get_affectedLonLatCounter

function get_affectedLonLatCounter!(xout::Array{Dict{Tuple{Int,Int}, Int}, 1}, CCobj::CCstats)
  if CCobj.LabelCube == zeros(Int, 1, 1, 1) error("empty CCobj.LabelCube") end
  #CCobj.affectedLonLatIdxCounts = Array(Dict{Tuple{Int,Int}, Int}, CCobj.numExtremes)
  #xout = CCobj.affectedLonLatIdxCounts
  for k = 1:size(CCobj.LabelCube, 3)
    for j = 1:size(CCobj.LabelCube, 2)
      for i = 1:size(CCobj.LabelCube, 1)
        l = CCobj.LabelCube[i,j,k]
        if l > 0
          if isassigned(xout, l)
            if haskey(xout[l], (i, j))
              xout[l][(i,j)] += 1
            else
              xout[l][(i,j)] = 1
            end
          else
            xout[l] = Dict((i, j) => 1) # lon -> i, lat -> j
          end
        end
      end
    end
  end
  for i=1:length(xout)
    isassigned(xout,i) || (xout[i]=Dict((1,1)=>1))
  end
  return xout
end

export get_affectedLonLatCounter!

function get_affectedLonLatTimeScoreValues(CCobj::CCstats, Scores::AbstractArray{<:AbstractFloat, 3})
  CCobj.affectedLonLatTimeScoreValue = Array{Dict{Tuple{Int,Int,Int}, Float64},1}(undef, CCobj.numExtremes)
  get_affectedLonLatTimeScoreValues!(CCobj.affectedLonLatTimeScoreValue, CCobj.LabelCube, Scores)
  return(CCobj.affectedLonLatTimeScoreValue)
end

export get_affectedLonLatTimeScoreValues

function get_affectedLonLatTimeScoreValues!(xout::Array{Dict{Tuple{Int,Int,Int}, T},1}, LabelCube::AbstractArray{<:Integer,3}, Scores::AbstractArray{T, 3}; tshift::Int = 0) where {T}
  if LabelCube == zeros(Int, 1, 1, 1) error("empty CCobj.LabelCube") end
  for k = 1:size(LabelCube, 3)
    for j = 1:size(LabelCube, 2)
      for i = 1:size(LabelCube, 1)
        l = LabelCube[i,j,k]
        if l > 0
          if k+tshift <= size(LabelCube, 3) && k+tshift > 0
            if isassigned(xout, l)
              if haskey(xout[l], (i, j, k))
                error("double defined $i $j $k")
              else
                xout[l][(i,j,k)] = Scores[i,j,k+tshift]
              end
            else
              xout[l] = Dict((i, j,k) => Scores[i,j,k+tshift]) # lon -> i, lat -> j, times -> k
            end
          end
        end
      end
    end
  end
  return xout
end

export get_affectedLonLatTimeScoreValues!

function get_affectedArea!(CCobj::CCstats)

  CCobj.affectedLonLatIdxCounts = get_affectedLonLatCounter(CCobj);

  CCobj.affectedPixels = zeros(Int, length(CCobj.affectedLonLatIdxCounts))
  CCobj.affectedArea = zeros(Float64, length(CCobj.affectedLonLatIdxCounts))
  for i = 1:length(CCobj.affectedLonLatIdxCounts)
    CCobj.affectedPixels[i] = length(CCobj.affectedLonLatIdxCounts[i])
    for lonlats = keys(CCobj.affectedLonLatIdxCounts[i])
      (lonidx, latidx) = lonlats
      CCobj.affectedArea[i] = CCobj.affectedArea[i] + CCobj.AreaLatitudes[latidx]
    end
  end
  return (CCobj.affectedArea, CCobj.affectedPixels)
end

function get_affectedAreaAndVolume!(CCobj::CCstats)

  CCobj.affectedLonLatIdxCounts = get_affectedLonLatCounter(CCobj);

  #@show CCobj.affectedLonLatIdxCounts
  timesteplength = median(CCobj.Times[2:end]-CCobj.Times[1:(end-1)])

  CCobj.affectedPixels = zeros(Int, length(CCobj.affectedLonLatIdxCounts))
  CCobj.affectedArea = zeros(Float64, length(CCobj.affectedLonLatIdxCounts))
  CCobj.affectedVolume = zeros(Float64, length(CCobj.affectedLonLatIdxCounts))
  CCobj.affectedVoxels = zeros(Int, length(CCobj.affectedLonLatIdxCounts))
  for i = 1:length(CCobj.affectedLonLatIdxCounts)
    CCobj.affectedPixels[i] = length(CCobj.affectedLonLatIdxCounts[i])
    for lonlats = keys(CCobj.affectedLonLatIdxCounts[i])
      (lonidx, latidx) = lonlats
      #print(" -- lonidx = $lonidx, latidx = $latidx")
      CCobj.affectedArea[i] = CCobj.affectedArea[i] + CCobj.AreaLatitudes[latidx]
      CCobj.affectedVoxels[i] = CCobj.affectedVoxels[i] + CCobj.affectedLonLatIdxCounts[i][lonlats]
      CCobj.affectedVolume[i] = CCobj.affectedVolume[i] + CCobj.AreaLatitudes[latidx] * CCobj.affectedLonLatIdxCounts[i][lonlats] * timesteplength
    end
  end
  return (CCobj.affectedArea, CCobj.affectedPixels, CCobj.affectedVoxels , CCobj.affectedVolume)
end

export get_affectedAreaAndVolume!

# get fill CCobj.AreaLatitudes wit the calculated areas
function getAreaLatitudes(CCobj::CCstats)
    CCobj.AreaLatitudes = zeros(Float64, length(CCobj.Latitudes))
    for i = 1:length(CCobj.Latitudes)
        CCobj.AreaLatitudes[i] = pixelArea(CCobj, i)
    end
    return(CCobj.AreaLatitudes)
end

function countVoxels_weighted(CCobj::CCstats)
  if CCobj.numExtremes == 0
    CCobj.numExtremes = maximum(CCobj.LabelCube)
  end
  lAr=zeros(Float64,CCobj.numExtremes)
  for t=1:size(CCobj.LabelCube, 3)
    for la=1:size(CCobj.LabelCube, 2)
      for lo=1:size(CCobj.LabelCube, 1)
        j=CCobj.LabelCube[lo, la, t]
        if j>0
          if t+1 > length(CCobj.Times)
            timesteplength = Float64(DateTime(year(CCobj.Times[t])+1,1,1)-CCobj.Times[t]) / (24 * 3600 * 1000)
          else
            timesteplength = Float64(CCobj.Times[t+1]-CCobj.Times[t]) / (24 * 3600 * 1000)
          end
          lAr[j] = lAr[j] + (CCobj.AreaLatitudes[la] * timesteplength)
        end
      end
    end
  end
  return lAr
end

function get_affectedVolume!(CCobj::CCstats)
  if CCobj.numExtremes == 0
    CCobj.numExtremes = maximum(CCobj.LabelCube)
  end
  timesteplength = median(CCobj.Times[2:end]-CCobj.Times[1:(end-1)])
  lAr=zeros(Float64,CCobj.numExtremes)
  for t=1:size(CCobj.LabelCube, 3)
    for la=1:size(CCobj.LabelCube, 2)
      for lo=1:size(CCobj.LabelCube, 1)
        j=CCobj.LabelCube[lo, la, t]
        if j>0
          lAr[j] = lAr[j] + (CCobj.AreaLatitudes[la] * timesteplength)
        end
      end
    end
  end
  CCobj.affectedVolume = lAr
  return CCobj.affectedVolume
end

export get_affectedVolume!

function any_loop(sel_labels::AbstractArray{tp, N}, testlabel::tp) where {tp, N}
    @inbounds for j = 1:length(sel_labels)
        if sel_labels[j] == testlabel
            return  true
        end
    end
    return false
end

# loop wrapper around any(sel_labels .!= testlabel)
function any_notloop(sel_labels::AbstractArray{tp, N}, testlabel::tp) where {tp, N}
    @inbounds for j = 1:length(sel_labels)
        if sel_labels[j] != testlabel
            return  true
        end
    end
    return false
end

# Set Bools to zero which are not within the sel_labels of the LabelCube
function clean_bools!(Bools::Array{<:Integer, N}, LabelCube::Array{<:Integer, N}, sel_labels::Array{<:Integer, 1}) where {N}
    for i = 1:length(Bools)
        if Bools[i] != 0
            if !any_loop(sel_labels, LabelCube[i])
                Bools[i] = 0
            end
        end
    end
    return(Bools)
end

# get a connection cube cleaned for the desired minimum of affected voxels (default) or the desired number of extremes
function getLabelCube!(Bools::Array{<:Integer, N}, ConnectionFootprint::Array{<:Integer, N}; mode::String = "Voxels", minAffectedVoxels::Int = 1, numExtremes::Int = 100) where {N}
    #ConnectionFootprint = trues(3,3,3) # consider all 3 * 3 * 3 = 27 connections
    # first run
    CC = label_components(Bools, ConnectionFootprint)
    affVoxels = countVoxels(CC)
    # second run to exclude low connectivity
    if mode == "Voxels"
        sel_labels = find(affVoxels .>= minAffectedVoxels)
    else
        perm = sortperm(affVoxels, rev = true)
        sel_labels = perm[1:numExtremes]
    end
    clean_bools!(Bools, CC, sel_labels)
    CC2 = label_components(Bools, ConnectionFootprint)
    return(Bools, CC2)
end

export getLabelCube!

function getLabelCube!(CCobj::CCstats)
    if CCobj.numExtremes == 0 && CCobj.minAffectedVoxels == 0 error("either numExtremes or minAffectedVoxels has to be != 0") end
    if CCobj.Bools == zeros(Int, 1,1,1) error("you forgot to write your Boolean cube to CCobj.Bools") end
    if all(CCobj.ConnectionFootprint .== 0)  error("you do not consider any connections in CCobj.ConnectionFootprint") end
    #ConnectionFootprint = trues(3,3,3) # consider all 3 * 3 * 3 = 27 connections
    # first run
    print("first run to compute connected components \n")
    @time CCobj.LabelCube = label_components(CCobj.Bools, CCobj.ConnectionFootprint)
    print("compute affected voxels \n")
    @time affVoxels = countVoxels(CCobj.LabelCube)
    print("get permutation \n")
    perm = sortperm(affVoxels, rev = true)
    # second run to exclude low connectivity
    if CCobj.minAffectedVoxels != 0
      print("select only events which are connected to more than $(CCobj.minAffectedVoxels) voxels \n")
        @time sel_labels = findall(affVoxels .>= CCobj.minAffectedVoxels)
    else
      print("select the $(CCobj.numExtremes) Events with the largest number of affected Voxels \n")
        @time begin
        sel_labels = perm[1:CCobj.numExtremes]
        end
    end
    CCobj.perm=perm[1:length(sel_labels)]
    print("remove the other events from Bools \n")
    @time clean_bools!(CCobj.Bools, CCobj.LabelCube, sel_labels)
    print("clean and sort label cube with cleaned bools \n")
    @time CCobj.LabelCube = label_components(CCobj.Bools, CCobj.ConnectionFootprint)
    #cleanAndSort_LabelCube!(CCobj, perm)
    CCobj.numExtremes = maximum(CCobj.LabelCube)
    return(CCobj.LabelCube)
end

export getLabelCube!

function cleanAndSort_LabelCube!(CCobj::CCstats, perm::Array{<:Integer, 1})
  #newlabelCube = similar(CCobj.LabelCube)
  for i = 1:length(CCobj.LabelCube)
    if CCobj.Bools[i] == 0 # Bools are already cleaned --> if zero label cube is also zero
      CCobj.LabelCube[i] = 0
    else
      CCobj.LabelCube[i] = findfirst(isequal(CCobj.LabelCube[i]), perm)
    end
  end
  return CCobj.LabelCube
end

# get minimum and maximum values of each dimension for each label of the LabelCube
function getBounds!(minbounds, maxbounds, LabelCube::Array{<:Integer, N}, numExtremes::Int = 0) where {N}
    if numExtremes == 0
     numExtremes = maximum(LabelCube)
    end
    for k = 1:size(LabelCube, 3)
        for j = 1:size(LabelCube, 2)
            for i = 1:size(LabelCube, 1)
                signat = LabelCube[i,j,k]
                if signat > 0
                    if minbounds[1, signat] > i  minbounds[1, signat] = i end
                    if minbounds[2, signat] > j  minbounds[2, signat] = j end
                    if minbounds[3, signat] > k  minbounds[3, signat] = k end
                    if maxbounds[1, signat] < i  maxbounds[1, signat] = i end
                    if maxbounds[2, signat] < j  maxbounds[2, signat] = j end
                    if maxbounds[3, signat] < k  maxbounds[3, signat] = k end
                end
            end
        end
    end
    return(minbounds, maxbounds)
end

function getBounds(LabelCube::Array{<:Integer, N}, numExtremes::Int = 0) where {N}
    if numExtremes == 0
     numExtremes = maximum(LabelCube)
    end
    minbounds=fill(maximum(size(LabelCube)) + 1,3,numExtremes)
    maxbounds=zeros(Int,3,numExtremes)
    getBounds!(minbounds, maxbounds, LabelCube, numExtremes)
    return(minbounds, maxbounds)
end

function getConCompStats(CCobj::CCstats)
  if CCobj.Times == [0.0] error("CCobj.Times is empty") end
  if minimum(CCobj.Latitudes) == 500.0 error("CCobj.Latitudes is undefined") end
  if minimum(CCobj.Longitudes) == 500.0 error("CCobj.Longitudes is undefined") end
  if CCobj.LabelCube == zeros(Int, 1,1,1) error("CCobj.LabelCube is empty") end

  CCobj.CosLatitudes = cos.(collect(CCobj.Latitudes) / 360.0 * 2pi)
  CCobj.AreaLatitudes = getAreaLatitudes(CCobj)

  CCobj.numExtremes = maximum(CCobj.LabelCube)

  #print("get affected voxels, volume, area, pixels \n")
  get_affectedAreaAndVolume!(CCobj)

  #print("get permutattion according to volume \n")
  CCobj.perm = sortperm(CCobj.affectedVolume, rev = true);
  CCobj.labels = collect(1:CCobj.numExtremes)

  #print("get minima and maxima of lon lat time \n")
  getLonLatTimebounds(CCobj)

  #print("get duration")
  TimeStepLength = median(CCobj.Times[2:end]-CCobj.Times[1:(end-1)])
  get_affectedTimes(CCobj, TimeStepLength)

  return CCobj
end

export getConCompStats

function getCentroids(CCobj::CCstats)
  if length(CCobj.affectedLonLatTimeScoreValue) == 0 error("CCobj.affectedLonLatTimeScoreValue is empty") end
  CCobj.LONcenters = zeros(Float64, CCobj.numExtremes)
  CCobj.LATcenters = zeros(Float64, CCobj.numExtremes)
  CCobj.Tcenters = zeros(Float64, CCobj.numExtremes)
  getCentroids!(CCobj.LONcenters, CCobj.LATcenters, CCobj.Tcenters, CCobj.affectedLonLatTimeScoreValue, CCobj)
  return CCobj.LONcenters, CCobj.LATcenters, CCobj.Tcenters
end

export getCentroids

function getCentroids!(LonCenters::Array{Float64, 1}, LatCenters::Array{Float64, 1}, TimeCenters::Array{Float64, 1}, LonLatTime::Array{Dict{Tuple{Int64,Int64,Int64},Float64},1}, CCobj::CCstats; ValueWeighting::Bool = true)
  sweights=zeros(Float64, CCobj.numExtremes)
  valweights=zeros(Float64, CCobj.numExtremes)
  for i = 1:length(LonLatTime)
    if isassigned(LonLatTime, i)
        for lonlattimeidx = keys(LonLatTime[i])
          if ValueWeighting
            val = LonLatTime[i][lonlattimeidx]
          else
            val = 1
          end
          innerstCentroids!(sweights, valweights, LatCenters, LonCenters, TimeCenters, CCobj, lonlattimeidx, val, i)
        end
    end
  end
  broadcast!(/ ,LatCenters, LatCenters, sweights)
  broadcast!(/, LonCenters, LonCenters, valweights)
  broadcast!(/ ,TimeCenters, TimeCenters, valweights)
  return LonCenters, LatCenters, TimeCenters
end

export getCentroids!

function innerstCentroids!(sweights, valweights, LatCenters, LonCenters, TimeCenters, CCobj, lonlattimeidx, val_orig, i)
  (lonidx, latidx, timeidx) = lonlattimeidx
  # weighting only works with positive values
  val = abs(val_orig)
  if lonidx < 1  lonidx = 1 end
  if latidx < 1  latidx = 1 end
  if timeidx < 1  timeidx = 1 end
  if lonidx > size(CCobj.LabelCube, 1)  lonidx = size(CCobj.LabelCube, 1) end
  if latidx > size(CCobj.LabelCube, 2)  latidx = size(CCobj.LabelCube, 2) end
  if timeidx > size(CCobj.LabelCube, 3)  timeidx = size(CCobj.LabelCube, 3) end
  LatCenters[i] = LatCenters[i] + CCobj.Latitudes[latidx]*CCobj.CosLatitudes[latidx]*val
  sweights[i]  += CCobj.CosLatitudes[latidx]*val
  LonCenters[i] = LonCenters[i] + CCobj.Longitudes[lonidx]*val
  TimeCenters[i] = TimeCenters[i] + CCobj.Times[timeidx]*val
  valweights[i] += val
  return sweights[i], valweights[i], LatCenters[i], LonCenters[i], TimeCenters[i]
end


text:
"""
`filelist(dir; join=true, sort=true)` return the list of files. For its keyword arguments, see `readdir`.
"""
function filelist(dir; join=true, sort=true)
    fs = readdir(dir; join=true, sort=sort) # files and directories
    if join # == true
        return fs[isfile.(fs)]
    else
        return basename.(fs[isfile.(joinpath.(dir,fs))])
    end
end

"""
`filelist(; join=true, sort=true)` returns the paths for all the files in the current directory.
"""
function filelist(;opts...)
    filelist(pwd(); opts...)
end

"""
`filelist(expr::Regex, dir; join=true, sort=true)` returns a vector of paths where the file matches `expr`.
"""
function filelist(expr::Regex, dir; opts...)
    fs = filelist(dir; opts...)
    target_id = occursin.(expr, basename.(fs))
    return fs[target_id]
end


text:
using SuyamaBayes
using Test

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

max_stars_repo_name: SosUts/SuyamaBayes

max_issues_repo_name: SosUts/SuyamaBayes

max_forks_repo_name: SosUts/SuyamaBayes


text:
"""
    {{ class }}({{ params }})

Lorem ipsum dolor sit amet.
"""
type {{ class }}

  {{fields}}

end


text:
## ################################################################
# constructors
###################################################################

"""
    calculateToeplitzKernel(shape, tr::Matrix{T}[; m = 4, σ = 2.0, window = :kaiser_bessel, fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.ESTIMATE), kwargs...])

Calculate the kernel for an implementation of the Gram matrix that utilizes its Toeplitz structure. The output is an array of twice the size of `shape`, as the Toeplitz trick requires an oversampling factor of 2 (cf. [Wajer, F. T. A. W., and K. P. Pruessmann. "Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories." Proc. Intl. Soc. Mag. Res. Med. 2001.](https://cds.ismrm.org/ismrm-2001/PDF3/0767.pdf)). The type of the kernel is `Complex{T}`, i.e. the complex of the k-space trajectory's type; for speed and memory efficiecy, call this function with `Float32.(tr)`, and the kernel will also be `Float32`.

# Required Arguments
- `shape::NTuple(Int)`: size of the image; e.g. `(256, 256)` for 2D imaging, or `(256,256,128)` for 3D imaging
- `tr::Matrix{T}`: non-Cartesian k-space trajectory in units revolutions/voxel, i.e. `tr[i] ∈ [-0.5, 0.5] ∀ i`. The matrix has the size `2 x Nsamples` for 2D imaging with a trajectory length `Nsamples`, and `3 x Nsamples` for 3D imaging.

# Optional Arguments:
- `m::Int`: nfft kernel size (used to calculate the Toeplitz kernel); `default = 4`
- `σ::Number`: nfft oversampling factor during the calculation of the Toeplitz kernel; `default = 2`
- `window::Symbol`: Window function of the nfft (c.f. [`getWindow`](@ref)); `default = :kaiser_bessel`
- `K::Int`: `default= 2000` # TODO: describe meaning of k

# Keyword Arguments:
- `fftplan`: plan for the final FFT of the Kernal from image to k-space. Therefore, it has to have twice the size of the original image. `default = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.ESTIMATE)`. If this constructor is used many times, it is worth to precompute the plan with `plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.MEASURE)` and reuse it.
- `kwargs`: Passed on to [`plan_fft!`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.plan_fft!) via `NFFTPlan`; can be used to modify the flags `flags=FFTW.ESTIMATE, timelimit=Inf`.

# Examples
```jldoctest; output = false, setup = :(using NFFT, Random; Random.seed!(1))
julia> Nx = 32;

julia> trj = Float32.(rand(2, 1000) .- 0.5);

julia> λ = calculateToeplitzKernel((Nx, Nx), trj)
64×64 Matrix{ComplexF32}:
 -432.209+297.833im    271.68-526.348im  …   3044.35+34.6322im
  990.057-244.569im  -582.825+473.084im      45.7631-87.8965im
  1922.01+223.525im   9248.63-452.04im       3649.78+108.941im
 -402.371-3.37619im   496.815+231.891im      322.646-329.089im
 -482.346-121.534im   559.756-106.981im      155.183+454.0im
   3293.8+194.388im  -3361.43+34.1272im  …   2672.16-526.853im
 -936.331+172.246im   4394.11-400.762im     -1121.04+160.219im
 -135.828+169.448im   3509.17+59.0678im      3883.84-501.913im
  395.143+158.638im   24.4377-387.153im       5731.3+173.827im
  925.902-117.765im    2935.3+346.28im      -1414.12-214.701im
         ⋮                               ⋱
  2239.72-295.883im   490.442+524.399im  …   2028.84-36.5824im
 -1108.11+227.146im   24.7403-455.661im     -549.699+105.319im
  1323.78+110.713im  -321.052+117.802im      651.944-443.178im
 -52.1597+288.0im    -326.042-516.516im       3619.1+44.4651im
  1180.56+73.3666im  -26.1233+155.148im     -869.065-405.832im
  3555.36+649.527im  -198.245-878.042im  …   1198.83-317.062im
 -368.958-177.954im  -360.343+406.469im     -1478.96-154.512im
  4861.29+38.9623im   6082.55-267.478im      2519.09+293.503im
  1022.55-185.869im   177.426+414.384im      3650.56-146.597im

julia> x = randn(ComplexF32, Nx, Nx);

julia> convolveToeplitzKernel!(x, λ)
32×32 Matrix{ComplexF32}:
  177.717-52.0493im   10.6059+20.7185im  …   746.131+330.005im
 -311.858+988.219im  -1216.83-1295.14im      410.732+751.925im
 -1082.27+872.968im   1023.97-1556.45im     -923.699+478.63im
 -999.035-525.161im  -1694.59-658.558im     -521.044+607.005im
 -342.999+1481.82im  -1729.18-2712.81im      56.5212+1394.81im
 -1187.65-1979.55im   1389.67+970.033im  …   2968.93-744.264im
 -666.533+485.511im  -1315.83+409.855im     -610.146+132.258im
  1159.57+64.6059im  -299.169-569.622im      663.802+396.827im
 -795.112-1464.63im   -462.43+2442.77im     -1622.72+1701.27im
 -1047.67+31.5578im  -1127.65-936.043im      474.071+797.911im
         ⋮                               ⋱
  327.345+2191.71im  -904.635+558.786im     -1449.32-276.721im
  -1047.2-71.362im    363.109+567.346im      -1974.0-2348.36im
 -1540.45+1661.77im   1175.75+1279.75im  …   1110.61+653.234im
 -526.832-435.297im  -265.021-2019.08im      68.5607-323.086im
 -1076.52-2719.16im  -477.005+2232.06im      -155.59-1275.66im
  1143.76-735.966im  -380.489+2485.78im      1812.17-261.191im
  1685.44-1243.29im   1911.16-1157.72im     -991.639-42.8214im
  1054.11-1282.41im   66.9358-588.991im  …  -952.238+1026.35im
 -417.276-273.367im  -946.698+1971.77im     -890.339-882.05im

```
"""
function calculateToeplitzKernel(shape, tr::AbstractMatrix{T}; m = 4, σ = 2.0, window = :kaiser_bessel, fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.ESTIMATE), kwargs...) where {T}

    shape_os = 2 .* shape

    p = plan_nfft(typeof(tr), tr, shape_os; m, σ, window, kwargs...)
    eigMat = adjoint(p) * OnesVector(Complex{T}, size(tr,2))
    return fftplan * fftshift(eigMat)
end

"""
    calculateToeplitzKernel!(f::Array{Complex{T}}, p::AbstractNFFTPlan, tr::Matrix{T}, fftplan)

Calculate the kernel for an implementation of the Gram matrix that utilizes its Toeplitz structure and writes it in-place in `f`, which has to be twice the size of the desired image matrix, as the Toeplitz trick requires an oversampling factor of 2 (cf. [Wajer, F. T. A. W., and K. P. Pruessmann. "Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories." Proc. Intl. Soc. Mag. Res. Med. 2001.](https://cds.ismrm.org/ismrm-2001/PDF3/0767.pdf)). The type of the kernel `f` has to be `Complex{T}`, i.e. the complex of the k-space trajectory's type; for speed and memory efficiecy, call this function with `Float32.(tr)`, and set the type of `f` accordingly.

# Required Arguments
- `f::Array{T}`: Array in which the kernel will be written.
- `p::AbstractNFFTPlan`: NFFTPlan with the same dimentions as `tr`, which will be overwritten in place.
- `tr::Matrix{T}`: non-Cartesian k-space trajectory in units revolutions/voxel, i.e. `tr[i] ∈ [-0.5, 0.5] ∀ i`. The matrix has the size `2 x Nsamples` for 2D imaging with a trajectory length `Nsamples`, and `3 x Nsamples` for 3D imaging.
- `fftplan`: plan for the final FFT of the Kernal from image to k-space. Therefore, it has to have twice the size of the original image. Calculate, e.g., with `fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.MEASURE)`, where `shape` is the size of the reconstructed image.

# Examples
```jldoctest; output = false, setup = :(using NFFT, Random; Random.seed!(1))
julia> using FFTW

julia> Nx = 32;

julia> trj = Float32.(rand(2, 1000) .- 0.5);

julia> p = plan_nfft(trj, (2Nx,2Nx))
NFFTPlan with 1000 sampling points for an input array of size(64, 64) and an output array of size(1000,) with dims 1:2

julia> fftplan = plan_fft(zeros(ComplexF32, (2Nx,2Nx)); flags=FFTW.MEASURE);

julia> λ = Array{ComplexF32}(undef, 2Nx, 2Nx);

julia> calculateToeplitzKernel!(λ, p, trj, fftplan);

julia> x = randn(ComplexF32, Nx, Nx);

julia> convolveToeplitzKernel!(x, λ);

```
"""
function calculateToeplitzKernel!(f::Array{Complex{T}}, p::AbstractNFFTPlan{T}, tr::Matrix{T}, fftplan) where T
    nodes!(p, tr)
    f = mul!(f, adjoint(p), OnesVector(Complex{T}, size(tr,2)))
    f2 = fftshift(f)
    mul!(f, fftplan, f2)
    return f
end

## ################################################################
# constructor for explicit/exact calculation of the Toeplitz kernel
# (slow)
###################################################################
function calculateToeplitzKernel_explicit(shape, tr::Matrix{T}, fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.ESTIMATE)) where T

    shape_os = 2 .* shape
    λ = Array{Complex{T}}(undef, shape_os)
    Threads.@threads for i ∈ CartesianIndices(λ)
        λ[i] = getMatrixElement(i, shape_os, tr)
    end
    return fftplan * fftshift(λ)
end

function getMatrixElement(idx, shape::Tuple, nodes::Matrix{T}) where T
    elem = zero(Complex{T})
    shape = T.(shape) # ensures the correct output type

    @fastmath @simd for i ∈ eachindex(view(nodes, 1, :))
        ϕ = zero(T)
        @inbounds for j ∈ eachindex(shape)
            ϕ += nodes[j,i] * (shape[j]/2 + 1 - idx[j])
        end
        elem += exp(-2im * π * ϕ)
    end
    return elem
end

## ################################################################
# apply Toeplitz kernel
###################################################################

"""
    convolveToeplitzKernel!(x::Array{T,N}, λ::Array{T,N}[, fftplan = plan_fft(λ; flags=FFTW.ESTIMATE), ifftplan = plan_ifft(λ; flags=FFTW.ESTIMATE), xOS1 = similar(λ), xOS2 = similar(λ)])

Convolves the image `x` with the Toeplitz kernel `λ` and overwrites `x` with the result. `x` is also returned for convenience. As this function is commonly applied many times, it is highly recommened to pre-allocate / pre-compute all optional arguments. By doing so, this entire function is non-allocating.

# Required Arguments
- `x::Array{T,N}`: Input image that will be overwritten with the result. `x` is a matrix (`N=2`) for 2D imaging and a 3D tensor (`N=3`) for 3D imaging. The type of the elments `T` must match the ones of `λ`.
- `λ::Array{T,N}`: Toeplitz kernel, which as to be the same type as `x`, but twice the size due to the required oversampling (cf. [`calculateToeplitzKernel`](@ref)).

# Optional, but highly recommended Arguments
- `fftplan`: plan for the oversampled FFT, i.e. it has to have twice the size of the original image. Calculate, e.g., with `fftplan = plan_fft(λ; flags=FFTW.MEASURE)`, where `shape` is the size of the reconstructed image.
- `ifftplan`: plan for the oversampled inverse FFT. Calculate, e.g., with `ifftplan = plan_ifft(λ; flags=FFTW.MEASURE)`, where `shape` is the size of the reconstructed image.
- `xOS1`: pre-allocated array of the size of `λ`. Pre-allocate with `xOS1 = similar(λ)`.
- `xOS2`: pre-allocated array of the size of `λ`. Pre-allocate with `xOS2 = similar(λ)`.

# Examples
```jldoctest; output = false, setup = :(using NFFT, Random; Random.seed!(1))
julia> using FFTW

julia> Nx = 32;

julia> trj = Float32.(rand(2, 1000) .- 0.5);

julia> λ = calculateToeplitzKernel((Nx, Nx), trj);

julia> xOS1 = similar(λ);

julia> xOS2 = similar(λ);

julia> fftplan = plan_fft(xOS1; flags=FFTW.MEASURE);

julia> ifftplan = plan_ifft(xOS1; flags=FFTW.MEASURE);

julia> x = randn(ComplexF32, Nx, Nx);

julia> convolveToeplitzKernel!(x, λ, fftplan, ifftplan, xOS1, xOS2)
32×32 Matrix{ComplexF32}:
  177.717-52.0493im   10.6059+20.7185im  …   746.131+330.005im
 -311.858+988.219im  -1216.83-1295.14im      410.732+751.925im
 -1082.27+872.968im   1023.97-1556.45im     -923.699+478.63im
 -999.035-525.161im  -1694.59-658.558im     -521.044+607.005im
 -342.999+1481.82im  -1729.18-2712.81im      56.5212+1394.81im
 -1187.65-1979.55im   1389.67+970.033im  …   2968.93-744.264im
 -666.533+485.511im  -1315.83+409.855im     -610.146+132.258im
  1159.57+64.6059im  -299.169-569.622im      663.802+396.827im
 -795.112-1464.63im   -462.43+2442.77im     -1622.72+1701.27im
 -1047.67+31.5578im  -1127.65-936.043im      474.071+797.911im
         ⋮                               ⋱
  327.345+2191.71im  -904.635+558.786im     -1449.32-276.721im
  -1047.2-71.362im    363.109+567.346im      -1974.0-2348.36im
 -1540.45+1661.77im   1175.75+1279.75im  …   1110.61+653.234im
 -526.832-435.297im  -265.021-2019.08im      68.5607-323.086im
 -1076.52-2719.16im  -477.005+2232.06im      -155.59-1275.66im
  1143.76-735.966im  -380.489+2485.78im      1812.17-261.191im
  1685.44-1243.29im   1911.16-1157.72im     -991.639-42.8214im
  1054.11-1282.41im   66.9358-588.991im  …  -952.238+1026.35im
 -417.276-273.367im  -946.698+1971.77im     -890.339-882.05im

```
"""
function convolveToeplitzKernel!(x::Array{T,N}, λ::Array{T,N},
    fftplan = plan_fft(λ; flags=FFTW.ESTIMATE),
    ifftplan = plan_ifft(λ; flags=FFTW.ESTIMATE),
    xOS1 = similar(λ),
    xOS2 = similar(λ)
    ) where {T,N}

    fill!(xOS1, 0)
    xOS1[CartesianIndices(x)] .= x
    mul!(xOS2, fftplan, xOS1)
    xOS2 .*= λ
    mul!(xOS1, ifftplan, xOS2)
    x .= @view xOS1[CartesianIndices(x)]
    return x
end

## ################################################################
# helper class
###################################################################
struct OnesVector{T} <: AbstractVector{T}
    elements::T
    length::Int
end

function OnesVector(T::Type, length::Int)
    return OnesVector(one(T), length)
end

function Base.size(A::OnesVector)
    return (A.length,)
end

function Base.length(A::OnesVector)
    return A.length
end

function Base.getindex(A::OnesVector, i::Int)
    return A.elements
end


text:
@testset "FaceValues" begin
for (func_interpol, quad_rule) in  (
                                    (Lagrange{1, RefCube, 1}(), QuadratureRule{0, RefCube}(2)),
                                    (Lagrange{1, RefCube, 2}(), QuadratureRule{0, RefCube}(2)),
                                    (Lagrange{2, RefCube, 1}(), QuadratureRule{1, RefCube}(2)),
                                    (Lagrange{2, RefCube, 2}(), QuadratureRule{1, RefCube}(2)),
                                    (Lagrange{2, RefTetrahedron, 1}(), QuadratureRule{1, RefTetrahedron}(2)),
                                    (Lagrange{2, RefTetrahedron, 2}(), QuadratureRule{1, RefTetrahedron}(2)),
                                    (Lagrange{3, RefCube, 1}(), QuadratureRule{2, RefCube}(2)),
                                    (Serendipity{2, RefCube, 2}(), QuadratureRule{1, RefCube}(2)),
                                    (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)),
                                    (Lagrange{3, RefTetrahedron, 2}(), QuadratureRule{2, RefTetrahedron}(2))
                                   )

    for fe_valtype in (FaceScalarValues, FaceVectorValues)
        fv = fe_valtype(quad_rule, func_interpol)
        ndim = JuAFEM.getdim(func_interpol)
        n_basefuncs = getnbasefunctions(func_interpol)

        fe_valtype == FaceScalarValues && @test getnbasefunctions(fv) == n_basefuncs
        fe_valtype == FaceVectorValues && @test getnbasefunctions(fv) == n_basefuncs * JuAFEM.getdim(func_interpol)

        xs, n = valid_coordinates_and_normals(func_interpol)
        for face in 1:getnfaces(func_interpol)
            reinit!(fv, xs, face)
            @test JuAFEM.getcurrentface(fv) == face

            # We test this by applying a given deformation gradient on all the nodes.
            # Since this is a linear deformation we should get back the exact values
            # from the interpolation.
            u = Vec{ndim, Float64}[zero(Tensor{1,ndim}) for i in 1:n_basefuncs]
            u_scal = zeros(n_basefuncs)
            H = rand(Tensor{2, ndim})
            V = rand(Tensor{1, ndim})
            for i in 1:n_basefuncs
                u[i] = H ⋅ xs[i]
                u_scal[i] = V ⋅ xs[i]
            end
            u_vector = reinterpret(Float64, u, (n_basefuncs*ndim,))

            for i in 1:length(getnquadpoints(fv))
                @test getnormal(fv, i) ≈ n[face]
                @test function_gradient(fv, i, u) ≈ H
                @test function_symmetric_gradient(fv, i, u) ≈ 0.5(H + H')
                @test function_divergence(fv, i, u) ≈ trace(H)
                function_value(fv, i, u)
                if isa(fv, FaceScalarValues)
                    @test function_gradient(fv, i, u_scal) ≈ V
                    function_value(fv, i, u_scal)
                elseif isa(fv, FaceVectorValues)
                    @test function_gradient(fv, i, u_vector) ≈ function_gradient(fv, i, u) ≈ H
                    @test function_value(fv, i, u_vector) ≈ function_value(fv, i, u)
                    @test function_divergence(fv, i, u_vector) ≈ function_divergence(fv, i, u) ≈ trace(H)
                end
            end

            # Test of volume
            vol = 0.0
            for i in 1:getnquadpoints(fv)
                vol += getdetJdV(fv,i)
            end
            x_face = xs[[JuAFEM.faces(func_interpol)[face]...]]
            @test vol ≈ calculate_volume(JuAFEM.getlowerdim(func_interpol), x_face)

            # Test quadrature rule after reinit! with ref. coords
            x = JuAFEM.reference_coordinates(func_interpol)
            reinit!(fv, x, face)
            vol = 0.0
            for i in 1:getnquadpoints(fv)
                vol += getdetJdV(fv, i)
            end
            @test vol ≈ reference_volume(func_interpol, face)

            # Test spatial coordinate (after reinit with ref.coords we should get back the quad_points)
            # TODO: Renable somehow after quad rule is no longer stored in FaceValues
            #for (i, qp_x) in enumerate(getpoints(quad_rule))
            #    @test spatial_coordinate(fv, i, x) ≈ qp_x
            #end

        end

        # test copy
        fvc = copy(fv)
        for fname in fieldnames(fv)
            @test typeof(fv) == typeof(fvc)
            if !isa(getfield(fv, fname), JuAFEM.ScalarWrapper)
                @test pointer(getfield(fv, fname)) != pointer(getfield(fvc, fname))
                @test getfield(fv, fname) == getfield(fvc, fname)
            end
        end
    end
end

end # of testset


text:
include("header.jl")
using SpecialFunctions
using Knet: reluback, sigmback, tanhback, invxback, eluback, seluback

@testset "unary" begin

    function frand(f,t,d...)
        r = rand(t,d...) .* t(0.5) .+ t(0.25)
        if in(f,(acosh,asec))
            return 1 ./ r
        else
            return r
        end
    end

    bcast(f)=(x->broadcast(f,x))

    unary_fns = Any[]
    for f in Knet.unary_ops
        if isa(f,Tuple); f=f[2]; end
        push!(unary_fns, eval(Meta.parse(f)))
    end

    skip_grads = [trigamma,lgamma]
    for f in unary_fns
        f in skip_grads && continue
        #@show f
        bf = bcast(f)
        for t in (Float32, Float64)
            #@show f,t
            sx = frand(f,t)
            @test isa(f(sx),t)
            @test gradcheck(f, sx)
            for n in (1,(1,1),2,(2,1),(1,2),(2,2))
                #@show f,t,n
                ax = frand(f,t,n)
                @test gradcheck(bf, ax)
                if gpu() >= 0
                    gx = KnetArray(ax)
                    cy = bf(ax)
                    gy = bf(gx)
                    @test isapprox(cy,Array(gy))
                    @test gradcheck(bf, gx)
                end
            end
        end
    end

    # Issue #456: 2nd derivative for MLP
    for trygpu in (false, true)
        trygpu && gpu() < 0 && continue
        (x,y,dy) = randn.((10,10,10))
        if trygpu; (x,y,dy) = KnetArray.((x,y,dy)); end
        (x,y,dy) = Param.((x,y,dy))
        for f in (relu,sigm,tanh,invx,selu,elu)
            f1(x) = f.(x); @test @gcheck f1(x)
            f1i(x,i) = f1(x)[i]; @test @gcheck f1i(x,1)
            g1i(x,i) = grad(f1i)(x,i); @test @gcheck g1i(x,1)
            g1ij(x,i,j) = g1i(x,i)[j]; @test @gcheck g1ij(x,1,1)
            h1ij(x,i,j) = grad(g1ij)(x,i,j); if h1ij(x,1,1) != nothing; @test @gcheck h1ij(x,1,1); end
        end
        @test @gcheck reluback.(dy,y)
        @test @gcheck sigmback.(dy,y)
        @test @gcheck tanhback.(dy,y)
        @test @gcheck invxback.(dy,y)
        @test @gcheck seluback.(dy,y)
        @test @gcheck  eluback.(dy,y)
    end
end

nothing


text:
function getargs(args::AbstractVector)
    if isempty(args)
        args = ["--help"]
    else
        if !('-' in args[1]) prepend!(args, Ref("-p")) end # if phase is first without -p
        if length(args) >= 2 && !("-p" in args || "--phase" in args) && !('-' in args[end-1]) # if phase is last without -p
            insert!(args, length(args), "-p")
        end
    end
    s = ArgParseSettings(
        exc_handler=exception_handler,
        add_version=true,
        version="v3.2.5",
        )
    @add_arg_table! s begin
        "--phase", "-p"
            help = "The phase image that should be unwrapped"
        "--magnitude", "-m"
            help = "The magnitude image (better unwrapping if specified)"
        "--output", "-o"
            help = "The output path or filename"
            default = "unwrapped.nii"
        "--echo-times", "-t"
            help = """The echo times required for temporal unwrapping
                specified in array or range syntax (eg. "[1.5,3.0]" or
                "3.5:3.5:14"). For identical echo times, "-t epi" can be
                used with the possibility to specify the echo time as
                e.g. "-t epi 5.3" (for B0 calculation)."""
            nargs = '+'
        "--mask", "-k"
            help = """nomask | qualitymask <threshold> | robustmask | <mask_file>.
                <threshold> for qualitymask in [0;1]"""
            default = ["robustmask"]
            nargs = '+'
        "--mask-unwrapped", "-u"
            help = """Apply the mask on the unwrapped result. If mask is
                "nomask", sets it to "robustmask"."""
            action = :store_true
        "--unwrap-echoes", "-e"
            help = "Load only the specified echoes from disk"
            default = [":"]
            nargs = '+'
        "--weights", "-w"
            help = """romeo | romeo2 | romeo3 | romeo4 | romeo6 |
                bestpath | <4d-weights-file> | <flags>.
                <flags> are up to 6 bits to activate individual weights
                (eg. "1010"). The weights are (1)phasecoherence
                (2)phasegradientcoherence (3)phaselinearity (4)magcoherence
                (5)magweight (6)magweight2"""
            default = "romeo"
        "--compute-B0", "-B"
            help = """Calculate combined B0 map in [Hz]. Phase offset
                correction might be necessary if not coil-combined with
                MCPC3Ds/ASPIRE."""
            action = :store_true
        "--phase-offset-correction", "--coil-combination"
            help = """on | off | bipolar.
                Applies the MCPC3Ds method to perform phase offset
                determination and removal (for multi-echo). This option also
                allows 5D input, where the 5th dimension is channels. Coil
                combination will be performed. "bipolar" removes eddy current
                artefacts (requires >= 3 echoes)."""
            default = "off"
            nargs = '?'
            constant = "on"
        "--individual-unwrapping", "-i"
            help = """Unwraps the echoes individually (not temporal).
                This might be necessary if there is large movement
                (timeseries) or phase-offset-correction is not
                applicable."""
            action = :store_true
        "--template"
            help = """Template echo that is spatially unwrapped and used for
                temporal unwrapping"""
            arg_type = Int
            default = 2
        "--no-mmap", "-N"
            help = """Deactivate memory mapping. Memory mapping might cause
                problems on network storage"""
            action = :store_true
        "--no-rescale"
            help = """Deactivate rescaling of input images. By default the
                input phase is rescaled to the range [-π;π]. This option
                allows inputting already unwrapped phase images without
                manually wrapping them first."""
            action = :store_true
        "--threshold"
            help = """<maximum number of wraps>.
                Threshold the unwrapped phase to the maximum number of wraps
                and sets exceeding values to 0"""
            arg_type = Float64
            default = Inf
        "--verbose", "-v"
            help = "verbose output messages"
            action = :store_true
        "--correct-global", "-g"
            help = """Phase is corrected to remove global n2π phase offset. The
                median of phase values (inside mask if given) is used to
                calculate the correction term"""
            action = :store_true
        "--write-quality", "-q"
            help = """Writes out the ROMEO quality map as a 3D image with one
                value per voxel"""
            action = :store_true
        "--write-quality-all", "-Q"
            help = """Writes out an individual quality map for each of the
                ROMEO weights."""
            action = :store_true
        "--max-seeds", "-s"
            help = """EXPERIMENTAL! Sets the maximum number of seeds for
                unwrapping. Higher values allow more seperated regions."""
            arg_type = Int
            default = 1
        "--merge-regions"
            help = """EXPERIMENTAL! Spatially merges neighboring regions after
                unwrapping."""
            action = :store_true
        "--correct-regions"
            help = """EXPERIMENTAL! Performed after merging. Brings the median
                of each region closest to 0 (mod 2π)."""
            action = :store_true
        "--wrap-addition"
            help = """[0;π] EXPERIMENTAL! Usually the true phase difference of
                neighboring voxels cannot exceed π to be able to unwrap
                them. This setting increases the limit and uses 'linear
                unwrapping' of 3 voxels in a line. Neighbors can have
                (π + wrap-addition) phase difference."""
            arg_type = Float64
            default = 0.0
        "--temporal-uncertain-unwrapping"
            help = """EXPERIMENTAL! Uses spatial unwrapping on voxels that have
                high uncertainty values after temporal unwrapping."""
            action = :store_true

    end
    return parse_args(args, s)
end

function exception_handler(settings::ArgParseSettings, err, err_code::Int=1)
    if err == ArgParseError("too many arguments")
        println(stderr,
            """wrong argument formatting!"""
        )
    end
    ArgParse.default_handler(settings, err, err_code)
end

function getechoes(settings, neco)
    echoes = eval(Meta.parse(join(settings["unwrap-echoes"], " ")))
    if echoes isa Int
        echoes = [echoes]
    elseif echoes isa Matrix
        echoes = echoes[:]
    end
    echoes = (1:neco)[echoes] # expands ":"
    if (length(echoes) == 1) echoes = echoes[1] end
    return echoes
end

function getTEs(settings, neco, echoes)
    if isempty(settings["echo-times"])
        if neco == 1 || length(echoes) == 1
            return 1
        else
            error("multi-echo data is used, but no echo times are given. Please specify the echo times using the -t option.")
        end
    end
    TEs = if settings["echo-times"][1] == "epi"
        ones(neco) .* if length(settings["echo-times"]) > 1; parse(Float64, settings["echo-times"][2]) else 1 end
    else
        eval(Meta.parse(join(settings["echo-times"], " ")))
    end
    if TEs isa Matrix
        TEs = TEs[:]
    end
    if length(TEs) == neco
        TEs = TEs[echoes]
    end
    return TEs
end

function parseweights(settings)
    if isfile(settings["weights"]) && splitext(settings["weights"])[2] != ""
        return UInt8.(niread(settings["weights"]))
    else
        try
            reform = "Bool[$(join(collect(settings["weights"]), ','))]"
            flags = falses(6)
            flags_tmp = eval(Meta.parse(reform))
            flags[1:length(flags_tmp)] = flags_tmp
            return flags
        catch
            return Symbol(settings["weights"])
        end
    end
end

function saveconfiguration(writedir, settings, args)
    writedir = abspath(writedir)
    open(joinpath(writedir, "settings_romeo.txt"), "w") do io
        for (fname, val) in settings
            if !(typeof(val) <: AbstractArray)
                println(io, "$fname: " * string(val))
            end
        end
        println(io, """Arguments: $(join(args, " "))""")
    end
end


text:
macro λ(ex)
    return parse_lambda(ex)
end

"""
    parse_lambda(ex)

Parse a lambda expression to `LegibleLambda` constructor expression.
"""
function parse_lambda(ex::Expr)
    if ex.head == :(->)
        ex_cut = ex |> (ex -> postwalk(cutlnn, ex)) |> (ex -> postwalk(cutblock, ex))
        name = (repr(ex_cut)[2:end])
        return quote
            LegibleLambda($(esc(ex)), $(QuoteNode(ex)), Base.@locals)
        end
    else
        throw("Must be called on a Lambda expression")
    end
end

struct LegibleLambda{F <: Function, T} <: Function
    λ::F
    ex::Expr
    vars::Dict{Symbol, T}
end

(f::LegibleLambda)(args...; kwargs...) = (f.λ)(args...; kwargs...)

function to_string(f::LegibleLambda)
    ex = replace_variable(f.ex, f.vars)

    ex_cut = ex |> (ex -> postwalk(cutlnn, ex)) |> (ex -> postwalk(cutblock, ex))
    return replace((repr(ex_cut |> pretty_kwargs)[2:end]), "->" => " -> ")
end

function Base.show(io::IO, f::LegibleLambda)
    print(io, to_string(f))
end

function Base.show(io::IO, ::MIME"text/plain", f::LegibleLambda)
    print(io, to_string(f))
end


text:
using Random

function is_sorted(a::Vector{Float64})
    for i = 1:length(a)-1
        if (a[i] > a[i + 1])
            return false
        end
    end
    return true
end

function bogo_sort(a::Vector{Float64})
    while(!is_sorted(a))
        shuffle!(a)
    end
end

function main()
    a = [1, 3, 2, 4]
    bogo_sort(a)
    println(a)
end

main()

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
using Test

using CSV
using DataFrames
using FlightMechanicsUtils

using FlightMechanicsSimulator

@testset "engine.jl" begin

    df = DataFrame(CSV.File("data/tgear.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.tgear(case.thtl)
        @test isapprox(rv1, case.tgear, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/pdot.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.pdot(case.p3, case.p1)
        @test isapprox(rv1, case.pdot, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/rtau.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.rtau(case.dp)
        @test isapprox(rv1, case.rtau, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/thrust.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.thrust(case.pow, case.alt * FT2M, case.rmach)
        @test isapprox(rv1 / (LB2KG * F16Stevens.GD * FT2M), case.thrust, atol = 1.0e-10)
    end

end


text:
import Pkg; Pkg.activate("."); Pkg.instantiate();
include("configs/config2.jl")
include("src/main.jl")
#MODEL
@show arrtype
M    = MACNetwork(o);
Knet.gc()
#FEATS
function benchmark(M::MACNetwork,o::Dict;N=10)
    B=64; L=25
    for i=1:N
        ids  = randperm(128)[1:B]
        xB   = arrtype(ones(Float32,1,B))
        xS   = arrtype(randn(Float32,14,14,1024,B))
        xQ   = [rand(1:84) for i=1:B*L]
        answers = [rand(1:28) for i=1:B]
        batchSizes = [B for i=1:L]
        xP   = nothing
        y    = @diff M(xQ,batchSizes,xS,xB,xP;answers=answers,p=o[:p],selfattn=o[:selfattn],gating=o[:gating])
    end
end

benchmark(M,o;N=100)

max_stars_repo_name: feyzaakyurek/MAC-Network

max_issues_repo_name: feyzaakyurek/MAC-Network

max_forks_repo_name: feyzaakyurek/MAC-Network


text:
using Test
using Mimi

include("../src/main_model.jl")
include("../src/compute_scc.jl")

m = getpage("RCP4.5 & SSP2", true, true, true)
run(m)
uds = [undiscount_scc(m, 2020), undiscount_scc(m, 2040), undiscount_scc(m, 2075)]
uds_compare = [0.8765, 0.5279, 0.2162] # From Excel PAGE-ICE
@test uds ≈ uds_compare rtol=1e-4

function checksccs(sccs, sccomp)
    sccobs = [quantile(sccs, [.05, .25, .5, .75, .95]); mean(sccs)]
    sccerror = abs.((sccobs - scccomp) ./ [10, 1, 1, 1, 10, 1])
    sccsd = std(sccs) / sqrt(1000)
    @test all(sccerror .< 3*sccsd) # average off by < 3x SD
end

m = getpage("RCP4.5 & SSP2", true, true, true)
sccs = compute_scc_mcs(m, 1000, year=2020)
checksccs(sccs, [37.031, 89.742, 161.868, 299.439, 847.020, 261.739])

# m = getpage("RCP8.5 & SSP5", true, true, true)
# sccs = compute_scc_mcs(m, 1000, year=2020)

m = getpage("RCP4.5 & SSP2", true, true, false)
sccs = compute_scc_mcs(m, 1000, year=2020)
checksccs(sccs, [53.903, 129.330, 217.337, 361.739, 886.998, 306.914])

# m = getpage("RCP8.5 & SSP5", true, true, false)
# sccs = compute_scc_mcs(m, 1000, year=2020)

end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using InfiniteLinearAlgebra, BlockArrays, ArrayLayouts, Test
import InfiniteLinearAlgebra: BlockTridiagonalToeplitzLayout, ul, adaptiveqr

@testset "∞-UL" begin
    @testset "Toeplitz" begin
        Δ = SymTridiagonal(Fill(-2,∞), Fill(1,∞))
        h = 0.01
        A = I - h*Δ
        U,L = ul(A)
        N = 10
        @test U[1:N,1:N+1]*L[1:N+1,1:N] ≈ A[1:N,1:N]

        u =  adaptiveqr(A) \ [1; zeros(∞)]
        v = L \ (U \ [1;zeros(∞)])
        @test u ≈ v
    end

    @testset "Symmetric" begin
        B = randn(2,2)
        A = randn(2,2) - 10I #; A = A + A'
        C = Matrix(B')

        J = mortar(Tridiagonal(Fill(C,∞), Fill(A,∞), Fill(B,∞)))

        @test MemoryLayout(J) isa BlockTridiagonalToeplitzLayout
        U,L = ul(J, Val(false))
        N = 10
        @test U[Block.(1:N),Block.(1:N+1)] * L[Block.(1:N+1),Block.(1:N)] ≈ J[Block.(1:N),Block.(1:N)]

        @test (J \ [1; zeros(∞)])[Block(1)] ≈ inv(L[Block(1,1)])[:,1]
    end

    @testset "Periodic Jacobi" begin
        e = 0.1
        B = [e 0; 1 e]
        A = [1 1; 1 -1.0]
        C = Matrix(B')
        J = mortar(Tridiagonal(Fill(C,∞), Fill(A,∞), Fill(B,∞))) - 10I
        U,L = ul(J, Val(false))
        N = 10;
        @test U[Block.(1:N),Block.(1:N+1)] * L[Block.(1:N+1),Block.(1:N)] ≈ J[Block.(1:N),Block.(1:N)]
    end
end


text:
@testset "Testing Point Constructors" begin
    @test Point(1,2) === Point(Vec(1.0,2.0))
    @test Point((1,2)) === Point(Vec(1.0,2.0))
    @test Point([1,2]) === Point(Vec(1.0,2.0))
end

@testset  "Adding a point to vector gives a point" begin
   a  =  Point(3, -2, 5)
   v =   Vec(-2, 3, 1)
   @test a + v == Point(1, 1, 6)
end
@testset  "Addition of point to vector is orderless" begin
    a  =  Point(3, -2, 5)
    v =   Vec(-2, 3, 1)
    @test a + v == v + a
 end

@testset "Subtracting two points gives a vector" begin
     p1  = Point(3, 2, 1)
     p2 =  Point(5, 6, 7)
    @test p1 - p2 === Vec(-2.0, -4.0, -6.0)
end

@testset "Subtracting a vector from a point gives a point" begin
     p = Point(3, 2, 1)
     v = Vec(5, 6, 7)
    @test p - v === Point(-2, -4, -6)
end

@testset "Subtracting a point from a vector gives a point" begin
    p = Point(3, 2, 1)
    v = Vec(5, 6, 7)
   @test v - p === Point(-2, -4, -6)
end


text:
@testset "MIDIEvent" begin
    # test one of each length of status bytes
    statusbytes = [MIDI.NOTEOFF, MIDI.CHANNELPRESSURE];

    @testset "it should test that MIDI events are read and written successfully when status is provided, and when running status is used" begin
        for status_ in statusbytes
            # Structure: dT (variable length), status byte (optional, based on running status), data
            # data length depends on the high nybble of the status byte
            # dT is a variable length value which is handled by readvariablelength and not readMIDIevent. We don't need to test it here, which
            # is why it's set to 0 in each case.
            input = vcat([0x00, status_], zeros(UInt8, MIDI.EVENTTYPETOLENGTH[status_]))
            input_no_status = vcat([0x00], zeros(UInt8, MIDI.EVENTTYPETOLENGTH[status_]))
            type = MIDI.MIDI_EVENTS_SPEC[status_]
            output = type(0, status_, zeros(Int, MIDI.EVENTTYPETOLENGTH[status_])...)

            result = MIDI.readMIDIevent(Int(input[1]), IOBuffer(input[2:length(input)]), UInt8(0))
            @test result.dT == output.dT && result.status == status_ && MIDI.encode(result) == MIDI.encode(output)

            result = MIDI.readMIDIevent(Int(input_no_status[1]), IOBuffer(input_no_status[2:length(input_no_status)]), status_)
            @test result.dT == output.dT && MIDI.encode(result) == MIDI.encode(output)

            input, output = output, input
            output_no_status = input_no_status

            buf = IOBuffer()
            MIDI.writeevent(buf, input, true)
            @test take!(buf) == output

            buf = IOBuffer()
            MIDI.writeevent(buf, input, false)
            @test take!(buf) == output_no_status
        end
    end
end


text:
function mouse_move(s::mjSim, window::GLFW.Window, xpos::Float64, ypos::Float64)
    # no buttons down: nothing to do
    if !s.button_left && !s.button_middle && !s.button_right
        return
    end

    # compute mouse displacement, save
    dx = xpos - s.lastx
    dy = ypos - s.lasty
    s.lastx = xpos
    s.lasty = ypos

    width, height = GLFW.GetWindowSize(window)

    mod_shift = GLFW.GetKey(window, GLFW.KEY_LEFT_SHIFT) || GLFW.GetKey(window, GLFW.KEY_RIGHT_SHIFT)

    # determine action based on mouse button
    if s.button_right
        action = mod_shift ? MJCore.mjMOUSE_MOVE_H : MJCore.mjMOUSE_MOVE_V
    elseif s.button_left
        action = mod_shift ? MJCore.mjMOUSE_ROTATE_H : MJCore.mjMOUSE_ROTATE_V
    else
        action = MJCore.mjMOUSE_ZOOM
    end

    # move perturb or camera
    if s.pert[].active != 0
        mjv_movePerturb(s.m.m, s.d.d, action,
                        dx / height, dy / height,
                        s.scn, s.pert)
    else
        mjv_moveCamera(s.m.m, action,
                        dx / height, dy / height,
                        s.scn, s.cam)
    end
end

# past data for double-click detection
function mouse_button(s::mjSim, window::GLFW.Window,
                      button::GLFW.MouseButton, act::GLFW.Action, mods::Int32)
    # update button state
    s.button_left = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_LEFT)
    s.button_middle = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_MIDDLE)
    s.button_right = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_RIGHT)

    # Alt: swap left and right
    if mods == GLFW.MOD_ALT
        tmp = s.button_left
        s.button_left = s.button_right
        s.button_right = tmp

        if button == GLFW.MOUSE_BUTTON_LEFT
            button = GLFW.MOUSE_BUTTON_RIGHT
        elseif button == GLFW.MOUSE_BUTTON_RIGHT
            button = GLFW.MOUSE_BUTTON_LEFT
        end
    end

    # update mouse position
    x, y = GLFW.GetCursorPos(window)
    s.lastx = x
    s.lasty = y

    # set perturbation
    newperturb = 0
    if act == GLFW.PRESS && mods == GLFW.MOD_CONTROL && s.pert[].select > 0
        # right: translate;  left: rotate
        if s.button_right
            newperturb = MJCore.mjPERT_TRANSLATE
        else    #if s.button_left
            newperturb = MJCore.mjPERT_ROTATE
        end
        # perturbation onset: reset reference
        #if newperturb>0
            s.pert[].active==0 && mjv_initPerturb(s.m.m, s.d.d, s.scn, s.pert)
        #end
    end
    s.pert[].active = newperturb

    # detect double-click (250 msec)
    if act == GLFW.PRESS && (time() - s.lastclicktm < 0.25) && (button == s.lastbutton)
        # determine selection mode
        if button == GLFW.MOUSE_BUTTON_LEFT
            selmode = 1
        elseif mods == GLFW.MOD_CONTROL
            selmode = 3; # CTRL + Right Click
        else
            selmode = 2; # Right Click
        end
        # get current window size
        width, height = GLFW.GetWindowSize(window)

        # find geom and 3D click point, get corresponding body
        selpnt = zeros(MVector{3,Float64})
        selgeom, selskin = Ref(Cint(0)), Ref(Cint(0))
        selbody = mjv_select(s.m.m, s.d.d, s.vopt,
                            width / height, x / width,
                            (height - y) / height,
                            s.scn, selpnt, selgeom, selskin)

        # set lookat point, start tracking if requested
        if selmode == 2 || selmode == 3
            # copy selpnt if geom clicked
            if selbody >= 0
                s.cam[].lookat = SVector{3,Float64}(selpnt...)
            end

            # switch to tracking camera
            if selmode == 3 && selbody >= 0
                s.cam[].type = MJCore.mjCAMERA_TRACKING
                s.cam[].trackbodyid = selbody
                s.cam[].fixedcamid = -1
            end
        else # set body selection
            if selbody >= 0
                # compute localpos
                tmp = selpnt - s.d.xpos[:,selbody+1]
                res = reshape(s.d.xmat[:,selbody+1], 3, 3)' * tmp
                s.pert[].localpos = SVector{3}(res)

                # record selection
                s.pert[].select = selbody
                s.pert[].skinselect = selskin[]
            else
                s.pert[].select = 0
                s.pert[].skinselect = -1
            end
        end

        # stop perturbation on select
        s.pert[].active = 0
    end
    # save info
    if act == GLFW.PRESS
        s.lastbutton = button
        s.lastclicktm = time()
    end
end

function scroll(s::mjSim, window::GLFW.Window, xoffset::Float64, yoffset::Float64)
    # scroll: emulate vertical mouse motion = 5% of window height
    mjv_moveCamera(s.m.m, MJCore.mjMOUSE_ZOOM, 0.0, -0.05 * yoffset, s.scn, s.cam)
end

function drop(window::GLFW.Window,
    count::Int, paths::String)
end


text:
function Worksheet(xf::XLSXFile, sheet_element::EzXML.Node)
    @assert EzXML.nodename(sheet_element) == "sheet"

    sheetId = parse(Int, sheet_element["sheetId"])
    relationship_id = sheet_element["r:id"]
    name = sheet_element["name"]

    target = "xl/" * get_relationship_target_by_id(xf.workbook, relationship_id)
    sheet_data = xf.data[target]

    return Worksheet(xf, sheetId, relationship_id, name, sheet_data)
end

isdate1904(ws::Worksheet) = isdate1904(ws.package)

"""
Retuns the dimension of this worksheet as a CellRange.
"""
function dimension(ws::Worksheet) :: CellRange
    xroot = EzXML.root(ws.data)
    @assert EzXML.nodename(xroot) == "worksheet" "Unicorn!"

    for dimension_element in EzXML.eachelement(xroot)
        if EzXML.nodename(dimension_element) == "dimension"
            ref_str = dimension_element["ref"]
            if is_valid_cellname(ref_str)
                return CellRange("$(ref_str):$(ref_str)")
            else
                return CellRange(ref_str)
            end
        end
    end

    error("Malformed Worksheet $(ws.name): dimension not found.")
end

"""
    getdata(sheet, ref)
    getdata(filepath, sheet, ref)

Returns a escalar or a matrix with values from a spreadsheet.
`ref` can be a cell reference or a range.

Example:

```julia
julia> v = XLSX.getdata("myfile.xlsx", "mysheet", "A1:B4")
```

Indexing in a `Worksheet` will dispatch to `getdata` method.
So the following example will have the same effect as the first example.

```julia
julia> f = XLSX.read("myfile.xlsx")

julia> sheet = f["mysheet"]

julia> v = sheet["A1:B4"]
```
"""
getdata(ws::Worksheet, single::CellRef) = celldata(ws, getcell(ws, single))

function getdata(ws::Worksheet, rng::CellRange) :: Array{Any,2}
    result = Array{Any, 2}(size(rng))
    fill!(result, Missings.missing)

    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    for sheetrow in eachrow(ws)
        if top <= sheetrow.row && sheetrow.row <= bottom
            for column in left:right
                cell = getcell(sheetrow, column)
                if !isempty(cell)
                    (r, c) = relative_cell_position(cell.ref, rng)
                    result[r, c] = celldata(ws, cell)
                end
            end
        end
    end

    return result
end

function getdata(ws::Worksheet, ref::AbstractString) :: Union{Array{Any,2}, Any}
    if is_valid_cellname(ref)
        return getdata(ws, CellRef(ref))
    elseif is_valid_cellrange(ref)
        return getdata(ws, CellRange(ref))
    else
        error("$ref is not a valid cell or range reference.")
    end
end

getdata(ws::Worksheet) = getdata(ws, dimension(ws))

Base.getindex(ws::Worksheet, r) = getdata(ws, r)
Base.getindex(ws::Worksheet, ::Colon) = getdata(ws)

Base.show(io::IO, ws::Worksheet) = print(io, "XLSX.Worksheet: \"$(ws.name)\". Dimension: $(dimension(ws)).")

"""
    getcell(sheet, ref)
    getcell(filepath, sheet, ref)

Returns an `AbstractCell` that represents a cell in the spreadsheet.

Example:

```julia
julia> sheet = XLSX.getsheet("myfile.xlsx", "mysheet")

julia> cell = XLSX.getcell(sheet, "A1")
```
"""
function getcell(ws::Worksheet, single::CellRef) :: AbstractCell

    for sheetrow in eachrow(ws)
        if row_number(sheetrow) == row_number(single)
            return getcell(sheetrow, column_number(single))
        end
    end

    return EmptyCell(single)
end

function getcell(ws::Worksheet, ref::AbstractString)
    if is_valid_cellname(ref)
        return getcell(ws, CellRef(ref))
    else
        error("$ref is not a valid cell reference.")
    end
end

"""
    getcellrange(sheet, rng)
    getcellrange(filepath, sheet, rng)

Returns a matrix with cells as `Array{AbstractCell, 2}`.
`rng` must be a valid cell range, as in `"A1:B2"`.

Example:

```julia
julia> XLSX.getcellrange("myfile.xlsx", "mysheet", "A1:B4")
```
"""
function getcellrange(ws::Worksheet, rng::CellRange) :: Array{AbstractCell,2}
    result = Array{AbstractCell, 2}(size(rng))
    for cellref in rng
        (r, c) = relative_cell_position(cellref, rng)
        result[r, c] = EmptyCell(cellref)
    end

    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    for sheetrow in eachrow(ws)
        if top <= sheetrow.row && sheetrow.row <= bottom
            for column in left:right
                cell = getcell(sheetrow, column)
                if !isempty(cell)
                    (r, c) = relative_cell_position(cell.ref, rng)
                    result[r, c] = cell
                end
            end
        end
    end

    return result
end

function getcellrange(ws::Worksheet, rng::AbstractString)
    if is_valid_cellrange(rng)
        return getcellrange(ws, CellRange(rng))
    else
        error("$rng is not a valid cell range.")
    end
end

"""
    gettable(sheet, [columns]; [first_row], [column_labels], [header], [infer_eltypes]) -> data, column_labels
    gettable(filepath, sheet, [columns]; [first_row], [column_labels], [header], [infer_eltypes]) -> data, column_labels

Returns tabular data from a spreadsheet as a tuple `(data, column_labels)`.
`data` is a vector of columns. `column_labels` is a vector of symbols.
Use this function to create a `DataFrame` from package `DataFrames.jl`.

Use `columns` argument to specify which columns to get.
For example, `columns="B:D"` will select columns `B`, `C` and `D`.
If `columns` is not given, the algorithm will find the first sequence
of consecutive non-empty cells.

Use `first_row` to indicate the first row from the table.
`first_row=5` will look for a table starting at sheet row `5`.
If `first_row` is not given, the algorithm will look for the first
non-empty row in the spreadsheet.

`header` is a `Bool` indicating if the first row is a header.
If `header=true` and `column_labels` is not specified, the column labels
for the table will be read from the first row of the table.
If `header=false` and `column_labels` is not specified, the algorithm
will generate column labels. The default value is `header=true`.

Use `column_labels` as a vector of symbols to specify names for the header of the table.

Use `infer_eltypes=true` to get `data` as a `Vector{Any}` of typed vectors.
The default value is `infer_eltypes=false`.

Example:

```julia
julia> using DataFrames, XLSX

julia> df = DataFrame(XLSX.gettable("myfile.xlsx", "mysheet")...)
```
"""
function gettable(sheet::Worksheet, cols::Union{ColumnRange, AbstractString}; first_row::Int=_find_first_row_with_data(sheet, convert(ColumnRange, cols).start), column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true, infer_eltypes::Bool=false)
    itr = TableRowIterator(sheet, cols; first_row=first_row, column_labels=column_labels, header=header)
    return gettable(itr; infer_eltypes=infer_eltypes)
end

function gettable(sheet::Worksheet; first_row::Int = 1, column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true, infer_eltypes::Bool=false)
    itr = TableRowIterator(sheet; first_row=first_row, column_labels=column_labels, header=header)
    return gettable(itr; infer_eltypes=infer_eltypes)
end

#
# Helper functions
#
getcell(filepath::AbstractString, sheet::Union{AbstractString, Int}, ref) = getcell( read(filepath)[sheet], ref )
getcellrange(filepath::AbstractString, sheet::Union{AbstractString, Int}, rng) = getcellrange( read(filepath)[sheet], rng )
getdata(filepath::AbstractString, sheet::Union{AbstractString, Int}, ref) = getdata( read(filepath)[sheet], ref )
gettable(filepath::AbstractString, sheet::Union{AbstractString, Int}, cols::Union{ColumnRange, AbstractString}; first_row::Int=_find_first_row_with_data(sheet, convert(ColumnRange, cols).start), column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true, infer_eltypes::Bool=false) = gettable( read(filepath)[sheet], cols; first_row=first_row, column_labels=column_labels, header=header, infer_eltypes=infer_eltypes )
gettable(filepath::AbstractString, sheet::Union{AbstractString, Int}; first_row::Int = 1, column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true, infer_eltypes::Bool=false) = gettable( read(filepath)[sheet]; first_row=first_row, column_labels=column_labels, header=header, infer_eltypes=infer_eltypes )


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct PhaseShiftingTransformer <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    r::Float64  # System per-unit value
    x::Float64  # System per-unit value
    primaryshunt::Float64
    tap::Float64
    α::Float64
    rate::Union{Nothing, Float64}
    internal::InfrastructureSystemsInternal
end

function PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, )
    PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, InfrastructureSystemsInternal())
end

function PhaseShiftingTransformer(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, )
    PhaseShiftingTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, α, rate, )
end

# Constructor for demo purposes; non-functional.

function PhaseShiftingTransformer(::Nothing)
    PhaseShiftingTransformer(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        tap=1.0,
        α=0.0,
        rate=0.0,
    )
end

"""Get PhaseShiftingTransformer name."""
get_name(value::PhaseShiftingTransformer) = value.name
"""Get PhaseShiftingTransformer available."""
get_available(value::PhaseShiftingTransformer) = value.available
"""Get PhaseShiftingTransformer activepower_flow."""
get_activepower_flow(value::PhaseShiftingTransformer) = value.activepower_flow
"""Get PhaseShiftingTransformer reactivepower_flow."""
get_reactivepower_flow(value::PhaseShiftingTransformer) = value.reactivepower_flow
"""Get PhaseShiftingTransformer arc."""
get_arc(value::PhaseShiftingTransformer) = value.arc
"""Get PhaseShiftingTransformer r."""
get_r(value::PhaseShiftingTransformer) = value.r
"""Get PhaseShiftingTransformer x."""
get_x(value::PhaseShiftingTransformer) = value.x
"""Get PhaseShiftingTransformer primaryshunt."""
get_primaryshunt(value::PhaseShiftingTransformer) = value.primaryshunt
"""Get PhaseShiftingTransformer tap."""
get_tap(value::PhaseShiftingTransformer) = value.tap
"""Get PhaseShiftingTransformer α."""
get_α(value::PhaseShiftingTransformer) = value.α
"""Get PhaseShiftingTransformer rate."""
get_rate(value::PhaseShiftingTransformer) = value.rate
"""Get PhaseShiftingTransformer internal."""
get_internal(value::PhaseShiftingTransformer) = value.internal


text:
@testset "Array" begin
    @testset "range" begin
        x = range(AFArray{Float32}, 1, 5)
        @test Array(x) == Array{Float32}(collect(1.0:5.0))
        @test typeof(x) == AFArray{Float32, 1}
        x = range(AFArray{Float32}, 3, 5)
        @test Array(x) == Array{Float32}(collect(3.0:7.0))
        @test typeof(x) == AFArray{Float32, 1}
        x = range(AFArray{Float32}, 2, 3, 10)
        @test Array(x) == collect(2.0:3:29.0)
        @test typeof(x) == AFArray{Float32, 1}

        @test sum(@inferred iota((2,3))) == 21
    end
end


text:
using MapTiles
using Test

@testset "Raster Base Maps" begin
    basemap = MapTiles.fetchmap(
        -1, # minlon
        48, # minlat
        3, # maxlon
        52, # maxlat
        4, # zoom
    )
    @test size(basemap.img) == (256,512)
    @test basemap.tilesize == (256,256)
    @test basemap.xmin == 7
    @test basemap.xmax == 8
    @test basemap.ymin == 5
    @test basemap.ymax == 5
    @test basemap.zoom == 4
end