text:
using EffectiveWaves, Plots

pyplot()

Maxtime=100.
T=Float64

spatial_dim = 3
medium = Acoustic(spatial_dim; ρ=1.0, c=1.0)

v = 1.0
s1 = Specie(
    Acoustic(spatial_dim; ρ=10.0, c=10.0), Sphere(1.0);
    volume_fraction = v * 3/4
);
s2 = Specie(
    Acoustic(spatial_dim; ρ=0.1, c=0.1), Sphere(1.0);
    volume_fraction = v * 1/4
);
species = [s1,s2]

k_low = ω ./ effective_medium(medium, species).c

## Swap a specie for the medium

    medium = Acoustic(spatial_dim; ρ=10.0, c=10.0)

    v = 1.0
    s2 = Specie(
        Acoustic(spatial_dim; ρ=0.1, c=0.1), Sphere(1.0);
        volume_fraction = v * 1/4
    );
    species = [s2]

    k_low = ω ./ effective_medium(medium, species).c

ω = 0.04187
basis_order = 1
tol = 1e-7
k = ω/medium.c

opts = Dict(
    :symmetry => PlanarAzimuthalSymmetry(),
    :tol => tol,
    # :num_wavenumbers => 4,
    # :mesh_size => 2.0, :mesh_points => 20,
    :basis_order => basis_order
)

detMM = dispersion_complex(ω, medium, species; opts...)

kps =  wavenumbers(ω, medium,  species; opts...)

k_low = ω ./ effective_medium(medium, species).c

# Explore around all the wavenumbers

inds = 1:1
x_max = maximum(abs.(real.(kps[inds])))
y_max = maximum(abs.(imag.(kps[inds])))

x = LinRange(-x_max,x_max,150)
y = LinRange(0.,y_max,100)

X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))

z_max = 2.0
Z0 = map( (x,y) -> abs(detMM(x + y*im)),X,Y);
Z = map( z -> (abs(z)> z_max) ? NaN : z, Z0);
contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M - high ω - low φ")
scatter!(kps[inds])

# Explore around one wavenumber

    k1 = abs(real(kps[1])) - imag(kps[1]) * im

    scale = 0.01
    x = LinRange(real(k1) - abs(real(k1)) * scale, real(k1) + abs(real(k1)) * scale,100)
    y = LinRange(imag(k1) - abs(imag(k1)) * scale, imag(k1) + abs(imag(k1)) * scale,100)

    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))

    Z0 = map( (x,y) -> abs(detMM(x + y*im)),X,Y);

    z_max = 2 * scale
    Z = map( z -> (abs(z)> z_max) ? NaN : z, Z0);
    contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M - high ω - low φ")
    scatter!([k1])

# Old code below

x = k0 .* LinRange(0.18,0.23,50)
y = k0 .* LinRange(0.01,0.17,50)

x = kx; y = ky;

X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map( (x,y) -> (z = detMM2([x,y]); (abs(z)> 0.5) ? NaN: z),X,Y)

contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M")

x = LinRange(3.,6.,40)
y = LinRange(-1.5,14.,40)

X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map( (x,y) -> (z = detMM2([x,y]); (abs(z)> 1.) ? NaN: z),X,Y)
# Z = map( (x,y) -> detMM2([x,y]),X,Y)

# contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M")
heatmap(x,y,Z, xlab = "Re k*", ylab = "Im k*", title="det M - high ω and almost dirichlet")

max_stars_repo_name: arturgower/EffectiveWaves


text:
module Test20211215

using PropertyDestructuring

@destructure const (; num, den) = 0x03//0x0f

end

ext: jl


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function runtests(name, isolate=true)
    if isolate
        mod_name = Symbol("TestMain_", replace(name, '/', '_'))
        m = eval(Main, :(module $mod_name end))
    else
        m = Main
    end
    eval(m, :(using Base.Test))
    @printf("     \033[1m*\033[0m \033[31m%-21s\033[0m", name)
    ex = quote
        @timed @testset $"$name" begin
            include($"$name.jl")
        end
    end
    res_and_time_data = eval(m, ex)
    rss = Sys.maxrss()
    @printf(" maxrss %7.2f MB\n", rss / 2^20)
    #res_and_time_data[1] is the testset
    passes,fails,error,broken,c_passes,c_fails,c_errors,c_broken = Base.Test.get_test_counts(res_and_time_data[1])
    if res_and_time_data[1].anynonpass == false
        res_and_time_data = (
                             (passes+c_passes,broken+c_broken),
                             res_and_time_data[2],
                             res_and_time_data[3],
                             res_and_time_data[4],
                             res_and_time_data[5])
    end
    vcat(collect(res_and_time_data), rss)
end

# looking in . messes things up badly
filter!(x->x!=".", LOAD_PATH)
nothing

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"""
    save_monitoring_results(filename::String, r::MonitoringResults; overwrite=false)
Store the given `MonitoringResults` (output of [`monitoring_stop`](@ref)) to disk as an HDF5 file with name `filename`.
"""
function save_monitoring_results(filename::String, r::MonitoringResults; overwrite=false)
    if isfile(filename) && !overwrite
        error(
            "Output file $filename already exists. If you want to overwrite, set `overwrite=true`.",
        )
    end

    h5open(filename, "w") do f
        f["times"] = r.times
        f["devices_strings"] = [str for (str, uuid) in r.devices]
        f["devices_uuids"] = [string(uuid) for (str, uuid) in r.devices]
        for key in keys(r.results)
            # store Vector{Vector{T}} as Matrix{T} where
            # different rows correspond to different devices
            data = reduce(hcat, r.results[key])
            f[string(key)] = data
        end
    end
    return nothing
end

"""
Given an HDF5 file created with [`save_monitoring_results`](@ref),
restore the saved monitoring results (i.e. output of [`monitoring_stop`](@ref)).
"""
function load_monitoring_results(filename::String)
    if !isfile(filename)
        error("Input file $filename doesn't exists.")
    end

    h5open(filename, "r") do f
        times = read(f["times"])
        dev_strs = read(f["devices_strings"])
        dev_uuids = read(f["devices_uuids"])
        devices = Tuple{String,UUID}[
            (dev_strs[i], UUID(dev_uuids[i])) for i in eachindex(dev_uuids)
        ]
        d = Dict{Symbol,Vector{Vector{Float64}}}()
        for key in keys(f)
            if key == "times" || key == "devices"
                continue
            end
            data = read(f[key])
            if data isa Matrix
                d[Symbol(key)] = [collect(col) for col in eachcol(data)]
            end
        end
        return MonitoringResults(times, devices, d)
    end
end


text:
abstract type MessagePassing <: GraphNet end

"""
    message(mp, x_i, x_j, e_ij)

Message function for message-passing scheme. This function can be overrided to dispatch to custom layers.
First argument should be message-passing layer, the rest of arguments can be `x_i`, `x_j` and `e_ij`.

# Arguments
- `mp`: message-passing layer.
- `x_i`: the feature of node `x_i`.
- `x_j`: the feature of neighbors of node `x_i`.
- `e_ij`: the feature of edge (`x_i`, `x_j`).
"""
@inline message(mp::T, x_i, x_j, e_ij) where {T<:MessagePassing} = x_j

"""
    update(mp, m, x)

Update function for message-passing scheme. This function can be overrided to dispatch to custom layers.
First argument should be message-passing layer, the rest of arguments can be `X` and `M`.

# Arguments
- `mp`: message-passing layer.
- `m`: the message aggregated from message function.
- `x`: the single node feature.
"""
@inline update(mp::T, m, x) where {T<:MessagePassing} = m

@inline function update_batch_edge(mp::T, adj, E::AbstractMatrix, X::AbstractMatrix) where {T<:MessagePassing}
    edge_idx = edge_index_table(adj)
    E_ = Vector[]
    for (i, js) = enumerate(adj)
        for j = js
            k = edge_idx[(i,j)]
            m = message(mp, get_feature(X, i), get_feature(X, j), get_feature(E, k))
            push!(E_, m)
        end
    end
    hcat(E_...)
end

@inline function update_batch_vertex(mp::T, M::AbstractMatrix, X::AbstractMatrix) where {T<:MessagePassing}
    X_ = Vector[]
    for i = 1:size(X,2)
        x = update(mp, get_feature(M, i), get_feature(X, i))
        push!(X_, x)
    end
    hcat(X_...)
end

@inline function aggregate_neighbors(mp::T, aggr::Symbol, M::AbstractMatrix, accu_edge, num_V, num_E) where {T<:MessagePassing}
    @assert !iszero(accu_edge) "accumulated edge must not be zero."
    cluster = generate_cluster(M, accu_edge, num_V, num_E)
    pool(aggr, cluster, M)
end

function propagate(mp::T, fg::FeaturedGraph, aggr::Symbol=:add) where {T<:MessagePassing}
    adj = adjacency_list(fg)
    num_V = nv(fg)
    accu_edge = accumulated_edges(adj)
    num_E = accu_edge[end]
    E = edge_feature(fg)
    X = node_feature(fg)

    E = update_batch_edge(mp, adj, E, X)

    M = aggregate_neighbors(mp, aggr, E, accu_edge, num_V, num_E)

    X = update_batch_vertex(mp, M, X)

    FeaturedGraph(graph(fg), X, E, zeros(0))
end


text:
using PyCall
import JSON

"""
Run an SCF in ABINIT starting from the input file `infile`
represented as a `abipy.abilab.AbinitInput` python object.
Write the result to the `output` directory in ETSF Nanoquanta format
and return the result as an `EtsfFolder` object.
"""
function run_abinit_scf(infile::PyObject, outdir)
    abilab = pyimport("abipy.abilab")
    flowtk = pyimport("abipy.flowtk")

    # Check ABINIT can be found
    abierrors = abilab.abicheck()
    abierrors != "" && error("abilab.abicheck() returned errors:\n$abierrors")

    # Create rundir
    rundir = joinpath(outdir, "abinit")
    if isdir(rundir)
        rm(rundir, recursive=true)
    end
    mkdir(rundir)

    # Adjust common infile settings:
    infile.set_vars(
        paral_kgb=0,    # Parallelization over k-Points and Bands
        iomode=3,       # Use NetCDF output
        istwfk="*1",    # Needed for extracting the wave function later
    )

    # Dump pspmap
    open(joinpath(outdir, "pspmap.json"), "w") do fp
        JSON.print(fp, infile.pspmap)
    end
    infile.extra = nothing  # Remove extra field (causes problems below)

    # Create flow and run it
    flow = flowtk.Flow(rundir, flowtk.TaskManager.from_user_config())
    work = flowtk.Work()
    scf_task = work.register_scf_task(infile)
    flow.register_work(work)
    flow.allocate()
    flow.make_scheduler().start()

    if !flow.all_ok
        @warn "Flow not all_ok ... check input and output files"
    end

    for file in collect(scf_task.out_files())
        if endswith(file, ".nc")
            cp(file, joinpath(outdir, basename(file)), force=true)
        end
    end

    EtsfFolder(outdir)
end

