text:
function MarkPerimetrInsideRectangle(r)
    counter = 0
    B = Array{Int, 1}(undef, 1000)
    B, counter = movetothecorner(r, West, B)
    findrectangle(r, Nord)
    marksideofrectangle(r, Nord)
    marksideofrectangle(r, Ost)
    marksideofrectangle(r, Sud)
    marksideofrectangle(r, West)
    marksideofrectangle(r, Nord)
    movements(r, West)
    movements(r, Sud)
    returnfromthecorner(r, Ost, B, counter)
end

function movements(r, side)
    while !isborder(r, side)
        move!(r, side)
    end
end

function movementsandcheckforborder(r, side)
    while !isborder(r, side)
        move!(r, side)
        if isborder(r, Nord)
            break
        end
    end
end

function marksideofrectangle(r, side)
    while isborder(r, side)
        move!(r, HorizonSide((Int(side) + 1) % 4))
        putmarker!(r)
    end
    move!(r, side)
    putmarker!(r)
end

function findrectangle(r, side)
    while !isborder(r, side)
        movementsandcheckforborder(r, HorizonSide((Int(side) + 3) % 4))
        if isborder(r, side)
            break
        end
        move!(r, side)
        movementsandcheckforborder(r, HorizonSide((Int(side) + 1) % 4))
        if isborder(r, side)
            break
        end
        if !isborder(r, side)
            move!(r, side)
        else
            break
        end
    end
end

function movetothecorner(r, side, B)
    counter = 1
    while ((!isborder(r, side)) || (!isborder(r, HorizonSide((Int(side) + 1) % 4))))
        if !isborder(r, side)
            move!(r, side)
            B[counter] = 1
            counter += 1
        end
        if !isborder(r, HorizonSide((Int(side) + 1) % 4))
            move!(r, HorizonSide((Int(side) + 1) % 4))
            B[counter] = 2
            counter += 1
        end
    end
    return B, counter
end

function returnfromthecorner(r, side, B, counter)
    for i in 1:(counter - 1)
        if B[counter - i] == 1
            move!(r, side)
        end
        if B[counter - i] == 2
            move!(r, HorizonSide((Int(side) + 1) % 4))
        end
    end
end


text:
export PoissonSparsePCG,solve!
using IterativeSolvers
using Preconditionners
using IncompleteLU

struct PoissonSparsePCG{Preco}
    spLxy::SparseMatrixCSC{Float64,Int64}
    p::Preco
    iluLxy::Any

    function PoissonSparsePCG(sp)
        tilu=@elapsed lu = ilu(sp, τ = 0.2)
        println("ILU factorization time:",tilu)
        new(sp,lu)
    end
end

function solve!(Xxy,Bxy,sp::PoissonSparsePCG)
    p = sp.iluLxy
    b=view(Bxy,1:length(Bxy))
    x=view(Xxy,1:length(Xxy))
    a,cgilulog=IterativeSolvers.cg!(x,sp.spLxy,b,Pl = p;log=true,tol = 1e-10,maxiter=1000)
    @show cgilulog
    cgilulog
end

ext: jl


text:
# TODO: Add tests that include conditions and iters etc. Make sure to test keywords for all types.

using AgentBased
using Base.Test
using HDF5

@testset "AgentBased.jl" begin

#==================================================================================================
TEST REPORTERS
==================================================================================================#

    floatreturner(x::Float64) = x

    r1 = Reporter("int", Int64, () -> 5)
    r2 = Reporter("float", Float64, floatreturner)

    @testset "Reporter call type stability" begin
        check_call(reporter::Reporter, args...) = reporter.call(args...)

        @test Base.return_types(check_call, Tuple{typeof(r1)}) == [Int64]
        @test Base.return_types(check_call, Tuple{typeof(r2), Float64}) == [Float64]
    end

