text:
using Distributed
addprocs(16)

using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using Random
using LinearAlgebra
using GLM
using CSV

# function for converting nongenetic covariates
function make_nongenetic_covariates()
	z = CSV.read("Covariate_polrGWAS.csv")
	z = [z[:, 3:6] z[:, 8] z[:, 12:end]]
	rename!(z, :x1 => :bmi)
	CSV.write("Covariate_Final.csv", z)
end

# runs a small example so cross validation code compiles
function run_iht_once_to_force_compile(d, l)
	n = 1000
	p = 10000
	k = 10

	#set random seed
	Random.seed!(2019)

	#construct x matrix and non genetic covariate (intercept)
	T = Float64
	x = simulate_random_snparray(n, p, undef)
	xbm = SnpBitMatrix{T}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(T, n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	true_b = zeros(T, p)
	true_b[1:k] .= collect(0.1:0.1:1.0)
	true_c = [T.(4.0)]
	shuffle!(true_b)
	correct_position = findall(!iszero, true_b)

	#simulate phenotypes (e.g. vector y)
	if d == Normal || d == Bernoulli || d == Poisson
	    prob = linkinv.(l, xbm * true_b)
	    clamp!(prob, -20, 20)
	    y = [rand(d(i)) for i in prob]
	    # prob = linkinv.(l, x * true_b + z * true_c)
	    # clamp!(prob, -20, 20)
	    # y = [rand(d(i)) for i in prob]
	    # k = k + 1
	elseif d == NegativeBinomial
	    nn = 10
	    μ = linkinv.(l, xbm * true_b)
	    # μ = linkinv.(l, x * true_b + z * true_c)
	    # k = k + 1
	    clamp!(μ, -20, 20)
	    prob = 1 ./ (1 .+ μ ./ nn)
	    y = [rand(d(nn, Float64(i))) for i in prob] #number of failtures before nn success occurs
	elseif d == Gamma
	    μ = linkinv.(l, xbm * true_b)
	    β = 1 ./ μ # here β is the rate parameter for gamma distribution
	    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
	end
	y .= T.(y)

	#specify path and folds
	path = collect(1:10)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run IHT
	cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, debias=false, parallel=true)
end

function run()
	d = Normal
	l = canonicallink(d())
	run_iht_once_to_force_compile(d, l)

	#Import covariates and phenotype. (data is in $SCRATCH)
	filedir = "/u/scratch/b/biona001/"
	x = SnpArray(filedir * "ukb.plink.filtered.bed")
	z = CSV.read(filedir * "Covariate_Final.csv")
	z = convert(Matrix{Float64}, z)
	standardize!(z)
	z = [ones(size(z, 1)) z]
	y = CSV.read(filedir * "Covariate_polrGWAS.csv")[:AveSBP] # average systolic blood pressure
	y = convert(Vector{Float64}, y)

	#set random seed
	Random.seed!(2019)

	#specify path and folds
	path = collect(1:5)
	num_folds = 3
	folds = rand(1:num_folds, size(x, 1))

	# run cross validation
	mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, debias=false, parallel=true)
	writedlm("mses_$(first(path))_to_$(last(path))", mses)
end


text:
using CosmoParticles
using FillArrays
using LazyArrays
using LinearAlgebra
using Rotations
using SortingAlgorithms
using StatsBase
using Test
using Unitful

const CP = CosmoParticles