"""
Run an SCF in ABINIT starting from a DFTK `Model` and some extra parameters.
Write the result to the `output` directory in ETSF Nanoquanta format
and return the `EtsfFolder` object.
"""
function run_abinit_scf(model::Model, outdir;
                        kgrid, abinitpseudos, Ecut, n_bands, tol=1e-6, kwargs...)
    abilab = pyimport("abipy.abilab")

    # Would be nice to generate the pseudofiles in ABINIT format on the fly,
    # but not for now. See the function psp10in  in 64_psp/m_psp_hgh.F90
    # for the parser, which is used in ABINIT for GTH and HGH files.
    structure = abilab.Structure.as_structure(pymatgen_structure(model))
    pseudos = pyimport("abipy.data").pseudos(abinitpseudos...)
    infile = abilab.AbinitInput(structure=structure, pseudos=pseudos)

    # Copy what can just be copied
    infile.set_kmesh(ngkpt=Vector{Int}(kgrid), shiftk=[0, 0, 0])
    infile.set_vars(
        ecut=Ecut,        # Hartree
        nband=n_bands,    # Number of bands
        tolvrs=tol,       # General tolerance settings
    )

    # Spin-polarization
    if model.spin_polarization == :spinless
        error("spin_polarization == spinless is not supported by abinit")
    elseif model.spin_polarization == :none
        infile.set_vars(nsppol=1, nspinor=1, nspden=1)
    elseif model.spin_polarization == :collinear
        infile.set_vars(nsppol=2, nspinor=1, nspden=2)
    elseif model.spin_polarization == :full
        infile.set_vars(nsppol=1, nspinor=2, nspden=4)
    end

    # Check all required terms are there
    if any(isnothing, indexin([Kinetic(), AtomicLocal(), AtomicNonlocal(),
                               Ewald(), PspCorrection(), Hartree()], model.term_types))
        error("run_abinit_scf only supports reduced Hartree-Fock or DFT models.")
    end

    # Parse XC term
    idcs_Xc = findall(t -> t isa Xc, model.term_types)
    @assert length(idcs_Xc) <= 1
    if isempty(idcs_Xc) || isempty(model.term_types[idcs_Xc[1]].functionals)
        infile.set_vars(ixc=0)  # Reduced HF
    else
        xcterm = model.term_types[idcs_Xc[1]]
        functionals = sort([func.identifier for func in xcterm.functionals])
        if functionals == [:lda_xc_teter93]
            infile.set_vars(ixc=1)
        elseif functionals == [:lda_c_vwn, :lda_x]
            infile.set_vars(ixc="-001007")
        elseif functionals == [:gga_c_pbe, :gga_x_pbe]
            infile.set_vars(ixc="-101130")  # Version implemented in libxc
            # infile.set_vars(ixc=11)       # Version implemented in ABINIT
            # NOTE: The results only agree to 7 digits or so on a few test calculations I did
        else
            error("Unknown functional combination: $functionals")
        end
    end

    # Spin-polarization
    if model.spin_polarization == :spinless
        error("spin_polarization == spinless is not supported by abinit")
    elseif model.spin_polarization == :none
        infile.set_vars(nsppol=1, nspinor=1, nspden=1)
    elseif model.spin_polarization == :collinear
        infile.set_vars(nsppol=2, nspinor=1, nspden=2)
    elseif model.spin_polarization == :full
        infile.set_vars(nsppol=1, nspinor=2, nspden=4)
    end

    # Parse occopt
    if isa(model.smearing, Smearing.None)
        infile.set_vars(occopt=1)
    elseif isa(model.smearing, Smearing.FermiDirac)
        infile.set_vars(occopt=3, tsmear=model.temperature)
    elseif isa(model.smearing, Smearing.MethfesselPaxton2)
        infile.set_vars(occopt=6, tsmear=model.temperature)
    elseif isa(model.smearing, Smearing.Gaussian)
        infile.set_vars(occopt=7, tsmear=model.temperature)
    else
        error("Smearing $(model.smearing) not implemented.")
    end

    # Add pspmap
    pspmap = Dict{Int, String}()
    for (element, positions) in model.atoms
        (element isa ElementPsp) && continue
        element.psp.identifier == "" && continue
        pspmap[element.Z] = element.psp.identifier
    end
    infile.pspmap = pspmap

    !isempty(kwargs) && infile.set_vars(;kwargs...)
    run_abinit_scf(infile, outdir)
end


text:
struct Trans2Covariate{S<:AbstractFloat, T<:AbstractFloat, U<:AbstractFloat} <: AbstractCovariate{U}
    name::String
    basecovariate1::AbstractCovariate{S}
    basecovariate2::AbstractCovariate{T}
    transform::Function
end

Base.length(var::Trans2Covariate) = length(var.basecovariate1)

function Trans2Covariate(name::String, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T},
                        transform::Function) where {S<:AbstractFloat} where {T<:AbstractFloat}
    U = promote_type(S, T)
    Trans2Covariate{S, T, U}(name, basecovariate1, basecovariate2, transform)
end

function Trans2Covariate(::Type{U}, name::String, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T},
                         transform::Function) where {U<:AbstractFloat} where {S<:AbstractFloat} where {T<:AbstractFloat}
    Trans2Covariate{S, T, U}(name, basecovariate1, basecovariate2, transform)
end

function slice(covariate::Trans2Covariate{S, T, U}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractFloat} where {T<:AbstractFloat} where {U<:AbstractFloat}
    basecov1 = covariate.basecovariate1
    basecov2 = covariate.basecovariate2
    f = covariate.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip(slice(basecov1, fromobs, toobs, slicelength), slice(basecov2, fromobs, toobs, slicelength))
    mapslice2(f, slices, slicelength, U)
end

function Base.map(covariate::Trans2Covariate, dataframe::AbstractDataFrame)
    Trans2Covariate(covariate.name, map(covariate.basecovariate1, dataframe), map(covariate.basecovariate1, dataframe),
                    covariate.transform)
end

function Base.broadcast(f, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T}) where {S<:AbstractFloat}  where {T<:AbstractFloat}
    if typeof(f(zero(S), zero(T))) <: AbstractFloat
        Trans2Covariate("$(getname(basecovariate1))$(f)$(getname(basecovariate2))", basecovariate1, basecovariate2, f)
    else
        Trans2CovBoolVariate{S, T}("$(getname(basecovariate1))$(f)$(getname(basecovariate2))", basecovariate1, basecovariate2, f)
    end
end

# function covariate(f, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T}) where {S<:AbstractFloat}  where {T<:AbstractFloat}
#     Trans2Covariate("$(f)($(getname(basecovariate1))$(getname(basecovariate2)))", basecovariate1, basecovariate2, f)
# end

# function covariate(f, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T}) where {S<:AbstractFloat}  where {T<:AbstractFloat}
#     Trans2Covariate("$(f)($(getname(basecovariate1))$(getname(basecovariate2)))", basecovariate1, basecovariate2, f)
# end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct InnerError <: SwaggerModel
    exceptiontype::Any # spec type: Union{ Nothing, String } # spec name: exceptiontype
    errordetail::Any # spec type: Union{ Nothing, String } # spec name: errordetail

    function InnerError(;exceptiontype=nothing, errordetail=nothing)
        o = new()
        validate_property(InnerError, Symbol("exceptiontype"), exceptiontype)
        setfield!(o, Symbol("exceptiontype"), exceptiontype)
        validate_property(InnerError, Symbol("errordetail"), errordetail)
        setfield!(o, Symbol("errordetail"), errordetail)
        o
    end
end # type InnerError

const _property_map_InnerError = Dict{Symbol,Symbol}(Symbol("exceptiontype")=>Symbol("exceptiontype"), Symbol("errordetail")=>Symbol("errordetail"))
const _property_types_InnerError = Dict{Symbol,String}(Symbol("exceptiontype")=>"String", Symbol("errordetail")=>"String")
Base.propertynames(::Type{ InnerError }) = collect(keys(_property_map_InnerError))
Swagger.property_type(::Type{ InnerError }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_InnerError[name]))}
Swagger.field_name(::Type{ InnerError }, property_name::Symbol) =  _property_map_InnerError[property_name]

function check_required(o::InnerError)
    true
end

function validate_property(::Type{ InnerError }, name::Symbol, val)
end


text:
include("copyFolder.jl")
using Base.Filesystem

ElementsSEB="/media/wang/ElementsSEB"
LinuxData_1T="/media/wang/LinuxData_1T"
Ex4= "/media/wang/ext4"
TwoTData = "/media/wang/2TData"

BackuRoot="backup"

srcRoot=Base.Filesystem.joinpath(  ElementsSEB ,BackuRoot )
tgrRoot=Base.Filesystem.joinpath(  TwoTData ,BackuRoot )

proCopyFolder( "$srcRoot","$tgrRoot")


text:
using Test
using BlockSystems
using ModelingToolkit
using ModelingToolkit: get_iv, get_eqs, get_states
using LinearAlgebra
using ModelingToolkit.Symbolics

@info "Testes of function_generation.jl"

@testset "function_generation.jl" begin
    @testset "transform algebraic equations" begin
        using BlockSystems: transform_algebraic_equations
        @parameters t a b
        @variables x(t) y(t)
        D = Differential(t)
        eqs = [D(x) ~ y,
               y ~ x,
               2*b + 8 ~ a + a,
               1 ~ x + y]
        eqs = transform_algebraic_equations(eqs)
        @test eqs == [D(x) ~ y,
                      0 ~ x - y,
                      0 ~ (a + a) - (2*b + 8),
                      0 ~ (x + y) - 1]
    end

    @testset "reorder of equations" begin
        using BlockSystems: reorder_by_states
        @parameters t
        @variables x(t) y(t) a(t) b(t)
        D = Differential(t)
        eqs = [D(x) ~ x, D(y) ~ y, 0 ~ x+y+a, 0 ~ x+y+b]
        @test reorder_by_states(eqs, [x,y,a,b]) == eqs[[1,2,3,4]]
        @test reorder_by_states(eqs, [y,x,a,b]) == eqs[[2,1,3,4]]
        @test reorder_by_states(eqs, [a,x,y,b]) == eqs[[3,1,2,4]]
        @test reorder_by_states(eqs, [x,a,y,b]) == eqs[[1,3,2,4]]
        @test reorder_by_states(eqs, [y,a,x,b]) == eqs[[2,3,1,4]]
        @test reorder_by_states(eqs, [a,y,x,b]) == eqs[[3,2,1,4]]
        @test reorder_by_states(eqs, [a,y,b,x]) == eqs[[3,2,4,1]]
        @test reorder_by_states(eqs, [y,a,b,x]) == eqs[[2,3,4,1]]
        @test reorder_by_states(eqs, [b,a,y,x]) == eqs[[3,4,2,1]]
        @test reorder_by_states(eqs, [a,b,y,x]) == eqs[[3,4,2,1]]
        @test reorder_by_states(eqs, [y,b,a,x]) == eqs[[2,3,4,1]]
        @test reorder_by_states(eqs, [b,y,a,x]) == eqs[[3,2,4,1]]
        @test reorder_by_states(eqs, [b,x,a,y]) == eqs[[3,1,4,2]]
        @test reorder_by_states(eqs, [x,b,a,y]) == eqs[[1,3,4,2]]
        @test reorder_by_states(eqs, [a,b,x,y]) == eqs[[3,4,1,2]]
        @test reorder_by_states(eqs, [b,a,x,y]) == eqs[[3,4,1,2]]
        @test reorder_by_states(eqs, [x,a,b,y]) == eqs[[1,3,4,2]]
        @test reorder_by_states(eqs, [a,x,b,y]) == eqs[[3,1,4,2]]
        @test reorder_by_states(eqs, [y,x,b,a]) == eqs[[2,1,3,4]]
        @test reorder_by_states(eqs, [x,y,b,a]) == eqs[[1,2,3,4]]
        @test reorder_by_states(eqs, [b,y,x,a]) == eqs[[3,2,1,4]]
        @test reorder_by_states(eqs, [y,b,x,a]) == eqs[[2,3,1,4]]
        @test reorder_by_states(eqs, [x,b,y,a]) == eqs[[1,3,2,4]]
        @test reorder_by_states(eqs, [b,x,y,a]) == eqs[[3,1,2,4]]

        eqs = [0 ~ x+y]
        @test reorder_by_states(eqs, [x]) == eqs
    end

    @testset "generate_massmatrix" begin
        using BlockSystems: generate_massmatrix
        @parameters t
        @variables x(t) y(t) a(t) b(t)
        D = Differential(t)
        eqs = [D(x) ~ x, D(y) ~ y, 0 ~ x+y+a, 0 ~ x+y+b]
        @test generate_massmatrix(eqs) == Diagonal([1,1,0,0])
        @test generate_massmatrix(eqs[[1,3,2,4]]) == Diagonal([1,0,1,0])
        eqs = [D(x) ~ x, D(y) ~ y]
        @test generate_massmatrix(eqs) == I
    end

    @testset "all_static" begin
        using BlockSystems: all_static
        @parameters t a b i1(t) i2(t)
        @variables o1(t) o2(t)
        D = Differential(t)
        eqs = [o1 ~ a*i1 + i2,
               o2 ~ b*i1 + i2]
        @test all_static(eqs)
        eqs = [D(o1) ~ a*i1 + i2,
               o2 ~ b*i1 + i2]
        @test !all_static(eqs)
        eqs = [o1 ~ a*i1 + i2 + o2,
               o2 ~ b*i1 + i2]
        @test !all_static(eqs)
    end

    @testset "test strip_iv" begin
        @parameters t a(t) b a
        @variables x(t) y(t) z
        using BlockSystems: strip_iv
        @test isequal(strip_iv(t, t), t)
        @test isequal(strip_iv(b, t), b)
        @test isequal(strip_iv(z, t), z)
        @test isequal(strip_iv(x, t), Sym{Real}(:x))
        @test isequal(strip_iv(y, t), Sym{Real}(:y))
        @test isequal(strip_iv(a, t), Sym{Real}(:a))
    end

    @testset "ode function generation" begin
        @parameters t a i(t)
        @variables x(t) y(t)
        D = Differential(t)
        iob = IOBlock([D(x) ~ a*i, D(y) ~ x], [i], [x, y])

        iof = generate_io_function(iob, warn=false);
        @test iof.massm == I
        @test isequal(iof.inputs, [Sym{Real}(:i)])
        @test isequal(iof.states, Sym{Real}.([:x, :y]))
        @test isequal(iof.params, [Sym{Real}(:a)])
        @test isequal(iof.params, [Sym{Real}(:a)])
        @test iof.g_ip === nothing
        @test iof.g_oop === nothing
        @test isempty(iof.rem_states)

        @test iof.f_oop([-1,1], [4], [-1], 0) == [-4, -1]
        a = Vector{Float64}(undef, 2)
        iof.f_ip(a, [-1,1], [4], [-1], 0)
        @test a == [-4, -1]

        iof = generate_io_function(iob, f_states=[y], warn=false)
        @test isequal(iof.states, Sym{Real}.([:y, :x]))
        iof = generate_io_function(iob, f_states=[iob.y], warn=false)
        @test isequal(iof.states, Sym{Real}.([:y, :x]))
    end

    @testset "static function generation" begin
        @parameters t a i1(t) i2(t)
        @variables x(t) y(t)
        D = Differential(t)
        iob = IOBlock([x ~ a*i1, y ~ i2], [i1, i2], [x, y])

        iof = generate_io_function(iob, warn=false);

        @test isequal(iof.inputs, Sym{Real}.([:i1, :i2]))
        @test isequal(iof.states, Sym{Real}.([:x, :y]))
        @test isequal(iof.params, [Sym{Real}(:a)])

        @test iof.f_oop([4, 1], [-1], 0) == [-4, 1]
        a = Vector{Float64}(undef, 2)
        iof.f_ip(a, [-1,1], [-1], 0)
        @test a == [1, 1]
    end

    @testset "parameter ordering" begin
        using ModelingToolkit.Symbolics: tosymbol

        @parameters t a b i1(t) i2(t)
        @variables o1(t) o2(t)
        D = Differential(t)

        iob = IOBlock([D(o1) ~ a*i1 + b, D(o2) ~ o1+i2], [i1, i2], [o1, o2])

        allequal(v1, v2) = all([isequal(e1, e2) for (e1, e2) ∈ zip(v1, v2)])

        gen = generate_io_function(iob, f_states=[o1, o2], warn=false)
        @test allequal(gen.states, strip_iv([o1, o2], t))

        gen = generate_io_function(iob, f_states=[o2, o1], warn=false)
        @test allequal(gen.states, strip_iv([o2, o1], t))

        gen = generate_io_function(iob, f_inputs=[i1, i2], warn=false)
        @test allequal(gen.inputs, strip_iv([i1, i2], t))
        gen = generate_io_function(iob, f_inputs=[i2, i1], warn=false)
        @test allequal(gen.inputs, strip_iv([i2, i1], t))

        gen = generate_io_function(iob, f_params=[a, b], warn=false)
        @test allequal(gen.params, strip_iv([a, b], t))
        gen = generate_io_function(iob, f_params=[b, a], warn=false)
        @test allequal(gen.params, strip_iv([b, a], t))

        # test with Symbols
        gen = generate_io_function(iob, f_states=[:o1, :o2], warn=false)
        @test allequal(gen.states, strip_iv([o1, o2], t))
        gen = generate_io_function(iob, f_states=[:o2, :o1], warn=false)
        @test allequal(gen.states, strip_iv([o2, o1], t))

        gen = generate_io_function(iob, f_inputs=[:i1, :i2], warn=false)
        @test allequal(gen.inputs, strip_iv([i1, i2], t))
        gen = generate_io_function(iob, f_inputs=[:i2, :i1], warn=false)
        @test allequal(gen.inputs, strip_iv([i2, i1], t))

        gen = generate_io_function(iob, f_params=[a, :b], warn=false)
        @test allequal(gen.params, strip_iv([a, b], t))
        gen = generate_io_function(iob, f_params=[:b, a], warn=false)
        @test allequal(gen.params, strip_iv([b, a], t))
    end
end


text:
# Autogenerated wrapper script for Blosc_jll for x86_64-apple-darwin
export libblosc

using Zlib_jll
using Zstd_jll
using Lz4_jll
JLLWrappers.@generate_wrapper_header("Blosc")
JLLWrappers.@declare_library_product(libblosc, "@rpath/libblosc.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Zstd_jll, Lz4_jll)
    JLLWrappers.@init_library_product(
        libblosc,
        "lib/libblosc.1.21.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Sample a particle from the prior model and then return it and the associated distance
"""
function sample_particle_distance(summaryStatistics::Any, sample_prior::Function, forward_model::Function, compute_distance::Function)
    proposal = sample_prior()
    simulatedData = forward_model(proposal)
    return proposal, compute_distance(simulatedData, summaryStatistics)
end

function find_particle(summaryStatistics::Any, previousStep::A, threshold::G, kernelbandwidth::G, sample_kernel::Function, forward_model::Function, compute_distance::Function) where {A <: AbcPmcStep, G <: Real}

    sampled = 0
    while true
        sampled += 1

        #Draw Sample
        proposal = StatsBase.sample(previousStep)
        proposal = sample_kernel(proposal, kernelbandwidth)

        #Simulate data
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Accept / reject
        if proposalDistance < threshold
            return proposal, proposalDistance, sampled
        end
    end
end

function find_particle(summaryStatistics::Any, previousStep::A, threshold::Array{G, 1}, kernelbandwidth::G, sample_kernel::Function, forward_model::Function, compute_distance::Function) where {A <: AbcPmcStep, G <: Real}

    sampled = 0
    while true
        sampled += 1

        #Draw Sample
        proposal = StatsBase.sample(previousStep)
        proposal = sample_kernel(proposal, kernelbandwidth)

        #Simulate data
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Accept / reject
        if all(proposalDistance .< threshold)
            return proposal, proposalDistance, sampled
        end
    end
end

function find_particle(summaryStatistics::Any, previousStep::A, threshold::G, kernelbandwidth::Array{G, 1}, sample_kernel::Function, forward_model::Function, compute_distance::Function) where {A <: AbcPmcStep, G <: Real}

    sampled = 0
    while true
        sampled += 1

        #Draw Sample
        proposal = StatsBase.sample(previousStep)
        proposal = sample_kernel(proposal, kernelbandwidth)

        #Simulate data
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Accept / reject
        if proposalDistance < threshold
            return proposal, proposalDistance, sampled
        end
    end
end

function find_particle(summaryStatistics::Any, previousStep::A, threshold::Array{G, 1}, kernelbandwidth::Array{G, 1}, sample_kernel::Function, forward_model::Function, compute_distance::Function) where {A <: AbcPmcStep, G <: Real}

    sampled = 0
    while true
        sampled += 1

        #Draw Sample
        proposal = StatsBase.sample(previousStep)
        proposal = sample_kernel(proposal, kernelbandwidth)

        #Simulate data
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Accept / reject
        if all(proposalDistance .< threshold)
            return proposal, proposalDistance, sampled
        end
    end
end

"""
Initialize ABC PMC Algorithm
"""
function ppmc_start(summaryStatistics::Any,
                   nparticles::Int64, ninitial::Int64,
                   sample_prior::Function, forward_model::Function,
                   compute_distance::Function, rank_distances::Function,
                   verbose::Bool = true)

    ##Draw first sample to determine typing
    proposal, dist = sample_particle_distance(summaryStatistics, sample_prior, forward_model, compute_distance)

    ##Initialize Arrays based on draws
    if typeof(proposal) <: Array
        particles = Array{eltype(proposal)}(ninitial, length(proposal))
    else
        particles = Array{eltype(proposal)}(ninitial)
    end

    if typeof(dist) <: Array
        distances = Array{eltype(dist)}(ninitial, length(dist))
    else
        distances = Array{eltype(dist)}(ninitial)
    end

    particles[1, :] = proposal
    distances[1, :] = dist

    ##Calculate initial samples
    if verbose println("Sampling Initial Particles and Computing Distances...") end
    np = nprocs()
    ii = 1
    @sync begin
        for proc = 1:np
            if proc != myid() || np == 1
                @async begin
                    while true
                        idx = (idx = ii; ii += 1; idx)
                        if idx > ninitial
                            break
                        end
                        particles[idx, :], distances[idx, :] = remotecall_fetch(sample_particle_distance, proc, summaryStatistics, sample_prior, forward_model, compute_distance)
                    end
                end
            end
        end
    end

    ##Rank Particles
    if verbose println("Selecting Particles...") end
    ranks = rank_distances(distances)
    sortedPermutation = sortperm(ranks)

    ##Subset Particles
    idxs = sortedPermutation[1:nparticles]
    if typeof(proposal) <: Array
        particles = particles[idxs, :]
    else
        particles = particles[idxs]
    end

    if typeof(dist) <: Array
        distances = distances[idxs, :]
        threshold = vec(maximum(distances, 1))
    else
        distances = distances[idxs]
        threshold = maximum(distances)
    end

    ##Return result
    weights = StatsBase.AnalyticWeights(fill(1.0 / nparticles, nparticles))
    nsampled = fill(floor(Int64, ninitial / nparticles), nparticles)
    return AbcPmc(particles, distances, weights, threshold, nsampled)
end

"""
Function to new step in ABC Population Monte Carlo Algorithm

# Arguments
* `previousStep::ABCPMCStep` object containing information on previous step
* `shrink_threshold::Function` rule for computing reduced acceptance threshold
* `transition_kernel::Function` transition kernel for adding noise to previously
 accepted parameter values
* `forward_model::Function` simulates new data given parameter value
* `compute_distance::Function` calculates a distance measure between output of
 `forawrd_model` and supplied summary statistics
* `density_kernel::Function` density value for transition kernel returning
likelihood for new parameter values given old parameter values
* `density_prior::Function` evaluates density of parameter value for prior function
* `verbose::Bool = true` default value is true, should phase of step be reported
# Value
ABCPMCStep object containing new posterior distribution, parameters and weights, as well
as computed distances, total number of parameters sampled, and acceptance threshold that
was used.
"""
function ppmc_step(previousStep::A, summaryStatistics::Any, nparticles::Int64,
                  kernel_bandwidth::Function, shrink_threshold::Function,
                  sample_kernel::Function, forward_model::Function,
                  compute_distance::Function, density_kernel::Function,
                  density_prior::Function, verbose::Bool = true) where A <: AbcPmcStep

    #Allocate variables
    newParticles = zeros(previousStep.particles)
    newDistances = zeros(previousStep.distances)
    nsampled = zeros(previousStep.nsampled)

    #Compute kernel variance / std
    if verbose println("Calculating Bandwidth...") end
    kernelbandwidth = kernel_bandwidth(previousStep)

    #Compute new threshold
    if verbose println("Calculating Threshold...") end
    threshold = shrink_threshold(previousStep)

    #Find new particles
    if verbose println("Sampling Particles...") end
    accepted = 1
    np = nprocs()
    @sync begin
        for proc = 1:np
            if proc != myid() || np == 1
                @async begin
                    while true
                        idx = (idx = accepted; accepted += 1; idx)
                        if idx > nparticles
                            break
                        end
                        newParticles[idx, :], newDistances[idx, :], nsampled[idx] = remotecall_fetch(find_particle, proc, summaryStatistics, previousStep, threshold, kernelbandwidth, sample_kernel, forward_model, compute_distance)
                    end
                end
            end
        end
    end

    if verbose  println("Calculating Weights...") end
    newWeights = kernel_weights(newParticles, previousStep, kernelbandwidth,
                                density_kernel, density_prior)

    #Return result as AbcPmcStep
    return AbcPmc(newParticles, newDistances, newWeights, threshold, nsampled)
end

"""
Population Monte Carlo Approximate Bayesian Computation Algorithm

# Arguments
* `summaryStatistics` data for posterior to be conditioned on
* `steps::Int64` number of iterations to run
* `nparticles::Int64` number of posterior samples to return for each iteration
* `ninitial::Int64` number of samples to base initial iterationon
* `sample_prior::Function` generates a random draw from the prior
* `density_prior::Function` pdf of prior
* `compute_distance::Function` compute distance metric between output of
  `forward_model` and provided summary statistics
* `forward_model::Function` function to simulate summary statistics according to new
parameter value
* `transition_kernel::Function` transition kernel or "noise" function applied to particles
* `density_kernel::Function` pdf functions for transition kernel
* `rank_distances::Function` method for ranking particles based on computed distances and
  choosing initial selection threshold
* `shrink_threshold::Function` method for shrinking threshold between iterations based on output
  from previous iteration
* `log::Bool = true` should progress be written to a log file
* `logFile::String = "log.txt"` file to write progress to, ignored if `log = false`
* `save::Bool = true` should results be written to a file
* `saveFile::String = "results.jld"` file to write results to
* `verbose::Bool = true` should progress be reported
"""
function pabc_pmc(summaryStatistics::Any, nsteps::Int64,
                 nparticles::Int64, ninitial::Int64,
                 sample_prior::Function, density_prior::Function,
                 sample_kernel::Function, density_kernel::Function,
                 forward_model::Function, compute_distance::Function,
                 rank_distances::Function, kernel_bandwidth::Function,
                 shrink_threshold::Function; verbose::Bool = true,
                 log::Bool = true, logFile::String = "log.txt",
                 save::Bool = true, saveFile::String = "results.jld")

    ##Write initial information is running log file
    if log
        startTime = now()
        initlog(logFile, startTime, ninitial, nparticles, nsteps)
    end

    ##Run first step
    results = [ppmc_start(summaryStatistics,
                                  nparticles, ninitial,
                                  sample_prior, forward_model,
                                  compute_distance, rank_distances,
                                  verbose)]
    ##Save results
    if verbose println("Saving Step...") end
    if save JLD.@save saveFile summaryStatistics results end

    ##Print to log file
    if verbose println(string("Step ", 1, " Complete")) end

    if log
        stepTime = now() - startTime
        totalTime = stepTime
        steplog(logFile, 1, totalTime, stepTime, results[1].acceptbw, ninitial)
    end

    #Subsequent steps
    for ii in 2:nsteps
        push!(results, ppmc_step(results[ii - 1], summaryStatistics, nparticles,
                                kernel_bandwidth, shrink_threshold, sample_kernel,
                                forward_model, compute_distance, density_kernel,
                                density_prior, verbose))

        ##Save results
        if verbose println("Saving Step...") end
        if save JLD.@save saveFile summaryStatistics results end

        #Print to log file
        if verbose println(string("Step ", ii, " Complete")) end

        if log
            stepTime = now() - startTime - totalTime
            totalTime += stepTime
            steplog(logFile, ii, totalTime, stepTime, results[ii].acceptbw, sum(results[ii].nsampled))
        end
    end
    return results
end

function pabc_pmc_warmstart(nsteps::Int64, nparticles::Int64,
                            sample_prior::Function, density_prior::Function,
                            sample_kernel::Function, density_kernel::Function,
                            forward_model::Function, compute_distance::Function,
                            kernel_bandwidth::Function, shrink_threshold::Function;
                            verbose::Bool = true,
                            log::Bool = true, logFile::String = "log.txt",
                            save::Bool = true, saveFile::String = "results.jld")

    ##Print to log file
    if verbose println("Resuming algorithm...") end

    ##Write initial information is running log file
    if log
        startTime = now()
        totalTime = now() - startTime
    end

    ##Save results
    summaryStatistics = JLD.load(saveFile, "summaryStatistics")
    results = JLD.load(saveFile, "results")

    ##Check Number of steps
    if nsteps <= length(results)
        error("$nsteps alreadty complete")
    end

    #Subsequent steps
    for ii in length(results) + 1:nsteps
        push!(results, ppmc_step(results[ii - 1], summaryStatistics, nparticles,
                                kernel_bandwidth, shrink_threshold, sample_kernel,
                                forward_model, compute_distance, density_kernel,
                                density_prior, verbose))

        ##Save results
        if verbose println("Saving Step...") end
        if save JLD.@save saveFile results end

        #Print to log file
        if verbose println(string("Step ", ii, " Complete")) end

        if log
            stepTime = now() - startTime - totalTime
            totalTime += stepTime
            steplog(logFile, ii, totalTime, stepTime, results[ii].acceptbw, sum(results[ii].nsampled))
        end
    end
    return results
end


text:
module lineq
using LinearAlgebra

function linear_index(i, j, sz)
    n_row = sz[1]
    return i + (j - 1) * n_row
end

function matrix_index(k, sz)
    n_row = sz[1]
    (j, r) = divrem(k, n_row)
    j = (r>0) + j
    i = k - (j-1) * n_row
    return (i, j)
end

struct JacobiOptions
    max_iter::Int
    n_iter_norm_check::Int
    delta_x_tol::Float64
    damping::Float64
end
function JacobiOptions()
    return JacobiOptions(100, 10, 1e-6, 1.0)
end
function JacobiOptions(;
    max_iter = 100,
    n_iter_norm_check = 10,
    delta_x_tol = 1e-15,
    damping = 1.0)

    return JacobiOptions(
        max_iter,
        n_iter_norm_check,
        delta_x_tol,
        damping)
end

function delta_x_norm_check(x_new, x_old, options)
    break_flag = false
    norm_dx = norm(x_new - x_old)
    if norm_dx <= options.delta_x_tol
        break_flag = true
    end
    return break_flag
end

function jacobi_iteration(A, b, options::JacobiOptions = JacobiOptions())

    n_dim = length(b)
    x_old = zeros(eltype(b), n_dim)
    x_new = zeros(eltype(b), n_dim)

    d_vec = diag(A)

    use_damping = options.damping == 1.0

    iter_norm_check = 0
    iter = 0
    for outer iter = 1:options.max_iter
        iter_norm_check += 1

        # stepping algorithm
        x_new = ( b + d_vec .* x_old - (A * x_old) ) ./ d_vec

        if use_damping
            x_new = x_new * options.damping + x_old * (1 - options.damping)
        end

        # check for stopping conditions
        if iter_norm_check >= options.n_iter_norm_check
            iter_norm_check = 0
            break_flag = delta_x_norm_check(x_new, x_old, options)
            if break_flag
                break
            end
        end

        # update x_old
        x_old = x_new
    end

    return (x_new, iter)
end

struct GaussSeidelOptions
    max_iter::Int
    n_iter_norm_check::Int
    delta_x_tol::Float64
    damping::Float64
end
function GaussSeidelOptions()
    return GaussSeidelOptions(100, 10, 1e-6, 1.0)
end
function GaussSeidelOptions(;
    max_iter = 100,
    n_iter_norm_check = 10,
    delta_x_tol = 1e-15,
    damping = 1.0)

    return GaussSeidelOptions(
        max_iter,
        n_iter_norm_check,
        delta_x_tol,
        damping)
end

function gauss_seidel_iteration(A, b, options::GaussSeidelOptions)

    n_dim = length(b)
    x_old = zeros(eltype(b), n_dim)
    x_new = zeros(eltype(b), n_dim)

    d_vec = diag(A)

    use_damping = options.damping == 1.0

    iter_norm_check = 0
    iter = 0
    for outer iter = 1:options.max_iter
        iter_norm_check += 1

        # stepping algorithm

        if use_damping
            x_new = x_new * options.damping + x_old * (1 - options.damping)
        end

        # check for stopping conditions
        if iter_norm_check >= options.n_iter_norm_check
            iter_norm_check = 0
            break_flag = delta_x_norm_check(x_new, x_old, options)
            if break_flag
                break
            end
        end

    end

    return (x_new, iter)
end

end # dev

module dtest
using ..lineq
using LinearAlgebra

function create_matrix(n_dim = 8)

    A = zeros(Float64, n_dim, n_dim)

    main_diag = [1; 4 * ones(n_dim-2); 4]
    side_diag = -ones(n_dim-1)
    sideside_diag = -ones(n_dim-3)

    main_diag_idx = diagind(A, 0)
    side_m1_diag_idx = diagind(A, 1)
    side_p1_diag_idx = diagind(A, -1)

    side_m3_diag_idx = diagind(A, 3)
    side_p3_diag_idx = diagind(A, -3)

    A[ main_diag_idx ] = main_diag
    A[ side_m1_diag_idx ] = side_diag
    A[ side_p1_diag_idx ] = side_diag
    A[ side_m3_diag_idx ] = sideside_diag
    A[ side_p3_diag_idx ] = sideside_diag

    return A
end

function create_vector(n_dim = 8)

    b = zeros(n_dim)
    b[1]   = -1.0
    b[end] = 2.0

    return b
end

n_dim = 8 # 100
A = create_matrix(n_dim)
b = create_vector(n_dim)

x_theo = A \ b

# jac_options = lineq.JacobiOptions(delta_x_tol = 1e-10, max_iter = 10000, damping = 1.0)
# (x_jac, iter) = lineq.jacobi_iteration(A, b, jac_options)

# println("Norm diff between theo and jac: $(norm(x_theo - x_jac))")
# println("Needed iterations: $(iter)")

gs_options = lineq.GaussSeidelOptions(delta_x_tol = 1e-10, max_iter = 100, damping = 1.0)
(x_gs, iter) = lineq.gauss_seidel_iteration(A, b, gs_options)

println("Norm diff between theo and gs: $(norm(x_theo - x_gs))")
println("Needed iterations: $(iter)")

D = diagm( diag(A) )
E = -LowerTriangular(A)
F = -UpperTriangular(A)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# ---
# title: Monogenic filters
# id: demo_phasecongmono
# cover: assets/phasecongmono.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

# Phase congruency marks all classes of features from steps to lines and is a dimensionless
# quantity that ranges from 0 to 1. This allows fixed thresholds to be used over wide
# classes of images.

using TestImages
using Images
using ImagePhaseCongruency

img = restrict(testimage("mandril_gray"))

(pc, or, ft, T) = phasecongmono(img)
nonmax = Images.thin_edges(pc, or)

mosaic(
    img,
    adjust_histogram(pc, LinearStretching()),
    nonmax,
    hysthresh(nonmax, 0.1, 0.2);
    nrow=2, rowmajor=true
)

# Images: 1) top left: original image 2) top right: phase congruency 3) bottom left:
# non-maximal suppression 4) bottom right: Hystersis thresholded

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "phasecongmono.png"), adjust_histogram(Gray.(pc), LinearStretching())) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
module ConvexHull

using JuMP

export homogeneous_system, read_ine, write_ine, read_ext, write_ext,
       double_description, canonicalize!,
       get_extrema, is_approx_included

importall Base

const ε = 10eps()

include("double_description.jl")
include("JuMP_interface.jl")
include("readers_writers.jl")

end


text:
# many I/O functions have UNIX names

touch("output.txt")
mkdir("docs")

# probably don't have permission
try
    touch("/output.txt")
    mkdir("/docs")
catch e
    warn(e)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    TravisCI(; config_file::Union{AbstractString, Nothing}="") -> TravisCI

Add `TravisCI` to a template's plugins to add a `.travis.yml` configuration file to
generated repositories, and an appropriate badge to the README.

# Keyword Arguments:
* `config_file::Union{AbstractString, Nothing}=""`: Path to a custom `.travis.yml`.
  If `nothing` is supplied, no file will be generated.
"""
struct TravisCI <: GenericPlugin
    gitignore::Vector{String}
    src::Union{String, Nothing}
    dest::String
    badges::Vector{Badge}
    view::Dict{String, Any}

    function TravisCI(; config_file::Union{AbstractString, Nothing}="")
        if config_file != nothing
            config_file = if isempty(config_file)
                config_file = joinpath(DEFAULTS_DIR, "travis.yml")
            elseif isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError("File $(abspath(config_file)) does not exist"))
            end
        end
        new(
            [],
            config_file,
            ".travis.yml",
            [
                Badge(
                    "Build Status",
                    "https://travis-ci.com/{{USER}}/{{PKGNAME}}.jl.svg?branch=master",
                    "https://travis-ci.com/{{USER}}/{{PKGNAME}}.jl",
                ),
            ],
            Dict{String, Any}(),
        )
    end
end

interactive(::Type{TravisCI}) = interactive(TravisCI; file="travis.yml")


text:
"""
    Multivector1{T <: Real} <: AbstractMultivector{T}

An immutable 2-tuple of real numbers that represents a member of a 1-dimensional multivector algebra.

Each `Multivector1` has the form
```math
    a+bA
```
where ``a`` and ``b`` are real (and of the same type), and ``A * A = 0``.
Here ``*`` is the wedge product. Note that
```math
    (a + bA) * (a + bA) = a^{2} + 2abA
```
which, in general, is not equal to zero.
"""
struct Multivector1{T <: Real} <: AbstractMultivector{T}
    l::T
    r::T

    Multivector1{U}(l::U, r::U) where {U <: Real} = new(l, r)
end

Multivector1(a::T, b::T) where {T <: Real} = Multivector1{T}(a, b)
Multivector1(a::Real, b::Real) = Multivector1(promote(a, b)...)
Multivector1(a::T) where {T <: Real} = Multivector1{T}(a, zero(T))

function show(io::IO, z::Multivector1)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", A: ")
    print(io, z.r)
    print(io, "]")
end

function real(z::Multivector1)
    z.l
end

function unreal(a::Real)
    zero(a)
end

function unreal(z::Multivector1)
    z.r
end

function zero(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(zero(z.l), zero(z.r))
end

function zero(::Type{Multivector1{T}}) where {T <: Real}
    Multivector1{T}(zero(T), zero(T))
end

function one(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(one(z.l), zero(z.r))
end

function one(::Type{Multivector1{T}}) where {T <: Real}
    Multivector1{T}(one(T), zero(T))
end

"""
    conj{T <: Real}(z::Multivector1{T})

The `Multivector1` conjugate. If ``z=a+bA``, then `conj(z)` gives
```math
    a-bA
```
This operation is an involution.
"""
function conj(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(z.l, -z.r)
end

"""
    cloak{T <: Real}(z::Multivector1{T})

The cloak conjugate changes the sign of even blades. If ``z=a+bA``, then `cloak(z)` gives
```math
    -a+bA
```
This operation is equivalent to `-conj(z)` and thus is also an involution.
"""
function cloak(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(-z.l, z.r)
end

"""
    dagger{T <: Real}(z::Multivector1{T})

The dagger conjugate changes the sign of odd blades. If ``z=a+bA``, then `dagger(z)` gives
```math
    a-bA
```
This operation is equivalent to `conj(z)` and thus is also an involution.
"""
function dagger(z::Multivector1)
    conj(z)
end

"""
    star{T <: Real}(z::Multivector1{T})

Returns the Hodge star conjugate. If ``z=a+bA``, then `star(z)` gives
```math
    b+aA
```
This operation is an involution.
"""
function star(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(z.r, z.l)
end

"""
    selfstar(z::Multivector1)

The self-star-conjugate part.
This operation is idempotent.
"""
function selfstar(z::Multivector1)
    (z + star(z)) / 2
end

"""
    antiselfstar(z::Multivector1)

The anti-self-star-conjugate part.
This operation is idempotent.
"""
function antiselfstar(z::Multivector1)
    (z - star(z)) / 2
end

function (+)(x::Multivector1, y::Multivector1)
    Multivector1(x.l + y.l, x.r + y.r)
end

function (+)(z::Multivector1, a::Real)
    Multivector1(z.l + a, z.r)
end

function (+)(a::Real, z::Multivector1)
    Multivector1(z.l + a, z.r)
end

function (-)(z::Multivector1{T}) where {T <: Real}
    Multivector1{T}(-z.l, -z.r)
end

function (-)(x::Multivector1, y::Multivector1)
    Multivector1(x.l - y.l, x.r - y.r)
end

function (-)(z::Multivector1, a::Real)
    Multivector1(z.l - a, z.r)
end

function (-)(a::Real, z::Multivector1)
    Multivector1(a - z.l, -z.r)
end

"""
    (*)(x::Multivector1, y::Multivector1)

Wedge product of 1-dimensional multivectors.
This operation is commutative and associative.
"""
function (*)(x::Multivector1, y::Multivector1)
    Multivector1(x.l * y.l, (x.r * y.l) + (y.r * x.l))
end

"""
    (*)(z::Multivector1, a::Real)
    (*)(a::Real, z::Multivector1)

Scaling and/or reflection of a `Multivector1` by a real number.
"""
function (*)(z::Multivector1, a::Real)
    Multivector1(z.l * a, z.r * a)
end

function (*)(a::Real, z::Multivector1)
    Multivector1(a * z.l, a * z.r)
end

function (/)(x::Multivector1, y::Multivector1)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    # Naive division algorithm
    Multivector1(x.l * y.l, (x.r * y.l) - (y.r * x.l)) / abs2(y)
end

function (\)(y::Multivector1, x::Multivector1)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    # Naive division algorithm
    Multivector1(x.l * y.l, (x.r * y.l) - (y.r * x.l)) / abs2(y)
end

function (/)(z::Multivector1, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Multivector1(z.l / a, z.r / a)
end

function (\)(a::Real, z::Multivector1)
    if iszero(a)
        error(ZeroDenominator)
    end

    Multivector1(a \ z.l, a \ z.r)
end

function random(::Type{Multivector1{T}}) where T <: Real
    Multivector1{T}(rand(T), rand(T))
end


text:
using DiffEqBase, DiffEqBiological, DiffEqJump, Plots, Statistics, DataFrames
gr()
fmt = :png

N = 256
h = 1 / N
rn = @empty_reaction_network
function getDiffNetwork!(rn,N)
    for i = 1:N
        addspecies!(rn, Symbol(:u, i))
    end
    addparam!(rn, :β)
    for i = 1:N
        (i < N) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i+1)=>1,))
        (i > 1) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i-1)=>1,))
    end
    rn
end
getDiffNetwork!(rn,N)
addjumps!(rn, build_regular_jumps=false, minimal_jumps=true)
rnpar = [1/(h*h)]
u0 = 10*ones(Int64, N)
tf = .01

methods = (Direct(),DirectFW(),SortingDirect(),NRM(),DirectCR(),RSSA())
shortlabels = [string(leg)[12:end-2] for leg in methods]
prob    = prob = DiscreteProblem(u0, (0.0, tf), rnpar)
ploth   = plot(reuse=false)
for (i,method) in enumerate(methods)
    println("Benchmarking method: ", method)
    jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
    plot!(ploth,sol.t,sol[Int(N//2),:],label=shortlabels[i], format=fmt)
end
plot!(ploth, title="Population at middle lattice site", xlabel="time",format=fmt)

function run_benchmark!(t, jump_prob, stepper)
    sol = solve(jump_prob, stepper)
    @inbounds for i in 1:length(t)
        t[i] = @elapsed (sol = solve(jump_prob, stepper))
    end
end

nsims = 50
benchmarks = Vector{Vector{Float64}}()
for method in methods
    jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
    stepper = SSAStepper()
    t = Vector{Float64}(undef,nsims)
    run_benchmark!(t, jump_prob, stepper)
    push!(benchmarks, t)
end

medtimes = Vector{Float64}(undef,length(methods))
stdtimes = Vector{Float64}(undef,length(methods))
avgtimes = Vector{Float64}(undef,length(methods))
for i in 1:length(methods)
    medtimes[i] = median(benchmarks[i])
    avgtimes[i] = mean(benchmarks[i])
    stdtimes[i] = std(benchmarks[i])
end

df = DataFrame(names=shortlabels,medtimes=medtimes,relmedtimes=(medtimes/medtimes[1]),
                avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)

sa = [string(round(mt,digits=4),"s") for mt in df.medtimes]
bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
ylabel!("median relative to Direct")
title!("256 Site 1D Diffusion CTRW")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
@testset "yen" begin
    img = testimage("cameraman")
    img2 = copy(img)
    img_binary = binarize(Yen(), img)

    #check original img is unchanged
    @test img == img2

    #check that the image only has ones or zeros
    non_zeros = findall(x -> x != 0.0 && x != 1.0, img_binary)
    @test length(non_zeros) == 0

    #check type of outputed img
    @test typeof(img_binary) == typeof(img)

    #check ones and zeros have been assigned to the correct side of the threshold
    max,maxpos=findmax(img)
    @test img_binary[maxpos] == 1
    min,minpos=findmin(img)
    @test img_binary[minpos] == 0
end


text:
α = 0.3
δ = 0.1
β = 0.9
k_l = 2
k_u = 0.001
n = 10
kgrid = collect(range(k_l, stop = k_u, length =n))

v = zeros(n)
tolerance = 0.001
imax = 1000
iter = 1
vnew = zeros(n)
while maximum(abs.(v-vnew) > tolerance) && iter < imax
v=vnew
c = zeros(n,n)
for i in 1:n
    for j in 1:n
    c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j]
        if c[i,j]<0
            c[i,j]=0
        end
    end
end
vnew = maximum(log.(c) .+ β*v' , dim = 2)
iter += 1
end

max_stars_repo_path: Dynamic Programming/Neoclassical.jl

max_issues_repo_path: Dynamic Programming/Neoclassical.jl

max_forks_repo_path: Dynamic Programming/Neoclassical.jl


text:
module EnsembleAnalysis

using DiffEqBase, StaticArrays, Statistics, RecursiveArrayTools

# Getters
get_timestep(sim,i) = (getindex(sol,i) for sol in sim)
get_timepoint(sim,t) = (sol(t) for sol in sim)
function componentwise_vectors_timestep(sim,i)
  arr = [get_timestep(sim,i)...]
  if typeof(arr[1]) <: AbstractArray
   return vecarr_to_vectors(VectorOfArray(arr))
  else
   return arr
  end
 end
 function componentwise_vectors_timepoint(sim,t)
   arr = [get_timepoint(sim,t)...]
   if typeof(arr[1]) <: AbstractArray
    return vecarr_to_vectors(VectorOfArray(arr))
   else
    return arr
   end
  end

# Timestep statistics
timestep_mean(sim,i) = componentwise_mean(get_timestep(sim,i))
timestep_mean(sim,::Colon) = timeseries_steps_mean(sim)
function timestep_median(sim,i)
  arr = componentwise_vectors_timestep(sim,i)
  if typeof(first(arr)) <: AbstractArray
    return reshape([median(x) for x in arr],size(sim[1][i])...)
  else
    return median(arr)
  end
end
timestep_median(sim,::Colon) = timeseries_steps_median(sim)
function timestep_quantile(sim,q,i)
  arr = componentwise_vectors_timestep(sim,i)
  if typeof(first(arr)) <: AbstractArray
    return reshape([quantile(x,q) for x in arr],size(sim[1][i])...)
  else
    return quantile(arr,q)
  end
end
timestep_quantile(sim,q,::Colon) = timeseries_steps_quantile(sim,q)
timestep_meanvar(sim,i) = componentwise_meanvar(get_timestep(sim,i))
timestep_meanvar(sim,::Colon) = timeseries_steps_meanvar(sim)
timestep_meancov(sim,i,j) = componentwise_meancov(get_timestep(sim,i),get_timestep(sim,j))
timestep_meancov(sim,::Colon,::Colon) = timeseries_steps_meancov(sim)
timestep_meancor(sim,i,j) = componentwise_meancor(get_timestep(sim,i),get_timestep(sim,j))
timestep_meancor(sim,::Colon,::Colon) = timeseries_steps_meancor(sim)
timestep_weighted_meancov(sim,W,i,j) = componentwise_weighted_meancov(get_timestep(sim,i),get_timestep(sim,j),W)
timestep_weighted_meancov(sim,W,::Colon,::Colon) = timeseries_steps_weighted_meancov(sim,W)

function timeseries_steps_mean(sim)
  DiffEqArray([timestep_mean(sim,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_median(sim)
  DiffEqArray([timestep_median(sim,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_quantile(sim,q)
  DiffEqArray([timestep_quantile(sim,q,i) for i in 1:length(sim[1])],sim[1].t)
end
function timeseries_steps_meanvar(sim)
  m,v = timestep_meanvar(sim,1)
  means = [m]
  vars = [v]
  for i in 2:length(sim[1])
    m,v = timestep_meanvar(sim,i)
    push!(means,m)
    push!(vars,v)
  end
  DiffEqArray(means,sim[1].t),DiffEqArray(vars,sim[1].t)
end
function timeseries_steps_meancov(sim)
  reshape([timestep_meancov(sim,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end
function timeseries_steps_meancor(sim)
  reshape([timestep_meancor(sim,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end
function timeseries_steps_weighted_meancov(sim,W)
  reshape([timestep_meancov(sim,W,i,j) for i in 1:length(sim[1]) for j in 1:length(sim[1])],length(sim[1]),length(sim[1]))
end

timepoint_mean(sim,t) = componentwise_mean(get_timepoint(sim,t))
function timepoint_median(sim,t)
  arr = componentwise_vectors_timepoint(sim,t)
  if typeof(first(arr)) <: AbstractArray
    return reshape([median(x) for x in arr],size(sim[1][1])...)
  else
    return median(arr)
  end
end
function timepoint_quantile(sim,q,t)
  arr = componentwise_vectors_timepoint(sim,t)
  if typeof(first(arr)) <: AbstractArray
    return reshape([quantile(x,q) for x in arr],size(sim[1][1])...)
  else
    return quantile(arr,q)
  end
end
timepoint_meanvar(sim,t) = componentwise_meanvar(get_timepoint(sim,t))
timepoint_meancov(sim,t1,t2) = componentwise_meancov(get_timepoint(sim,t1),get_timepoint(sim,t2))
timepoint_meancor(sim,t1,t2) = componentwise_meancor(get_timepoint(sim,t1),get_timepoint(sim,t2))
timepoint_weighted_meancov(sim,W,t1,t2) = componentwise_weighted_meancov(get_timepoint(sim,t1),get_timepoint(sim,t2),W)

function DiffEqBase.EnsembleSummary(sim::DiffEqBase.AbstractEnsembleSolution{T,N},
                         t=sim[1].t;quantiles=[0.05,0.95]) where {T,N}
  if typeof(sim[1]) <: DiffEqBase.DESolution
    m,v = timeseries_point_meanvar(sim,t)
    qlow = timeseries_point_quantile(sim,quantiles[1],t)
    qhigh = timeseries_point_quantile(sim,quantiles[2],t)
  else
    m,v = timeseries_steps_meanvar(sim)
    qlow = timeseries_steps_quantile(sim,quantiles[1])
    qhigh = timeseries_steps_quantile(sim,quantiles[2])
  end

  trajectories = length(sim)
  DiffEqBase.EnsembleSummary{T,N,typeof(t),typeof(m),typeof(v),typeof(qlow),typeof(qhigh)}(t,m,v,qlow,qhigh,trajectories,sim.elapsedTime,sim.converged)
end

function timeseries_point_mean(sim,ts)
  DiffEqArray([timepoint_mean(sim,t) for t in ts],ts)
end
function timeseries_point_median(sim,ts)
  DiffEqArray([timepoint_median(sim,t) for t in ts],ts)
end
function timeseries_point_quantile(sim,q,ts)
  DiffEqArray([timepoint_quantile(sim,q,t) for t in ts],ts)
end
function timeseries_point_meanvar(sim,ts)
  m,v = timepoint_meanvar(sim,first(ts))
  means = [m]
  vars = [v]
  for t in Iterators.drop(ts,1)
    m,v = timepoint_meanvar(sim,t)
    push!(means,m)
    push!(vars,v)
  end
  DiffEqArray(means,ts),DiffEqArray(vars,ts)
end
function timeseries_point_meancov(sim,ts1,ts2)
  reshape([timepoint_meancov(sim,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end
function timeseries_point_meancor(sim,ts1,ts2)
  reshape([timepoint_meancor(sim,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end
function timeseries_point_weighted_meancov(sim,W,ts1,ts2)
  reshape([timepoint_meancov(sim,W,t1,t2) for t1 in ts1 for t2 in ts2],length(ts1),length(ts2))
end

function componentwise_mean(A)
  x0 = first(A)
  n = 0
  mean = zero(x0)
  for x in A
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      mean .+= x
    else
      mean += x
    end
  end
  if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
    mean ./= n
  else
    mean /= n
  end
  mean
end

# Welford algorithm
# https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
function componentwise_meanvar(A;bessel=true)
  x0 = first(A)
  n = 0
  mean = zero(x0)
  M2 = zero(x0)
  delta = zero(x0)
  delta2 = zero(x0)
  for x in A
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      delta .= x .- mean
      mean .+= delta./n
      delta2 .= x .- mean
      M2 .+= delta.*delta2
    else
      delta = x .- mean
      mean += delta./n
      delta2 = x .- mean
      M2 += delta.*delta2
    end
  end
  if n < 2
    return NaN
  else
    if bessel
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        M2 .= M2 ./ (n .- 1)
      else
        M2 = M2 ./ (n .- 1)
      end
    else
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        M2 .= M2 ./ n
      else
        M2 = M2 ./ n
      end
    end
    return mean,M2
  end
end

function componentwise_meancov(A,B;bessel=true)
  x0 = first(A)
  y0 = first(B)
  n = 0
  meanx = zero(x0)
  meany = zero(y0)
  C = zero(x0)
  dx = zero(x0)
  for (x,y) in zip(A,B)
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      dx .= x .- meanx
      meanx .+= dx./n
      meany .+= (y.-meany)./n
      C .+= dx .* (y .- meany)
    else
      dx = x .- meanx
      meanx += dx./n
      meany += (y.-meany)./n
      C += dx .* (y .- meany)
    end
  end
  if n < 2
    return NaN
  else
    if bessel
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (n .- 1)
      else
        C = C ./ (n .- 1)
      end
    else
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ n
      else
        C = C ./ n
      end
    end
    return meanx,meany,C
  end
end

function componentwise_meancor(A,B;bessel=true)
  mx,my,cov = componentwise_meancov(A,B;bessel=bessel)
  mx,vx = componentwise_meanvar(A;bessel=bessel)
  my,vy = componentwise_meanvar(B;bessel=bessel)
  if typeof(vx) <: AbstractArray
    vx .= sqrt.(vx)
    vy .= sqrt.(vy)
  else
    vx = sqrt.(vx)
    vy = sqrt.(vy)
  end
  mx,my,cov./(vx.*vy)
end

function componentwise_weighted_meancov(A,B,W;weight_type=:reliability)
  x0 = first(A)
  y0 = first(B)
  w0 = first(W)
  n = 0
  meanx = zero(x0)
  meany = zero(y0)
  wsum = zero(w0)
  wsum2 = zero(w0)
  C = zero(x0)
  dx = zero(x0)
  for (x,y,w) in zip(A,B,W)
    n += 1
    if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
      wsum .+= w
      wsum2 .+= w.*w
      dx .= x .- meanx
      meanx .+= (w ./ wsum) .* dx
      meany .+= (w ./ wsum) .* (y .- meany)
      C .+= w .* dx .* (y .- meany)
    else
      wsum += w
      wsum2 += w.*w
      dx = x .- meanx
      meanx += (w ./ wsum) .* dx
      meany += (w ./ wsum) .* (y .- meany)
      C += w .* dx .* (y .- meany)
    end
  end
  if n < 2
    return NaN
  else
    if weight_type == :population
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ wsum
      else
        C = C ./ wsum
      end
    elseif weight_type == :reliability
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (wsum .- wsum2 ./ wsum)
      else
        C = C ./ (wsum .- wsum2 ./ wsum)
      end
    elseif weight_type == :frequency
      if typeof(x0) <: AbstractArray && !(typeof(x0) <: SArray)
        C .= C ./ (wsum .- 1)
      else
        C = C ./ (wsum .- 1)
      end
    else
      error("The weight_type which was chosen is not allowed.")
    end
    return meanx,meany,C
  end
end

export get_timestep, get_timepoint, apply_timestep, apply_timepoint,
       componentwise_vectors_timestep, componentwise_vectors_timepoint

export componentwise_mean, componentwise_meanvar

export timestep_mean, timestep_median, timestep_quantile, timestep_meanvar,
       timestep_meancov, timestep_meancor, timestep_weighted_meancov

export timeseries_steps_mean, timeseries_steps_median, timeseries_steps_quantile,
       timeseries_steps_meanvar, timeseries_steps_meancov,
       timeseries_steps_meancor, timeseries_steps_weighted_meancov

export timepoint_mean, timepoint_median, timepoint_quantile,
       timepoint_meanvar, timepoint_meancov,
       timepoint_meancor, timepoint_weighted_meancov

export timeseries_point_mean, timeseries_point_median, timeseries_point_quantile,
       timeseries_point_meanvar, timeseries_point_meancov,
       timeseries_point_meancor, timeseries_point_weighted_meancov

end


text:
using Statistics
import .Plots; const pl = Plots

export plot_gp_intervals

"""
    plot_gp_intervals(gp::GaussianProcess, x::AbstractMatrix,
                      y::Union{AbstractMatrix,Nothing} = nothing; input_dim = 1,
                      output_dim = 1)

Creates a simple 95% uncertainty interval plot for the GP at a set of ordered
data points ``x``. Optionally also plots a set of data points ``y``. For
multivariateand vector GPs, plots the dimensions specified by
`input_dim` and `output_dim`.
"""
function plot_gp_intervals(gp::GaussianProcess, x::AbstractMatrix, y::Union{AbstractMatrix,Nothing} = nothing; input_dim = 1, output_dim = 1)
  f = gp(x)[output_dim,:,:]
  n_samples = min(size(f,ndims(f)),32)

  x_d = x[input_dim,:]
  m = mean(f;dims=2) |> vec
  u = mapslices(v->quantile(v,0.975), f; dims=2) |> vec
  u_2sd = mean(f;dims=2) + mapslices(v->1.96*sqrt(var(v)), f; dims=2) |> vec
  l = mapslices(v->quantile(v,0.025), f; dims=2) |> vec
  l_2sd = mean(f;dims=2) - mapslices(v->1.96*sqrt(var(v)), f; dims=2) |> vec

  p = pl.plot(size=(800,600))

  y isa Nothing ? () : pl.scatter!(p, x_d, y[output_dim,:], label=nothing)
  pl.plot!(p, x_d, m, ribbon = (abs.(l.-m),abs.(u.-m)), label=nothing)
  for i in 1:n_samples
    pl.plot!(p, x_d, f[:,i], color=:gray, label=nothing)
  end

  p
end


text:
function mk_cv_tbl(dts, dvs, gts, gvs, nts, nvs, G, trait, erc)
    # target is each country in turn
    rst = zeros(6, 7)           # 6 = w/nw by (DGN); 7 = D G N DG DN GN DGN

    @info "Prediction of one country with one of 3 coutries"
    # D → D
    rst[1, 1], rst[4, 1] = cv_G(dts, dvs, G, trait, erc)
    # G → D
    rst[1, 2], rst[4, 2] = cv_G(gts, dvs, G, trait, erc)
    # N → D
    rst[1, 3], rst[4, 3] = cv_G(nts, dvs, G, trait, erc)
    # D → G
    rst[2, 1], rst[5, 1] = cv_G(dts, gvs, G, trait, erc)
    # G → G
    rst[2, 2], rst[5, 2] = cv_G(gts, gvs, G, trait, erc)
    # N → G
    rst[2, 3], rst[5, 3] = cv_G(nts, gvs, G, trait, erc)
    # D → G
    rst[3, 1], rst[6, 1] = cv_G(dts, nvs, G, trait, erc)
    # G → G
    rst[3, 2], rst[6, 2] = cv_G(gts, nvs, G, trait, erc)
    # N → G
    rst[3, 3], rst[6, 3] = cv_G(nts, nvs, G, trait, erc)

    @info "Prediction of one country with two of 3 coutries"
    dgts = vcat(dts, gts)
    dnts = vcat(dts, nts)
    gnts = vcat(gts, nts)
    # D + G → D
    rst[1, 4], rst[4, 4] = cv_G(dgts, dvs, G, trait, erc)
    # D + N → D
    rst[1, 5], rst[4, 5] = cv_G(dnts, dvs, G, trait, erc)
    # G + N → D
    rst[1, 6], rst[4, 6] = cv_G(gnts, dvs, G, trait, erc)
    # D + G → G
    rst[2, 4], rst[5, 4] = cv_G(dgts, gvs, G, trait, erc)
    # D + N → G
    rst[2, 5], rst[5, 5] = cv_G(dnts, gvs, G, trait, erc)
    # G + N → G
    rst[2, 6], rst[5, 6] = cv_G(gnts, gvs, G, trait, erc)
    # D + G → N
    rst[3, 4], rst[6, 4] = cv_G(dgts, nvs, G, trait, erc)
    # D + N → N
    rst[3, 5], rst[6, 5] = cv_G(dnts, nvs, G, trait, erc)
    # G + N → N
    rst[3, 6], rst[6, 6] = cv_G(gnts, nvs, G, trait, erc)
    @info "Prediction of one country with all 3 coutries"
    dgnts = vcat(dts, gts, nts)
    # D + G + N → D
    rst[1, 7], rst[4, 7] = cv_G(dgnts, dvs, G, trait, erc)
    # D + G + N → G
    rst[2, 7], rst[5, 7] = cv_G(dgnts, gvs, G, trait, erc)
    # D + G + N → N
    rst[3, 7], rst[6, 7] = cv_G(dgnts, nvs, G, trait, erc)

    round.(rst, sigdigits = 2)
end

"""
Test all the combinations to predict a 3rd country with
- 3 x 1→1
- 3 x 2→2
- 1 x 3→1
## DRP
- drp-training-2021-04-08.jld: with Norwegian DRP updated
- drp-training-2021-04-21.jld: with german data updated
"""
function cv_21_04_22(drp)
    @load "$dat_dir/jld/cv-setup.jld" dcs gcs ncs
    @load "$dat_dir/jld/$drp"         dts gts nts
    @load "$dat_dir/jld/ebv-all.jld"  dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"     G A ID

    begin                       # The milk trait
        dt, dv = mk_milk_dt(dts, dcs, dbv)
        gt, gv = mk_milk_dt(gts, gcs, gbv)
        nt, nv = mk_milk_dt(nts, ncs, nbv)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Milk, :Milk_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end

    begin                       # the SCS trait
        dt, dv = mk_scs_dt(dts, dcs, dbv)
        gt, gv = mk_scs_dt(gts, gcs, gbv)
        nt, nv = mk_scs_dt(nts, ncs, nbv)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :SCS, :SCS_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
end

"""
    function cv_21_04_22_2(drp)
---
Include the validation set also
"""
function cv_21_04_22_2(drp)
    @load "$dat_dir/jld/cv-setup.jld" dcs gcs ncs
    @load "$dat_dir/jld/$drp"         dts gts nts
    @load "$dat_dir/jld/ebv-all.jld"  dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"     G A ID
    G += 0.01I

    begin                       # The milk trait
        dt, dv = mk_milk_dt_2(dts, dcs, dbv)
        gt, gv = mk_milk_dt_2(gts, gcs, gbv)
        nt, nv = mk_milk_dt_2(nts, ncs, nbv)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Milk, :Milk_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end

    begin                       # the SCS trait
        dt, dv = mk_scs_dt(dts, dcs, dbv)
        gt, gv = mk_scs_dt(gts, gcs, gbv)
        nt, nv = mk_scs_dt(nts, ncs, nbv)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :SCS, :SCS_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
end

function mk_milk_dt_2(ts, cs, bv)
    df1 = select(innerjoin(cs, ts, on = :ID), :ix, :Milk, :Milk_ERC)
    dft = filter(row -> row.Milk < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, :Milk)
    return dft, dfv
end

function mk_scs_dt_2(ts, cs, bv)
    df1 = select(innerjoin(cs, ts, on = :ID), :ix, :SCS, :SCS_ERC)
    dft = filter(row -> row.SCS < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, :SCS)
    return dft, dfv
end


text:
include("utils.jl")

@testset "SymmetricMatrices" begin
    M = SymmetricMatrices(3, ℝ)
    A = [1 2 3; 4 5 6; 7 8 9]
    A_sym = [1 2 3; 2 5 -1; 3 -1 9]
    A_sym2 = [1 2 3; 2 5 -1; 3 -1 9]
    B_sym = [1 2 3; 2 5 1; 3 1 -1]
    M_complex = SymmetricMatrices(3, ℂ)
    @test repr(M_complex) == "SymmetricMatrices(3, ℂ)"
    @test Manifolds.allocation_promotion_function(M_complex, get_vector, ()) === complex
    C = [1 1 -im; 1 2 -im; im im -1]
    D = [1 0; 0 1]
    X = zeros(3, 3)
    @testset "Real Symmetric Matrices Basics" begin
        @test repr(M) == "SymmetricMatrices(3, ℝ)"
        @test representation_size(M) == (3, 3)
        @test base_manifold(M) === M
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3,3},ℝ}
        @test check_point(M, B_sym) === nothing
        @test_throws DomainError is_point(M, A, true)
        @test_throws DomainError is_point(M, C, true)
        @test_throws DomainError is_point(M, D, true)
        @test check_vector(M, B_sym, B_sym) === nothing
        @test_throws DomainError is_vector(M, B_sym, A, true)
        @test_throws DomainError is_vector(M, A, B_sym, true)
        @test_throws DomainError is_vector(M, B_sym, D, true)
        @test_throws DomainError is_vector(M, B_sym, 1 * im * zero_vector(M, B_sym), true)
        @test manifold_dimension(M) == 6
        @test manifold_dimension(M_complex) == 9
        @test A_sym2 == project!(M, A_sym, A_sym)
        @test A_sym2 == project(M, A_sym, A_sym)
        A_sym3 = similar(A_sym)
        embed!(M, A_sym3, A_sym)
        A_sym4 = embed(M, A_sym)
        @test A_sym3 == A_sym
        @test A_sym4 == A_sym
    end
    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{3,3,Float64,9})

    bases = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    for T in types
        pts = [convert(T, A_sym), convert(T, B_sym), convert(T, X)]
        @testset "Type $T" begin
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_reverse_diff=isa(T, Vector),
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                basis_types_vecs=(
                    DiagonalizingOrthonormalBasis(log(M, pts[1], pts[2])),
                    bases...,
                ),
                basis_types_to_from=bases,
                is_tangent_atol_multiplier=1,
            )
            test_manifold(
                M_complex,
                pts,
                test_injectivity_radius=false,
                test_reverse_diff=isa(T, Vector),
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                basis_types_vecs=(DefaultOrthonormalBasis(ℂ),),
                basis_types_to_from=(DefaultOrthonormalBasis(ℂ),),
                is_tangent_atol_multiplier=1,
            )
            @test isapprox(-pts[1], exp(M, pts[1], log(M, pts[1], -pts[1])))
        end # testset type $T
    end # for
end # test SymmetricMatrices


text:
mvs = []
rectangle = [Ost, Nord, West, Sud]
cardinal_directions = [Nord, West, Sud, Ost]

function inverse(Side::HorizonSide)
    return HorizonSide((mod(Int(Side)+2, 4)))
end

function MoveIfPossible(r::Robot, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
    end
end

function snake_dif!(r::Robot)
    while !isborder(r, Nord) || !isborder(r, Ost)
        while !isborder(r, Ost)
            move!(r, Ost)
            if isborder(r, Nord)
                break
            end
        end
        if !isborder(r, Nord)
            move!(r, Nord)
        else
            break
        end
        while !isborder(r, West)
            move!(r, West)
            if isborder(r, Nord)
                break
            end
        end
        if !isborder(r, Nord)
            move!(r, Nord)
        else
            break
        end
    end
end

function go_to_corner!(r::Robot, mvs::AbstractArray)
    while isborder(r, Sud) == false || isborder(r, West) == false
        moveIfPossible!(r, mvs, Sud)
        moveIfPossible!(r, mvs, West)
    end
end

function moveIfPossible!(r::Robot, mvs::AbstractArray, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
        push!(mvs, side)
    end
end

function return_back!(r::Robot, mvs::AbstractArray)
    while length(mvs) > 0
        last_move = pop!(mvs)
        move!(r, inverse(last_move))
    end
end

function go_to_corner_Sud_West!(r::Robot)
    while !isborder(r, Sud) || !isborder(r, West)
        while !isborder(r, Sud)
            move!(r, Sud)
        end
        while !isborder(r, West)
            move!(r, West)
        end
    end
end

function main(r::Robot)
    go_to_corner!(r, mvs)
    snake_dif!(r)
    while isborder(r, Nord)
        move!(r, West)
    end
    move!(r, Ost)
    for l in [1, 2, 3, 4]
        while isborder(r, cardinal_directions[l])
            putmarker!(r)
            move!(r, rectangle[l])
        end
        move!(r, cardinal_directions[l])
    end
    go_to_corner_Sud_West!(r)
    return_back!(r, mvs)
end


text:
# Converts from cartesian coordinates to generalized
# Jacobi Keplerian coordinates, and then converts to orbital
# elements:

using DelimitedFiles
using LinearAlgebra

if ~@isdefined YEAR
  const YEAR  = 365.242
  const GNEWT = 39.4845/YEAR^2
end

include("../../../src/NbodyGradient/src/init_nbody.jl")

function amatrix(mass::Array{T,1},n_body) where {T <: Real}
# Set up "A" matrix (Hamers & Portegies-Zwart 2016) which transforms from
# cartesian coordinates to Keplerian orbits (we are using opposite sign
# convention of HPZ16, for example, r_1 = R_2-R_1).
# Mass vector is "mass".  Returns A matrix, "amat":
amat = zeros(T,n_body,n_body)
# Sum of masses for each Keplerian:
mkep = zeros(T,n_body-1)
# Get the indices:
indices = get_indices_planetary(n_body)
# Fill in the A matrix & compute the Keplerian elements:
for i=1:n_body-1
  # Sums of masses for two components of Keplerian:
  m1 = zero(T)
  m2 = zero(T)
  for j=1:n_body
    if indices[i,j] == 1
      m1 += mass[j]
    end
    if indices[i,j] == -1
      m2 += mass[j]
    end
  end
  mkep[i] = m1+m2
  # Compute Kepler problem: r is a vector of positions of "body" 2 with respect to "body" 1; rdot is velocity vector
  # For now set inclination to Inclination = pi/2 and longitude of nodes to Omega = pi:
  # Fill in the A matrix
  for j=1:n_body
    if indices[i,j] == 1
      amat[i,j] = -mass[j]/m1
    end
    if indices[i,j] == -1
      amat[i,j] =  mass[j]/m2
    end
  end
end
mtot = sum(mass)
# Final row is for center-of-mass of system:
for j=1:n_body
  amat[n_body,j] = mass[j]/mtot
end
return amat,mkep
end

function cartesian_to_elements(t::T,M::T,r::Array{T,1},rdot::Array{T,1}) where {T <: Real}
# Takes 3-d vectors r and rdot, and converts
# these to [a,lambda,e,omega] (assuming y=0):
# Equations (2.126-2.130) in M&D:
r2 = dot(r,r)
rabs = sqrt(r2)
v2 = dot(rdot,rdot)
rrdot = dot(r,rdot)
h = cross(r,rdot)
h2 = dot(h,h)
habs = sqrt(h2)
dotr = sign(rrdot)*sqrt(abs(v2-h2/r2))
# Now, compute elements:
semi = 1/(2/rabs-v2/(GNEWT*M))
ecc2 = 1-h2/(GNEWT*M*semi)
ecc  = sqrt(ecc2)
sinompf = -r[3]/rabs
#cosompf = -r[1]/rabs  # Do I want to assum \Omega = pi? [ ]
cosompf = r[1]/rabs  # Do I want to assum \Omega = pi? [ ]
ompf = atan(sinompf,cosompf)
sinf = semi*(1-ecc2)*dotr/(habs*ecc)
cosf = (semi*(1-ecc2)/rabs-1)/ecc
omega = atan(sinompf*cosf-cosompf*sinf,cosompf*cosf+sinompf*sinf)
# Eccentric anomaly (2.41): (does this have ambiguity? [ ])
#ean  = acos((1-rabs/semi)/ecc)
# Eccentric anomaly (2.46):
ean = 2 * atan(sqrt(1-ecc)*sinf,sqrt(1+ecc)*(1+cosf))
# Mean anomaly:
man  = ean - ecc*sin(ean)
# Now, we can compute lambda:
lam = man + omega
return [semi,lam,ecc,omega]
end

# Function which takes file name, reads in data, and computes
# the orbital elements:
function compute_elements(fname::String,masses::Array{T,1},n::Int64) where {T <: Real}
data = readdlm(fname)
# Fill out position and velocity arrays:
t = data[:,1]
nt = length(t)
x = reshape(data[:,2:n*3+1],nt,3,n)
v = reshape(data[:,n*3+2:end],nt,3,n)
# First, convert from Cartesian coordinates to
# Keplerian coordinates:
rkepler    = zeros(T,nt,n,3)
rdotkepler = zeros(T,nt,n,3)
# Construct the A matrix:
amat,mkep = amatrix(masses,n)
# Loop over time, and convert each set of cartesian coordinates
# to Keplerian coordinates:
for it=1:nt
  for i=1:n,j=1:3,k=1:n
    rkepler[it,i,j] += amat[i,k]*x[it,j,k]
    rdotkepler[it,i,j] += amat[i,k]*v[it,j,k]
  end
end
# Okay, now for conversion of Keplerian coordinates to
# orbital elements.  We will assume plane-parallel,
# and will compute [a,\lambda,e,\omega] at each time step:
elements = zeros(T,nt,4,n)
for it=1:nt
  # Compute orbital elements:
  for i=1:n-1
    elements[it,:,i] = cartesian_to_elements(t[it],mkep[i],rkepler[it,i,:],rdotkepler[it,i,:])
    if it > 1
      while elements[it,2,i] < (-pi/2 + elements[it-1,2,i])
        elements[it,2,i] += 2pi
      end
      while elements[it,2,i] > ( pi/2 + elements[it-1,2,i])
        elements[it,2,i] -= 2pi
      end
    end
  end
end
# Return array of elements as a function of time for each Keplerian
# in the Jacobi hierarchy:
return t,elements
end


text:
using Overseer

abstract type TComp <: ComponentData end

@component struct T1 <: TComp end
@component struct T2 end
@component struct T3 <: TComp end
@component struct T4 end

struct TSys <: System end

function Overseer.prepare(::TSys, m::AbstractLedger)
    if isempty(entities(m))
        Entity(m, T4())
    end
end

Overseer.requested_components(::TSys) = (T1, T2, T3, T4)

function Overseer.update(::TSys, m::AbstractLedger)
    t1 = m[T1]
    t2 = m[T2]
    t3 = m[T3]
    t4 = m[T4]
    for e in @entities_in(t1 && t2 && t3)
        t4[e] = T4()
    end
end

m = Ledger(Stage(:default, [TSys()]))

Entity(m, T1(), T2())
Entity(m, T2(), T3())
for i = 1:10
    Entity(m, T1(), T2(), T3())
end

@test length(valid_entities(m)) == 13

@test length(m[T1]) == 11

update(m)

@test length(m[T4]) == 11

empty!(m[T4])
update(stage(m, :default), m)
@test length(m[T4]) == 10

@test length(m[Entity(4)]) == 4

delete!(m, Entity(4))
@test !isempty(Overseer.free_entities(m))
@test length(filter(x->x==Entity(0), m.entities)) == 1

Entity(m, T3())
@test m[T3][Entity(4)] == T3()
@test isempty(Overseer.free_entities(m))

@test length(m[T4]) == 9

for i = 5:10
    schedule_delete!(m, Entity(i))
end
delete_scheduled!(m)

@test length(m[T4]) == 3

# Ensure delete_scheduled!() does nothing when no entities are scheduled for
# deletion, but some of the old ids have been recycled
for _ = 5:10
    Entity(m, T4())
end
delete_scheduled!(m)
@test length(m[T4]) == 9

empty!(m)
@test isempty(m.entities)
@test isempty(m.stages)
@test isempty(m.components)

push!(m, Stage(:default, [TSys()]))

@test length(m.components) == 8

struct TSys2 <: System end

push!(m, :default, TSys())

@test length(last(stage(m, :default))) == 2

insert!(m, :default, 1, TSys2())

@test last(stage(m, :default))[1] == TSys2()

insert!(m, 1, Stage(:test, [TSys(), TSys2()]))

@test first(m.stages[1]) == :test

@test eltype(m[T4]) == T4

prepare(m)
@test !isempty(entities(m))
@test singleton(m, T4) == T4()

struct SmallSys <: System end

Overseer.requested_components(::SmallSys) = (T1, T3)

m2 = Ledger(Stage(:default, [SmallSys()]))

@test m2.components[2] === Overseer.EMPTY_COMPONENT

e = Entity(m2)
m2[e] = T2()

empty_entities!(m2)
@test isempty(m2.entities)
@test !isempty(m2.components)

@test length(components(m2, TComp)) == 2

empty!(m)
for i=1:10
    e1 = Entity(m,
                Test1(length(entities(m))),
                Test2(length(entities(m))),
                Test3(length(entities(m))))

    e1 = Entity(m,
                Test1(length(entities(m))),
                Test2(length(entities(m))))

    e2 = Entity(m,
                Test2(length(entities(m))),
                )
end

t2 = group(m, Test1, Test2)
t1 = group(m, Test1, Test2, Test3)

@test_throws ArgumentError group(m, Test1, Test3)

@test length(groups(m)) == 1
@test t1.child == t2
@test t1.child.parent == t1 == t2.parent

@test length(t1) == 10
@test length(t2) == 20

@test m[Test1].indices.packed[1:length(t1)] == m[Test2].indices.packed[1:length(t1)] == m[Test3].indices.packed[1:length(t1)] == t1.indices.packed[1:length(t1)]
@test m[Test1].indices.packed[1:length(t2)] == m[Test2].indices.packed[1:length(t2)]  == t2.indices.packed[1:length(t2)]

test_entity = last(entities(m))
t2[test_entity] = Test1()
@test length(t1) == 10
t1[test_entity] = Test3()
@test length(t1) == 11

@test t1[test_entity] == (m[Test1][test_entity], m[Test2][test_entity], m[Test3][test_entity])

tot = 0
for e in @entities_in(group(m, Test1, Test2))
    global tot += 1
end
@test tot == length(group(m, Test1, Test2))

test_entity = entities(m)[end-1]
m[test_entity] = Test3(test_entity.id)

@test m[Test1].indices.packed[1:length(t1)] == m[Test2].indices.packed[1:length(t1)] == m[Test3].indices.packed[1:length(t1)] == t1.indices.packed[1:length(t1)]
@test m[Test1].indices.packed[1:length(t2)] == m[Test2].indices.packed[1:length(t2)]  == t2.indices.packed[1:length(t2)]

@test t1[test_entity] == (m[Test1][test_entity], m[Test2][test_entity], m[Test3][test_entity])

# tot = 0
# for e in @entities_in(group(m, Test2, Test3))
#     global tot += 1
# end
# @test tot == 8 == length(group(m, Test2, Test3))

# tot = 0
# for e in @entities_in(group(m, Test2, Test3) && group(m, Test1, Test2))
#     global tot += 1
# end

# tot2 = 0
# for e in @entities_in(m[Test1] && m[Test2] && m[Test3])
#     global tot2 += 1
# end
# @test tot == tot2

# remove_group!(m, Test2, Test3)
# @test length(groups(m)) == 2

# tg = create_group!(m, Test1, Test2; ordered=true)
# @test length(groups(m)) == 2

# @test groups(m)[1] isa Overseer.OrderedGroup

# tg = group(m, Test1, Test2, Test3)
# beforelen = length(tg)
# m[Test3].shared[m[Test3].data[beforelen+1]] != Test3(5)

# m[Entity(1)] = Test3(5)
# @test length(tg) == beforelen+1
# @test m[Test3].shared[m[Test3].data[length(tg)]] == Test3(5)
# @test m[Test2].data[length(tg)] == Test2(0)


text:
struct ZGroup{S<:AbstractStore}
    storage::S
    path::String
    arrays::Dict{String, ZArray}
    groups::Dict{String, ZGroup}
    attrs::Dict
    writeable::Bool
end

# path can also be a SubString{String}
ZGroup(storage, path::AbstractString, arrays, groups, attrs, writeable) =
    ZGroup(storage, String(path), arrays, groups, attrs, writeable)

zname(g::ZGroup) = zname(g.path)

#Open an existing ZGroup
function ZGroup(s::T,mode="r",path="") where T <: AbstractStore
  arrays = Dict{String, ZArray}()
  groups = Dict{String, ZGroup}()

  for d in subdirs(s,path)
    dshort = split(d,'/')[end]
    m = zopen_noerr(s,mode,path=_concatpath(path,dshort))
    if isa(m, ZArray)
      arrays[dshort] = m
    elseif isa(m, ZGroup)
      groups[dshort] = m
    end
  end
  attrs = getattrs(s,path)
  startswith(path,"/") && error("Paths should never start with a leading '/'")
  ZGroup(s, path, arrays, groups, attrs,mode=="w")
end

"""
    zopen_noerr(AbstractStore, mode = "r"; consolidated = false)

Works like `zopen` with the single difference that no error is thrown when
the path or store does not point to a valid zarr array or group, but nothing
is returned instead.
"""
function zopen_noerr(s::AbstractStore, mode="r"; consolidated = false, path="", lru = 0)
    consolidated && isinitialized(s,".zmetadata") && return zopen(ConsolidatedStore(s, path), mode, path=path)
    lru !== 0 && return zopen(LRUStore(s,maxsize=lru),mode,path=path)
    if is_zarray(s, path)
        return ZArray(s,mode,path)
    elseif is_zgroup(s,path)
        return ZGroup(s,mode,path)
    else
        return nothing
    end
end

function Base.show(io::IO, g::ZGroup)
    print(io, "ZarrGroup at ", g.storage, " and path ", g.path)
    !isempty(g.arrays) && print(io, "\nVariables: ", map(i -> string(zname(i), " "), values(g.arrays))...)
    !isempty(g.groups) && print(io, "\nGroups: ", map(i -> string(zname(i), " "), values(g.groups))...)
    nothing
end
Base.haskey(g::ZGroup,k)= haskey(g.groups,k) || haskey(g.arrays,k)

function Base.getindex(g::ZGroup, k)
    if haskey(g.groups, k)
        return g.groups[k]
    elseif haskey(g.arrays, k)
       return g.arrays[k]
    else
       throw(KeyError("Zarr Dataset does not contain $k"))
    end
end

"""
    zopen(s::AbstractStore, mode="r"; consolidated = false, path = "", lru = 0)

Opens a zarr Array or Group at Store `s`. If `consolidated` is set to "true",
Zarr will search for a consolidated metadata field as created by the python zarr
`consolidate_metadata` function. This can substantially speed up metadata parsing
of large zarr groups. Setting `lru` to a value > 0 means that chunks that have been
accessed before will be cached and consecutive reads will happen from the cache.
Here, `lru` denotes the number of chunks that remain in memory.
"""
function zopen(s::AbstractStore, mode="r"; consolidated = false, path = "", lru = 0)
    # add interfaces to Stores later
    r = zopen_noerr(s,mode; consolidated=consolidated, path=path, lru=lru)
    if r === nothing
        throw(ArgumentError("Specified store $s in path $(path) is neither a ZArray nor a ZGroup"))
    else
        return r
    end
end

"""
    zopen(p::String, mode="r")

Open a zarr Array or group at disc path p.
"""
function zopen(s::String, mode="r"; kwargs...)
  store, path = storefromstring(s,false)
  zopen(store, mode; path=path, kwargs...)
end

function storefromstring(s, create=true)
  for (r,t) in storageregexlist
    if match(r,s) !== nothing
      return storefromstring(t,s,create)
    end
  end
  if create || isdir(s)
    return DirectoryStore(s), ""
  else
    throw(ArgumentError("Path $s is not a directory."))
  end
end

"""
    zgroup(s::AbstractStore; attrs=Dict())

Create a new zgroup in the store `s`
"""
function zgroup(s::AbstractStore, path::String=""; attrs=Dict())
    d = Dict("zarr_format"=>2)
    isemptysub(s, path) || error("Store is not empty")
    b = IOBuffer()
    JSON.print(b,d)
    s[path,".zgroup"]=take!(b)
    writeattrs(s,path,attrs)
    ZGroup(s, path, Dict{String,ZArray}(), Dict{String,ZGroup}(), attrs,true)
end

zgroup(s::String;kwargs...)=zgroup(storefromstring(s, true)...;kwargs...)

"Create a subgroup of the group g"
function zgroup(g::ZGroup, name; attrs=Dict())
  g.writeable || throw(IOError("Zarr group is not writeable. Please re-open in write mode to create an array"))
  g.groups[name] = zgroup(g.storage,_concatpath(g.path,name),attrs=attrs)
end

"Create a new subarray of the group g"
function zcreate(::Type{T},g::ZGroup, name::String, addargs...; kwargs...) where T
  g.writeable || throw(IOError("Zarr group is not writeable. Please re-open in write mode to create an array"))

  z = zcreate(T, g.storage, addargs...; path = _concatpath(g.path,name), kwargs...)
  g.arrays[name] = z
  return z
end

HTTP.serve(s::Union{ZArray,ZGroup}, args...; kwargs...) = HTTP.serve(s.storage, s.path, args...; kwargs...)
function consolidate_metadata(z::Union{ZArray,ZGroup})
  z.writeable || throw(Base.IOError("Zarr group is not writeable. Please re-open in write mode to create an array",0))
  consolidate_metadata(z.storage,z.path)
end


text:
"""
A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension
"""
function unit_indices(N::Int) #create unit CartesianIndex for each dimension
    out = Vector{CartesianIndex{N}}(undef, N)
    null = zeros(Int64, N)
    for i in 1:N
        unit_i = copy(null)
        unit_i[i] = 1
        out[i] = CartesianIndex(Tuple(unit_i))
    end
    Tuple(out)
end

add_dims(A::AbstractArray, n::Int) = cat(ndims(A) + n, A)

""
perpindex(A, dim::Integer) = A[setdiff(1:length(A), dim)]

"""
the size of A perpendicular to dim
"""
perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)]


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

using Base.Dates, Requires

module Timestamps

using Base.Dates, Requires

import Base: getindex, show, convert, sum, prod, mean, var, std, maximum, minimum, merge

export Timestamp, TimestampArray

include("timestamp.jl")
include("conversion.jl")
include("operators.jl")
include("arraymethods.jl")
include("utils.jl")

end