#==================================================================================================
TEST COLLECTORS
==================================================================================================#

    c1 = Collector(r1, r2)

    @testset "Collector calls type stability" begin
        check_call_1(coll::Collector, args...) = coll.calls[1](args...)
        check_call_2(coll::Collector, args...) = coll.calls[2](args...)
        check_call_both(coll::Collector, args...) = (coll.calls[1](), coll.calls[2](args...))

        @test Base.return_types(check_call_1, Tuple{typeof(c1)}) == [Int64]
        @test Base.return_types(check_call_2, Tuple{typeof(c1), Float64}) == [Float64]
        @test Base.return_types(check_call_both, Tuple{typeof(c1), Float64}) ==
            [Tuple{Int64,Float64}]
    end

#==================================================================================================
TEST TYPEDBUFFERS
==================================================================================================#

    t1 = AgentBased.TypedBuffer(c1, 3)

    @testset "TypedBuffer internals" begin
        @test t1.size == 0
        @test t1.maxsize == 3
        @test length(t1.data) == 2
        @test all(length(t1.data[i]) == 3 for i in length(t1.data))
    end

    @testset "TypedBuffer data type stability" begin
        function check_data(tb::AgentBased.TypedBuffer)
            foo = tb.data[1][1]
            bar = tb.data[2][1]
            return foo, bar
        end

        @test Base.return_types(check_data, Tuple{typeof(t1)}) == [Tuple{Int64,Float64}]
        @test Base.return_types(AgentBased.TypedBuffer, Tuple{typeof(c1), Int}) ==
            [AgentBased.TypedBuffer{Tuple{Vector{Int64},Vector{Float64}}}]
    end

#==================================================================================================
TEST HDF5WRITER
==================================================================================================#

    h1 = HDF5Writer("test.h5", "test", c1)
    h1 = HDF5Writer("test.h5", "test", c1; chunksz = 3, channels_per_worker = 1, writemode = "w")

    @testset "HDF5Writer internals" begin
        coll_type = typeof(h1).parameters[1]
        buff_type = typeof(h1).parameters[2]
        @test all(buff_type.parameters[i].parameters[1] == coll_type.parameters[1].parameters[i]
                for i in 1:length(buff_type.parameters))
    end

    @testset "HDF5Writer data type stability" begin
        function check_data(writer::HDF5Writer)
            foo = writer.localbuffer.data[1][1]
            bar = writer.localbuffer.data[2][1]
            return foo, bar
        end

        @test Base.return_types(check_data, Tuple{typeof(h1)}) == [Tuple{Int64,Float64}]
    end

    @testset "HDF5Writer calls type stability" begin
        check_call_1(writer::HDF5Writer, args...) = writer.collector.calls[1](args...)
        check_call_2(writer::HDF5Writer, args...) = writer.collector.calls[2](args...)
        check_call_both(writer::HDF5Writer, args...) = (writer.collector.calls[1](),
                                                        writer.collector.calls[2](args...))

        @test Base.return_types(check_call_1, Tuple{typeof(h1)}) == [Int64]
        @test Base.return_types(check_call_2, Tuple{typeof(h1), Float64}) == [Float64]
        @test Base.return_types(check_call_both, Tuple{typeof(h1), Float64}) ==
            [Tuple{Int64,Float64}]
    end

#==================================================================================================
TEST COLLECTDATA, ENSUREROOM, AND FLUSHDATA
==================================================================================================#

    @testset "HDF5Writer data collection" begin
        test_iter() = 1:4
        test_r1 = Reporter("iter", Int64, (i) -> i)
        test_r2 = Reporter("float", Float64, (i) -> convert(Float64, i))
        test_c = Collector(test_r1, test_r2; iter = test_iter)
        test_w = HDF5Writer("test.h5", "test", test_c; chunksz = 3, channels_per_worker = 1)

        @sync begin
            @async begin
                collectdata(test_w)
                AgentBased.flushdata(test_w)
            end
            @async begin
                take!(test_w.remotequeue) ==
                    AgentBased.TypedBuffer{Tuple{Vector{Int64},Vector{Float64}}}(
                        3, 3, ([1, 2, 3], [1.0, 2.0, 3.0]))
                take!(test_w.remotequeue) ==
                    AgentBased.TypedBuffer{Tuple{Vector{Int64},Vector{Float64}}}(
                        1, 3, ([4, 2, 3], [4.0, 2.0, 3.0]))
            end
        end
    end

#==================================================================================================
TEST RUNBATCH
==================================================================================================#

#==================================================================================================
TEST INSERTHDF5
==================================================================================================#

end


text:
using MDCDL
using LinearAlgebra
using Images, TestImages
using Plots

using Base.Filesystem
using Dates

function main()
    ########## Configurations #########
    # choose NSOLT type: (Rnsolt | Cnsolt)
    Nsolt = Rnsolt
    # TP := eltype(Nsolt) (<:AbstractFloat)
    TP = Float64
    # decimation factor: (<:NTuple{D,Int} where D is #dims)
    df = (2,2)
    # polyphase order: (<:NTuple{D,Int} where D)
    ord = (2,2)
    # number of channels: (<:Union{Integer,Tuple{Int,Int}} for Rnsolt)
    #                     (<:Integer for Cnsolt)
    nch = 8
    # number of tree level (<: Integer)
    level = 3

    # size of minibatches (<:NTuple{D,Int})
    szx = (16,16)
    # number of minibatches (<:Integer)
    nSubData = 16

    # path of log files (do nothing if isa(logdir, Nothing))
    # logdir = joinpath(@__DIR__, "results", Dates.format(now(), "yyyymmdd_HH_MM_SS_sss"))
    logdir = nothing
    # save minibatches?
    do_save_trainingset = false
    do_export_atoms = false

    # options for sparse coding
    sparsity = 0.25
    sparsecoder = SparseCoders.IHT
    sparsecoder_options = (
        iterations = 400,
        nonzeros = trunc(Int, sparsity * prod(szx)), filter_domain=:convolution,)
    # options for dictionary update
    optimizer = Optimizers.Steepest
    optimizer_options = (
        rate = 1e-4,
        iterations=10,)

    # general options of dictionary learning
    options = ( epochs  = 1000,
                verbose = :standard, # :none, :standard, :specified, :loquacious
                sparsecoder = sparsecoder,
                sparsecoder_options = sparsecoder_options,
                optimizer = optimizer,
                optimizer_options = optimizer_options,
                logdir = logdir,)
    ####################################
    logdir != nothing && !isdir(logdir) && mkpath(logdir)

    # original image
    orgImg = testimage("cameraman")

    # generate minibatches of orgImage
    trainingSet = map(1:nSubData) do nsd
        pos = rand.(UnitRange.(0,size(orgImg) .- szx))
        orgImg[UnitRange.(1 .+ pos, szx .+ pos)...]
    end

    if logdir != nothing && do_save_trainingset
        datadir = joinpath(logdir, "data")
        !isdir(datadir) && mkpath(datadir)
        map(idx->save(joinpath(datadir, "$idx.png"), trainingSet[idx]), 1:nSubData)
    end

    # create NSOLT instance
    nsolt = Nsolt(TP, df, ord, nch)
    # set random orthonormal matrices to the initial matrices.
    MDCDL.rand!(nsolt, isPropMat = false, isPropAng = false, isSymmetry = false)
    istype1(nsolt) && MDCDL.vm1constraint!(nsolt)

    msnsolt = Multiscale([ deepcopy(nsolt) for l in 1:level ]...)
    # rand!(msnsolt.filterbanks[level], isSymmetry = false, isPropAng=false, isPropMat = false)

    # dictionary learning
    MDCDL.train!(msnsolt, trainingSet; options...)

    if logdir != nothing && do_export_atoms
        for idx = 1:length(msnsolt)
            png(plot(msnsolt.filterbanks[idx]), joinpath(logdir, string("nsolt_",idx,".png")))
        end
    end

    return msnsolt
end

main()

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using Tables, Formatting

include("../general/examples.jl")

for tblname in tblnames

  fname = "./examples/$(string(tblname)).txt"
  f = open(fname, "w")
  println("opened: $fname")
  rendercells(f, tocellarray(eval(tblname)))
  close(f)
  println("closed: $fname")
end


text:
# main.jl

# (c) 2020 Johannes Boehm <johannes.boehm@sciencespo.fr>

using Revise
using DataFrames

using Random
using Distances

using Test

include("src/JacobiAlgorithm.jl")

# Locations
rng = MersenneTwister(777)
nbz = 3
z_df = DataFrame(h = rand(rng, nbz), v = rand(rng, nbz))

# Supply
nbi = 200
rng = MersenneTwister(778)
i_df = DataFrame(h = rand(rng, nbi),
    v = rand(rng, nbi),
    τ = rand(rng, nbi),
    λ = 10.0 * rand(rng, nbi)
    )

# Demand
nbj = 500
rng = MersenneTwister(779)
j_df = DataFrame(h = rand(rng, nbj),
    v = rand(rng, nbj),
    σ = 1.0 ./ rand(rng, nbj),
    ϵ = 20.0 .* rand(rng, nbj),
    η = rand(rng, nbj)
    )

avg_speed_drive = 25.0
avg_speed_walk = 4.0

T_iz = pairwise(Euclidean(), Matrix(i_df[:,[:h,:v]])', Matrix(z_df[:,[:h,:v]])' ) ./ avg_speed_drive
T_jz = pairwise(Euclidean(), Matrix(j_df[:,[:h,:v]])', Matrix(z_df[:,[:h,:v]])' ) ./ avg_speed_walk

function utility_firms(p::R, τ::R, λ::R, T::R) where R<:Real
    return p^(one(R) - τ) - λ * T
end
function utility_consumers(p::R, σ::R, ϵ::R, T::R) where R<:Real
    expo = one(R) - one(R) / σ
    return ( p^expo + (ϵ*T)^expo )^(one(R)/expo)
end

# it would be more "Julian" to have a separate type for firms

function demand_z(prices_z::Vector{R}, consumers_j::DataFrame, T::Matrix{R}) where R<:Real
    number_of_markets = size(prices_z, 1)
    number_of_consumers = size(consumers_j, 1)
    d = zeros(R, number_of_markets)
    # construct matrix of utilities
    u_jz = zeros(R, number_of_consumers, number_of_markets)
    for z = 1:number_of_markets
        for j = 1:number_of_consumers
            u_jz[j,z] = utility_consumers(prices_z[z],consumers_j[j,:σ],consumers_j[j,:ϵ], T[j,z])
        end
    end
    # for each firm i, increment the supply vector at the position where it has the highest value
    for j=1:number_of_consumers
        min_index = -1
        min = consumers_j[j,:η]
        for z = 1:number_of_markets
            if u_jz[j,z] < min
                min = u_jz[j,z]
                min_index = z
            end
        end
        if min_index != -1
            d[min_index] += 1
        else
            # outside option is better
        end
    end
    return d
end
function supply_z(prices_z::Vector{R}, firms_i::DataFrame, T::Matrix{R}) where R<:Real
    number_of_firms = size(firms_i, 1)
    number_of_markets = size(prices_z, 1)
    s = zeros(R, number_of_markets)
    # construct matrix of utilities
    u_iz = zeros(R, number_of_firms, number_of_markets)
    for z = 1:number_of_markets
        for i = 1:number_of_firms
            u_iz[i,z] = utility_firms(prices_z[z],firms_i[i,:τ],firms_i[i,:λ],T[i,z])
        end
    end
    # for each firm i, increment the supply vector at the position where it has the highest value
    for i=1:number_of_firms
        max_index = -1
        max = zero(R)
        for z = 1:number_of_markets
            if u_iz[i,z] > max
                max = u_iz[i,z]
                max_index = z
            end
        end
        s[max_index] += 1
    end
    return s
end

function ssmooth_z(prices_z::Vector{R}, firms_i::DataFrame, T::Matrix{R}, smoothpar::R) where R<:Real
    number_of_firms = size(firms_i, 1)
    number_of_markets = size(prices_z, 1)
    s = zeros(R, number_of_markets)
    # construct matrix of utilities
    u_iz = zeros(R, number_of_firms, number_of_markets)
    for z = 1:number_of_markets
        for i = 1:number_of_firms
            u_iz[i,z] = utility_firms(prices_z[z],firms_i[i,:τ],firms_i[i,:λ],T[i,z])
        end
    end
    # calculate smooth max
    u_i0 = zeros(R, number_of_firms)
    max_i = [max(zero(R), maximum([u_iz[i,z] for z=1:number_of_markets])) for i=1:number_of_firms]
    utilde_i = max_i .- smoothpar.* log.( exp.( (u_i0 .- max_i)./smoothpar ) .+ sum( exp.( (u_iz .- repeat(max_i,1,number_of_markets))./smoothpar ), dims=2 ) )
    return sum( exp.( ( u_iz .- repeat(utilde_i,1,number_of_markets)) ./ smoothpar) , dims=1)
end
function dsmooth_z(prices_z::Vector{R}, consumers_j::DataFrame, T::Matrix{R}, smoothpar::R) where R<:Real
    # here I denote the cost by u_jz
    number_of_markets = size(prices_z, 1)
    number_of_consumers = size(consumers_j, 1)
    d = zeros(R, number_of_markets)
    # construct matrix of utilities
    u_jz = zeros(R, number_of_consumers, number_of_markets)
    for z = 1:number_of_markets
        for j = 1:number_of_consumers
            u_jz[j,z] = utility_consumers(prices_z[z],consumers_j[j,:σ],consumers_j[j,:ϵ], T[j,z])
        end
    end
    # calculate smooth max
    u_j0 = [consumers_j[j,:η] for j=1:number_of_consumers]
    min_j = [min(u_j0[j], minimum([u_jz[j,z] for z=1:number_of_markets])) for j=1:number_of_consumers]
    utilde_j = min_j .- smoothpar.* log.( exp.( (-u_j0 .+ min_j)./smoothpar ) .+ sum( exp.( (-u_jz .+ repeat(min_j,1,number_of_markets))./smoothpar ), dims=2 ) )
    return sum( exp.( ( -u_jz .+ repeat(utilde_j,1,number_of_markets)) ./ smoothpar) , dims=1)
end
# excess supply function, smoothed
function esmooth_z(prices_z::Vector{R}, firms_i::DataFrame, T_firms::Matrix{R}, consumers_j::DataFrame, T_consumers::Matrix{R}, smoothpar::R) where R<:Real
    return ssmooth_z(prices_z, firms_i, T_firms, smoothpar) .- dsmooth_z(prices_z, consumers_j, T_consumers, smoothpar)
end

# show the validity of the approximations

# exact
supply_z(ones(Float64,3), i_df, T_iz)
# approximation
ssmooth_z(ones(Float64,3), i_df, T_iz, 0.0000001)
# test it
@test maximum(abs.(supply_z(ones(Float64,3), i_df, T_iz) .- ssmooth_z(ones(Float64,3), i_df, T_iz, 0.0000001))) ≈ 0.0 atol= 1e04

# exact
demand_z(ones(Float64,3).*0.1, j_df, T_jz)
# approximation
dsmooth_z(ones(Float64,3).*0.1, j_df, T_jz, 0.0000001)
# test it
@test maximum(abs.(demand_z(ones(Float64,3).*0.1, j_df, T_jz) .- dsmooth_z(ones(Float64,3).*0.1, j_df, T_jz, 0.0000001))) ≈ 0.0 atol= 1e04

include("src/JacobiAlgorithm.jl")

# define the problem
e = JacobiAlgorithm.EquilibriumProblem(nbz, p->esmooth_z(p,i_df, T_iz, j_df, T_jz, 1e-7), 0.0, 1e5, zeros(Float64,nbz), zeros(Float64,nbz) )

# try out excess supply fct
e.e_z(zeros(Float64,3)) .- e.q_z'

# try out coordinate update
p_z = ones(Float64, 3).*0.1
JacobiAlgorithm.coordinate_update(e,1,p_z)
JacobiAlgorithm.coordinate_update(e,2,p_z)
JacobiAlgorithm.coordinate_update(e,3,p_z)

# try out the first iterations of Jacobi algorithm
p_z = ones(Float64, 3).*0.1
for i = 1:5
    global p_z
    p_z , _ = JacobiAlgorithm.fJ_z(e,p_z)
    @show p_z
end

# solve the problem using the jacobi algorithm: about 4.2 seconds
JacobiAlgorithm.solve(e, steptol=1e-12, valtol=1e-5, output=2)

# threaded version: < 2 seconds
JacobiAlgorithm.solve_threaded(e, steptol=1e-12, valtol=1e-5, output=2)


text:
@testset "dimension reducing methods" begin
    a = [1 2; 3 4]
    dimz = (X((143, 145)), Y((-38, -36)))
    da = DimensionalArray(a, dimz)

    @test sum(da; dims=X) == sum(da; dims=1)
    @test sum(da; dims=Y()) == sum(da; dims=2)
    @test dims(sum(da; dims=Y())) == (X(LinRange(143.0, 145.0, 2)), Y(LinRange(-38.0, -38.0, 1)))
    @test prod(da; dims=X) == [3 8]
    @test prod(da; dims=Y()) == [2 12]'
    @test dims(prod(da; dims=X())) == (X(LinRange(143.0, 143.0, 1)), Y(LinRange(-38.0, -36.0, 2)))
    @test maximum(da; dims=X) == [3 4]
    @test maximum(da; dims=Y()) == [2 4]'
    @test minimum(da; dims=X) == [1 2]
    @test minimum(da; dims=Y()) == [1 3]'
    @test dims(minimum(da; dims=X())) == (X(LinRange(143.0, 143.0, 1)), Y(LinRange(-38.0, -36.0, 2)))
    @test mean(da; dims=X) == [2.0 3.0]
    @test mean(da; dims=Y()) == [1.5 3.5]'
    @test dims(mean(da; dims=Y())) == (X(LinRange(143.0, 145.0, 2)), Y(LinRange(-38.0, -38.0, 1)))
    @test reduce(+, da; dims=X) == [4 6]
    @test reduce(+, da; dims=Y()) == [3 7]'
    @test dims(reduce(+, da; dims=Y())) == (X(LinRange(143.0, 145.0, 2)), Y(LinRange(-38.0, -38.0, 1)))
    @test mapreduce(x-> x > 3, +, da; dims=X) == [0 1]
    @test mapreduce(x-> x > 3, +, da; dims=Y()) == [0 1]'
    @test dims(mapreduce(x-> x > 3, +, da; dims=Y())) == (X(LinRange(143.0, 145.0, 2)), Y(LinRange(-38.0, -38.0, 1)))
    @test std(da; dims=X()) == [1.4142135623730951 1.4142135623730951]
    @test std(da; dims=Y()) == [0.7071067811865476 0.7071067811865476]'
    @test var(da; dims=X()) == [2.0 2.0]
    @test var(da; dims=Y()) == [0.5 0.5]'
    @test dims(var(da; dims=Y())) == (X(LinRange(143.0, 145.0, 2)), Y(LinRange(-38.0, -38.0, 1)))
    a = [1 2 3; 4 5 6]
    da = DimensionalArray(a, dimz)
    @test median(da; dims=Y()) == [2.0 5.0]'
    @test median(da; dims=X()) == [2.5 3.5 4.5]

    @test dropdims(da[X(1:1)]; dims=X) == [1, 2, 3]
    @test dropdims(da[2:2, 1:1]; dims=(X(), Y()))[] == 4
    @test typeof(dropdims(da[2:2, 1:1]; dims=(X(), Y()))) <: DimensionalArray{Int,0,Tuple{}}
    # TODO: test refdims after dropdims
end

if VERSION > v"1.1-"
    @testset "iteration methods" begin
        a = [1 2 3 4
             3 4 5 6
             5 6 7 8]
        # eachslice
        da = DimensionalArray(a, (Y(10:30), Time(1:4)))
        @test [mean(s) for s in eachslice(da; dims=Time)] == [3.0, 4.0, 5.0, 6.0]
        slices = [s .* 2 for s in eachslice(da; dims=Y)]
        @test slices[1] == [2, 4, 6, 8]
        @test slices[2] == [6, 8, 10, 12]
        @test slices[3] == [10, 12, 14, 16]
        dims(slices[1]) == (Time(1.0:1.0:4.0),)
        slices = [s .* 2 for s in eachslice(da; dims=Time)]
        @test slices[1] == [2, 6, 10]
        dims(slices[1]) == (Y(10.0:10.0:30.0),)
    end
end

@testset "dimension reordering methods" begin
    da = DimensionalArray(zeros(5, 4), (Y(10:20), X(1:4)))
    tda = transpose(da)
    @test dims(tda) == (X(LinRange(1.0, 4.0, 4)), Y(LinRange(10.0, 20.0, 5)))
    @test size(tda) == (4, 5)
    ada = adjoint(da)
    @test dims(ada) == (X(LinRange(1.0, 4.0, 4)), Y(LinRange(10.0, 20.0, 5)))
    @test size(ada) == (4, 5)

    dsp = permutedims(da)
    @test parent(dsp) == permutedims(parent(da))
    @test dims(dsp) == reverse(dims(da))
    da = DimensionalArray(ones(5, 2, 4), (Y(10:20), Time(10:11), X(1:4)))
    dsp = permutedims(da, [3, 1, 2])

    @test permutedims(da, [X, Y, Time]) == permutedims(da, (X, Y, Time))
    @test permutedims(da, [X(), Y(), Time()]) == permutedims(da, (X(), Y(), Time()))
    dsp = permutedims(da, (X(), Y(), Time()))
    @test dims(dsp) == (X(LinRange(1.0, 4.0, 4)), Y(LinRange(10.0, 20.0, 5)), Time(LinRange(10.0, 11.0, 2)))
end

@testset "dimension mirroring methods" begin
    # Need to think about dims for these, currently (Y, Y) etc.
    # But you can't index (Y, Y) with dims as you get the
    # first Y both times. It will plot correctly at least.
    a = rand(5, 4)
    da = DimensionalArray(a, (Y(10:20), X(1:4)))

    cvda = cov(da; dims=X)
    @test cvda == cov(a; dims=2)
    @test dims(cvda) == (X(LinRange(1.0, 4.0, 4)), X(LinRange(1.0, 4.0, 4)))
    crda = cor(da; dims=Y)
    @test crda == cor(a; dims=1)
    @test dims(crda) == (Y(LinRange(10.0, 20.0, 5)), Y(LinRange(10.0, 20.0, 5)))
end

@testset "mapslices" begin
    a = [1 2 3 4
         3 4 5 6
         5 6 7 8]
    da = DimensionalArray(a, (Y(10:30), Time(1:4)))
    ms = mapslices(sum, da; dims=Y)
    @test ms == [9 12 15 18]
    @test dims(ms) == (Time(LinRange(1.0, 4.0, 4)),)
    @test refdims(ms) == (Y(10.0),)
    ms = mapslices(sum, da; dims=Time)
    @test parent(ms) == [10 18 26]'
    @test dims(ms) == (Y(LinRange(10.0, 30.0, 3)),)
    @test refdims(ms) == (Time(1.0),)
end

# These need fixes in base. kwargs are ::Integer so we can't add methods
# or dispatch on AbstractDimension in underscore _methods
accumulate
cumsum
cumprod


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    𝓒(x)

Return closure of `x`.
"""
𝓒(x) = x ./ sum(x)

"""
    Composition(partscomps)
    Composition(parts, comps)
    Composition(part₁=comp₁, part₂=part₂, ...)
    Composition(comps)
    Composition(comp₁, comp₂, ...)

A D-part composition as defined by Aitchison 1986.

## Examples

A 2-part composition with parts `a=0.1` and `b=0.8`:

```julia
julia> Composition(a=0.2, b=0.8)
julia> Composition((a=0.2, b=0.8))
julia> Composition((:a, :b), (0.2, 0.8))
```

When the names of the parts are not specified, the
constructor uses default names `part1`, `part2`,
..., `partD`:

```
julia> Composition(0.1, 0.8)
julia> Composition((0.1, 0.8))
```
"""
struct Composition{D,PARTS}
  data::NamedTuple{PARTS,NTuple{D,Union{Float64,Missing}}}
end

Composition(data::NamedTuple) =
  Composition{length(data),keys(data)}(data)

Composition(; data...) = Composition((; data...))

Composition(parts::NTuple, comps) =
  Composition((; zip(parts, Tuple(comps))...))

Composition(comps) =
  Composition(ntuple(i->Symbol("w$i"), length(comps)), comps)

Composition(comp::Real, comps...) = Composition((comp, comps...))

"""
    parts(c)

Parts in the composition `c`.
"""
parts(::Composition{D,PARTS}) where {D,PARTS} = PARTS

"""
    components(c)

Components in the composition `c`.
"""
components(c::Composition) = getfield(c, :data) |> values |> SVector

Base.getproperty(c::Composition, p::Symbol) = getproperty(getfield(c, :data), p)

# -------------
# VECTOR SPACE
# -------------

+(c₁::Composition{D,PARTS}, c₂::Composition{D,PARTS}) where {D,PARTS} =
  Composition(PARTS, 𝓒(components(c₁) .* components(c₂)))

-(c::Composition) = Composition(parts(c), 𝓒(1 ./ components(c)))

-(c₁::Composition, c₂::Composition) = c₁ + -c₂

*(λ::Real, c::Composition) = Composition(parts(c), 𝓒(components(c).^λ))

==(c₁::Composition, c₂::Composition) =
  parts(c₁) == parts(c₂) && 𝓒(components(c₁)) ≈ 𝓒(components(c₂))

⋅(c₁::Composition{D}, c₂::Composition{D}) where {D} = begin
  x, y = components(c₁), components(c₂)
  sum(log(x[i]/x[j])*log(y[i]/y[j]) for j=1:D for i=j+1:D) / D
end

norm(c::Composition) = √(c⋅c)

distance(c₁::Composition, c₂::Composition) = norm(c₁ - c₂)

# --------------------
# RANDOM COMPOSITIONS
# --------------------

"""
    rand(Composition{D}, n=1)

Generates `D`-part composition at random according to a
balanced Dirichlet distribution.
"""
function rand(rng::Random.AbstractRNG,
              ::Random.SamplerType{<:Composition{D}}) where {D}
  α = fill(1.0, D)
  d = Dirichlet(α)
  x = rand(rng, d)
  Composition(x)
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, c::Composition)
  w = [(@sprintf "%.03f" w) for w in components(c)]
  show(io, join(w, " : "))
end

function Base.show(io::IO, mime::MIME"text/plain",
                   c::Composition{D,PARTS}) where {D,PARTS}
  w = components(c)
  x = Vector{Float64}()
  p = Vector{Symbol}()
  m = Vector{Symbol}()
  for i in 1:D
    if ismissing(w[i])
      push!(m, PARTS[i])
    else
      push!(p, PARTS[i])
      push!(x, w[i])
    end
  end
  plt = barplot(p, x, title="$D-part composition")
  isempty(m) || annotate!(plt, :t, "missing: $(join(m,", "))")
  show(io, mime, plt)
end