@testset "CosmoParticles.jl" begin

    @testset "Particles" begin
        p = Particles(:dm)

        p.id = sample(1:1000, 100; replace=false)
        p.pos = rand(3, 100)

        # different id and pos for copied Particles
        pcid = rand(100)
        pcpos = rand(3, 100)

        # Particle specific
        pc = Particles(:dm, :id => p.id, :pos => p.pos)
        @test pc == p

        pc = copy(p)
        @test pc == p
        @test pc !== p
        @test pc.props == p.props
        @test pc.props !== p.props

        pc = empty(p)
        @test pc.type === :dm
        @test isempty(pc.props)
        @test isempty(pc)

        @test CP.particle_name(p) == "Particles"
        @test issetequal(propertynames(p), [:type, :id, :pos])
        @test issetequal(propertynames(p; private=true), [:type, :props, :id, :pos])

        io = IOBuffer()
        show(io, "text/plain", p)
        @test String(take!(io)) == "dm: 100 Particles\n id pos"

        # test AbstractParticles with basic implementation Particles
        @test CP.get_props(p) === p.props

        @test getproperty(p, :id) === p.props[:id]
        @test getproperty(p, :type) === p.type === :dm
        @test getproperty(p, :props) === p.props
        @test p.id === p.props[:id]

        pc = copy(p)
        setproperty!(p, :pos, pcpos)
        @test p.pos === pcpos
        pc.id = pcid
        @test pc.id === pcid
        @test p.id !== pcid

        @test_throws ErrorException("""
            setfield!: immutable struct of type Particles cannot be changed""") p.type = :gas

        @test p[:id] === p.props[:id]

        pc = copy(p)
        p[:id] = pcid
        @test p.id === pcid

        @test keys(p) == keys(p.props)
        @test haskey(p, :id)
        @test !haskey(p, :vel)
        @test !haskey(p, :type)

        @test values(p) == values(p.props)

        @test issetequal(propertynames(p), [:type, :id, :pos])

        # test default implementations for AbstractParticles
        @test issetequal(Base.@invoke(propertynames(p::AbstractParticles)), [:id, :pos])

        pc = deepcopy(p)
        empty!(pc)
        @test isempty(pc.props)
        @test isempty(pc)

        @test !isempty(p)

        copy!(pc, p)
        @test pc !== p && pc == p && !isempty(pc)

        @test CP.particle_number(p) == 100
        @test CP.particle_number(Particles(:gas)) == 0
        @test CP.particle_number(Particles(:gas, Dict{Symbol,Any}(:mass => 3))) == 0

        @test isnothing(Base.@invoke(CP.particle_name(p::AbstractParticles)))

        io = IOBuffer()
        CP.show_properties(io, "text/plain", p)
        @test String(take!(io)) == "100 Particles\n id pos"
    end

    @testset "AllParticles" begin
        dm = Particles(:dm)
        dm.id = sample(1:1000, 100; replace=false)
        dm.pos = rand(3, 100)
        dm.mass = rand()
        dm.test = rand()

        gas = Particles(:gas)
        gas.id = sample(1001:2000, 100; replace=false)
        gas.pos = rand(3, 100)
        gas.mass = rand(100)
        gas.temp = rand(100)
        gas.zs = rand(2, 100)
        gas.test = rand(2, 100)

        pcoll = ParticleCollection(dm, gas)
        p = AllParticles(pcoll)

        @test CP.particle_collection(p) === pcoll
        # @test get_props(p) === p.props # not part of the API

        if p.id[1] == dm.id[1]
            @test p.id == vcat(dm.id, gas.id)
            @test p.pos == hcat(dm.pos, gas.pos)
            @test p.mass == vcat(fill(dm.mass, length(dm.id)), gas.mass)
            @test isequal.(p.temp, vcat(fill(missing, length(dm.id)), gas.temp)) |> all
            @test isequal.(p.zs, hcat(fill(missing, 2, length(dm.id)), gas.zs)) |> all
            @test p.test == hcat(fill(dm.test, 2, length(dm.id)), gas.test)
        else
            @test p.id == vcat(gas.id, dm.id)
            @test p.pos == hcat(gas.pos, dm.pos)
            @test p.mass == vcat(gas.mass, fill(dm.mass, length(dm.id)))
            @test isequal.(p.temp, vcat(gas.temp, fill(missing, length(dm.id)))) |> all
            @test isequal.(p.zs, hcat(gas.zs, fill(missing, 2, length(dm.id)))) |> all
            @test p.test == hcat(gas.test, fill(dm.test, 2, length(dm.id)))
        end

        @test_throws ErrorException p.id = rand(100)
        @test_throws KeyError p.vel

        @test issetequal(keys(p), [:id, :mass, :pos, :temp, :zs, :test])
        @test keys(p) == keys(AllParticles(pcoll))

        @test haskey(p, :id) && haskey(AllParticles(pcoll), :id)
        @test !haskey(p, :vel) && !haskey(AllParticles(pcoll), :vel)

        @test all(val === p[key] for (key, val) in zip(keys(p), values(p)))
        @test all(val === p[key] for (key, val) in pairs(p))

        @test_throws ErrorException copy(p)
        @test_throws ErrorException copy!(p, AllParticles(pcoll))
        @test_throws ErrorException copy!(p, dm)

        pc = Particles(:all)
        copy!(pc, p)
        collect(values(p)) # materializes properties in Dict
        @test pc.id == p.id
        @test pc.mass == p.mass
        @test pc.pos == p.pos
        @test all(isequal.(pc.temp, p.temp))
        @test all(isequal.(pc.zs, p.zs))
        @test pc.test == p.test

        copy!(pc, p, [:id, :mass, :temp])
        @test haskey.((pc,), [:id, :mass, :temp]) |> all
        @test !haskey(pc, :pos)
        @test isa.([pc.id, pc.mass, pc.temp], Array) |> all

        @test Particles(p, [:id, :mass, :temp]) == pc

        @test_throws ErrorException empty(p)
        @test_throws ErrorException empty!(p)

        @test !isempty(p)
        @test isempty(ParticleCollection() |> AllParticles)
        @test isempty(ParticleCollection(Particles(:dm), Particles(:gas)) |> AllParticles)

        @test p == AllParticles(pcoll)

        @test CP.particle_name(p) == "Particles"
        @test CP.particle_number(p) == 200

        list = [:id, :mass, :pos, :temp, :zs, :test]
        @test issetequal(propertynames(p), list)
        @test issetequal(propertynames(p; private=true), [fieldnames(AllParticles) |> collect; list])

        io = IOBuffer()
        show(io, "text/plain", p)
        @test String(take!(io)) == "all: 200 Particles\n id mass pos temp test zs"

        @test CP.particle_number(AllParticles(ParticleCollection())) == 0

        mask = isodd.(p.id)
        ind = findall(mask)
        @test_throws ErrorException p[mask]
        @test_throws ErrorException p[ind]

        # cover case where the dimension is not found from a vector or matrix
        # and where only one particle subtype exists
        pos = rand(3, 100)
        mass = rand()
        pc = ParticleCollection(Particles(:dm, :pos => pos, :mass => mass))
        ap = pc.all
        @test ap.pos == pos
        @test ap.mass == fill(mass, 100)
    end

    @testset "Particle Collection" begin
        dm = Particles(:dm)
        dm.id = sample(1:1000, 100; replace=false)
        dm.pos = rand(3, 100)

        gas = Particles(:gas)
        gas.id = sample(1001:2000, 100; replace=false)
        gas.pos = rand(3, 100)

        @testset "ParticleCollection" begin
            pc = ParticleCollection()
            @test pc isa ParticleCollection{Particles}
            @test pc.particles == Dict{Symbol,Particles}()

            @test ParticleCollection(Particles) == pc

            pc = ParticleCollection(:dm => dm, :gas => gas)
            @test pc isa ParticleCollection{Particles}
            @test pc.particles == Dict{Symbol,Particles}(:dm => dm, :gas => gas)

            @test ParticleCollection(dm, gas) == pc

            @test redshift(pc) == 0

            pc = ParticleCollection(dm)
            pcc = copy(pc)
            @test pcc.particles[:dm] === pc.particles[:dm]
            pcc.particles[:dm] = gas
            @test pcc.particles[:dm] !== pc.particles[:dm]

            @test empty(pc) == ParticleCollection()

            pc = ParticleCollection(dm, gas)
            io = IOBuffer()
            show(io, "text/plain", pc)
            @test String(take!(io)) ==
                  "ParticleCollection\ndm: 100 Particles\n id pos\ngas: 100 Particles\n id pos"
        end

        @testset "RedshiftParticleCollection" begin
            z = 0.5

            pc = RedshiftParticleCollection(0.5)
            @test pc isa RedshiftParticleCollection{Particles}
            @test pc.particles == Dict{Symbol,Particles}()
            @test pc.z == z

            @test RedshiftParticleCollection(Particles, z) == pc

            pc = RedshiftParticleCollection(z, :dm => dm, :gas => gas)
            @test pc isa RedshiftParticleCollection{Particles}
            @test pc.particles == Dict{Symbol,Particles}(:dm => dm, :gas => gas)

            @test RedshiftParticleCollection(z, dm, gas) == pc

            @test redshift(pc) == z

            pc = RedshiftParticleCollection(z, dm)
            pcc = copy(pc)
            @test pcc.particles[:dm] === pc.particles[:dm]
            pcc.particles[:dm] = gas
            @test pcc.particles[:dm] !== pc.particles[:dm]

            @test empty(pc) == RedshiftParticleCollection(z)

            pc = RedshiftParticleCollection(z, dm, gas)
            io = IOBuffer()
            show(io, "text/plain", pc)
            @test String(take!(io)) ==
                  "ParticleCollection at z = $z\ndm: 100 Particles\n id pos\ngas: 100 Particles\n id pos"

            @test issetequal(propertynames(pc), [:z, :all, :dm, :gas])
            @test issetequal(propertynames(pc; private=true), [:z, :all, :particles, :dm, :gas])

            rpc = RedshiftParticleCollection(0, dm, gas)
            pc = ParticleCollection(dm, gas)
            @test rpc == pc && pc == rpc
            rpc = RedshiftParticleCollection(z, dm, gas)
            @test rpc != pc && pc != rpc
        end

        @testset "AbstractParticleCollection" begin
            pc = ParticleCollection(dm, gas)

            @test CP.get_particles(pc) === pc.particles

            @test getproperty(pc, :dm) === pc.particles[:dm]
            @test getproperty(pc, :particles) === pc.particles
            @test pc.dm === pc.particles[:dm]
            @test pc[:dm] === pc.particles[:dm]

            @test_throws ErrorException(
                "setfield!: immutable struct of type ParticleCollection cannot be changed",
            ) pc.particles = Dict{Symbol,Particles}()
            pc.dm = gas
            @test pc[:dm] === gas

            pc[:dm] = dm
            @test pc.dm === dm

            @test pc.all == AllParticles(pc)
            @test_throws ErrorException pc.all = dm
            @test_throws KeyError pc[:all]

            @test keys(pc) == keys(pc.particles)
            @test haskey(pc, :dm)
            @test !haskey(pc, :stars)
            @test !haskey(pc, :particles)

            @test values(pc) == values(pc.particles)

            @test issetequal(propertynames(pc), [:all, :dm, :gas])

            pcc = deepcopy(pc)
            empty!(pcc)
            @test isempty(pcc.particles)
            @test isempty(pcc)

            @test !isempty(pc)

            copy!(pcc, pc)
            @test pcc !== pc && pcc == pc && !isempty(pcc)

            io = IOBuffer()
            CP.show_properties(io, "text/plain", pc)
            @test String(take!(io)) == "dm: 100 Particles\n id pos\ngas: 100 Particles\n id pos"
        end
    end

    @testset "Utils" begin
        @testset "Index particles" begin
            p = Particles(:dm)
            p.id = sample(1:1000, 100; replace=false)
            p.pos = rand(3, 100)
            p.mass = 4u"kg"

            mask = rand(100) .> 0.5
            ind = findall(mask)

            @test CP._applyind(3, mask) == 3
            @test CP._applyind(p.mass, mask) == p.mass
            @test CP._applyind(p.id, mask) == p.id[mask]
            @test CP._applyind(p.pos, mask) == p.pos[:, mask]
            @test CP._applyind(3, ind) == 3
            @test CP._applyind(p.mass, ind) == p.mass
            @test CP._applyind(p.id, ind) == p.id[ind]
            @test CP._applyind(p.pos, ind) == p.pos[:, ind]

            a = rand(45)
            b = rand(55)
            c = ApplyArray(vcat, a, b)
            @test CP._applyind(c, mask) == Array(c)[mask]
            @test CP._applyind(c, ind) == Array(c)[ind]

            a = rand(3, 45)
            b = rand(3, 55)
            c = ApplyArray(hcat, a, b)
            @test CP._applyind(c, mask) == Array(c)[:, mask]
            @test CP._applyind(c, ind) == Array(c)[:, ind]

            pc = CP.applyind(p, mask)
            @test pc.id == p.id[mask]
            @test pc.pos == p.pos[:, mask]
            @test pc.mass == p.mass

            @test CP.applyind(p, ind) == pc
            @test p[ind] == pc
            @test p[mask] == pc

            pcc = deepcopy(p)
            CP.applyind!(pcc, mask)
            @test pcc == pc

            pcc = deepcopy(p)
            CP.applyind!(pcc, ind)
            @test pcc == pc

            pc = CP.applyind(p, mask; affect=[:pos, :mass, :foo])
            @test !haskey(pc, :id)
            @test !haskey(pc, :foo)
            @test pc.pos == p.pos[:, mask]
            @test pc.mass == p.mass

            pc = ParticleCollection(p)
            ap = pc.all
            @test_throws ErrorException CP.applyind(ap, ind)
            @test_throws ErrorException CP.applyind!(ap, ind)
            @test_throws ErrorException CP.applyind(ap, mask)
            @test_throws ErrorException CP.applyind!(ap, mask)
        end

        @testset "Remove particles" begin
            p = Particles(:dm)
            p.id = sample(1:1000, 100; replace=false)
            p.pos = rand(3, 100)
            p.mass = 4u"kg"

            mask = rand(100) .> 0.5
            ind = findall(mask)

            @test CP._removeind(3, ind) == 3
            @test CP._removeind(p.mass, ind) == p.mass
            @test CP._removeind(p.id, ind) == deleteat!(copy(p.id), ind)
            posremoved = CP._removeind(p.pos, ind)
            @test posremoved[1, :] == deleteat!(p.pos[1, :], ind)
            @test posremoved[2, :] == deleteat!(p.pos[2, :], ind)
            @test posremoved[3, :] == deleteat!(p.pos[3, :], ind)

            a = rand(45)
            b = rand(55)
            c = ApplyArray(vcat, a, b)
            @test CP._removeind(c, ind) == deleteat!(Array(c), ind)

            a = rand(3, 45)
            b = rand(3, 55)
            c = ApplyArray(hcat, a, b)
            cremoved = CP._removeind(c, ind)
            carray = Array(c)
            @test cremoved[1, :] == deleteat!(carray[1, :], ind)
            @test cremoved[2, :] == deleteat!(carray[2, :], ind)
            @test cremoved[3, :] == deleteat!(carray[3, :], ind)

            pc = CP.applyind(p, .~mask)

            @test CP.removeind(p, ind) == pc

            pcc = deepcopy(p)
            CP.removeind!(pcc, ind)
            @test pcc == pc

            affect = [:pos, :mass, :foo]
            pc = CP.removeind(p, ind; affect)
            @test pc == CP.applyind(p, .~mask; affect)

            @test CP.deleteat(p, ind) == CP.removeind(p, ind)
            pc = deepcopy(p)
            @test CP.deleteat!(pc, ind) == CP.removeind(p, ind)

            pc = ParticleCollection(p)
            ap = pc.all
            @test_throws Exception CP.removeind(ap, ind)
            @test_throws Exception CP.removeind!(ap, ind)
        end

        @testset "Find all in" begin
            a = sample(1:10000, 6000; replace=false)
            set = sample(1:10000, 100; replace=false)
            inboth = intersect(a, set)

            ind = CP.findall_in(a, set)
            @test length(ind) == length(inboth)
            @test issetequal(a[ind], inboth)
            @test CP.findall_in(a, Set(set)) == ind

            sort!(a)
            sort!(set)
            ind = CP.findall_in(a, set)
            @test length(ind) == length(inboth)
            @test issetequal(a[ind], inboth)
            @test CP.findall_in(a, Set(set)) == ind
            @test CP.findall_in_sorted(a, set) == ind

            alarge = collect(1:1_200_000)
            alarge2 = collect(200_000:1_600_000)
            setlarge = sample(1:1_200_000, 100_000; replace=false) |> sort!
            inboth = intersect(alarge, setlarge)
            inboth2 = intersect(alarge2, setlarge)

            ind = CP.findall_in(alarge, setlarge)
            @test length(ind) == length(inboth)
            @test issetequal(alarge[ind], inboth)
            @test CP.findall_in(alarge, Set(setlarge)) == ind

            ind = CP.findall_in(alarge2, setlarge)
            @test length(ind) == length(inboth2)
            @test issetequal(alarge2[ind], inboth2)
            @test CP.findall_in(alarge2, Set(setlarge)) == ind

            # test empty arrays
            @test CP.findall_in_sorted([], []) |> isempty
            @test CP.findall_in_sorted(a, []) |> isempty
            @test CP.findall_in_sorted([], set) |> isempty

            @test CP.findall_in([], Set([])) |> isempty
            @test CP.findall_in(a, Set([])) |> isempty
            @test CP.findall_in([], Set(set)) |> isempty
        end

        @testset "Unitful operations" begin
            a = rand(100)
            au = a * u"m"
            af = convert(Vector{Float32}, a)
            afu = af * u"m"
            b = rand()

            @test CP.product_preserve_type(a, b) == a .* b
            @test CP.product_preserve_type(au, b) == au .* b
            @test CP.product_preserve_type(af, b) == af .* convert(Float32, b)
            @test CP.product_preserve_type(afu, b) == afu .* convert(Float32, b)

            ac = copy(a)
            CP.product_preserve_type!(ac, b)
            @test ac == a .* b

            auc = copy(au)
            CP.product_preserve_type!(auc, b)
            @test auc == au .* b

            afc = copy(af)
            CP.product_preserve_type!(afc, b)
            @test afc == af .* convert(Float32, b)

            afuc = copy(afu)
            CP.product_preserve_type!(afuc, b)
            @test afuc == afu .* convert(Float32, b)
        end

        @testset "Lazy ustrip" begin
            @test CP.ustrip_lazy(3) == 3

            a = rand(4)
            b = rand(3)
            au = a * u"m"
            bu = b * u"m"
            @test CP.ustrip_lazy(a) === a
            @test CP.ustrip_lazy(au) == a
            @test CP.ustrip_lazy(Diagonal(au)) == Diagonal(a)
            @test CP.ustrip_lazy(Bidiagonal(au, bu, :U)) == Bidiagonal(a, b, :U)
            @test CP.ustrip_lazy(Tridiagonal(bu, au, bu)) == Tridiagonal(b, a, b)
            @test CP.ustrip_lazy(SymTridiagonal(au, bu)) == SymTridiagonal(a, b)

            auview = @view au[1:3]
            @test CP.ustrip_lazy(auview) == a[1:3]

            c = vcat(a, b)
            cu = ApplyArray(vcat, au, bu)
            @test CP.ustrip_lazy(cu) == c

            c = fill(2.5, 4)
            cu = Fill(2.5u"m", 4)
            @test CP.ustrip_lazy(cu) == c

            auf = rand(Float32, 4) * u"m"
            au_strip = CP.ustrip_lazy(u"km", au)
            auf_strip = CP.ustrip_lazy(u"km", auf)
            @test au_strip ≈ ustrip(au ./ 1000)
            @test eltype(au_strip) === Float64
            @test auf_strip ≈ ustrip(auf ./ 1000)
            @test eltype(auf_strip) === Float32

            auc = copy(au)
            aufc = copy(auf)
            auc_strip = CP.ustrip_lazy!(u"km", auc)
            aufc_strip = CP.ustrip_lazy!(u"km", aufc)
            @test auc_strip == au_strip
            @test aufc_strip == auf_strip
            @test CP.ustrip_lazy(auc) == au_strip
            @test CP.ustrip_lazy(aufc) == auf_strip

            auc = copy(au)
            aufc = copy(auf)
            auc_conv = CP.uconvert_lazy!(u"km", auc)
            aufc_conv = CP.uconvert_lazy!(u"km", aufc)
            @test auc_conv ≈ uconvert.(u"km", au)
            @test aufc_conv ≈ uconvert.(u"km", auf)
            @test CP.ustrip_lazy(auc) == CP.ustrip_lazy(auc_conv)
            @test CP.ustrip_lazy(aufc) == CP.ustrip_lazy(aufc_conv)

            auc = copy(au)
            aufc = copy(auf)
            auc_conv = CP.uconvert_lazy!(u"m", auc)
            aufc_conv = CP.uconvert_lazy!(u"m", aufc)
            @test auc_conv ≈ uconvert.(u"m", au)
            @test aufc_conv ≈ uconvert.(u"m", auf)
            @test CP.ustrip_lazy(auc) == CP.ustrip_lazy(auc_conv)
            @test CP.ustrip_lazy(aufc) == CP.ustrip_lazy(aufc_conv)

            au = rand(4) * u"°C"
            auf = rand(Float32, 4) * u"°C"
            auc = copy(au)
            aufc = copy(auf)
            auc_conv = CP.uconvert_lazy!(u"K", auc)
            aufc_conv = CP.uconvert_lazy!(u"K", aufc)
            @test auc_conv ≈ uconvert.(u"K", au)
            @test aufc_conv ≈ uconvert.(u"K", auf)
            @test CP.ustrip_lazy(auc) == CP.ustrip_lazy(auc_conv)
            @test CP.ustrip_lazy(aufc) == CP.ustrip_lazy(aufc_conv)
        end
    end

    @testset "Transformations" begin
        @testset "Rotation" begin
            rotmat = rand(RotMatrix{3})
            a = rand(3, 100)
            au = a * u"m"

            arot = rotmat * a
            aurot = rotmat * au

            @test CP.matrix_rotate(a, rotmat) == arot
            @test CP.matrix_rotate(au, rotmat) == aurot

            ac = copy(a)
            CP.matrix_rotate!(ac, rotmat)
            @test ac == arot

            auc = copy(au)
            CP.matrix_rotate!(auc, rotmat)
            @test auc == aurot

            ac = similar(a)
            CP.matrix_rotate!(ac, a, rotmat)
            @test ac == arot

            auc = similar(au)
            CP.matrix_rotate!(auc, au, rotmat)
            @test auc == aurot

            @testset "Particle Rotation" begin
                p = Particles(:dm)
                p.pos = copy(au)
                p.vel = copy(a)

                pc = rotate(p, rotmat, :pos)
                @test p.pos == au
                @test pc.pos == aurot
                @test p.vel === pc.vel == a

                pc = rotate(p, rotmat)
                @test p.pos == au
                @test pc.pos == aurot
                @test p.vel == a
                @test pc.vel == arot

                pc = deepcopy(p)
                rotate!(pc, rotmat, :vel)
                @test pc.pos == au
                @test pc.vel == arot

                pc = deepcopy(p)
                rotate!(pc, rotmat)
                @test pc.pos == aurot
                @test pc.vel == arot
            end

            @testset "Particle Collection Rotation" begin
                dm = Particles(:dm)
                dm.pos = copy(au)
                dm.vel = copy(a)
                gas = Particles(:gas)
                gas.pos = copy(au)
                gas.vel = copy(a)

                pc = ParticleCollection(dm, gas)

                pcc = rotate(pc, rotmat, :pos)
                @test pc.dm.pos == pc.gas.pos == au
                @test pcc.dm.pos == pcc.gas.pos == aurot
                @test pc.dm.vel === pcc.dm.vel == a

                pcc = rotate(pc, rotmat)
                @test pc.dm.pos == pc.gas.pos == au
                @test pcc.dm.pos == pcc.gas.pos == aurot
                @test pc.dm.vel == pc.gas.vel == a
                @test pcc.dm.vel == pcc.gas.vel == arot

                pcc = deepcopy(pc)
                rotate!(pcc, rotmat, :vel)
                @test pcc.dm.pos == pcc.gas.pos == au
                @test pcc.dm.vel == pcc.gas.vel == arot

                pcc = deepcopy(pc)
                rotate!(pcc, rotmat)
                @test pcc.dm.pos == pcc.gas.pos == aurot
                @test pcc.dm.vel == pcc.gas.vel == arot

                pccc = deepcopy(pc)
                ap = pccc.all
                @test_throws ErrorException rotate(ap, rotmat)

                rotate!(ap, rotmat)
                @test pccc == pcc
            end
        end

        @testset "Translations" begin
            da = rand(3)
            dau = rand(3) * u"m"
            a = rand(3, 100)
            au = a * u"m"

            atrans = a .+ da
            autrans = au .+ dau

            @testset "Particle Translations" begin
                p = Particles(:dm)
                p.pos = copy(au)
                p.vel = copy(a)

                pc = translate(p, dau)
                @test p.pos == au
                @test pc.pos == autrans
                @test p.vel === pc.vel == a

                pc = translate(p, da, :vel)
                @test p.pos === pc.pos == au
                @test p.vel == a
                @test pc.vel == atrans

                pc = deepcopy(p)
                translate!(pc, dau)
                @test pc.pos == autrans
                @test pc.vel == a

                pc = deepcopy(p)
                translate!(pc, da, :vel)
                @test pc.pos == au
                @test pc.vel == atrans

                # periodic translations
                period, period_half = 10, 5
                CP._translate_periodic(1, 2, period) == 1
                CP._translate_periodic(8, 2, period, period_half) == -2
                CP._translate_periodic(2, 8, period, period_half) == 12

                o = oneunit(eltype(dau))

                pc = translate_periodic(p, dau, o)
                @test all(dau .- 0.5o .≤ minimum(pc.pos; dims=2)) &&
                      all(maximum(pc.pos; dims=2) .≤ dau .+ 0.5o)

                pc = deepcopy(p)
                translate_periodic!(pc, dau, o)
                @test all(dau .- 0.5o .≤ minimum(pc.pos; dims=2)) &&
                      all(maximum(pc.pos; dims=2) .≤ dau .+ 0.5o)

                pcc = translate_periodic_to_center(p, dau, o)
                @test pcc == translate(pc, -dau)

                pcc = deepcopy(p)
                translate_periodic_to_center!(pcc, dau, o)
                @test pcc == translate(pc, -dau)
            end

            @testset "Particle Collection Translations" begin
                dm = Particles(:dm)
                dm.pos = copy(au)
                dm.vel = copy(a)
                gas = Particles(:gas)
                gas.pos = copy(au)
                gas.vel = copy(a)

                pc = ParticleCollection(dm, gas)

                pcc = translate(pc, dau)
                @test pc.dm.pos == pc.gas.pos == au
                @test pcc.dm.pos == pcc.gas.pos == autrans
                @test pc.dm.vel === pcc.dm.vel == a

                pcc = translate(pc, da, :vel)
                @test pc.dm.pos === pcc.dm.pos == au
                @test pc.dm.vel == pc.gas.vel == a
                @test pcc.dm.vel == pcc.gas.vel == atrans

                pcc = deepcopy(pc)
                translate!(pcc, dau)
                @test pcc.dm.pos == pcc.gas.pos == autrans
                @test pcc.dm.vel == pcc.gas.vel == a

                pcc = deepcopy(pc)
                translate!(pcc, da, :vel)
                @test pcc.dm.pos == pcc.gas.pos == au
                @test pcc.dm.vel == pcc.dm.vel == atrans

                pccc = deepcopy(pc)
                ap = pccc.all
                @test_throws ErrorException translate(ap, da)

                translate!(ap, dau)
                @test pccc == translate(pc, dau)

                pccc = deepcopy(pc)
                ap = pccc.all
                translate!(ap, da, :vel)
                @test pccc == pcc

                o = oneunit(eltype(dau))

                pcc = translate_periodic(pc, dau, o)
                @test pcc.dm == translate_periodic(pc.dm, dau, o)
                @test pcc.gas == translate_periodic(pc.gas, dau, o)

                pccc = deepcopy(pc)
                translate_periodic!(pccc, dau, o)
                @test pccc == pcc

                pcc = translate_periodic_to_center(pc, dau, o)
                @test pcc.dm == translate_periodic_to_center(pc.dm, dau, o)
                @test pcc.gas == translate_periodic_to_center(pc.gas, dau, o)

                pccc = deepcopy(pc)
                translate_periodic_to_center!(pccc, dau, o)
                @test pccc == pcc
            end
        end

        @testset "Comoving" begin
            z = 1.5
            a = rand(Float32, 3, 100)
            au = a * u"m"
            n = rand(Float32, 100)
            nu = n * u"m^-3"
            m = copy(n)

            propexp = [(:pos, 1), (:vel, 1), (:n, -3)]

            acom = a .* (1 + z)
            aucom = au .* (1 + z)
            ncom = n .* (1 + z)^-3
            nucom = nu .* (1 + z)^-3

            aphys = a ./ (1 + z)
            auphys = au ./ (1 + z)
            nphys = n ./ (1 + z)^-3
            nuphys = nu ./ (1 + z)^-3

            @testset "Particle Comoving" begin
                p = Particles(:dm)
                p.pos = a
                p.n = n
                p.mass = m

                @test a .* CP.factor_to_comoving(z, 1) ≈ acom
                @test n .* CP.factor_to_comoving(z, -3) ≈ ncom
                @test a .* CP.factor_to_physical(z, 1) ≈ aphys
                @test n .* CP.factor_to_physical(z, -3) ≈ nphys

                pc = to_comoving(p, z; propexp)
                @test pc.pos ≈ acom
                @test pc.n ≈ ncom
                @test pc.mass ≈ m

                pcc = deepcopy(p)
                to_comoving!(pcc, z; propexp)
                @test pcc == pc

                pc = to_physical(p, z; propexp)
                @test pc.pos ≈ aphys
                @test pc.n ≈ nphys
                @test pc.mass ≈ m

                pcc = deepcopy(p)
                to_physical!(pcc, z; propexp)
                @test pcc == pc

                p = Particles(:dm)
                p.pos = au
                p.n = nu
                p.mass = m

                pc = to_comoving(p, z; propexp)
                @test pc.pos ≈ aucom
                @test pc.n ≈ nucom
                @test pc.mass ≈ m

                pcc = deepcopy(p)
                to_comoving!(pcc, z; propexp)
                @test pcc == pc

                pc = to_physical(p, z; propexp)
                @test pc.pos ≈ auphys
                @test pc.n ≈ nuphys
                @test pc.mass ≈ m

                pcc = deepcopy(p)
                to_physical!(pcc, z; propexp)
                @test pcc == pc
            end

            @testset "Particle Collection Comoving" begin
                dm = Particles(:dm)
                dm.pos = copy(au)
                dm.n = copy(n)
                dm.mass = copy(m)
                gas = Particles(:gas)
                gas.pos = copy(a)
                gas.n = copy(nu)
                gas.mass = copy(m)

                pc = RedshiftParticleCollection(z, dm, gas)

                pcc = to_comoving(pc; propexp)
                @test pcc.dm == to_comoving(dm, z; propexp)
                @test pcc.gas == to_comoving(gas, z; propexp)

                pccc = deepcopy(pc)
                to_comoving!(pccc; propexp)
                @test pccc == pcc

                pcc = to_physical(pc; propexp)
                @test pcc.dm == to_physical(dm, z; propexp)
                @test pcc.gas == to_physical(gas, z; propexp)

                pccc = deepcopy(pc)
                to_physical!(pccc; propexp)
                @test pccc == pcc

                pccc = deepcopy(pc)
                ap = pccc.all
                @test_throws ErrorException to_comoving(ap, z; propexp)
                @test_throws ErrorException to_physical(ap, z; propexp)

                to_comoving!(ap, z; propexp)
                @test pccc == to_comoving(pc; propexp)

                pccc = deepcopy(pc)
                ap = pccc.all
                to_physical!(ap, z; propexp)
                @test pccc == to_physical(pc; propexp)
            end
        end
    end

    @testset "Misc Operations" begin
        ids = sample(1:1000, 100; replace=false)
        a = rand(3, 100)
        bu = rand(100) * u"kg"

        p = Particles(:dm)
        p.id = ids
        p.pos = copy(a)
        p.mass = copy(bu)

        @testset "Particle Sorting" begin
            pc = sort(p, :id)
            @test issorted(pc.id)
            @test p.id == ids
            ind = searchsortedfirst(pc.id, p.id[1])
            @test pc.pos[:, ind] == p.pos[:, 1]
            @test pc.mass[ind] == p.mass[1]

            pc = sort(p, :mass; affect=[:mass, :pos], alg=RadixSort)
            @test issorted(pc.mass)
            @test haskey(pc, :pos)
            @test !haskey(pc, :id)

            pc = sort(p, :mass; affect=[:pos], alg=RadixSort)
            @test !haskey(pc, :mass)

            pc = deepcopy(p)
            sort!(pc, :id; alg=RadixSort)
            @test issorted(pc.id)
            ind = searchsortedfirst(pc.id, p.id[1])
            @test pc.pos[:, ind] == p.pos[:, 1]
            @test pc.mass[ind] == p.mass[1]
        end

        @testset "Particle Collection Sorting" begin
            dm = p
            gas = Particles(:gas, deepcopy(p.props))
            pc = ParticleCollection(dm, gas)

            pcc = sort(pc, :id)
            @test pc != pcc
            @test pcc.dm == sort(dm, :id)
            @test pcc.gas == sort(gas, :id)

            affect = [:mass, :pos]
            pcc = sort(pc, :mass; affect, alg=RadixSort)
            @test pcc.dm == sort(dm, :mass; affect)
            @test pcc.gas == sort(gas, :mass; affect)

            affect = [(:dm, [:mass, :pos]), (:gas, [:mass])]
            pcc = sort(pc, :mass; affect, alg=QuickSort)
            @test pcc.dm == sort(dm, :mass; affect=[:mass, :pos])
            @test pcc.gas == sort(gas, :mass; affect=[:mass])

            pcc = deepcopy(pc)
            sort!(pcc, :id; alg=RadixSort)
            @test pcc == sort(pc, :id)

            ap = pc.all
            @test_throws ErrorException sort(ap, :id; affect)
            @test_throws ErrorException sort!(ap, :id; alg=RadixSort)
        end

        @testset "Particle Filtering" begin
            massmin = 0.5u"kg"
            pc = filter(p -> p.mass .> massmin, p)
            @test all(pc.mass .> massmin)
            @test p.mass == bu
            ind = findfirst(>(massmin), p.mass)
            @test pc.pos[:, 1] == p.pos[:, ind]
            @test pc.id[1] == p.id[ind]

            @test filter(p -> findall(p.mass .> massmin), p) == pc

            pc = filter(p -> p.mass .> massmin, p; affect=[:id, :mass])
            @test all(pc.mass .> massmin)
            @test haskey(pc, :id)
            @test !haskey(pc, :pos)

            pc = deepcopy(p)
            filter!(p -> p.mass .> massmin, pc)
            @test all(pc.mass .> massmin)
            @test p.mass == bu
            ind = findfirst(>(massmin), p.mass)
            @test pc.pos[:, 1] == p.pos[:, ind]
            @test pc.id[1] == p.id[ind]

            ids_wanted = sample(1:1000, 100; replace=false)

            pc = filter(p; ids=ids_wanted)
            @test all(in.(pc.id, (ids_wanted,)))
            @test !any(in.(setdiff(ids, pc.id), (ids_wanted,)))
            ind = findfirst(in(ids_wanted), p.id)
            @test pc.pos[:, 1] == p.pos[:, ind]
            @test pc.mass[1] == p.mass[ind]

            @test filter(p; ids=Set(ids_wanted)) == pc

            pc = filter(p; ids=ids_wanted, affect=[:id, :mass])
            @test all(in.(pc.id, (ids_wanted,)))
            @test !any(in.(setdiff(ids, pc.id), (ids_wanted,)))
            @test haskey(pc, :mass)
            @test !haskey(pc, :pos)

            pc = deepcopy(p)
            filter!(pc; ids=ids_wanted)
            @test all(in.(pc.id, (ids_wanted,)))
            @test !any(in.(setdiff(ids, pc.id), (ids_wanted,)))
            ind = findfirst(in(ids_wanted), p.id)
            @test pc.pos[:, 1] == p.pos[:, ind]
            @test pc.mass[1] == p.mass[ind]

            @test filter!(deepcopy(p); ids=Set(ids_wanted)) == pc
        end

        @testset "Particle Collection Filtering" begin
            dm = p
            gas = Particles(:gas, deepcopy(p.props))
            pc = ParticleCollection(dm, gas)

            massmin = 0.5u"kg"
            f = p -> p.mass .> massmin
            pcc = filter(f, pc)
            @test pcc.dm == filter(f, dm)

            @test filter(p -> findall(p.mass .> massmin), pc) == pcc

            affect = [:id, :mass]
            pcc = filter(f, pc; affect)
            @test pcc.dm == filter(f, dm; affect)

            affect = [(:dm, [:id, :mass])]
            pcc = filter(f, pc; affect)
            @test pcc.dm == filter(f, dm; affect=[:id, :mass])
            @test !haskey(pcc, :gas)

            pcc = deepcopy(pc)
            filter!(f, pcc)
            @test pcc == filter(f, pc)

            ids_wanted = sample(1:1000, 100; replace=false)

            pcc = filter(pc; ids=ids_wanted)
            @test pcc.dm == filter(dm; ids=ids_wanted)
            @test filter(pc; ids=Set(ids_wanted)) == pcc

            affect = [:id, :mass]
            pcc = filter(pc; ids=ids_wanted, affect)
            @test pcc.dm == filter(dm; ids=ids_wanted, affect)

            affect = [(:gas, [:id])]
            pcc = filter(pc; ids=ids_wanted, affect)
            @test !haskey(pcc, :dm)
            @test pcc.gas == filter(gas; ids=ids_wanted, affect=[:id])

            pcc = deepcopy(pc)
            filter!(pcc; ids=ids_wanted)
            @test pcc == filter(pcc; ids=ids_wanted)

            @test filter!(deepcopy(pc); ids=Set(ids_wanted)) == pcc

            ap = pc.all
            @test_throws ErrorException filter(f, ap; affect)
            @test_throws ErrorException filter!(f, ap)
            @test_throws ErrorException filter(ap; ids=ids_wanted)
            @test_throws ErrorException filter!(ap; ids=ids_wanted)
        end

        @testset "Particle Removing" begin
            ids_remove = sample(1:1000, 100; replace=false)
            ids_keep = setdiff(1:1000, ids_remove)

            pc_check = filter(p; ids=ids_keep)

            @test delete(p; ids=ids_remove) == pc_check
            @test delete(p; ids=Set(ids_remove)) == pc_check

            pc = delete(p; ids=ids_remove, affect=[:id, :mass])
            @test pc == filter(p; ids=ids_keep, affect=[:id, :mass])

            pc = deepcopy(p)
            delete!(pc; ids=ids_remove)
            @test pc == pc_check

            @test delete!(deepcopy(p); ids=Set(ids_remove)) == pc_check
        end

        @testset "Particle Collection Removing" begin
            dm = p
            gas = Particles(:gas, deepcopy(p.props))
            pc = ParticleCollection(dm, gas)

            ids_remove = sample(1:1000, 100; replace=false)

            pcc = delete(pc; ids=ids_remove)
            @test pcc.dm == delete(dm; ids=ids_remove)
            @test delete(pc; ids=Set(ids_remove)) == pcc

            affect = [:id, :mass]
            pcc = delete(pc; ids=ids_remove, affect)
            @test pcc.dm == delete(dm; ids=ids_remove, affect)

            affect = [(:gas, [:id])]
            pcc = delete(pc; ids=ids_remove, affect)
            @test !haskey(pcc, :dm)
            @test pcc.gas == delete(gas; ids=ids_remove, affect=[:id])

            pcc = deepcopy(pc)
            delete!(pcc; ids=ids_remove)
            @test pcc == delete(pcc; ids=ids_remove)

            @test delete!(deepcopy(pc); ids=Set(ids_remove)) == pcc

            ap = pc.all
            @test_throws Exception delete(ap; ids=ids_remove)
            @test_throws Exception delete!(ap; ids=ids_remove)
        end
    end

    @testset "Geometry" begin
        # create square/cube/hypercube of positions between 0 and 1
        pos2 = reshape(Iterators.product(0:0.1:1, 0:0.1:1) |> collect, :)
        pos2 = reduce(hcat, getindex.(pos2, i) for i in eachindex(pos2[1]))'

        pos3 = reshape(Iterators.product(0:0.1:1, 0:0.1:1, 0:0.1:1) |> collect, :)
        pos3 = reduce(hcat, getindex.(pos3, i) for i in eachindex(pos3[1]))'

        pos4 = reshape(Iterators.product(0:0.1:1, 0:0.1:1, 0:0.1:1, 0:0.1:1) |> collect, :)
        pos4 = reduce(hcat, getindex.(pos4, i) for i in eachindex(pos4[1]))'

        @testset "Hyperrectangle" begin
            hrect = CosmoHyperrectangle(0.5rand(4), 0.5(1 .+ rand(4)))
            @test hrect isa CosmoHyperrectangle{Float64,4}

            hrect = CosmoHyperrectangle([0, 0, 0] .// 100, [15, 25, 20] .// 100)
            @test hrect isa CosmoHyperrectangle{Float64,3}
            @test hrect isa CosmoCuboid{Float64}

            # construct hyperrectangles from center and side lengths
            @test CosmoHyperrectangle([75, 125, 100] .// 1000, ([15, 25, 20] .// 100)...) == hrect
            @test CosmoCuboid([75, 125, 100] .// 1000, ([15, 25, 20] .// 100)...) == hrect
            @test CosmoRectangle([75, 125] .// 1000, ([15, 25] .// 100)...) ==
                  CosmoRectangle([0, 0] .// 100, [15, 25] .// 100)

            # test constructors
            @test CosmoHyperrectangle([0, 0, 0], [0.15, 0.25, 0.20]) isa CosmoCuboid{Float64}
            @test CosmoHyperrectangle{Float64,3}([0, 0, 0.0], [0.15, 0.25, 0.20]) isa CosmoCuboid{Float64}
            @test CosmoHyperrectangle{Float64,3}([0, 0, 0], [0.15, 0.25, 1 // 10]) isa CosmoCuboid{Float64}
            @test_throws AssertionError CosmoHyperrectangle{Float64,4}([0, 0, 0.0], [0.15, 0.25, 0.20])
            @test CosmoCuboid([0, 0, 0], [0.15, 0.25, 0.20]) isa CosmoCuboid{Float64}
            @test_throws AssertionError CosmoCuboid([0, 0], [0.15, 0.25])
            @test CosmoRectangle([0, 0], [0.25, 0.20]) isa CosmoRectangle{Float64}
            @test_throws AssertionError CosmoRectangle([0, 0, 0], [0.15, 0.25, 0.20])
            @test_throws AssertionError CosmoRectangle([0, 0], [0.15, 0.25, 0.20])
            @test CosmoCuboid{Float64}([0, 0, 0.0], [0.15, 0.25, 0.20]) isa CosmoCuboid{Float64}
            @test CosmoCuboid{Float32}([0, 0, 0], [0.15, 0.25, 1 // 10]) isa CosmoCuboid{Float32}
            @test CosmoRectangle{Float64}([0, 0.0], [0.25, 0.20]) isa CosmoRectangle{Float64}
            @test CosmoRectangle{Float64}([0, 0], [0.25, 0.20]) isa CosmoRectangle{Float64}

            @test CosmoHyperrectangle([0, 0], 0.15, 0.25) isa CosmoRectangle{Float64}
            @test CosmoHyperrectangle([0, 0, 0], 0.15, 0.25, 0.20) isa CosmoCuboid{Float64}
            @test CosmoHyperrectangle([0, 0, 0, 0], 0.15, 0.25, 0.20, 0.54) isa CosmoHyperrectangle{Float64,4}
            @test_throws Exception CosmoHyperrectangle([0, 0], 0.15, 0.25, 0.54)

            r = 5 // 2
            @test CosmoHypercube([0, 0, 0, 0], r) == CosmoHyperrectangle([-r, -r, -r, -r], [r, r, r, r])
            @test CosmoCube([0, 0, 0], r) == CosmoCuboid([-r, -r, -r], [r, r, r])
            @test CosmoSquare([0, 0], r) == CosmoRectangle([-r, -r], [r, r])

            cuboid = CosmoCuboid([0, 0, 0] .// 100, [15, 25, 19] .// 100)
            @test CP.geometry_enclosing_center(cuboid) ≈ [15, 25, 19] .// 200 |> float
            @test CP.geometry_enclosing_center(CosmoCube(zeros(Int, 3), 4 // 1)) == zeros(3)
            @test all(CP.geometry_enclosing_corners(cuboid) .≈ (zeros(3), [15, 25, 19] .// 100 |> float))

            pos3in = pos3[:, CP.mask_in(pos3, cuboid)]
            @test all(0 .≤ pos3in[1, :] .≤ 0.15) &&
                  all(0 .≤ pos3in[2, :] .≤ 0.25) &&
                  all(0 .≤ pos3in[3, :] .≤ 0.19)

            rect = CosmoRectangle([0, 0] .// 100, [25, 19] .// 100)
            pos2in = pos2[:, CP.mask_in(pos2, rect)]
            @test all(0 .≤ pos2in[1, :] .≤ 0.25) && all(0 .≤ pos2in[2, :] .≤ 0.19)

            hrect = CosmoHyperrectangle([0, 0, 0, 0] .// 100, [15, 25, 19, 35] .// 100)
            pos4in = pos4[:, CP.mask_in(pos4, hrect)]
            @test all(0 .≤ pos4in[1, :] .≤ 0.15) &&
                  all(0 .≤ pos4in[2, :] .≤ 0.25) &&
                  all(0 .≤ pos4in[3, :] .≤ 0.19) &&
                  all(0 .≤ pos4in[4, :] .≤ 0.35)

            @test_throws AssertionError CP.mask_in(pos3, hrect)
        end

        @testset "Hypersphere" begin
            hsphere = CosmoHypersphere(0.5rand(4), 0.5)
            @test hsphere isa CosmoHypersphere{Float64,4}

            hsphere = CosmoHypersphere([0, 0, 0] .// 100, 15 // 100)
            @test hsphere isa CosmoHypersphere{Float64,3}
            @test hsphere isa CosmoSphere{Float64}

            @test hsphere == CosmoHypersphere([0, 0, 0], 15 // 100)

            # test constructors
            @test CosmoHypersphere([0.15, 0.25, 0.20], 2) isa CosmoSphere{Float64}
            @test CosmoHypersphere{Float64,3}([0.15, 0.25, 0.20], 0.2) isa CosmoSphere{Float64}
            @test CosmoHypersphere{Float32,3}([0.15, 0.25, 0.20], 0.2) isa CosmoSphere{Float32}
            @test_throws AssertionError CosmoHypersphere{Float64,4}([0.15, 0.25, 0.20], 0.2)
            @test CosmoSphere([0.15, 0.25, 0.20], 0.2) isa CosmoSphere{Float64}
            @test_throws AssertionError CosmoSphere([0.15, 0.25], 0.3)
            @test CosmoCircle([0.25, 0.20], 0.2) isa CosmoCircle{Float64}
            @test_throws AssertionError CosmoCircle([0.15, 0.25, 0.20], 0.2)
            @test CosmoSphere{Float64}([0.15, 0.25, 0.20], 0.2) isa CosmoSphere{Float64}
            @test CosmoSphere{Float32}([0.15, 0.25, 0.20], 1) isa CosmoSphere{Float32}
            @test CosmoCircle{Float64}([0.25, 0.20], 0.2) isa CosmoCircle{Float64}

            r = 12 // 100

            center = [10, 20, 30] .// 100
            sphere = CosmoSphere(center, r)
            @test CP.geometry_enclosing_center(sphere) ≈ center |> float
            @test all(
                CP.geometry_enclosing_corners(sphere) .≈
                ([-2, 8, 18] .// 100 |> float, [22, 32, 42] .// 100 |> float),
            )

            pos3in = pos3[:, CP.mask_in(pos3, sphere)]
            @test all(sum(abs2, pos3in .- center; dims=1) .≤ r^2)
            pos3notin = pos3[:, .~CP.mask_in(pos3, sphere)]
            @test all(sum(abs2, pos3notin .- center; dims=1) .> r^2)

            center = [20, 30] .// 100
            circle = CosmoCircle(center, r)
            pos2in = pos2[:, CP.mask_in(pos2, circle)]
            @test all(sum(abs2, pos2in .- center; dims=1) .≤ r^2)

            center = [10, 20, 30, 40] .// 100
            hsphere = CosmoHypersphere(center, r)
            pos4in = pos4[:, CP.mask_in(pos4, hsphere)]
            @test all(sum(abs2, pos4in .- center; dims=1) .≤ r^2)

            sphere = CosmoSphere(zeros(3), r)
            @test CosmoSphere(r) == sphere
            pos3in = pos3[:, CP.mask_in(pos3, sphere)]
            @test all(sum(abs2, pos3in; dims=1) .≤ r^2)

            circle = CosmoCircle(zeros(2), r)
            @test CosmoCircle(r) == circle
            pos2in = pos2[:, CP.mask_in(pos2, circle)]
            @test all(sum(abs2, pos2in; dims=1) .≤ r^2)

            hsphere = CosmoHypersphere(zeros(4), r)
            @test CosmoHypersphere(4, r) == hsphere
            pos4in = pos4[:, CP.mask_in(pos4, hsphere)]
            @test all(sum(abs2, pos4in; dims=1) .≤ r^2)

            @test_throws AssertionError CP.mask_in(pos3, hsphere)
        end

        @testset "Cylinder" begin
            scyl = CosmoStandingCylinder(0.5rand(3), 0.5, 0.5)
            @test scyl isa CosmoStandingCylinder{Float64}

            scyl = CosmoStandingCylinder([1, 2, 3] .// 10, 2 // 10, 1 // 10)
            @test scyl isa CosmoStandingCylinder{Float64}

            @test CosmoStandingCylinder(rand(3), 1, 1 // 10) isa CosmoStandingCylinder{Float64}
            @test CosmoStandingCylinder{Float32}(rand(3), 1, 1 // 10) isa CosmoStandingCylinder{Float32}
            @test_throws AssertionError CosmoStandingCylinder(rand(4), 1, 1 // 10)

            center = [1, 2, 3] .// 10
            h = 24 // 100
            r = 15 // 100
            scyl = CosmoStandingCylinder(center, h, r)

            @test CP.geometry_enclosing_center(scyl) ≈ center |> float
            @test all(
                CP.geometry_enclosing_corners(scyl) .≈
                ([-5, 5, 18] .// 100 |> float, [25, 35, 42] .// 100 |> float),
            )

            pos3in = pos3[:, CP.mask_in(pos3, scyl)]
            @test @views all(sum(abs2, pos3in[1:2, :] .- center[1:2]; dims=1) .≤ r^2)
            @test @views all(18 // 100 .≤ pos3in[3, :] .≤ 42 // 100)
            pos3notin = pos3[:, .~CP.mask_in(pos3, scyl)]
            @test @views all(
                vec(sum(abs2, pos3notin[1:2, :] .- center[1:2]; dims=1) .> r^2) .|
                (18 // 100 .> pos3notin[3, :]) .|
                (pos3notin[3, :] .> 42 // 100),
            )
            @test_throws AssertionError CP.mask_in(pos2, scyl)

            cyl = CosmoCylinder(0.5rand(3), 0.5rand(3), 0.5)
            @test cyl isa CosmoCylinder{Float64}

            cyl = CosmoCylinder([1, 2, 3] .// 10, [2, 3, 4] .// 10, 1 // 10)
            @test cyl isa CosmoCylinder{Float64}

            @test CosmoCylinder(rand(3), rand(3), 1 // 10) isa CosmoCylinder{Float64}
            @test CosmoCylinder(rand(3), [0, 0, 0], 1 // 10) isa CosmoCylinder{Float64}
            @test CosmoCylinder{Float32}(rand(3), [0, 0, 0], 1 // 10) isa CosmoCylinder{Float32}
            @test_throws AssertionError CosmoCylinder(rand(4), rand(3), 0.1)
            @test_throws AssertionError CosmoCylinder(rand(3), rand(2), 0.1)

            center = center
            startpos = [10, 20, 18] .// 100
            endpos = [10, 20, 42] .// 100
            r = r
            cyl = CosmoCylinder(startpos, endpos, r)

            @test CP.geometry_enclosing_center(cyl) ≈ center |> float
            @test all(CP.geometry_enclosing_corners(cyl) .≈ ([-5, 5, 18] .// 100, [25, 35, 42] .// 100))

            # right mask for standing cylinder
            @test scyl == cyl
            @test CP.mask_in(pos3, scyl) == CP.mask_in(pos3, cyl)
            @test_throws AssertionError CP.mask_in(pos2, cyl)

            @test CosmoStandingCylinder(cyl) == scyl
            @test CosmoCylinder(scyl) == cyl

            # same result switching random start and end points of the cylinder
            p1 = rand(3)
            p2 = rand(3)
            r = 0.15
            c1 = CosmoCylinder(p1, p2, r)
            c2 = CosmoCylinder(p2, p1, r)
            @test CP.geometry_enclosing_center(c1) ≈ CP.geometry_enclosing_center(c2)
            @test all(CP.geometry_enclosing_corners(c1) .≈ CP.geometry_enclosing_corners(c2))
            @test CP.mask_in(pos3, c1) == CP.mask_in(pos3, c2)

            @test_throws AssertionError CosmoStandingCylinder(c1)

            scyl = CosmoStandingCylinder([1, 2, 3] .// 10, 2 // 10, 1 // 10)
            @test scyl == CosmoStandingCylinder([1, 2, 3] .// 10, 2 // 10, 1 // 10)
            @test scyl ≈ CosmoCylinder([1, 2, 2] .// 10, [1, 2, 4] .// 10, 1 // 10)
            @test CosmoCylinder([1, 2, 2] .// 10, [1, 2, 4] .// 10, 1 // 10) ≈ scyl
        end

        @testset "Filtering" begin
            p = Particles(:dm)
            p.pos = pos3
            p.mass = rand(size(pos3, 2))

            pu = Particles(:dm)
            pu.pos = pos3 .* u"m"
            pu.mass = rand(size(pos3, 2))

            center = rand(3)
            sph = CosmoSphere(center, 0.2)
            sphu = CosmoSphere(center * u"m", 0.2u"m")
            mask = CP.mask_in(pos3, sph)

            @test filter(p, sph) == p[mask]
            @test filter(pu, sphu) == pu[mask]

            pc = filter(p, sph; affect=[:pos])
            @test pc.pos == p.pos[:, mask]
            @test !haskey(pc, :mass)

            pc = deepcopy(p)
            filter!(pc, sph)
            @test pc == p[mask]

            # Particle collection
            dm = p
            gas = Particles(:gas, deepcopy(p.props))
            pc = ParticleCollection(dm, gas)

            pcc = filter(pc, sph)
            @test pcc.dm == dm[mask]

            pcc = filter(pc, sph; affect=[:pos])
            @test pcc.dm == filter(dm, sph; affect=[:pos])

            pcc = filter(pc, sph; affect=[(:dm, [:pos, :mass])])
            @test pcc.dm == filter(dm, sph; affect=[:pos, :mass])
            @test !haskey(pcc, :gas)

            pcc = deepcopy(pc)
            filter!(pcc, sph)
            @test pcc == filter(pc, sph)

            ap = pc.all
            @test_throws ErrorException filter(ap, sph)
            @test_throws ErrorException filter!(ap, sph)
        end
    end

    @testset "Properties" begin
        @testset "Norm" begin
            a2 = rand(2, 100)
            a3 = rand(3, 100)
            a4 = rand(4, 100)

            n2 = norm.(eachcol(a2))
            n3 = norm.(eachcol(a3))
            n4 = norm.(eachcol(a4))

            @test colnorm(a2) ≈ n2
            @test colnorm(a3) ≈ n3
            @test colnorm(a4) ≈ n4

            @test colnorm(a2 * u"m") ≈ n2 * u"m"
            @test colnorm(a3 * u"m") ≈ n3 * u"m"
            @test colnorm(a4 * u"m") ≈ n4 * u"m"

            @test colnorm2(a2) ≈ n2 .^ 2
            @test colnorm2(a3) ≈ n3 .^ 2
            @test colnorm2(a4) ≈ n4 .^ 2

            @test colnorm2(a2 * u"m") ≈ (n2 * u"m") .^ 2
            @test colnorm2(a3 * u"m") ≈ (n3 * u"m") .^ 2
            @test colnorm2(a4 * u"m") ≈ (n4 * u"m") .^ 2

            am2 = similar(a2, Union{Float64,Missing})
            am3 = similar(a3, Union{Float64,Missing})
            am4 = similar(a4, Union{Float64,Missing})
            am2 .= a2
            am3 .= a3
            am4 .= a4
            am2[:, 1] .= missing
            am3[:, 1] .= missing
            am4[:, 1] .= missing

            @test all(colnorm(am2) .≈ n2) |> ismissing
            @test all(colnorm(am3) .≈ n3) |> ismissing
            @test all(colnorm(am4) .≈ n4) |> ismissing

            @test all(colnorm(am2 * u"m") .≈ n2 * u"m") |> ismissing
            @test all(colnorm(am3 * u"m") .≈ n3 * u"m") |> ismissing
            @test all(colnorm(am4 * u"m") .≈ n4 * u"m") |> ismissing

            @test all(colnorm2(am2) .≈ n2 .^ 2) |> ismissing
            @test all(colnorm2(am3) .≈ n3 .^ 2) |> ismissing
            @test all(colnorm2(am4) .≈ n4 .^ 2) |> ismissing

            @test all(colnorm2(am2 * u"m") .≈ (n2 * u"m") .^ 2) |> ismissing
            @test all(colnorm2(am3 * u"m") .≈ (n3 * u"m") .^ 2) |> ismissing
            @test all(colnorm2(am4 * u"m") .≈ (n4 * u"m") .^ 2) |> ismissing
        end

        @testset "Norm around Origin" begin
            a2 = rand(2, 100)
            a3 = rand(3, 100)
            a4 = rand(4, 100)

            o2 = rand(2)
            o3 = rand(3)
            o4 = rand(4)

            n2 = norm.(eachcol(a2 .- o2))
            n3 = norm.(eachcol(a3 .- o3))
            n4 = norm.(eachcol(a4 .- o4))

            @test colnorm(a2, o2) ≈ n2
            @test colnorm(a3, o3) ≈ n3
            @test colnorm(a4, o4) ≈ n4

            @test colnorm(a2 * u"m", o2 * u"m") ≈ n2 * u"m"
            @test colnorm(a3 * u"m", o3 * u"m") ≈ n3 * u"m"
            @test colnorm(a4 * u"m", o4 * u"m") ≈ n4 * u"m"

            @test colnorm2(a2, o2) ≈ n2 .^ 2
            @test colnorm2(a3, o3) ≈ n3 .^ 2
            @test colnorm2(a4, o4) ≈ n4 .^ 2

            @test colnorm2(a2 * u"m", o2 * u"m") ≈ (n2 * u"m") .^ 2
            @test colnorm2(a3 * u"m", o3 * u"m") ≈ (n3 * u"m") .^ 2
            @test colnorm2(a4 * u"m", o4 * u"m") ≈ (n4 * u"m") .^ 2

            am2 = similar(a2, Union{Float64,Missing})
            am3 = similar(a3, Union{Float64,Missing})
            am4 = similar(a4, Union{Float64,Missing})
            am2 .= a2
            am3 .= a3
            am4 .= a4
            am2[:, 1] .= missing
            am3[:, 1] .= missing
            am4[:, 1] .= missing

            @test all(colnorm(am2, o2) .≈ n2) |> ismissing
            @test all(colnorm(am3, o3) .≈ n3) |> ismissing
            @test all(colnorm(am4, o4) .≈ n4) |> ismissing

            @test all(colnorm(am2 * u"m", o2 * u"m") .≈ n2 * u"m") |> ismissing
            @test all(colnorm(am3 * u"m", o3 * u"m") .≈ n3 * u"m") |> ismissing
            @test all(colnorm(am4 * u"m", o4 * u"m") .≈ n4 * u"m") |> ismissing

            @test all(colnorm2(am2, o2) .≈ n2 .^ 2) |> ismissing
            @test all(colnorm2(am3, o3) .≈ n3 .^ 2) |> ismissing
            @test all(colnorm2(am4, o4) .≈ n4 .^ 2) |> ismissing

            @test all(colnorm2(am2 * u"m", o2 * u"m") .≈ (n2 * u"m") .^ 2) |> ismissing
            @test all(colnorm2(am3 * u"m", o3 * u"m") .≈ (n3 * u"m") .^ 2) |> ismissing
            @test all(colnorm2(am4 * u"m", o4 * u"m") .≈ (n4 * u"m") .^ 2) |> ismissing
        end

        @testset "Dot Product" begin
            a2 = rand(2, 100)
            a3 = rand(3, 100)
            a4 = rand(4, 100)
            b2 = rand(2, 100)
            b3 = rand(3, 100)
            b4 = rand(4, 100)

            c2 = reduce(hcat, [dot(ai, bi) for (ai, bi) in zip(eachcol(a2), eachcol(b2))]) |> vec
            c3 = reduce(hcat, [dot(ai, bi) for (ai, bi) in zip(eachcol(a3), eachcol(b3))]) |> vec
            c4 = reduce(hcat, [dot(ai, bi) for (ai, bi) in zip(eachcol(a4), eachcol(b4))]) |> vec

            for (a, b, c) in zip([a2, a3, a4], [b2, b3, b4], [c2, c3, c4])
                @test coldot(a, b) ≈ c
                @test coldot(a * u"m", b) ≈ c * u"m"
                @test coldot(a, b * u"m") ≈ c * u"m"
                @test coldot(a * u"m", b * u"m") ≈ c * u"m^2"

                am = similar(a, Union{Float64,Missing})
                am .= a
                am[:, 1] .= missing
                bm = similar(b, Union{Float64,Missing})
                bm .= b
                bm[:, 1] .= missing

                @test all(coldot(am, b) .≈ c) |> ismissing
                @test all(coldot(a * u"m", bm) .≈ c * u"m") |> ismissing
                @test all(coldot(a, bm * u"m") .≈ c * u"m") |> ismissing
                @test all(coldot(am * u"m", bm) .≈ c * u"m") |> ismissing
                @test all(coldot(am * u"m", bm * u"m") .≈ c * u"m^2") |> ismissing
            end

            am = similar(a3, Union{Float64,Nothing})
            am .= a3
            @test coldot(am, b3) ≈ c3
        end

        @testset "Cross Product" begin
            a = rand(3, 100)
            b = rand(3, 100)
            au = a * u"m"
            bu = b * u"m"

            c = reduce(hcat, [cross(ai, bi) for (ai, bi) in zip(eachcol(a), eachcol(b))])

            @test colcross(a, b) ≈ c
            @test colcross(au, b) ≈ c * u"m"
            @test colcross(a, bu) ≈ c * u"m"
            @test colcross(au, bu) ≈ c * u"m^2"

            am = similar(a, Union{Float64,Missing})
            am .= a
            am[:, 1] .= missing
            bm = similar(b, Union{Float64,Missing})
            bm .= b
            bm[:, 1] .= missing

            @test all(colcross(am, b) .≈ c) |> ismissing
            @test all(colcross(au, bm) .≈ c * u"m") |> ismissing
            @test all(colcross(a, bm * u"m") .≈ c * u"m") |> ismissing
            @test all(colcross(am * u"m", bm) .≈ c * u"m") |> ismissing
            @test all(colcross(am * u"m", bm * u"m") .≈ c * u"m^2") |> ismissing

            am = similar(a, Union{Float64,Nothing})
            am .= a
            @test colcross(am, b) ≈ c
        end

        @testset "Angular momentum" begin
            a = rand(3, 100)
            b = rand(3, 100)
            au = a * u"m"
            bu = b * u"m"

            for m in [rand(), rand(100)]
                mu = m * u"kg"

                j = m' .* colcross(a, b)

                @test angmom(a, b, m) ≈ j
                @test angmom(au, b, m) ≈ j * u"m"
                @test angmom(a, bu, m) ≈ j * u"m"
                @test angmom(au, bu, m) ≈ j * u"m^2"
                @test angmom(a, b, mu) ≈ j * u"kg"
                @test angmom(au, b, mu) ≈ j * u"m*kg"
                @test angmom(a, bu, mu) ≈ j * u"m*kg"
                @test angmom(au, bu, mu) ≈ j * u"m^2*kg"

                for func in [angmomtot, angmomtot_stable]
                    @test func(a, b, m) ≈ sum(j; dims=2)
                    @test func(au, b, m) ≈ sum(j * u"m"; dims=2)
                    @test func(a, bu, m) ≈ sum(j * u"m"; dims=2)
                    @test func(au, bu, m) ≈ sum(j * u"m^2"; dims=2)
                    @test func(a, b, mu) ≈ sum(j * u"kg"; dims=2)
                    @test func(au, b, mu) ≈ sum(j * u"m*kg"; dims=2)
                    @test func(a, bu, mu) ≈ sum(j * u"m*kg"; dims=2)
                    @test func(au, bu, mu) ≈ sum(j * u"m^2*kg"; dims=2)
                end

                origin = rand(3)
                velorigin = rand(3)
                jorigin = m' .* colcross(a .- origin, b)
                @test angmom(au, bu, mu; origin=origin * u"m") ≈ jorigin * u"m^2*kg"
                @test angmomtot(au, bu, mu; origin=origin * u"m") ≈ sum(jorigin * u"m^2*kg"; dims=2)
                jorigin = m' .* colcross(a, b .- velorigin)
                @test angmom(au, bu, mu; velorigin=velorigin * u"m") ≈ jorigin * u"m^2*kg"
                @test angmomtot(au, bu, mu; velorigin=velorigin * u"m") ≈ sum(jorigin * u"m^2*kg"; dims=2)
                @test angmomtot_stable(au, bu, mu; velorigin=velorigin * u"m") ≈
                      sum(jorigin * u"m^2*kg"; dims=2)
                jorigin = m' .* colcross(a .- origin, b .- velorigin)
                @test angmom(au, bu, mu; origin=origin * u"m", velorigin=velorigin * u"m") ≈
                      jorigin * u"m^2*kg"
                @test angmomtot(au, bu, mu; origin=origin * u"m", velorigin=velorigin * u"m") ≈
                      sum(jorigin * u"m^2*kg"; dims=2)
                @test angmomtot_stable(au, bu, mu; origin=origin * u"m", velorigin=velorigin * u"m") ≈
                      sum(jorigin * u"m^2*kg"; dims=2)

                am = similar(a, Union{Float64,Missing})
                am .= a
                am[:, 1] .= missing
                bm = similar(b, Union{Float64,Missing})
                bm .= b
                bm[:, 1] .= missing

                @test all(angmom(am, b, m) .≈ j) |> ismissing
                @test all(angmom(au, bm, m) .≈ j * u"m") |> ismissing
                @test all(angmom(a, bm * u"m", m) .≈ j * u"m") |> ismissing
                @test all(angmom(am * u"m", bm, m) .≈ j * u"m") |> ismissing
                @test all(angmom(am * u"m", bm * u"m", m) .≈ j * u"m^2") |> ismissing
                @test all(angmom(am, b, mu) .≈ j * u"kg") |> ismissing
                @test all(angmom(au, bm, mu) .≈ j * u"m*kg") |> ismissing
                @test all(angmom(a, bm * u"m", mu) .≈ j * u"m*kg") |> ismissing
                @test all(angmom(am * u"m", bm, mu) .≈ j * u"m*kg") |> ismissing
                @test all(angmom(am * u"m", bm * u"m", mu) .≈ j * u"m^2*kg") |> ismissing

                am = similar(a, Union{Float64,Nothing})
                am .= a
                @test angmom(am, b, m) ≈ j
            end
        end

        @testset "Particles properties" begin
            dm = Particles(:dm, :pos => rand(3, 100), :vel => rand(Float32, 3, 100), :mass => 2)
            gas = Particles(
                :gas,
                :pos => rand(3, 100),
                :vel => rand(Float32, 3, 100),
                :mass => rand(100),
                :temp => rand(100),
                :test => 1,
                :mass2 => Fill(2, 100),
            )

            dmu = deepcopy(dm)
            gasu = deepcopy(gas)
            dmu.pos *= u"m"
            dmu.vel *= u"km/s"
            dmu.mass *= u"kg"
            gasu.pos *= u"m"
            gasu.vel *= u"km/s"
            gasu.mass *= u"kg"
            gasu.mass2 *= u"kg"
            gasu.temp *= u"K"

            for (dm, gas) in zip([dm, dmu], [gas, gasu])
                cp = ParticleCollection(dm, gas)
                ap = cp.all
                p = Particles(ap)

                @test meanpos(dm) == mean(dm.pos; dims=2)
                @test meanpos(dm; massweighted=false) == mean(dm.pos; dims=2)
                @test meanpos(gas) ≈ sum(gas.pos .* gas.mass'; dims=2) / sum(gas.mass)
                @test meanpos(gas; massweighted=false) == mean(gas.pos; dims=2)
                @test meanpos(ap) ≈ meanpos(p)

                @test meanprop(dm, :pos) == mean(dm.pos; dims=2)
                @test meanprop(gas, :pos; massprop=:mass2) == mean(gas.pos; dims=2)
                @test meanprop(gas, :temp) ≈ sum(gas.temp .* gas.mass) / sum(gas.mass)
                @test meanprop(gas, :temp; massweighted=false) == mean(gas.temp)
                @test meanprop(gas, :temp; massprop=:mass2) == mean(gas.temp)
                @test meanprop(dm, :mass) == dm.mass
                @test meanprop(dm, :mass; massweighted=false) == dm.mass
                @test meanprop(gas, :test) == gas.test
                @test meanprop(gas, :mass2) == gas.mass2[1]
                @test meanprop(gas, :mass2; massweighted=false) == gas.mass2[1]
                @test meanprop(gas, :mass2; massprop=:mass2) == gas.mass2[1]
                @test meanprop(ap, :mass) ≈ meanvel(p, :mass)

                @test meanvel(dm) == mean(dm.vel; dims=2)
                @test meanvel(dm; massweighted=false) == mean(dm.vel; dims=2)
                @test meanvel(gas) ≈ sum(gas.vel .* gas.mass'; dims=2) / sum(gas.mass)
                @test meanvel(gas; massweighted=false) == mean(gas.vel; dims=2)
                @test meanvel(ap) ≈ meanvel(p)

                @test sumprop(dm, :pos) == sum(dm.pos; dims=2)
                @test sumprop(gas, :pos) == sum(gas.pos; dims=2)
                @test sumprop(dm, :mass) == CP.particle_number(dm) * dm.mass
                @test sumprop(gas, :mass2) == sum(gas.mass2)
                @test sumprop(ap, :pos) == sumprop(p, :pos)

                @test angmom(dm; angmomprop=:j) == angmom(dm.pos, dm.vel, dm.mass)
                @test angmomtot(dm; angmomprop=:j) == angmomtot(dm.pos, dm.vel, dm.mass)
                @test angmomtot_stable(dm; angmomprop=:j) == angmomtot_stable(dm.pos, dm.vel, dm.mass)
                dm.j = angmom(dm)
                @test angmom(dm; angmomprop=:j) ≈ angmom(dm.pos, dm.vel, dm.mass)
                @test angmomtot(dm; angmomprop=:j) ≈ angmomtot(dm.pos, dm.vel, dm.mass)
                @test angmomtot_stable(dm; angmomprop=:j) ≈ angmomtot_stable(dm.pos, dm.vel, dm.mass)
                @test angmom(gas; angmomprop=:j) == angmom(gas.pos, gas.vel, gas.mass)
                @test angmomtot(gas; angmomprop=:j) == angmomtot(gas.pos, gas.vel, gas.mass)
                @test angmomtot_stable(gas; angmomprop=:j) == angmomtot_stable(gas.pos, gas.vel, gas.mass)

                origin = 2 * dm.pos[:, 1]
                velorigin = 2 * dm.vel[:, 1]

                for func in [angmom, angmomtot, angmomtot_stable]
                    @test func(gas; origin) == func(gas.pos, gas.vel, gas.mass; origin)
                    @test func(gas; velorigin) == func(gas.pos, gas.vel, gas.mass; velorigin)
                    @test func(gas; origin, velorigin) == func(gas.pos, gas.vel, gas.mass; origin, velorigin)
                end
            end
        end
    end
end

ext: jl


text:
"""
```julia
function lda(Z[, x]; conf, xc_func! = xc_lda(LDA_X, LDA_C_VWN), Vex = r -> -Z/r, δn=1e-8, maxit=100, xmax=25, Α=1e5)
```
Solve Kohn-Sham DFT Self-Consistent Field equations for an atom using
[local density approximation (LDA)](https://en.wikipedia.org/wiki/Local-density_approximation)

On input:

* `Z`        - nuclear charge
* `x`        - grid, xᵢ = log(rᵢ), default is 401 points -35:0.1:5
* `conf`     - electronic configuration, i.e. `c"[He] 2s2 2p4"`, by default `conf = atom[Z]`
* `xc_func!` - LDA exchange-correlation functional, default is Slater exchange + VWN correlation
* `Vex`      - external potential, Coulomb -Z/r by default
* `δn`       - density convergence criterion
* `maxit`    - maximum permissible number of SCF iterations
* `xmax`     - grid is enlarged to this value for the Poisson equation solution
* `Α`        - parameter for the convex matrix pencil formulation

On output the function returns tuple: `(energy  = (total = E, kinetic = Ek, hartree = Eh, xc = Exc, potential = Ep), density = ρ, orbitals = ψ, iterations = Int(it))`
where

*  `lda(...).energy.total` - total energy, a.u.
*  `lda(...).density`      - electron density ρ(x) = 1/4π * ∑nᵢ*ψᵢ²(x)
*  `lda(...).orbitals`     - orbitals ψᵢ(x), energy levels εᵢ, azimuthal and radial quantum numbers l, nᵣ, and populations nᵢ (as listed in `conf`).

# Example
[Hooke atom](https://en.wikipedia.org/wiki/Hooke's_atom)

```jldoctest
julia> lda(2, Vex = r -> 1/8 * r^2).energy.total
[ Info: A grid of 401 points x = -35.0000:0.1000:5.0000 is used
[ Info: Starting SCF procedure with convergence threshold δn ≤ 1.000000e-08
[ Info: cycle           energy          δn
[ Info:   0           2.103849      1.89556590
[ Info:   1           2.013080      0.35778036
[ Info:   2           2.023521      0.03964354
[ Info:   3           2.026081      0.00360748
[ Info:   4           2.026212      0.00040386
[ Info:   5           2.026228      0.00004527
[ Info:   6           2.026229      0.00000518
[ Info:   7           2.026229      0.00000060
[ Info:   8           2.026229      0.00000007
[ Info:   9           2.026229      0.00000001
[ Info: SCF has converged after 9 iterations
┌ Info: RESULTS SUMMARY:
│       ELECTRON KINETIC               0.627459
│       ELECTRON-ELECTRON              1.022579
│       EXCHANGE-CORRELATION          -0.523773
│       ELECTRON-NUCLEAR               0.899965
│       TOTAL ENERGY                   2.026229
└       VIRIAL RATIO                  -2.229260
2.0262293420299144
```
"""
function lda(Z::Real,
             x::AbstractRange = -35:0.1:5;
          conf::Array{Array{T,1},1} = atom[Z],
            xc::XC_Functional = xc_func(),
           Vex::Function = r -> -Z / r,
            δn::Real = sqrt(eps()),
         maxit::Integer = 100,
          xmax::Real = 25,
             Α::Real = 1e5) where T <: Real

  xc_func! = xc.xc_func

  μ = 1    # electron mass, a.u.
  β = 0.8  # initial value of the density mixing parameter

  # radial grid
  r, r², sqr, n, dx = exp.(x), exp.(2x), exp.(x/2), length(x), step(x)
  @info @sprintf("A grid of %3i points x = %5.4f:%5.4f:%5.4f is used", n, x[1], dx, x[end])

  ρᵢₙ, ρₒᵤₜ, prev_ρᵢₙ, prev_ρₒᵤₜ, V, vp, vh, vxc, εxc, ε = [similar(r) for _ = 1:10]

  Q = Nₑ(conf)

  # Thomas-Fermi density for neutral atom with Q electrons
  @. ρᵢₙ = TF(r, Q) * r
  @. vp  = Vex(r)

  # Indices of the main diagonal
  D = 1:n+1:n*n

  H = -1/2μ*laplacian(x); K = H[D]
  S = copy(H)

  # Large grid for the Poisson equation
  xl = first(x):step(x):max(xmax, last(x))
  nl = length(xl)
  rmin, rmax = exp(first(xl)), exp(last(xl))

  # Differential operator of the Poisson equation
  Δ = laplacian(xl) - 1/4*I

  @info @sprintf("Starting SCF procedure with convergence threshold δn ≤ %e", δn)
  @info "cycle\t\tenergy\t\tδn"

  Eₜₒₜ = 0.0; Δρ = 0.0; it = 0
  for i = 0:maxit
    # Poisson equation solution is Hartree potential, vh
    v = Δ \ vcat(-4π .* ρᵢₙ .* sqr .* r, zeros(nl - n))

    # Apply boundary conditions at r → 0 and r → ∞
    c₁ =  (√rmin * v[1] - √rmax * v[end] + Q * (1 - rmin)) / (rmax - rmin)
    c₂ = -(rmin  * √rmin * v[1] - √rmax * rmin * v[end] + Q * rmin * (1 - rmax)) / (rmax - rmin)
    @. @views vh = v[1:n] + c₁ * sqr + c₂ / sqr

    # Change variable vh(x) → vh(r)
    vh ./= sqr

    # Exchange-correlation potential and energy density
    xc_func!(ρᵢₙ ./ r, vxc, εxc)

    # Kohn-Sham potential
    @. V = (vp + vh + vxc) * r²

    # Schrödinger equation
    ∑nᵢεᵢ = 0.0; ρₒᵤₜ .= 0.0

    for l=0:length(conf)-1
      @. H[D] = K    + V + 1/2μ * (l + 1/2)^2
      @. S[D] = H[D] + Α * r²
      θ, y = eigen(Symmetric(H), Symmetric(S))
      @. ε = Α * θ / (1 - θ)
      for (nᵣ, nᵢ) in enumerate(conf[l+1])
        @views y[:, nᵣ] ./= sqrt(∫(dx, y[:, nᵣ] .^ 2 .* r²))
        @views ρₒᵤₜ .+= nᵢ / 4π * y[:, nᵣ] .^ 2
        ∑nᵢεᵢ += nᵢ * ε[nᵣ]
      end
    end

    # DFT total energy
    Eₜₒₜ = ∑nᵢεᵢ + 4π * ∫(dx, ρₒᵤₜ .* (-1/2 * vh .- vxc .+ εxc) .* r²)

    Δρ = 4π * ∫(dx, abs.(ρₒᵤₜ .- ρᵢₙ) .* r²)
    @info @sprintf "%3i\t%14.6f\t%14.8f" i Eₜₒₜ Δρ
    if Δρ < δn
      it = i
      break
    end

    # Convergence acceleration scheme
    if i > 0
      dρ = ∫(dx, ρₒᵤₜ .- prev_ρₒᵤₜ) / ∫(dx, ρᵢₙ .- prev_ρᵢₙ)
      β = min(max(1 / (1 - dρ), 0.1), 0.9)
    end
    @. prev_ρᵢₙ  = ρᵢₙ
    @. prev_ρₒᵤₜ = ρₒᵤₜ
    @. ρᵢₙ = (1 - β) * ρᵢₙ + β * ρₒᵤₜ
  end

  if Δρ > δn
    @warn "SCF does not converged after $maxit iterations!"
  else
    @info "SCF has converged after $it iterations"
  end

  # Wavefunctions (orbitals)
  @. V = vp + vh + vxc
  ∑ε, ρ, ψ = radial_shr_eq(V, x, conf = conf, μ = 1, Α = 1e5)

  # Total energy components:
  Ek  = ∑ε - 4π * ∫(dx, ρ .* V .* r²)
  Eh  = 2π * ∫(dx, ρ .* vh .* r²)
  Exc = 4π * ∫(dx, ρ .* εxc .* r²)
  Ep  = 4π * ∫(dx, ρ .* vp .* r²)
  E   = Ek + Ep + Eh + Exc

  # Virial ratio
  VR  = (Eh + Exc + Ep) / -Ek

  @info @sprintf("RESULTS SUMMARY:
      ELECTRON KINETIC       %16.6f
      ELECTRON-ELECTRON      %16.6f
      EXCHANGE-CORRELATION   %16.6f
      ELECTRON-NUCLEAR       %16.6f
      TOTAL ENERGY           %16.6f
      VIRIAL RATIO           %16.6f\n", Ek, Eh, Exc, Ep, E, VR)

  return (energy  = (total = E, kinetic = Ek, hartree = Eh, xc = Exc, potential = Ep),
          density = ρ, orbitals = ψ, iterations = Int(it))
end


text:
m = [01 02 03 04 05 06;
     19 20 21 22 23 24;
     07 08 09 10 11 12;
     13 14 15 16 17 18;
     25 26 27 28 29 30]

spp = ["sp1", "sp2", "sp3", "sp4", "sp5"]
spp_scrambled = ["sp1", "sp4", "sp2", "sp3", "sp5"]
dir = SpeciesDirectory(spp)

@testset "SpeciesDataMatrix constructor" begin
    dm = SpeciesDataMatrix{Int}(m, dir, spp_scrambled)
    show(dm)
    m2 = [01 02 03 04 05 06;
          07 08 09 10 11 12;
          13 14 15 16 17 18;
          19 20 21 22 23 24;
          25 26 27 28 29 30]

    @test dm.data == m2

    addspecies!(dir, "sp6")

    @test_throws ArgumentError SpeciesDataMatrix{Int}(m, dir)

    dm = SpeciesDataMatrix{Float64}(dir, 12)
    @test size(dm.data) == (6, 12)

    dm = SpeciesDataMatrix{Float64}(12.3, dir, 12)
    @test size(dm.data) == (6, 12)
    @test dm.data[3, 3] == 12.3
end

m = [01 02 03 04 05 06;
     19 20 21 22 23 24;
     07 08 09 10 11 12;
     13 14 15 16 17 18;
     25 26 27 28 29 30]

dir = SpeciesDirectory(spp)
dm = SpeciesDataMatrix{Int}(m, dir, spp_scrambled)

@testset "getindex and setindex!" begin
    @test size(dm) == (5, 6)
    @test firstindex(dm) == 1
    @test firstindex(dm, 1) == 1
    @test lastindex(dm) == 30
    @test lastindex(dm, 2) == 6
    @test dm[3, 5] == 17
    @test dm["sp4"] == [19, 20, 21, 22, 23, 24]
    @test dm[["sp1", "sp2"]] == [01 02 03 04 05 06;
        07 08 09 10 11 12]
    @test dm[["sp1", "sp2"], 1:3] == [01 02 03;
        07 08 09]
    @test dm["sp4", 3] == 21
    @test dm["sp4", 2:end] == [20, 21, 22, 23, 24]
    dm["sp2"] = [0, 0, 0, 0, 0, 0]
    @test dm["sp2"] == [0, 0, 0, 0, 0, 0]
    dm[1, 4] = 0
    @test dm[1, 4] == 0
    dm["sp4", 1:3] = [0, 0, 0]
    @test dm["sp4", 1:3] == [0, 0, 0]
    dm[["sp4", "sp5"], 1:3] = [33 33 33; 11 11 11]
    @test dm[["sp4", "sp5"], 1:3] == [33 33 33; 11 11 11]
    dm["sp5"] = [2, 2, 2, 2, 2, 2]
    @test dm["sp5"] == [2, 2, 2, 2, 2, 2]
end

@testset "I/O" begin
    write("tmp.txt", dm)

    dm2 = read_species_data("tmp.txt", Int)

    @test dm.dir.list == dm2.dir.list
    @test dm.data == dm2.data

    dm3 = read_species_data("tmp.txt", Int, dir)
    @test dm.dir.dict == dm3.dir.dict

    rm("tmp.txt")
end


text:
using EXPRES-StellarSignals
using Test

@testset "EXPRES-StellarSignals.jl" begin
    # Write your tests here.
end


text:
using Combinatorics
export importTechFile
export writedotfile
export generateCombinations

mutable struct MasterTech
    name::String
    inputs::Array
    outputs::Array
    functional_group::String
    appscore::Float64
    inrel::String
    outrel::String
    transC
    transC_reliability
    MasterTech() = new()
end

"""
This function reads a .csv file with technology and relationships
and returns a tuple of an array with sources and an array of all technologies
"""
function importTechFile(techFile::String; sourceGroup::String="U",
                        sourceAddGroup::String="Uadd", sinkGroup::String="D")

    techTable = readdlm(techFile, ';')

    ## Load the table into suitable Data Structure
    techList = MasterTech[]
    for i in 2:size(techTable,2)

        ## Copy some information
        currentTech = MasterTech()
        currentTech.name = techTable[1,i]
        currentTech.functional_group = techTable[2,i]
        currentTech.appscore = parse(Float64, string(techTable[6,i]))
        currentTech.inputs = Array{String, 1}(0)
        currentTech.outputs = Array{String, 1}(0)
        currentTech.inrel = techTable[4,i]
        currentTech.outrel = techTable[5,i]

        ## Separate Inputs and outputs.
        tmp = split(techTable[3,i], " ")
        filter!(x -> x != "", tmp)
        j = 1
        while (j <= length(tmp)) && (tmp[j] != "->")
            tmp[j] = replace(tmp[j], ",", "")
            push!(currentTech.inputs, tmp[j])
            j += 1
        end
        j += 1
        while j <= length(tmp)
            tmp[j] = replace(tmp[j], ",", "")
            push!(currentTech.outputs, tmp[j])
            j += 1
        end

        if length(intersect(currentTech.inputs, currentTech.outputs)) > 0
            error("Same input and output product(s) for Tech $(techTable[1,i]) (column $i)!
  Product(s): $(intersect(currentTech.inputs, currentTech.outputs))")
        end

        ## read transfer coefficients

        function make_transC_substance_dict(s::String)
            d = Dict{Product, Float64}()
            tmp = [split(replace(x, " ", ""), '=') for x in split(s, ',')]
            for t in tmp
                d[Product(t[1])] = parse(Float64, t[2])
            end
            return d
        end

        transC = Dict{String, Dict{Product, Float64}}()
        transC_rel = Dict{String, Float64}()
        for k in 1:NSUBSTANCE
            transC[SUBSTANCE_NAMES[k]] = make_transC_substance_dict(String(techTable[5+2*k,i]))
            transC_rel[SUBSTANCE_NAMES[k]] = convert(Float64, techTable[6+2*k,i])
        end
        currentTech.transC = transC
        currentTech.transC_reliability = transC_rel

        for sub in SUBSTANCE_NAMES
            if '>' in currentTech.outrel   # is transport technolgy
                needed_keys = Set([Product("x"), Product("airloss"), Product("soilloss"), Product("waterloss")])
            else
                if length(currentTech.outputs) > 0 # is not sink
                    out_tmp = [replace(p, "transported", "") for p in currentTech.outputs]
                    needed_keys = Set(vcat(Product.(out_tmp), [Product("airloss"), Product("soilloss"), Product("waterloss")]))
                else                # is sink
                    needed_keys = Set([Product("recovered"), Product("airloss"), Product("soilloss"), Product("waterloss")])
                end
            end
            if needed_keys != Set(keys(currentTech.transC[sub]))
                diff_element = setdiff(needed_keys, Set(keys(currentTech.transC[sub])))
                error("Transfere coefficients for Tech $(techTable[1,i]) (column $i) do not match for $sub.
  No match for: $([String(e.name) for e in diff_element])")
            end
        end

        ## Add to Array
        push!(techList, currentTech)

    end

    ## help function to generate two Techs with only "transported" or "not transported" output products
    function make2(t::Tech)

        if t.functional_group != Symbol(sinkGroup)      # of t is not a sink
            outs = t.outputs
            out_transported = filter(p -> contains(string(p.name), "transported"), outs)
            out_NOTtransported = filter(p -> !contains(string(p.name), "transported"), outs)

            ## produce array of two Techs
            techs = Tech[Tech(t.inputs,
                              out_NOTtransported,
                              t.name,
                              t.functional_group,
                              t.appscore,
                              t.n_inputs,
                              t.transC,
                              t.transC_reliability),
                         Tech(t.inputs,
                              out_transported,
                              t.name*"_trans",
                              t.functional_group,
                              t.appscore,
                              t.n_inputs,
                              t.transC,
                              t.transC_reliability)]

            return filter!(t -> length(t.outputs)>0, techs)
        else                    # if sink just return it as is
            [t]
        end
    end

    ## Now generate subtechs
    subTechList = AbstractTech[]
    for currentTech in techList
        ## Case 1: Copy | Copy
        if (currentTech.inrel == "NA" || currentTech.inrel == "AND") && (currentTech.outrel == "NA" || currentTech.outrel == "AND")

            newSubTech = Tech(currentTech.inputs,
                              currentTech.outputs,
                              currentTech.name,
                              currentTech.functional_group,
                              currentTech.appscore,
                              currentTech.transC,
                              currentTech.transC_reliability)

            append!(subTechList, make2(newSubTech))

            ## Case 2: Copy | GenVar
        elseif (currentTech.inrel == "NA" || currentTech.inrel == "AND") && (currentTech.outrel == "OR")
            error(" The relationship of output products cannot be defined as 'OR'!")

            ## Case 3: Copy | Gen1
        elseif (currentTech.inrel == "NA" || currentTech.inrel == "AND") && (currentTech.outrel == "XOR")
            error(" The relationship of output products cannot be defined as 'XOR'!")

            ## Case 4: GenVar | Hierarchy
        elseif currentTech.inrel == "OR" && (currentTech.outrel != "NA" && currentTech.outrel != "AND")

            ## Prepare (get the hierarchy)
            hira = split(currentTech.outrel, " > ")
            if length(hira) != length(currentTech.inputs)
                error("Case 4: Lenght of Hierarchy is expected to be the same as number of Input Products")
            end
            c = 1

            ## This bit gets the possible combinations of in-products
            for i in 1:length(currentTech.inputs)
                possibleCombis = collect(combinations(currentTech.inputs, i))
                for subTechIn in possibleCombis
                    ## This bit checks for the most dominant in-product to determine output
                    for candidateOut in hira
                        if in(candidateOut, subTechIn)

                            ## This bit generates the SubTech
                            subTechName = join([currentTech.name, c], "_")
                            subTechOut = join(["transported", candidateOut], "")
                            newSubTech = Tech(subTechIn,
                                              [subTechOut],
                                              subTechName,
                                              currentTech.functional_group,
                                              currentTech.appscore,
                                              currentTech.transC,
                                              currentTech.transC_reliability)
                            push!(subTechList, newSubTech)

                            c += 1
                            break
                        end
                    end
                end
            end

            ## Case 5: Gen1 | Select
        elseif currentTech.inrel == "XOR" && (currentTech.outrel != "NA" && currentTech.outrel != "AND")
            for i in 1:length(currentTech.inputs)
                subTechName = join([currentTech.name, i], "_")
                subTechIn = currentTech.inputs[i]
                subTechOut = join(["transported", subTechIn], "")
                newSubTech = Tech([subTechIn],
                                  [subTechOut],
                                  subTechName,
                                  currentTech.functional_group,
                                  currentTech.appscore,
                                  currentTech.transC,
                                  currentTech.transC_reliability)
                push!(subTechList, newSubTech)
            end

            ## Case 6: GenVar | Copy
        elseif currentTech.inrel == "OR" && (currentTech.outrel == "AND" || currentTech.outrel == "NA")
            c = 1
            for i in 1:length(currentTech.inputs)
                possibleCombis = collect(combinations(currentTech.inputs, i))
                for subTechIn in possibleCombis
                    subTechName = join([currentTech.name, c], "_")
                    newSubTech = Tech(subTechIn,
                                      currentTech.outputs,
                                      subTechName,
                                      currentTech.functional_group,
                                      currentTech.appscore,
                                      currentTech.transC,
                                      currentTech.transC_reliability)
                    append!(subTechList, make2(newSubTech))
                    c += 1
                end
            end

            ## Case 7: Gen1 | Copy
        elseif currentTech.inrel == "XOR" && (currentTech.outrel == "AND" || currentTech.outrel == "NA")
            for i in 1:length(currentTech.inputs)
                subTechIn = currentTech.inputs[i]
                subTechName = join([currentTech.name, i], "_")
                newSubTech = Tech([subTechIn],
                                  currentTech.outputs,
                                  subTechName,
                                  currentTech.functional_group,
                                  currentTech.appscore,
                                  currentTech.transC,
                                  currentTech.transC_reliability)
                append!(subTechList, make2(newSubTech))
            end

## Case 8: Other syntax
else
print(currentTech.inrel)
print(currentTech.outrel)
error("Case 8: Unrecognized Syntax Inrel.")
end
end

## Filter: for all Techs that are not in c_group:
##         either ALL inputs and outputs are "transported" or ALL inputs and outputs are "NOT transported".
function islegal(t::Tech)
    inp  = [string(s.name) for s in t.inputs]
    outp = [string(s.name) for s in t.outputs]

    inp_trans  = [contains(i, "transported") for i in inp]
    outp_trans = [contains(i, "transported") for i in outp]

    ##        is transport              |         everything is transported  | nothing is transported
    (t.functional_group == Symbol("C")) | (all(inp_trans) & all(outp_trans)) | (all(.!(inp_trans)) & all(.!(outp_trans)))
end

subTechFiltered = filter(islegal, subTechList)

# Test if all technologies in the input file exist at least in one form in the final tech lists
for t in techList
    if !any(contains(x.name, t.name) for x in subTechFiltered)
        error("Tech '$(t.name)' is not imported! Check csv file carefully!")
    end
end

## separate sources, sourcesAdd, and Technologies
sources = filter(t -> t.functional_group == Symbol(sourceGroup), subTechFiltered)
sourcesAdd = filter(t -> t.functional_group == Symbol(sourceAddGroup), subTechFiltered)
techs = filter(t -> (t.functional_group != Symbol(sourceGroup) && t.functional_group != Symbol(sourceAddGroup)), subTechFiltered)

return sources, sourcesAdd, techs

end

"""
     This function generates possible combinations of one source and all sourcesAdds
     returns all possible combinations.
    """
function generateCombinations{T <: AbstractTech}(source::T, sourcesAdd::Array{T})

    src_comb = Array{Tech}[]
    push!(src_comb, [source])

    for s_pos in Combinatorics.combinations(sourcesAdd)
        push!(s_pos, source)
        push!(src_comb, s_pos)
    end

    return src_comb
end

## ---------------------------------
## write dot file for visualisation with graphviz

"""
 Writes a DOT file of a `System`. The resulting file can be visualized with GraphViz, e,g.:
 ```
 dot -Tpng file.dot -o graph.png
 ```
 ## Arguments
 nogroup    Array of functional groups which should not be grouped in the plot

 """
function writedotfile(sys::System, file::String, no_group::Array{String}=["S", "C", "T"], options::String="")
    make_legal(name::String) = replace(name, " :: ", "")

    open(file, "w") do f
        println(f, "digraph system {")
        println(f, "rankdir=LR;")
        println(f, "node[style=filled colorscheme=pastel15];") # accent5
        if options!=""
            println(f, "$(options);")
        end

        if haskey(sys.properties, "ID")
            println(f, "label=\"ID: $(sys.properties["ID"])\";")
            println(f, "labelfontsize=22.0;")
            println(f, "labelloc=\"top\";")
            println(f, "labeljust=left")
        end

        ## define colors for function groups
        fgroups = sort(unique(t.functional_group for t in sys.techs))
        colors = Dict(fgroups[i] => mod(i,5)+1 for i in 1:length(fgroups))

        colors = Dict(:U => "# F15A31", :S => "# F99D34", :C => "# C1C430", :T => "# 70BF54", :D => "# 00B6CD")

        ## define nodes
        for t in vcat(sys.techs...)
            label = "$(t.name)\n"
            label = label * "($(t.functional_group))"
            println(f, replace("$(make_legal(t.name)) [shape=box, fillcolor=\"$(get(colors, t.functional_group, "# 999999"))\" label=\"$label\"];", ".", "_"))
        end
        ## edges
        for c in sys.connections
            println(f, replace("$(make_legal(c[2].name)) -> $(make_legal(c[3].name)) [label=\"$(c[1].name)\"];", ".", "_"))
        end

        no_group = [Symbol(x) for x in no_group]
        ## group according to functional groups
        for fg in filter(x -> !(x in no_group), fgroups)
            names = [t.name for t in sys.techs if t.functional_group==fg]
            names = map(n -> make_legal(replace(n, ".", "_")), names)
            println(f, "{ rank=same $(join(names, ' ')) }")
        end

        println(f, "}")
    end
end


text:
# This code provides helpers for loading data from the HIPPARCOS GAIA Catalog of Accelerations.
#

using FITSIO
using Tables
using Measurements

"""
    ProperMotionAnomHGCA(;gaia_id=1234)

Load proper motion anomaly data from the HIPPARCOS-GAIA Catalog of Accelerations (Brandt et al)
for a star with catalog id `gaia_id`.
The resulting velocities are in mas/yr and have the long term trend between HIPPARCOS and GAIA
already subtracted out. e.g. we would expect 0 pma if there is no companion.
"""
function ProperMotionAnomHGCA(;gaia_id,catalog=(datadep"HGCA_eDR3")*"/HGCA_vEDR3.fits")

    ## Load the Hipparcos-GAIA catalog of accelerations
    hgca = FITS(catalog) do hdu
        Tables.columntable(hdu[2])
    end

    idx = findfirst(==(gaia_id), hgca.gaia_source_id)

    # Proper motion anomaly
    # The difference between the ~instant proper motion measured by GAIA compared to the
    # long term trend between Hipparcos and GAIA
    Δμ_gaia_ra = (hgca.pmra_gaia[idx] ± hgca.pmra_gaia_error[idx]) - (hgca.pmra_hg[idx] ± hgca.pmra_hg_error[idx])
    Δμ_gaia_dec = (hgca.pmdec_gaia[idx] ± hgca.pmdec_gaia_error[idx]) - (hgca.pmdec_hg[idx] ± hgca.pmdec_hg_error[idx])

    Δμ_hip_ra = (hgca.pmra_hip[idx] ± hgca.pmra_hip_error[idx]) - (hgca.pmra_hg[idx] ± hgca.pmra_hg_error[idx])
    Δμ_hip_dec = (hgca.pmdec_hip[idx] ± hgca.pmdec_hip_error[idx]) - (hgca.pmdec_hg[idx] ± hgca.pmdec_hg_error[idx])

    return ProperMotionAnom(
        # Hipparcos epoch
        (;
            ra_epoch=years2mjd(hgca.epoch_ra_hip[idx]),
            dec_epoch=years2mjd(hgca.epoch_dec_hip[idx]),

            pm_ra=Measurements.value(Δμ_hip_ra),
            σ_pm_ra=Measurements.uncertainty(Δμ_hip_ra),

            pm_dec=Measurements.value(Δμ_hip_dec),
            σ_pm_dec=Measurements.uncertainty(Δμ_hip_dec),
        ),
        # GAIA epoch
        (;
            ra_epoch=years2mjd(hgca.epoch_ra_gaia[idx]),
            dec_epoch=years2mjd(hgca.epoch_dec_gaia[idx]),

            pm_ra=Measurements.value(Δμ_gaia_ra),
            σ_pm_ra=Measurements.uncertainty(Δμ_gaia_ra),

            pm_dec=Measurements.value(Δμ_gaia_dec),
            σ_pm_dec=Measurements.uncertainty(Δμ_gaia_dec),

        ),
    )
end
export ProperMotionAnomHGCA

"""
    gaia_plx(gaia_id=12123)

Get a distribution (TruncatedNormal) of parallax distance in mas of a source with
GAIA catalog id `gaia_id`.
"""
function gaia_plx(;gaia_id,catalog=(datadep"HGCA_eDR3")*"/HGCA_vEDR3.fits")

    # Load the Hipparcos-GAIA catalog of accelerations as a basic column table
    hgca = FITS(catalog) do hdu
        Tables.columntable(hdu[2])
    end

    idx = findfirst(==(gaia_id), hgca.gaia_source_id)
    return TruncatedNormal(hgca.parallax_gaia[idx,], hgca.parallax_gaia_error[idx,], 0, Inf)
end
export gaia_plx


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestOntologyDBs
using Test

using Catlab
using SemanticFlowGraphs

# Local file
############

# Load concepts.
db = OntologyDB()
@test !has_concept(db, "foo")
@test_throws OntologyError concept(db, "foo")
load_ontology_file(db, joinpath(@__DIR__, "data", "foobar.json"))
@test has_concept(db, "foo")
@test concept(db, "foo") isa Monocl.Ob
@test concept(db, "bar-from-foo") isa Monocl.Hom

# Concept accessors.
@test concepts(db) isa Presentation
@test concepts(db, ["foo", "bar-from-foo"]) ==
  [ concept(db, "foo"), concept(db, "bar-from-foo") ]

# Remote database
#################

# Load all concepts.
db = OntologyDB()
@test !has_concept(db, "model")
load_concepts(db)
@test has_concept(db, "model")
@test concept(db, "model") isa Monocl.Ob
@test concept(db, "fit") isa Monocl.Hom

# Load single annotation.
df_id = AnnotationID("python", "pandas", "data-frame")
@test !has_annotation(db, df_id)
@test_throws OntologyError annotation(db, df_id)
@test isa(load_annotation(db, df_id), ObAnnotation)
@test isa(annotation(db, df_id), ObAnnotation)
@test isa(annotation_document(db, df_id), AbstractDict)
@test has_annotation(db, df_id)
@test annotation(db, df_id) == annotation(db, "python/pandas/data-frame")
@test annotation(db, df_id) == annotation(db, "annotation/python/pandas/data-frame")

bad_id = AnnotationID("python", "pandas", "xxx")
@test_throws OntologyError load_annotation(db, bad_id)

# Load all annotations in a package.
series_id = AnnotationID("python", "pandas", "series")
ndarray_id = AnnotationID("python", "numpy", "ndarray")
@test_throws OntologyError annotation(db, series_id)
@test_throws OntologyError annotation(db, ndarray_id)
load_annotations(db, language="python", package="pandas")
@test_throws OntologyError annotation(db, ndarray_id)

note = annotation(db, series_id)
@test isa(note, Annotation)
@test isa(note.definition, Monocl.Ob)

note = annotation(db, AnnotationID("python", "pandas", "read-table"))
@test isa(note, Annotation)
@test isa(note.definition, Monocl.Hom)

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph