text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo
using Plots
include("methods/gradient_bundles.jl")
include("methods/least_squares.jl")
include("methods/utilities.jl")

# ## cone
contact_type = :linear
## contact_type = :nonlinear

# ## scale system for nice plots
mech = get_mechanism(:box2D,
    timestep=0.1,
    g=-1.0,
    friction_coefficient=1.0,
    contact=true,
    contact_type=contact_type);

# ## inputs
Fs = 0.01:0.01:0.2
Fsref = 0.001:0.001:0.2

# ## options
undercut = 1.0

mode = :friction
## mode = :impact
if mode == :friction && contact_type == :nonlinear
    filename = "nonlinear_friction"
elseif mode == :friction && contact_type == :linear
    filename = "linear_friction"
else
    filename = "impact"
end

plt = plot(layout=(3,1), size=(500,800), legend=:topleft)

# ## hard contact
for btol = 1e-10
    res = [box2D_dojo(mech, F, btol=btol, mode=mode) for F = Fsref]
    write_csv(["F", "x", "grad"], [[Fsref[i]; res[i]...] for i=1:length(Fsref)], "data/$(filename)_dojo_"*scn(btol)[2:end]*".csv")
    x = [r[1] for r in res]
    ∇x = [r[2] for r in res]
    plot!(plt[1,1], Fsref,  x, xlabel="F", ylabel="x", linewidth=2.0,
        label="true dynamics", color=:black)
    plot!(plt[2,1], Fsref, ∇x, xlabel="F", ylabel="∇x", linewidth=2.0,
        label="true dynamics", color=:black)#ylims=(-0.1, 1.5))
    plot!(plt[3,1], Fsref, ∇x, xlabel="F", ylabel="∇x", linewidth=2.0,
        label="true dynamics", color=:black)#, ylims=(-0.1, 1.5))
    display(plt)
end

# ## smooth contact
for btol in [1e-4, 1e-5, 1e-6, 1e-7, 1e-8] .* undercut
    res = [box2D_dojo(mech, F, btol=btol, mode=mode, undercut=undercut) for F = Fsref]
    write_csv(["F", "x", "grad"], [[Fsref[i]; res[i]...] for i=1:length(Fsref)], "data/$(filename)_dojo_"*scn(btol)[2:end]*".csv")
    x = [r[1] for r in res]
    ∇x = [r[2] for r in res]
    plot!(plt[1,1], Fsref,  x, xlabel="F", ylabel="x", linewidth=3.0, linestyle=:dot,
        label="Dojo btol="*scn(btol))
    plot!(plt[2,1], Fsref, ∇x, xlabel="F", ylabel="∇x", linewidth=3.0, linestyle=:dot,
        label="Dojo btol="*scn(btol)[2:end])
    display(plt)
end

# ## gradient bundles
for Σ in [1e-2, 1e-3, 1e-4]
    res = [box2D_gradientbundle(mech, F, N=500, Σ=Σ*I, mode=mode) for F = Fs]
    write_csv(["F", "x", "gb0", "gb1"], [[Fs[i]; res[i]...] for i=1:length(Fs)], "data/$(filename)_gb_"*scn(Σ)[2:end]*".csv")
    x = [r[1] for r in res]
    ∇x0 = [r[2] for r in res]
    ∇x1 = [r[3] for r in res]
    # plot!(plt[1,1], Fs,  x, xlabel="F", ylabel="x", linewidth=3.0, linestyle=:dot)
    plot!(plt[3,1], Fs, ∇x0, xlabel="F", ylabel="∇x", linewidth=6.0,
        linestyle=:dot,  label="0th order G.B. Σ="*scn(Σ)[2:end])
    plot!(plt[3,1], Fs, ∇x1, xlabel="F", ylabel="∇x", linewidth=3.0,
        linestyle=:dash, label="1st order G.B. Σ="*scn(Σ)[2:end])
    display(plt)
end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
hist(y, bins=50, density=true)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Salary", fontsize=14); ylabel("Density", fontsize=14)

edfit = D.fit_mle(D.Exponential, y)
xx = range(minimum(y), 2500, length=100)
yy = pdf.(edfit, xx)
plot(xx, yy, lw=3, label="Exponential distribution fit")

legend(fontsize=12)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g2.svg")) # hide

ext: jl


text:
using JuMP
using ParameterJuMP
using Polyhedra
mutable struct BalancedComplexPolytope{T, CT, VT <: AbstractVector{CT}, D<:Polyhedra.FullDim}
    d::D
    points::Vector{VT}
    optimizer_constructor
    model::Union{Nothing, MOI.ModelLike}
    sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    real_z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
    imag_z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
    t_0::Union{Nothing, MOI.VariableIndex}
    function BalancedComplexPolytope{T, CT, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        optimizer_constructor=nothing) where {T, CT, VT, D}
        new{T, CT, VT, D}(Polyhedra.FullDim_convert(D, d),
                          Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing, nothing)
    end
end
function BalancedComplexPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    VT = eltype(points)
    CT = Polyhedra.coefficient_type(VT)
    return BalancedComplexPolytope{real(CT), CT, VT, typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::BalancedComplexPolytope) = v.d

function _build_model(brp::BalancedComplexPolytope{T}) where T
    @assert brp.model === nothing
    n = length(brp.points)
    brp.model = _instantiate(brp.optimizer_constructor, T)
    # Real part of `t`
    α = Vector{MOI.VariableIndex}(undef, n)
    # Imaginary part of `t`
    β = Vector{MOI.VariableIndex}(undef, n)
    q = Vector{MOI.VariableIndex}(undef, n)
    # |t| ≤ q
    for i in 1:n
        vars, con = MOI.add_constrained_variables(brp.model, MOI.SecondOrderCone(3))
        q[i], α[i], β[i] = vars
    end
    # sum |t| ≤ sum q ≤ 1
    brp.sum_con = MOI.add_constraint(brp.model, _sum(q, T), MOI.LessThan(one(T)))
    rp = [real.(p) for p in brp.points]
    ip = [imag.(p) for p in brp.points]
    real_hull = MOI.Utilities.operate!.(
        -, T,
        _convex_combination(α, rp, brp.d),
        _convex_combination(β, ip, brp.d))
    imag_hull = MOI.Utilities.operate!.(
        +, T,
        _convex_combination(α, ip, brp.d),
        _convex_combination(β, rp, brp.d))
    return real_hull, imag_hull
end

function _parametrized_model(brp::BalancedComplexPolytope, v)
    real_hull, imag_hull = _build_model(brp)
    brp.real_z = [MOI.add_constraint(brp.model, real_hull[i], MOI.EqualTo(real(v[i]))) for i in 1:brp.d]
    brp.imag_z = [MOI.add_constraint(brp.model, imag_hull[i], MOI.EqualTo(imag(v[i]))) for i in 1:brp.d]
end

function _ratio_model(brp::BalancedComplexPolytope{T}, v) where T
    real_hull, imag_hull = _build_model(brp)
    brp.t_0 = MOI.add_variable(brp.model)
    for i in 1:brp.d
        push!(real_hull[i].terms, MOI.ScalarAffineTerm(-real(v[i]), brp.t_0))
        push!(imag_hull[i].terms, MOI.ScalarAffineTerm(-imag(v[i]), brp.t_0))
    end
    brp.real_z = [MOI.add_constraint(brp.model, real_hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
    brp.imag_z = [MOI.add_constraint(brp.model, imag_hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
end

function _update_model(p::BalancedComplexPolytope{T}, v) where T
    vars, con = MOI.add_constrained_variables(p.model, MOI.SecondOrderCone(3))
    q, α, β = vars
    MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(q, one(T)))
    for i in 1:p.d
        MOI.modify(p.model, p.real_z[i], MOI.ScalarCoefficientChange(α, real(v[i])))
        MOI.modify(p.model, p.imag_z[i], MOI.ScalarCoefficientChange(β, imag(v[i])))
    end
end

function _fix(bcp::BalancedComplexPolytope, v)
    for i in 1:bcp.d
        if bcp.t_0 === nothing
            MOI.set(bcp.model, MOI.ConstraintSet(), bcp.real_z[i], MOI.EqualTo(real(v[i])))
            MOI.set(bcp.model, MOI.ConstraintSet(), bcp.imag_z[i], MOI.EqualTo(imag(v[i])))
        else
            MOI.modify(bcp.model, bcp.real_z[i], MOI.ScalarCoefficientChange(bcp.t_0, -real(v[i])))
            MOI.modify(bcp.model, bcp.imag_z[i], MOI.ScalarCoefficientChange(bcp.t_0, -imag(v[i])))
        end
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: dlm
using AWS.Compat
using AWS.UUIDs
"""
    CreateLifecyclePolicy()

Creates a policy to manage the lifecycle of the specified AWS resources. You can create up to 100 lifecycle policies.

# Required Parameters
- `Description`: A description of the lifecycle policy. The characters ^[0-9A-Za-z _-]+ are supported.
- `ExecutionRoleArn`: The Amazon Resource Name (ARN) of the IAM role used to run the operations specified by the lifecycle policy.
- `PolicyDetails`: The configuration details of the lifecycle policy.
- `State`: The desired activation state of the lifecycle policy after creation.

# Optional Parameters
- `Tags`: The tags to apply to the lifecycle policy during creation.
"""

create_lifecycle_policy(Description, ExecutionRoleArn, PolicyDetails, State; aws_config::AWSConfig=global_aws_config()) = dlm("POST", "/policies", Dict{String, Any}("Description"=>Description, "ExecutionRoleArn"=>ExecutionRoleArn, "PolicyDetails"=>PolicyDetails, "State"=>State); aws_config=aws_config)
create_lifecycle_policy(Description, ExecutionRoleArn, PolicyDetails, State, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("POST", "/policies", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Description"=>Description, "ExecutionRoleArn"=>ExecutionRoleArn, "PolicyDetails"=>PolicyDetails, "State"=>State), args)); aws_config=aws_config)

"""
    DeleteLifecyclePolicy()

Deletes the specified lifecycle policy and halts the automated operations that the policy specified.

# Required Parameters
- `policyId`: The identifier of the lifecycle policy.

"""

delete_lifecycle_policy(policyId; aws_config::AWSConfig=global_aws_config()) = dlm("DELETE", "/policies/$(policyId)/"; aws_config=aws_config)
delete_lifecycle_policy(policyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("DELETE", "/policies/$(policyId)/", args; aws_config=aws_config)

"""
    GetLifecyclePolicies()

Gets summary information about all or the specified data lifecycle policies. To get complete information about a policy, use GetLifecyclePolicy.

# Optional Parameters
- `policyIds`: The identifiers of the data lifecycle policies.
- `resourceTypes`: The resource type.
- `state`: The activation state.
- `tagsToAdd`: The tags to add to objects created by the policy. Tags are strings in the format key=value. These user-defined tags are added in addition to the AWS-added lifecycle tags.
- `targetTags`: The target tag for a policy. Tags are strings in the format key=value.
"""

get_lifecycle_policies(; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/policies"; aws_config=aws_config)
get_lifecycle_policies(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/policies", args; aws_config=aws_config)

"""
    GetLifecyclePolicy()

Gets detailed information about the specified lifecycle policy.

# Required Parameters
- `policyId`: The identifier of the lifecycle policy.

"""

get_lifecycle_policy(policyId; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/policies/$(policyId)/"; aws_config=aws_config)
get_lifecycle_policy(policyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/policies/$(policyId)/", args; aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags for the specified resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    TagResource()

Adds the specified tags to the specified resource.

# Required Parameters
- `Tags`: One or more tags.
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.

"""

tag_resource(Tags, resourceArn; aws_config::AWSConfig=global_aws_config()) = dlm("POST", "/tags/$(resourceArn)", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag_resource(Tags, resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes the specified tags from the specified resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.
- `tagKeys`: The tag keys.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = dlm("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateLifecyclePolicy()

Updates the specified lifecycle policy.

# Required Parameters
- `policyId`: The identifier of the lifecycle policy.

# Optional Parameters
- `Description`: A description of the lifecycle policy.
- `ExecutionRoleArn`: The Amazon Resource Name (ARN) of the IAM role used to run the operations specified by the lifecycle policy.
- `PolicyDetails`: The configuration of the lifecycle policy. You cannot update the policy type or the resource type.
- `State`: The desired activation state of the lifecycle policy after creation.
"""

update_lifecycle_policy(policyId; aws_config::AWSConfig=global_aws_config()) = dlm("PATCH", "/policies/$(policyId)"; aws_config=aws_config)
update_lifecycle_policy(policyId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dlm("PATCH", "/policies/$(policyId)", args; aws_config=aws_config)


text:
"""
    DirectionalDEAModel
An data structure representing a directional distance function DEA model.
"""
struct DirectionalDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    rts::Symbol
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
end

"""
    deaddf(X, Y, Gx, Gy)
Compute data envelopment analysis directional distance function model for inputs
`X` and outputs `Y`, using directions `Gx` and `Gy`.

# Optional Arguments
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaddf(X, Y, ones(size(X)), ones(size(Y)))
Directional DF DEA Model
DMUs = 11; Inputs = 2; Outputs = 1
Returns to Scale = CRS
────────────────
      efficiency
────────────────
1   -3.43053e-16
2    3.21996
3    2.12169
4    0.0
5    6.73567
6    1.94595
7    0.0
8    3.63586
9    1.83784
10  10.2311
11   0.0
────────────────
```
"""
function deaddf(X::Matrix, Y::Matrix, Gx::Matrix, Gy::Matrix; rts::Symbol = :CRS, slack = true, Xref::Matrix = X, Yref::Matrix = Y)::DirectionalDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    nGx, mGx = size(Gx)
    nGy, sGy = size(Gy)

    nrefx, mref = size(Xref)
    nrefy, sref = size(Yref)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if nrefx != nrefy
        error("number of observations is different in inputs reference set and ouputs reference set")
    end
    if m != mref
        error("number of inputs in evaluation set and reference set is different")
    end
    if s != sref
        error("number of outputs in evaluation set and reference set is different")
    end
    if size(Gx) != size(X)
        error("size of inputs should be equal to size of inputs direction")
    end
    if size(Gy) != size(Y)
        error("size of outputs should be equal to size of outputs direction")
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Directions to use
        Gx0 = Gx[i,:]
        Gy0 = Gy[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(GLPK.Optimizer))
        @variable(deamodel, eff)
        @variable(deamodel, lambda[1:nref] >= 0)

        @objective(deamodel, Max, eff)

        @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= x0[j] - eff * Gx0[j])
        @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j] + eff * Gy0[j])

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            error("Invalid returns to scale $rts. Returns to scale should be :CRS or :VRS")
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Compute slacks
    if slack == true

        # Get first-stage efficient X and Y
        Xeff = X .- effi .* Gx
        Yeff = Y .+ effi .* Gy

        # Use additive model with radial efficient X and Y to get slacks
        radialSlacks = deaadd(Xeff, Yeff, :Ones, rts = rts, Xref = Xref, Yref = Yref)
        slackX = slacks(radialSlacks, :X)
        slackY = slacks(radialSlacks, :Y)
    else
        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    return DirectionalDEAModel(n, m, s, rts, effi, slackX, slackY, lambdaeff)

end

function deaddf(X::Vector, Y::Matrix, Gx::Vector, Gy::Matrix; rts::Symbol = :CRS, slack = true, Xref::Vector = X, Yref::Matrix = Y)::DirectionalDEAModel
    X = X[:,:]
    Xref = Xref[:,:]
    Gx = Gx[:,:]
    return deaddf(X, Y, Gx, Gy, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function deaddf(X::Matrix, Y::Vector, Gx::Matrix, Gy::Vector; rts::Symbol = :CRS, slack = true, Xref::Matrix = X, Yref::Vector = Y)::DirectionalDEAModel
    Y = Y[:,:]
    Yref = Yref[:,:]
    Gy = Gy[:,:]
    return deaddf(X, Y, Gx, Gy, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function deaddf(X::Vector, Y::Vector, Gx::Vector, Gy::Vector; rts::Symbol = :CRS, slack = true, Xref::Vector = X, Yref::Vector = Y)::DirectionalDEAModel
    X = X[:,:]
    Xref = Xref[:,:]
    Gx = Gx[:,:]
    Y = Y[:,:]
    Yref = Yref[:,:]
    Gy = Gy[:,:]
    return deaddf(X, Y, Gx, Gy, rts = rts, slack = slack, Xref = Xref, Yref = Yref)
end

function Base.show(io::IO, x::DirectionalDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Directional DF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ; ["slackY$i" for i in 1:s ]], ["$i" for i in 1:n]))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], ["$i" for i in 1:n]))
        end
    end

end


text:
# from Posit to Float64,32,16
Float64(x::Posit8) = ccall((:convertP8ToDouble, SoftPositPath), Float64, (Posit8,),x)
Float64(x::Posit16) = ccall((:convertP16ToDouble, SoftPositPath), Float64, (Posit16,),x)
Float64(x::Posit32) = ccall((:convertP32ToDouble, SoftPositPath), Float64, (Posit32,),x)

# conversion from PX2 to Float64,32,16
Float64(x::Posit8_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit8_2,),x)
Float64(x::Posit16_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit16_2,),x)
Float64(x::Posit24_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit24_2,),x)

# conversion from PX1 to Float64,32,16
Float64(x::Posit8_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit8_1,),x)
Float64(x::Posit16_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit16_1,),x)
Float64(x::Posit24_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit24_1,),x)

# conversion to Float32,16
Float32(x::AbstractPosit) = Float32(Float64(x))
Float16(x::AbstractPosit) = Float16(Float64(x))


text:
# Copyright 2007, Timothy A. Davis, Univ. of Florida
# Copyright 2017, Andrey Prokopenko, ORNL

# Adapted from Tim Davis' meshnd code
function meshsparse(G, stencil = -1)
"""MESHSPARSE convert a 2D or 3D mesh into a sparse matrix matrix.

Example:
A = meshsparse(G)
A = meshsparse(G, 5)            # 2D 5-point stencil (default for 2D case)
A = meshsparse(G, 9)            # 2D 9-point stencil
A = meshsparse(G, 7)            # 3D 7-point stencil (default for 3D case)
A = meshsparse(G, 27)           # 3D 27-point stencil
A = meshsparse(G, stencil)      # user-provided stencil

To create a sparse matrix for an m-by-n 2D mesh or m-by-n-by-k 3D mesh, use

A = meshsparse(meshnd(m, n))
A = meshsparse(meshnd(m, n, k))

G is an m-by-n-by-k matrix, with entries numbered 1 to m*n*k (with k=1 for
the 2D case). The entries in G can appear in any order, but no duplicate
entries can appear. That is sort(G[:])' must equal 1:m*n*k. A is returned as
a sparse matrix with m*n*k rows and columns whose pattern depends on the
stencil. The number of nonzeros in most rows/columns of A is equal to the
number of points in the stencil. For examples on how to specify your own
stencil, see the contents of meshsparse.m.
"""

if stencil == -1
    (m, n, k) = size(G)
    if k == 1
        stencil = 5     # 2D default is a 5-point stencil
    else
        stencil = 7     # 3D default is a 7-point stencil
    end
end

if length(stencil) == 1
    # create the stencil

    if (stencil == 5)
        # 5-point stencil (2D)
        stencil = [
                   -1  0    # north
                    1  0    # south
                    0  1    # east
                    0 -1 ]  # west

    elseif stencil == 9
        # 9-point stencil (2D)
        stencil = [
            -1  0       # north
             1  0       # south
             0  1       # east
             0 -1       # west
            -1 -1       # north-west
            -1  1       # north-east
             1 -1       # south-west
             1  1 ]     # south-east

    elseif (stencil == 7)
        # 7-point stencil (3D)
        stencil = [
            -1  0  0    # north
             1  0  0    # south
             0  1  0    # east
             0 -1  0    # west
             0  0 -1    # up
             0  0  1 ]  # down

    elseif stencil == 27
        # 27-point stencil (3D)
        stencil = zeros(26, 3)
        t = 0
        for i = -1:1
            for j = -1:1
                for k = -1:1
                    if ~(i == 0 & j == 0 & k == 0)
                        t = t + 1
                        stencil[t,:] = [i j k]
                    end
                end
            end
        end
    end
end

stencil = trunc(stencil)
(npoints, d) = size(stencil)
if d == 2
    # append zeros onto a 2D stencil to make it "3D"
    stencil = [stencil zeros(Int64, npoints, 1)]
end
(npoints, d) = size(stencil)
if d != 3
    error("invalid stencil")
end

(m, n, k) = size(G)
i1 = 1:m
j1 = 1:n
k1 = 1:k

Ti = zeros(Int64, npoints*m*n*k, 1)
Tj = zeros(Int64, npoints*m*n*k, 1)
nz = 0

for point = 1:npoints
    # find the overlapping rows of G
    idelta = stencil[point,1]
    i2 = i1 + idelta
    ki = find(ii -> (ii >= 1 && ii <= m), i2)

    # find the overlapping columns of G
    jdelta = stencil[point,2]
    j2 = j1 + jdelta
    kj = find(jj -> (jj >= 1 && jj <= n), j2)

    # find the overlapping slices of G
    kdelta = stencil[point,3]
    k2 = k1 + kdelta
    kk = find(kk -> (kk >=1 && kk <= k), k2)

    # find the nodes in G the shifted G that touch
    g2 = G[i2[ki], j2[kj], k2[kk]]     # shifted mesh
    g1 = G[i1[ki], j1[kj], k1[kk]]     # unshifted mesh

    # place the edges in the triplet list
    e = length(g1)
    Ti[(nz+1):(nz+e)] = g1[:]
    Tj[(nz+1):(nz+e)] = g2[:]
    nz = nz + e
end

# convert the triplets into a sparse matrix
Ti = Ti[1:nz]
Tj = Tj[1:nz]

A = npoints * speye(m*n*k) - sparse(Ti, Tj, 1, m*n*k, m*n*k)

return A

end # function

max_stars_repo_name: aprokop/ca-preconditioners

max_issues_repo_name: aprokop/ca-preconditioners

max_forks_repo_name: aprokop/ca-preconditioners


text:
using Test
using GaussQuadrature
using ClimateMachine.VTK: writemesh

@testset "VTK" begin
    for dim in 1:3
        N = 5
        nelem = 3
        T = Float64

        Nq = N + 1
        (r, _) = GaussQuadrature.legendre(T, Nq, GaussQuadrature.both)
        s = dim > 1 ? r : [0]
        t = dim > 2 ? r : [0]
        Nqi = length(r)
        Nqj = length(s)
        Nqk = length(t)
        Np = Nqi * Nqj * Nqk

        x1 = Array{T, 4}(undef, Nqi, Nqj, Nqk, nelem)
        x2 = Array{T, 4}(undef, Nqi, Nqj, Nqk, nelem)
        x3 = Array{T, 4}(undef, Nqi, Nqj, Nqk, nelem)

        for e in 1:nelem, k in 1:Nqk, j in 1:Nqj, i in 1:Nqi
            xoffset = nelem + 1 - 2e
            x1[i, j, k, e], x2[i, j, k, e], x3[i, j, k, e] =
                r[i] - xoffset, s[j], t[k]
        end

        if dim == 1
            x1 = x1 .^ 3
        elseif dim == 2
            x1, x2 = x1 + sin.(π * x2) / 5, x2 + exp.(-x1 .^ 2)
        else
            x1, x2, x3 = x1 + sin.(π * x2) / 5,
            x2 + exp.(-hypot.(x1, x3) .^ 2),
            x3 + sin.(π * x1) / 5
        end
        d = exp.(sin.(hypot.(x1, x2, x3)))
        s = copy(d)

        if dim == 1
            @test "test$(dim)d.vtu" ==
                  writemesh("test$(dim)d", x1; fields = (("d", d), ("s", s)))[1]
            @test "test$(dim)d.vtu" == writemesh(
                "test$(dim)d",
                x1;
                x2 = x2,
                fields = (("d", d), ("s", s)),
            )[1]
            @test "test$(dim)d.vtu" == writemesh(
                "test$(dim)d",
                x1;
                x2 = x2,
                x3 = x3,
                fields = (("d", d), ("s", s)),
            )[1]
        elseif dim == 2
            @test "test$(dim)d.vtu" == writemesh(
                "test$(dim)d",
                x1,
                x2;
                fields = (("d", d), ("s", s)),
            )[1]
            @test "test$(dim)d.vtu" == writemesh(
                "test$(dim)d",
                x1,
                x2;
                x3 = x3,
                fields = (("d", d), ("s", s)),
            )[1]
        elseif dim == 3
            @test "test$(dim)d.vtu" == writemesh(
                "test$(dim)d",
                x1,
                x2,
                x3;
                fields = (("d", d), ("s", s)),
            )[1]
        end
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 19602f12-7fc4-11eb-33c6-c7fc690e6da2
begin
	using ApproxFun, SpecialFunctions, LinearAlgebra, BenchmarkTools
	import ClassicalOrthogonalPolynomials
end

# ╔═╡ 3497f8a0-7fc4-11eb-09a4-99ede0577168
begin
	using PyPlot
	PyPlot.svg(true)
	function with_pyplot(f::Function)
		f()
		fig = gcf()
		close(fig)
		return fig
	end
end

# ╔═╡ 77a708ca-7fd3-11eb-2c99-3919a290433f
using TimerOutputs

# ╔═╡ 20713f94-7fc4-11eb-051e-fd71664ec440
Tcheb(n,x) = ClassicalOrthogonalPolynomials.chebyshevt(n, x)

# ╔═╡ 2902c6b4-7fc4-11eb-3e3d-171afb2a489c
begin
	S = Chebyshev()
	G = Fun(x -> Tcheb(128, x), S)
	I = ApproxFun.I
end

# ╔═╡ 74d93d2a-7fc4-11eb-0767-f5b765f1f70e
begin
	α = 0.
	ν = 40.
	c = 40
	c2 = Fun(t -> (1+t)^2, S)
	c1 = Fun(t -> (2α + 3)*(1+t), S)
	c0 = Fun(t -> c^2 * (1+t)^2 + (α + 1)^2 - ν^2, S)
	b1 = Fun(t -> (1+t), S)
	b0 = Fun(t -> 2α + 1, S)

	bf = BasisFunctional(c/(ν+1) >= 20 ? Integer(ceil(c)) : 0, space(G))
end;

# ╔═╡ 5f850bf2-7fc4-11eb-0b68-437c354291e2
A = [
	c2 * 𝒟^2 + c1 * 𝒟 + c0 * I    0.;
	b1 * 𝒟 + b0 * I               I;
	bf                            0;
	0                             bf
];

# ╔═╡ 5f6fa354-7fc4-11eb-3206-dbaa9c57c1fc
b = [-G;
	  0.;
	  0.;
	  0.]

# ╔═╡ 9933d9f4-7fc5-11eb-0472-975929d365b8
F, H = \(A, b, tolerance=1maximum(G.coefficients)*eps())

# ╔═╡ ee258fc0-7fc5-11eb-16fe-7512390173e1
integral(t) = (1+t)^(α+1) * (
	c * (1+t) * (1/2) * (
		besselj(ν-1, c*(1+t)) - besselj(1+ν, c*(1+t))
		) * F(t) + besselj(ν, c*(1+t)) * (α * F(t) + H(t)))

# ╔═╡ a1e75c1c-7fc6-11eb-0b03-ffb83c581358
result = integral(1) - integral(-1)

# ╔═╡ 47101788-803c-11eb-37cc-ef405ac24eed

# ╔═╡ 60a6f6dc-7fc7-11eb-39db-4d3f30fc79aa
begin
	ref = big"-5.331844823289785705118221686111276809121403304e-7"
	"abserr", (result - ref)
end

# ╔═╡ 7ec6de8e-7fc7-11eb-0e82-ab2d6d5cbef6
# @btime \($A, $b, tolerance=10maximum($G.coefficients)*eps())

# ╔═╡ 88cdb6d0-7fd3-11eb-2508-69fe5d96f67a
const to = TimerOutput()

# ╔═╡ fbaa9bac-7fd1-11eb-18a4-d36066740213
function compute_moments(ν, c, α=0.0)

	S = Chebyshev()
	I = ApproxFun.I

	c2 = Fun(t -> (1+t)^2, S)
	c1 = Fun(t -> (2α + 3)*(1+t), S)
	c0 = Fun(t -> c^2 * (1+t)^2 + (α + 1)^2 - ν^2, S)
	b1 = Fun(t -> (1+t), S)
	b0 = Fun(t -> 2α + 1, S)

	bf = BasisFunctional(c/(ν+1) >= 20 ? Integer(ceil(c)) : 0, S)

	A = [
		c2 * 𝒟^2 + c1 * 𝒟 + c0 * I    0.;
		b1 * 𝒟 + b0 * I               I;
		bf                            0.;
		0.                             bf
	];

	vals = Float64[]
	qrA = qr(A)  # solve A once
	for n in 0:128
		G = Fun(x -> Tcheb(n, x), S)
		b = [-G;
		  0.;
		  0.;
		  0.
		]
		F, H = \(qrA, b, tolerance=10maximum(G.coefficients)*eps())

		t = 1.0
		moment = (1+t)^(α+1) * (
			c * (1+t) * (1/2) * (
			besselj(ν-1, c*(1+t)) - besselj(1+ν, c*(1+t))) *
			F(t) + besselj(ν, c*(1+t)) * (α * F(t) + H(t)))
		push!(vals, moment)
	end
	return vals
end

# ╔═╡ 69731f42-7fd2-11eb-0c1f-575407f9ee64
@time compute_moments(100.0, 40)

# ╔═╡ b75114de-7fd3-11eb-224e-71b3c667824e
begin
	reset_timer!(to)
	compute_moments(500.0, 300.0)
	to
end

# ╔═╡ b739d9d6-7fd3-11eb-37cc-e7ac8c81c6d9

# ╔═╡ b71f9fa8-7fd3-11eb-3049-05b0c265a291

# ╔═╡ 2700ed36-7fcf-11eb-19c7-59fe0b88eea8
# function do_keller(ν, c, α=0.0)

# 	S = Chebyshev()
# 	G = Fun(x -> chebyshevt(4, x), S)
# 	I = ApproxFun.I

# 	c2 = Fun(t -> (1+t)^2, S)
# 	c1 = Fun(t -> (2α + 3)*(1+t), S)
# 	c0 = Fun(t -> c^2 * (1+t)^2 + (α + 1)^2 - ν^2, S)
# 	b1 = Fun(t -> (1+t), S)
# 	b0 = Fun(t -> 2α + 1, S)

# 	bf = BasisFunctional(c/(ν+1) >= 20 ? Integer(ceil(c)) : 0, space(G))

# 	A = [
# 		c2 * 𝒟^2 + c1 * 𝒟 + c0 * I    0.;
# 		b1 * 𝒟 + b0 * I               I;
# 		bf                            0;
# 		0                             bf
# 	];

# 	b = [-G;
# 	  0.;
# 	  0.;
# 	  0.
# 	]
# 	@time F, H = \(A, b, tolerance=10maximum(G.coefficients)*eps())
# 	t = 1.0

# 	return (1+t)^(α+1) * (
# 	c * (1+t) * (1/2) * (
# 		besselj(ν-1, c*(1+t)) - besselj(1+ν, c*(1+t))
# 		) * F(t) + besselj(ν, c*(1+t)) * (α * F(t) + H(t)))
# end

# ╔═╡ 8ce980b8-7fcf-11eb-1fbb-fbd0693d81d6
# @time do_keller(ν, c, α, S, G, I, c2, c1, c0, b1, b0, bf, A)

# ╔═╡ 32f414aa-7fd0-11eb-0aa0-efbd6ae3e8af
 # do_keller(3000.0, 1e-1)

# ╔═╡ 65e8471e-7fd0-11eb-0c7a-1550ff5a2094

# ╔═╡ b453481e-7fc5-11eb-2f5e-bb09402f0b88
# with_pyplot() do
# 	# axvline(r)
# 	plot(H.coefficients, "-")
# 	yscale("symlog", linthresh=1e-16)
# end

# ╔═╡ 36cf695a-7fc4-11eb-2968-61eb0cc280bd
# with_pyplot() do
# 	figure(figsize=(10,4))
# 	plot()
# end

# ╔═╡ Cell order:
# ╠═19602f12-7fc4-11eb-33c6-c7fc690e6da2
# ╟─3497f8a0-7fc4-11eb-09a4-99ede0577168
# ╠═20713f94-7fc4-11eb-051e-fd71664ec440
# ╠═2902c6b4-7fc4-11eb-3e3d-171afb2a489c
# ╠═74d93d2a-7fc4-11eb-0767-f5b765f1f70e
# ╠═5f850bf2-7fc4-11eb-0b68-437c354291e2
# ╠═5f6fa354-7fc4-11eb-3206-dbaa9c57c1fc
# ╠═9933d9f4-7fc5-11eb-0472-975929d365b8
# ╠═ee258fc0-7fc5-11eb-16fe-7512390173e1
# ╠═a1e75c1c-7fc6-11eb-0b03-ffb83c581358
# ╠═47101788-803c-11eb-37cc-ef405ac24eed
# ╠═60a6f6dc-7fc7-11eb-39db-4d3f30fc79aa
# ╠═7ec6de8e-7fc7-11eb-0e82-ab2d6d5cbef6
# ╠═77a708ca-7fd3-11eb-2c99-3919a290433f
# ╠═88cdb6d0-7fd3-11eb-2508-69fe5d96f67a
# ╠═fbaa9bac-7fd1-11eb-18a4-d36066740213
# ╠═69731f42-7fd2-11eb-0c1f-575407f9ee64
# ╠═b75114de-7fd3-11eb-224e-71b3c667824e
# ╠═b739d9d6-7fd3-11eb-37cc-e7ac8c81c6d9
# ╠═b71f9fa8-7fd3-11eb-3049-05b0c265a291
# ╠═2700ed36-7fcf-11eb-19c7-59fe0b88eea8
# ╠═8ce980b8-7fcf-11eb-1fbb-fbd0693d81d6
# ╠═32f414aa-7fd0-11eb-0aa0-efbd6ae3e8af
# ╠═65e8471e-7fd0-11eb-0c7a-1550ff5a2094
# ╠═b453481e-7fc5-11eb-2f5e-bb09402f0b88
# ╠═36cf695a-7fc4-11eb-2968-61eb0cc280bd


text:
for T in (:Double64, :ComplexD64)
  @eval begin
    Base.isnan(x::$T) = isnan(Hi(x))
    Base.isinf(x::$T) = isinf(Hi(x))
    Base.isfinite(x::$T) = isfinite(Hi(x))
    Base.iszero(x::$T) = iszero(Hi(x))
    Base.isone(x::$T) = isone(Hi(x)) && iszero(Lo(x))
    Base.isinteger(x::$T) = isinteger(Hi(x)) && isinteger(Lo(x))
  end
end

for F in (:iseven, :isodd)
  for T in (:Float64, :Double64)
    @eval begin

      function Base.$F(x::$T)
          !isinteger(x) && return false
          absx = abs(x)
          if absx < typemax(Int64)
            $F(Int64(x))
          elseif absx < typemax(Int128)
             $F(Int128(x))
          else
             $F(BigInt(x))
          end
      end

    end
  end
end


text:
module GridsTests

using Test
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.CellData

using Gridap.Geometry: GridMock

trian = GridMock()
test_grid(trian)

q1i = Point(0.5,0.5)
np1 = 4
q1 = fill(q1i,np1)

q2i = Point(0.25,0.25)
np2 = 3
q2 = fill(q2i,np2)

q = CompressedArray([q1,q2],get_cell_type(trian))

cell_map = get_cell_map(trian)

x1i = Point(0.5, 0.5)
x2i = Point(1.25, 0.25)
x3i = Point(1.75, 0.5)
x4i = Point(0.5, 1.5)
x5i = Point(1.5, 1.5)
x1 = fill(x1i,np1)
x2 = fill(x2i,np2)
x3 = fill(x3i,np2)
x4 = fill(x4i,np1)
x5 = fill(x5i,np1)
x = [x1,x2,x3,x4,x5]

@test all(lazy_map(test_field,cell_map,q,x))
_x = lazy_map(evaluate,cell_map,q)
test_array(_x,x)

# from LagrangianRefFE

quad8 = LagrangianRefFE(Float64,QUAD,2)
grid = Grid(quad8)
@test num_nodes(grid) == num_nodes(quad8)

# from Polytope

grid = Grid(ReferenceFE{2},WEDGE)
@test num_cells(grid) == 5
@test num_cell_dims(grid) == 2
@test num_point_dims(grid) == 3

grid = Grid(ReferenceFE{3},WEDGE)
@test num_cells(grid) == 1
@test num_cell_dims(grid) == 3
@test num_point_dims(grid) == 3

## get low dim grid
#
#grid = Grid(ReferenceFE{1},trian)

# Test linear grid

grid = compute_reference_grid(TRI,5)
model = UnstructuredDiscreteModel(grid)
f_to_c = get_faces(get_grid_topology(model),1,2)
num_boundary_facets =  num_facets(TRI)*5

@test count(isequal(1),map(length,f_to_c)) == num_boundary_facets

Ω_tri = Triangulation(model)
dΩ = Measure(Ω_tri,2)

test_grid(grid)

@test sum( ∫(1)dΩ ) ≈ 1/2

grid = compute_reference_grid(TET,5)
model = UnstructuredDiscreteModel(grid)
f_to_c = get_faces(get_grid_topology(model),2,3)
num_boundary_facets = num_facets(TET)*num_cells(Ω_tri)

@test count(isequal(1),map(length,f_to_c)) == num_boundary_facets

Ω_tet = Triangulation(model)
dΩ = Measure(Ω_tet,2)

test_grid(grid)

@test sum( ∫(1)dΩ ) ≈ 1/6

end # module


text:
#==========================================================================================#

# TYPE

mutable struct OLS <: MLE

    sample::Microdata
    β::Vector{Float64}
    V::Matrix{Float64}

    OLS() = new()
end

#==========================================================================================#

# CONSTRUCTOR

function OLS(MD::Microdata)
    obj        = OLS()
    obj.sample = MD
    return obj
end

#==========================================================================================#

# ESTIMATION

function _fit!(obj::OLS, ::UnitWeights)
    y     = getvector(obj, :response)
    x     = getmatrix(obj, :control)
    obj.β = x \ y
end

function _fit!(obj::OLS, w::AbstractWeights)
    y     = getvector(obj, :response)
    x     = getmatrix(obj, :control)
    v     = Diagonal(w) * x
    obj.β =  (v' * x) \ (v' * y)
end

#==========================================================================================#

# SCORE (MOMENT CONDITIONS)

score(obj::OLS) = Diagonal(residuals(obj)) * getmatrix(obj, :control)

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

jacobian(obj::OLS, ::UnitWeights)      = - crossprod(getmatrix(obj, :control))
jacobian(obj::OLS, w::AbstractWeights) = - crossprod(getmatrix(obj, :control), w)

# HOMOSCEDASTIC VARIANCE MATRIX

function _vcov!(obj::OLS, corr::Homoscedastic, w::UnitWeights)
    σ²    = sum(abs2, residuals(obj)) / dof_residual(obj)
    obj.V = lmul!(-σ², inv(jacobian(obj, w)))
end

function _vcov!(obj::OLS, corr::Homoscedastic, w::AbstractWeights)
    σ²    = sum(abs2.(residuals(obj)), w) / dof_residual(obj)
    obj.V = lmul!(-σ², inv(jacobian(obj, w)))
end

#==========================================================================================#

# LINEAR PREDICTOR

function predict(obj::OLS, MD::Microdata)
    (getnames(obj, :control) != getnames(MD, :control)) && throw("missing variables")
    getmatrix(MD, :control) * obj.β
end

# FITTED VALUES

fitted(obj::OLS, MD::Microdata) = predict(obj, MD)

# DERIVATIVE OF FITTED VALUES

jacobexp(obj::OLS) = copy(getmatrix(obj, :control))

#==========================================================================================#

# UTILITIES

coefnames(obj::OLS) = getnames(obj, :control)
mtitle(obj::OLS)    = "OLS"

# COEFFICIENT OF DETERMINATION

adjr2(obj::OLS) = 1.0 - (1.0 - r2(obj)) * (nobs(obj) - 1) / dof_residual(obj)
r2(obj::OLS)    = _r2(obj, getweights(obj))

function _r2(obj::OLS, ::UnitWeights)
    y   = response(obj)
    ŷ   = fitted(obj)
    rss = sum(abs2, y .- ŷ)
    tss = sum(abs2, y .- mean(y))
    return 1.0 - rss / tss
end

function _r2(obj::OLS, w::AbstractWeights)
    y   = response(obj)
    μ   = mean(y, w)
    ŷ   = fitted(obj)
    rss = sum(abs2.(y .- ŷ), w)
    tss = sum(abs2.(y .- μ), w)
    return 1.0 - rss / tss
end

# LIKELIHOOD FUNCTION

function _loglikelihood(obj::OLS, w::AbstractWeights)
    n  = sum(w)
    r  = residuals(obj)
    r .= abs2.(r)
    σ² = sum(r, w) / n
    return - (log(σ²) + log2π + 1.0) * n / 2.0
end

# LIKELIHOOD FUNCTION UNDER NULL MODEL

function _nullloglikelihood(obj::OLS, w::AbstractWeights)
    n  = sum(w)
    r  = copy(getvector(obj, :response))
    μ  = mean(r, w)
    r .= abs2.(r .- μ)
    σ² = sum(r, w) / n
    return - (log(σ²) + log2π + 1.0) * n / 2.0
end

# DEVIANCE

function _deviance(obj::OLS, w::AbstractWeights)
    r  = residuals(obj)
    r .= abs2.(r)
    return sum(r, w)
end

# DEVIANCE UNDER NULL MODEL

function _nulldeviance(obj::OLS, w::AbstractWeights)
    r  = copy(getvector(obj, :response))
    μ  = mean(r, w)
    r .= abs2.(r .- μ)
    return sum(r)
end


text:
let
    problem = BabyPOMDP()

    solver = RandomSolver(rng=MersenneTwister(1))

    policy = solve(solver, problem)

    sim = RolloutSimulator(max_steps=10, rng=MersenneTwister(1))

    r = simulate(sim, problem, policy, updater(policy), initialstate_distribution(problem))

    @test isapprox(r, -27.27829, atol=1e-3)
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
using Query
using FileIO

"""
Create a parameter `component`_`name` with the given value,
and connect parameter `name` within `component` to this distinct global parameter.
"""
function setdistinctparameter(m::Model, component::Symbol, name::Symbol, value)
    globalname = Symbol(string(component, '_', name))
    set_param!(m, component, name, globalname, value)
end

"""
Load raw RV output into reformat_RV_outputs
"""
function load_RV(filename::String, RVname::String;
                    output_path::String=joinpath(@__DIR__, "../../output/"),
                    time_filter::Int=2200,
                    region_filter::String="LatAmerica")

    df = DataFrame(load(joinpath(output_path, "$filename.csv")))
    cols = names(df)

    # apply filters if necessary, currently the function supports a time filter
    # of a single time value and a region filter of a single region
    if in("time", cols)

        if in("region", cols)
            filtered_df = df |> @query(i, begin
                @where i.time == time_filter
                @where i.region == region_filter
                @select i
            end) |> DataFrame

        else
            filtered_df = df |> @query(i, begin
                @where i.time == time_filter
                @select i
            end) |> DataFrame
        end

        return filtered_df[!, Symbol(RVname)]

    # no filters applied
    else
        return df[!, Symbol(RVname)]
    end

end


text:
"""
This struct contains all the data for creating a market project
    for price forcasts using CEM and endogeneous Economic Dispatch.
"""

mutable struct MarketProject

    name::String                         # name of project
    project_type::String                 # is the project of storage type
    tech_type::String                    # technology type
    fixed_cost::Float64                   # $/unit/investment period
    queue_cost::Vector{Float64}           # queue cost
    marginal_energy_cost::Float64         # $/MW/hour
    marginal_reserveup_cost::Float64      # $/MW/hour
    marginal_reservedown_cost::Float64    # $/MW/hour
    expansion_cost::Vector{Float64}        # $/unit
    discount_rate::Float64                # discount rate
    min_gen::Float64                      # MW/unit
    max_gen::Float64                      # MW/unit
    min_input::Float64                    # MW/unit
    max_input::Float64                    # MW/unit
    efficiency_in::Float64               # input efficiency
    efficiency_out::Float64              # output efficiency
    min_storage::Float64                  # minimum storage capacity
    max_storage::Float64                  # maximum storage capacity
    init_storage::Float64                 # initial storage
    availability::Array{Float64, 2}      # Hourly availability factor
    derating_factor::Float64              # Derating factor
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}    # MW/unit/hour
    max_reserveup::Float64                # MW/unit/hour
    max_reservedown::Float64              # MW/unit/hour
    existing_units::Float64               # existing units
    units_in_queue::Vector{Float64}        # units in queue
    build_lead_time::Int64                 # build lead time
    remaining_build_time::Int64            # remaining build time
    max_new_options::Int                   # maximum units
    base_cost_units::Int64                 # number of units that can be built without capital cost multiplier
    capex_years::Int64                    # capital cost recovery years
    life_time::Int64                      # life_time
    remaining_life::Int64                # Remaining life_time
    capacity_eligible::Bool              # eligible for capacity market participation
    rec_eligible::Bool                   # eligible for rps compliance
    zone::String                         # project zone
    ownedby::Vector{String}              # onned by which investors

    function MarketProject(
        name::String,
        project_type::String,
        tech_type::String,
        fixed_cost::Number,
        queue_cost::Vector{<: Number},
        marginal_energy_cost::Number,
        marginal_reserveup_cost::Number,
        marginal_reservedown_cost::Number,
        expansion_cost::Vector{Float64},
        discount_rate::Number,
        min_gen::Number,
        max_gen::Number,
        min_input::Number,
        max_input::Number,
        efficiency_in::Number,
        efficiency_out::Number,
        min_storage::Number,
        max_storage::Number,
        init_storage::Number,
        availability::Array{<: Number, 2},
        derating_factor::Number,
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}},
        max_reserveup::Number,
        max_reservedown::Number,
        existing_units::Number,
        units_in_queue::Vector{<: Number},
        build_lead_time::Number,
        remaining_build_time::Number,
        max_new_options::Number,
        base_cost_units::Number,
        capex_years::Number,
        life_time::Number,
        remaining_life::Int64,
        capacity_eligible::Bool,
        rec_eligible::Bool,
        zone::String,
        ownedby::Vector{String})

        @assert fixed_cost >= 0
        @assert marginal_energy_cost >= 0
        @assert marginal_reserveup_cost >= 0
        @assert marginal_reservedown_cost >= 0
        @assert all(expansion_cost .>= 0)
        @assert all(queue_cost .>= 0)
        @assert discount_rate >= 0
        @assert min_gen >= 0
        @assert max_gen >= 0
        @assert min_input >= 0
        @assert max_input >= 0
        @assert efficiency_in >= 0
        @assert efficiency_out >= 0
        @assert min_storage >= 0
        @assert max_storage >= 0
        @assert init_storage >= 0
        @assert all(availability .>= 0)
        @assert derating_factor >= 0
        @assert existing_units >= 0
        @assert all(units_in_queue .>= 0)
        @assert build_lead_time >= 0
        @assert remaining_build_time >= 0
        @assert base_cost_units >= 0
        @assert max_new_options >= base_cost_units
        @assert capex_years >- 0
        @assert life_time >= 0
        @assert remaining_life >= 0

        new(name,
            project_type,
            tech_type,
            fixed_cost,
            queue_cost,
            marginal_energy_cost,
            marginal_reserveup_cost,
            marginal_reservedown_cost,
            expansion_cost,
            discount_rate,
            min_gen,
            max_gen,
            min_input,
            max_input,
            efficiency_in,
            efficiency_out,
            min_storage,
            max_storage,
            init_storage,
            availability,
            derating_factor,
            ramp_limits,
            max_reserveup,
            max_reservedown,
            existing_units,
            units_in_queue,
            build_lead_time,
            remaining_build_time,
            max_new_options,
            base_cost_units,
            capex_years,
            life_time,
            remaining_life,
            rec_eligible,
            capacity_eligible,
            zone,
            ownedby
            )

    end

end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
#using Revise
using CutCell
include("useful_routines.jl")

@test allapprox(CutCell.reference_bottom_face_midpoint(),[0.,-1.])
@test allapprox(CutCell.reference_right_face_midpoint(),[1.,0.])
@test allapprox(CutCell.reference_top_face_midpoint(),[0.,1.])
@test allapprox(CutCell.reference_left_face_midpoint(),[-1.,0.])

refmidpoints = CutCell.reference_face_midpoints()
cellmap = CutCell.CellMap([0.,0.],[1.,1.])
spmidpoints = cellmap.(refmidpoints)
@test length(spmidpoints) == 4
@test allapprox(spmidpoints[1],[0.5,0.])
@test allapprox(spmidpoints[2],[1.,0.5])
@test allapprox(spmidpoints[3],[0.5,1.])
@test allapprox(spmidpoints[4],[0.,0.5])

p = [1. 2. 3.]

bp = CutCell.extend_to_face(p,1)
testbp = vcat(p,-ones(3)')
@test allapprox(bp,testbp)

rp = CutCell.extend_to_face(p,2)
testrp = vcat(ones(3)',p)
@test allapprox(rp,testrp)

tp = CutCell.extend_to_face(p,3)
testtp = vcat(p,ones(3)')
@test allapprox(tp,testtp)

lp = CutCell.extend_to_face(p,4)
testlp = vcat(-ones(3)',p)
@test allapprox(testlp,lp)

points = [1.0 2.0 3.0]
bp = vcat(points,-ones(3)')
rp = vcat(ones(3)',points)
tp = vcat(points,ones(3)')
lp = vcat(-ones(3)',points)

@test allapprox(CutCell.extend_to_face(points,1),bp)
@test allapprox(CutCell.extend_to_face(points,2),rp)
@test allapprox(CutCell.extend_to_face(points,3),tp)
@test allapprox(CutCell.extend_to_face(points,4),lp)

quad1d = tensor_product_quadrature(1,4)
points = quad1d.points
bp = vcat(points,-ones(4)')
rp = vcat(ones(4)',points)
tp = vcat(points,ones(4)')
lp = vcat(-ones(4)',points)
bq = CutCell.extend_to_face(quad1d,1)
rq = CutCell.extend_to_face(quad1d,2)
tq = CutCell.extend_to_face(quad1d,3)
lq = CutCell.extend_to_face(quad1d,4)
@test allapprox(bq.points,bp)
@test allapprox(rq.points,rp)
@test allapprox(tq.points,tp)
@test allapprox(lq.points,lp)
@test allapprox(bq.weights,quad1d.weights)
@test allapprox(rq.weights,quad1d.weights)
@test allapprox(tq.weights,quad1d.weights)
@test allapprox(lq.weights,quad1d.weights)

facequads = CutCell.face_quadratures(4)
@test allapprox(facequads[1].points,bp)
@test allapprox(facequads[2].points,rp)
@test allapprox(facequads[3].points,tp)
@test allapprox(facequads[4].points,lp)
@test allapprox(facequads[1].weights,quad1d.weights)
@test allapprox(facequads[2].weights,quad1d.weights)
@test allapprox(facequads[3].weights,quad1d.weights)
@test allapprox(facequads[4].weights,quad1d.weights)

polyorder = 1
numqp = 2
basis = TensorProductBasis(2, polyorder)
levelset = InterpolatingPolynomial(1, basis)
nf = CutCell.number_of_basis_functions(basis)
mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [2, 1], nf)
nodalcoordinates = CutCell.nodal_coordinates(mesh)

normal = [1.0, 0.0]
x0 = [0.5, 0.0]
levelsetcoeffs = plane_distance_function(nodalcoordinates, normal, x0)

cutmesh = CutCell.CutMesh(levelset,levelsetcoeffs,mesh)
facequads = CutCell.FaceQuadratures(levelset,levelsetcoeffs,cutmesh,numqp)
idx = [5,9,6,10,7,11,8,12,1,0,2,0,3,0,4,0]
testfacetoquad = reshape(idx,2,4,2)
@test allequal(testfacetoquad,facequads.facetoquad)


text:
using StatsBase
using Base.Test

@test xlogx(0) === 0.0
@test_approx_eq xlogx(2) 2.0 * log(2.0)

@test xlogy(0, 1) === 0.0
@test_approx_eq xlogy(2, 3) 2.0 * log(3.0)

@test_approx_eq logistic(2) 1.0 / (1.0 + exp(-2.0))
@test_approx_eq logit(0.5) 0.0
@test_approx_eq logit(logistic(2)) 2.0

@test_approx_eq softplus(2.0) log(1.0 + exp(2.0))
@test_approx_eq softplus(-2.0) log(1.0 + exp(-2.0))
@test_approx_eq softplus(10000) 10000.0
@test_approx_eq softplus(-10000) 0.0

@test_approx_eq invsoftplus(softplus(2.0)) 2.0
@test_approx_eq invsoftplus(softplus(-2.0)) -2.0

@test_approx_eq logsumexp(2.0, 3.0) log(exp(2.0) + exp(3.0))
@test_approx_eq logsumexp(10002, 10003) 10000 + logsumexp(2.0, 3.0)

@test_approx_eq logsumexp([1.0, 2.0, 3.0]) 3.40760596444438
@test_approx_eq logsumexp([1.0, 2.0, 3.0] .+ 1000.) 1003.40760596444438

x = [1.0, 2.0, 3.0]
r = exp(x) ./ sum(exp(x))
@test_approx_eq softmax([1.0, 2.0, 3.0]) r
softmax!(x)
@test_approx_eq x r


text:
# ImmutableRef base-class:
abstract type ImmutableRef{T} <: Ref{T} end
Base.setindex!(r::ImmutableRef, val, idcs...) = error("setindex! object of type $(typeof(r)) is immutable")

struct ImmutableRefValue{T} <: ImmutableRef{T}
    _val::T
end
Base.getindex(r::ImmutableRefValue) = r._val

ImmutableRef(x::T) where {T} = ImmutableRefValue{T}(x)

# DelayFunction is a type for specifying a delay object's computation;
# once the computation is peformed, the ref that holds the DelayFunction
# object will point to the result of the computation instead
struct _DelayFunction
    _f::Function
    _mux::Threads.SpinLock
    function _DelayFunction(f::Function)
        new(f, Threads.SpinLock())
    end
end

"""
    Delay{T}

Delays objects are ImmutableRef objects that store a function of zero arguments;
when the delay object is dereferenced (`delay[]`) then the function is called and
its return value is saved in the delay for all subsequent dereferences. Delays
are thread-safe as long as their functions are also thread-safe.

# Examples
```jtdoctest
# Define a function that will be calculated in the delay:
julia> function test_delay()
           println("Calculating...")
           10
       end
test_delay (generic function with 1 method)

julia> d = Delay(test_delay)
Delay{Any}(<unrealized>)

julia> d[]
Calculating...
10

julia> d[]
10

julia> d
Delay{Any}(10)

# We can also specify a type for the delay:
julia> d = Delay{Integer}(test_delay)
Delay{Integer}(<unrealized>)

julia> d[]
Calculating...
10

julia> d = Delay{String}(test_delay)
Delay{String}(<unrealized>)

julia> d[]
MethodError: Cannot `convert` an object of type Int64 to an object of type String
```
"""
mutable struct Delay{T} <: ImmutableRef{T}
    _val::Union{T,_DelayFunction}
    function Delay{T}(df::_DelayFunction) where {T}
        if T <: _DelayFunction
            error("Delay{$T} cannot reference _DelayFunction")
        end
        new(df)
    end
end
Delay{T}(f::Function) where {T} = Delay{T}(_DelayFunction(f))
Delay(f::Function) = Delay{Any}(_DelayFunction(f))
# we don't want people poking around in Delay objects:
Base.getproperty(delay::Delay, sym::Symbol) = error("type $(typeof(delay)) has no public field $sym")
Base.setproperty!(delay::Delay, sym::Symbol, val) = error("type $(typeof(delay)) has no public field $sym")
Base.propertynames(delay::Delay) = ()
# The accessor, which may run the delay:
function Base.getindex(d::Delay{T}) where {T}
    val = getfield(d, :_val)
    if isa(val, _DelayFunction)
        mux = val._mux
        try
            lock(mux)
            # make sure someone else didn't calculate before we locked:
            lockedval = getfield(d, :_val)
            if lockedval === val
                # it hasn't been memoized yet and we have the lock:
                val = val._f()
                tval::T = val
                if isa(tval, _DelayFunction)
                    error("getindex Delay function returned a value with type _DelayFunction")
                end
                setfield!(d, :_val, tval)
            else
                val = lockedval
            end
        finally
            unlock(mux)
        end
    end
    val
end
Base.fetch(d::Delay) = d[]
"""
    isrealized(delay::Delay)
Yields true if the given delay object has already been calculated and
memoized and false otherwise.
"""
isrealized(delay::Delay) = !isa(getfield(delay, :_val), _DelayFunction)
function Base.show(io::IO, ::MIME"text/plain", d::Delay{T}) where {T}
    val = getfield(d, :_val)
    if isa(val, _DelayFunction)
        print(io, "$(typeof(d))(<unrealized>)")
    else
        print(io, "$(typeof(d))($val)")
    end
end

# Promises

# Generic box type to denote an undelivered promise
struct _PromiseState{T}
    _val::Ref{T}
    _cond::Threads.Condition

    function _PromiseState{T}() where {T}
        new(Ref{T}(), Threads.Condition())
    end
end
"""
    Promise{T}
Promise objects are repositories for values that may or may not yet have
been calculated. To wait on a promise or to retrieve the value once it
has been delivered, simply access its value: `promise[]`. To deliver a
value to a promise, use `deliver(promise, value)`. Promises are safe to
use across multiple threads as long as its fields are not changed.
"""
mutable struct Promise{T} <: ImmutableRef{T}
    _state::Union{T,_PromiseState}
    function Promise{T}() where {T}
        if T <: _PromiseState
            error("Promise object cannot be initialized with _PromiseState type")
        end
        new(_PromiseState{T}())
    end
end
# we don't want people poking around in Promise objects:
Base.getproperty(promise::Promise, sym::Symbol) = error("type $(typeof(promise)) has no public field $sym")
Base.setproperty!(promise::Promise, sym::Symbol, val) =
    error("type $(typeof(promise)) has no public field $sym")
Base.propertynames(promise::Promise) = ()
function Base.getindex(promise::Promise{T}) where {T}
    st = getfield(promise, :_state)
    if isa(st, _PromiseState)
        c = st._cond
        lock(c)
        try
            # make sure nobody else got there before we locked the mutex
            if st === getfield(promise, :_state)
                # okay, we just need to wait on the condition...
                wait(c)
                # when the condition is finished, promise._state will be the
                # delivered value, so we don't need to do anything else
            end
            st = getfield(promise, :_state)
        finally
            unlock(c)
        end
    end
    st
end
Base.fetch(p::Promise) = d[]
"""
    deliver(promise::Promise{T}, val::T)
Delivers the given value to the given promise then returns val. If the
promise has already been realized, then an error is thrown.
"""
function deliver(promise::Promise{T}, val::T) where {T}
    st = getfield(promise, :_state)
    if isa(st, _PromiseState)
        c = st._cond
        lock(c)
        try
            if st === getfield(promise, :_state)
                setfield!(promise, :_state, val)
                notify(c, val)
            else
                error("deliver cannot deliver value to a realized promise")
            end
        finally
            unlock(c)
        end
    else
        error("deliver cannot deliver value to a realized promise")
    end
    val
end
"""
    isrealized(promise::Promise{T})
Yields true if the given promise object has been realized and false otherwise.
"""
isrealized(promise::Promise{T}) where {T} = !isa(getfield(promise, :_state), _PromiseState)

function Base.show(io::IO, ::MIME"text/plain", p::Promise{T}) where {T}
    st = getfield(p, :_state)
    if isa(st, _PromiseState)
        print(io, "$(typeof(p))(<unrealized>)")
    else
        print(io, "$(typeof(p))($st)")
    end
end

max_issues_repo_name: noahbenson/FunctionalUtilities

max_forks_repo_name: noahbenson/FunctionalUtilities


text:
using Pkg
Pkg.activate(pwd())

# # Introduction to Flux

include(joinpath(pwd(), "lecture_11", "utilities.jl"))

using RDatasets
using Random

Random.seed!(666)
iris = dataset("datasets", "iris")

#+

X = Matrix(iris[:, 1:4])
y = iris.Species

X_train, y_train, X_test, y_test, classes = prepare_data(X', y; dims=2)

# ## Creating the network

using Flux

n_hidden = 5
m = Chain(
    Dense(size(X_train,1), n_hidden, relu),
    Dense(n_hidden, size(y_train,1), identity),
    softmax,
)

#+

m(X_train)
params(m[2])[2] .= [-1;0;1]

## Training the network

using Flux: crossentropy

L(x, y) = crossentropy(m(x), y)
L(X_train, y_train)

#+

ps = params(m)
grad = gradient(() -> L(X_train, y_train), params(X_train))
size(grad[X_train])

#+

using Plots

opt = Descent(0.1)
max_iter = 250

acc_test = zeros(max_iter)
for i in 1:max_iter
    gs = gradient(() -> L(X_train, y_train), ps)
    Flux.Optimise.update!(opt, ps, gs)
    acc_test[i] = accuracy(X_test, y_test)
end

plot(acc_test, xlabel="Iteration", ylabel="Test accuracy", label="", ylim=(-0.01,1.01))


text:
#===============================================================================
# Some practice using Julia control flow structures
#
# Author: Carl Cortright
# Date: 12/27/2016
#
===============================================================================#

# Compound expressions
z = begin a=1; b=2; a+b end
x = (a=1; b=2; a+b)

# Conditionals
x=2
y=3
if x < y
  println("x is less than y")
end

# Note: Variables introducted in if statements will stay in scope after execution

#
# Tests a basic if-statement
#
function controlTest(x,y)
  # Note: julia returns the return value of the last executed statement from an if
  z = if x < y
    println("x is less than y")
  elseif x > y
    println("y is greater than x")
  else
    println("x equals y")
  end
  println(z)
end
controlTest(2,1)

# Note: Julia only takes true or false as arguments to conditional statements
# Note: Julia has short-circuit evaluation for && and ||

# While loops
i = 1
while i<=5
  println(i)
  i += 1
end

# For loops
# 1:5 is a range
for i = 1:5
  println(i)
end

for i in range(1,5)
  println(i)
end

# A cool way to write that ^
for i ∈ range(1,5)
  println(i)
end

# Testing a break statement
while i < 20
  if i === 10
    break
  end
  i += 1
end
println(i)

for i = 1:10
   if i % 3 != 0
    #  continue moves us to the next iteration
     continue
   end
   println(i)
 end

# Julia will nest loops over multiple ranges
for i=1:2, j=2:8
  println("i=$i j=$j")
end

# Playing with exceptions
type weirderror <: Exception end
f(x) = x >= 0 ? exp(-x) : throw(DomainError)
f(1)

# Error stops execution immediately
fuzzy_sqrt(x) = x >= 0 ? sqrt(x) : error("x must be positive")
fuzzy_sqrt(1)
warn("we are plaing with errors")
info("here's some extra info to help you diagnose what's happening")

function g(x)
  try
    sqrt(x)
  catch
    sqrt(complex(x,0))
  end
end
println(g(1))
println(g(-1))

# Try-catch statements will also save the type of error to a variable.
try
  sqrt(-1)
catch y
  if isa(y, DomainError)
    println("Found the problem")
  end
end

# Julia also includes finally blocks
try
  sqrt(-1)
catch
  println("handling the exception")
finally
  println("executing the finally")
end

# Using tasks to solve the producer/consumer problem
function producer()
  produce("start")
  for i in 1:4
    produce(i)
  end
  produce("end")
end
p = Task(producer)
c = ""
while c != "end"
  c = consume(p)
  println(c)
end
# A task is also an interable object
for x in Task(producer)
  println(x)
end

# This is all pretty powerful, but in reality, both produce and consume use
# yeildto() to avoid race conditions and implement monitor-like behavior

max_issues_repo_name: carlcortright/Playground


text:
using PyPlot
x = range(0, stop=1, length=50)
plot(x, sin.(2x).*exp.(-x/3))
savefig(joinpath(@__DIR__, "output", "script2.svg"))

max_stars_repo_name: thanktua-spp/website

max_issues_repo_name: thanktua-spp/website

max_forks_repo_name: thanktua-spp/website


text:
module TestUtil

using Test
using AMLPipelineBase
using AMLPipelineBase.Utils
using DataFrames: nrow
using Random

function test_utils()
  data = getiris()

  catnum=find_catnum_columns(data)
  @test catnum == ([5],[1,2,3,4])

  (r,l)=holdout(nrow(data),0.9)
  @test length.((r,l)) == (15,135)

  kf=kfold(nrow(data),10)
  @test size.(kf) |> unique == [(135,)]

  Y=data.Species |> collect
  @test score(:accuracy,Y,Y) == 100.0
  @test score(:accuracy,Y,reverse(Y)) |> round == 33.0

  @test (infer_eltype.(eachcol(data)) .<: [Float64,Float64,Float64,Float64,AbstractString]) |> sum == 5

  dc = Dict(Dict(:a=>:c,:d=>Dict(:e=>:f)))
  @test nested_dict_to_tuples(dc) == Set(Any[([:d, :e], :f), ([:a], :c)])

  dx = deepcopy(dc)
  nested_dict_set!(dx,[:d,:e],:h)
  @test dx[:d][:e] == :h

  @test nested_dict_merge(dx,dc) == dc
  @test mergedict(dc,dx) == dx

  x = [2,200,332,201,missing,50,missing,6,missing,7,8,109]
  @test skipmean(x) |> round == 102.0
  @test skipmedian(x) == 50.0
  @test skipstd(x) |> round == 118.0
end
@testset "Utils" begin
  test_utils()
end

function test_metric()
   data = getiris()
   X=data[:,2:end]
   Y=data[:,1] |> collect
   rf = RandomForest()
   catf = CatFeatureSelector()
   numf = NumFeatureSelector()
   ohe = OneHotEncoder()

   pl = (catf |> ohe) + (numf) |> rf

   rmse(X,Y) = score(:rmse,X,Y)
   @test crossvalidate(pl,X,Y,rmse,10,false).mean < 0.5
   @test crossvalidate(pl,X,Y;metric=rmse,verbose=false).mean < 0.5

   mse(X,Y) = score(:mse,X,Y)
   @test crossvalidate(pl,X,Y,mse,10,false).mean < 0.2
   @test crossvalidate(pl,X,Y;metric=mse,verbose=false).mean < 0.2

   mae(X,Y) = score(:mae,X,Y)
   @test crossvalidate(pl,X,Y,mae,10,false).mean < 0.2
   @test crossvalidate(pl,X,Y;metric=mae,verbose=false).mean < 0.2

   @test crossvalidate(pl,X,Y,mae,10,false).mean < 0.2
   @test crossvalidate(pl,X,Y;metric=mae,verbose=false).mean < 0.2

   tx,ty,tsx,tsy = train_test_split(data[:,1:4],data.Species |> collect)
   fit!(pl,tx,ty)
   score(:accuracy,transform!(pl,tsx),tsy) > 0.50
   pl1 = fit(pl,tx,ty)
   score(:accuracy,transform(pl1,tsx),tsy) > 0.50
end
@testset "Utils Metrics" begin
  test_metric()
end

end


text:
# 3d.jl 3D arrays

export jim

# 3D (or higher), built on array version
function jim(z::AbstractArray{T} ; kwargs...) where {T <: Number}
    zz = reshape(z, size(z,1), size(z,2), :)
    out = [@view zz[:,:,i] for i in 1:size(zz,3)]
    jim(out ; kwargs...)
end


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control

using  ModiaBase
using  Unitful
using  Measurements
import MonteCarloMeasurements
using  DataStructures: OrderedDict, OrderedSet
using  DataFrames

export SimulationModel, measurementToString, get_lastValue

"""
    baseType(T)

Return the base type of a type T.

# Examples
```
baseType(Float32)                # Float32
baseType(Measurement{Float64})   # Float64
```
"""
baseType(::Type{T})                                           where {T}                  = T
baseType(::Type{Measurements.Measurement{T}})                 where {T<:AbstractFloat}   = T
baseType(::Type{MonteCarloMeasurements.Particles{T,N}})       where {T<:AbstractFloat,N} = T
baseType(::Type{MonteCarloMeasurements.StaticParticles{T,N}}) where {T<:AbstractFloat,N} = T

"""
    str = measurementToString(v)

Return variable `v::Measurements.Measurement{FloatType}` or a vector of such variables
in form of a string will the full number of significant digits.
"""
measurementToString(v::Measurements.Measurement{FloatType}) where {FloatType} =
       string(Measurements.value(v)) * " ± " * string(Measurements.uncertainty(v))

function measurementToString(v::Vector{Measurements.Measurement{FloatType}}) where {FloatType}
    str = string(typeof(v[1])) * "["
    for (i,vi) in enumerate(v)
        if i > 1
            str = str * ", "
        end
        str = str * measurementToString(vi)
    end
    str = str * "]"
    return str
end

"""
    simulationModel = SimulationModel{FloatType}(
            modelModule, modelName, getDerivatives!, equationInfo, x_startValues,
            parameters, variableNames;
            modelModule = nothing,
            vSolvedWithInitValuesAndUnit::OrderedDict{String,Any}(),
            vEliminated::Vector{Int}=Int[],
            vProperty::Vector{Int}=Int[],
            var_name::Function = v->nothing)

# Arguments

- `modelModule`: Module in which `@instantiateModel` is invoked (it is used for `Core.eval(modelModule, ...)`),
  that is evaluation of expressions in the environment of the user.
- `modelName::String`: Name of the model
- `getDerivatives::Function`: Function that is used to evaluate the model equations,
  typically generated with [`TinyModia.generate_getDerivatives!`].
- `equationInfo::ModiaBase.EquationInfo`: Information about the states and the equations.
- `x_startValues`:: Deprecated (is no longer used).
- `parameters`: A hierarchical NamedTuple of (key, value) pairs defining the parameter and init/start values.
- variableNames: A vector of variable names. A name can be a Symbol or a String.
"""
mutable struct SimulationModel{FloatType}
    modelModule
    modelName::String
    getDerivatives!::Function
    equationInfo::ModiaBase.EquationInfo
    linearEquations::Vector{ModiaBase.LinearEquations{FloatType}}
    variables::OrderedDict{String,Int}      # Dictionary of variables and their result indices (negated alias has negativ index)
    zeroVariables::OrderedSet{String}       # Set of variables that are identically to zero
    vSolvedWithInitValuesAndUnit::OrderedDict{String,Any}   # Dictionary of (names, init values with units) for all explicitly solved variables with init-values defined
    p::AbstractVector                       # Parameter and init/start values
                                            # (parameter values are used in getDerivatives!,
                                            #  init/start values are extracted and stored in x_start)
    separateObjects::OrderedDict{Int,Any}   # Dictionary of separate objects
    storeResult::Bool
    isInitial::Bool
    isTerminal::Bool
    time::FloatType
    nGetDerivatives::Int                    # Number of getDerivatives! calls
    x_start::Vector{FloatType}
    der_x::Vector{FloatType}
    result::Vector{Tuple}
    algorithmType::DataType                  # Type of integration algorithm (used in default-heading of plot)

    function SimulationModel{FloatType}(modelModule, modelName, getDerivatives!, equationInfo, x_startValues,
                                        parameters, variableNames;
                                        vSolvedWithInitValuesAndUnit::AbstractDict = OrderedDict{String,Any}(),
                                        vEliminated::Vector{Int} = Int[],
                                        vProperty::Vector{Int}   = Int[],
                                        var_name::Function       = v -> nothing) where {FloatType}

        # Construct result dictionaries
        variables = OrderedDict{String,Int}()
        zeroVariables = OrderedSet{String}()
        vSolvedWithInitValuesAndUnit2 = OrderedDict{String,Any}( [(string(key),vSolvedWithInitValuesAndUnit[key]) for key in keys(vSolvedWithInitValuesAndUnit)] )
            # Store variables
            for (i, name) in enumerate(variableNames)
                variables[string(name)] = i
            end

            # Store eliminated variables
            for v in vEliminated
                name = var_name(v)
                if ModiaBase.isZero(vProperty, v)
                    push!(zeroVariables, name)
                elseif ModiaBase.isAlias(vProperty, v)
                    name2 = var_name( ModiaBase.alias(vProperty, v) )
                    variables[name] = variables[name2]
                else # negated alias
                    name2 = var_name( ModiaBase.negAlias(vProperty, v) )
                    variables[name] = -variables[name2]
                end
            end

        # Construct parameter values that are copied into the code
        parameterValues = [eval(p) for p in values(parameters)]

        # Construct data structure for linear equations
        linearEquations = ModiaBase.LinearEquations{FloatType}[]
        for leq in equationInfo.linearEquations
            push!(linearEquations, ModiaBase.LinearEquations{FloatType}(leq...))
        end

        # Construct dictionary for separate objects
        separateObjects = OrderedDict{Int,Any}()

        # Initialize execution flags
        isInitial       = true
        isTerminal      = false
        storeResult     = false
        nGetDerivatives = 0

        new(modelModule, modelName, getDerivatives!, equationInfo, linearEquations, variables, zeroVariables,
            vSolvedWithInitValuesAndUnit2, parameterValues,
            separateObjects, storeResult, isInitial, isTerminal, convert(FloatType, 0), nGetDerivatives,
            zeros(FloatType,0), zeros(FloatType,0), Tuple[])
    end
end

"""
    floatType = getFloatType(simulationModel::SimulationModel)

Return the floating point type with which `simulationModel` is parameterized
(for example returns: `Float64, Float32, DoubleFloat, Measurements.Measurement{Float64}`).
"""
getFloatType(m::SimulationModel{FloatType}) where {FloatType} = FloatType

"""
    hasParticles(value)

Return true, if `value` is of type `MonteCarloMeasurements.StaticParticles` or
`MonteCarloMeasurements.Particles`.
"""
hasParticles(value) = typeof(value) <: MonteCarloMeasurements.StaticParticles ||
                      typeof(value) <: MonteCarloMeasurements.Particles

"""
    get_value(obj::NamedTuple, name::String)

Return the value identified by `name` from the potentially hierarchically
`NamedTuple obj`. If `name` is not in `obj`, the function returns `missing`.

# Examples
```julia
s1 = (a = 1, b = 2, c = 3)
s2 = (v1 = s1, v2 = (d = 4, e = 5))
s3 = (v3 = s2, v4 = s1)

@show get_value(s3, "v3.v1.b")   # returns 2
@show get_value(s3, "v3.v2.e")   # returns 5
@show get_value(s3, "v3.v1.e")   # returns missing
```
"""
function get_value(obj::NamedTuple, name::String)
    if length(name) == 0 || length(obj) == 0
        return missing
    end
    j = findnext('.', name, 1)
    if isnothing(j)
        key = Symbol(name)
        return haskey(obj,key) ? obj[key] : missing
    elseif j == 1
        return missing
    else
        key = Symbol(name[1:j-1])
        if haskey(obj,key) && typeof(obj[key]) <: NamedTuple && length(name) > j
            get_value(obj[key], name[j+1:end])
        else
            return missing
        end
    end
end

"""
    appendName(path::String, name::Symbol)

Return `path` appended with `.` and `string(name)`.
"""
appendName(path, key) = path == "" ? string(key) : path * "." * string(key)

"""
    get_names(obj::NamedTuple)

Return `Vector{String}` containing all the names present in `obj`

# Examples
```julia
s1 = (a = 1, b = 2, c = 3)
s2 = (v1 = s1, v2 = (d = 4, e = 5))
s3 = (v3 = s2, v4 = s1)

@show get_names(s3)
```
"""
function get_names(obj::NamedTuple)
    names = String[]
    get_names!(obj, names, "")
    return names
end
function get_names!(obj::NamedTuple, names::Vector{String}, path::String)::Nothing
    for (key,value) in zip(keys(obj), obj)
        name = appendName(path, key)
        if typeof(value) <: NamedTuple
            get_names!(value, names, name)
        else
            push!(names, name)
        end
    end
    return nothing
end

"""
    get_lastValue(model::SimulationModel, name::String; unit=true)

Return the last stored value of variable `name` from `model`.
If `unit=true` return the value with its unit, otherwise with stripped unit.

If `name` is not known or no result values yet available, an info message is printed
and the function returns `nothing`.
"""
function get_lastValue(m::SimulationModel, name::String; unit=true)
    if haskey(m.variables, name)
        if length(m.result) == 0
            @info "get_lastValue(model,\"$name\"): No results yet available."
            return nothing
        end

        resIndex = m.variables[name]
        negAlias = false
        if resIndex < 0
            resIndex = -resIndex
            negAlias = true
        end
        value = m.result[end][resIndex]
        if negAlias
            value = -value
        end
    elseif name in m.zeroVariables
        # Unit missing (needs to be fixed)
        value = 0.0
    else
        value = get_value(m.p[1], name)
        if ismissing(value)
            @info "get_lastValue: $name is not known and is ignored."
            return nothing;
        end
    end

    return unit ? value : ustrip(value)
end

get_xe(x, xe_info) = xe_info.length == 1 ? x[xe_info.startIndex] : x[xe_info.startIndex:xe_info.startIndex + xe_info.length-1]

function set_xe!(x, xe_info, value)::Nothing
    if xe_info.length == 1
        x[xe_info.startIndex] = value
    else
        x[xe_info.startIndex:xe_info.startIndex + xe_info.length-1] = value
    end
    return nothing
end

"""
    success = init!(simulationModel, startTime, tolerance, merge, log, logParameters, logStates)

Initialize `simulationModel::SimulationModel` at `startTime`. In particular:

- Empty result data structure.

- Merge parameter and init/start values into simulationModel.

- Construct x_start.

- Call simulationModel.getDerivatives! once with isInitial = true to
  compute and store all variables in the result data structure at `startTime`
  and initialize simulationModel.linearEquations.

- Check whether explicitly solved variables that have init-values defined,
  have the required value after initialization (-> otherwise error).

If initialization is successful return true, otherwise false.
"""
function init!(m::SimulationModel, startTime, tolerance, merge,
               log::Bool, logParameters::Bool, logStates::Bool)::Bool
    empty!(m.result)

	# Apply updates from merge Map and log parameters
    if isnothing(merge)
        if logParameters
            parameters = m.p[1]
            @showModel parameters
        end
    else
        m.p = [recursiveMerge(m.p[1], merge)]
        if logParameters
            updatedParameters = m.p[1]
            @showModel updatedParameters
        end
    end

    # Re-initialize dictionary of separate objects
    empty!(m.separateObjects)

    # Update startIndex, length in x_info and determine x_start
    FloatType = getFloatType(m)
    x_start = []
    startIndex = 1
    for xe_info in m.equationInfo.x_info
        xe_info.startIndex = startIndex

        # Determine init/start value in m.p[1]
        xe_value = get_value(m.p[1], xe_info.x_name)
        if ismissing(xe_value)
            error("Missing start/init value ", xe_info.x_name)
        end
        if hasParticles(xe_value)
            len = 1
            push!(x_start, xe_value)
        else
            len = length(xe_value)
            push!(x_start, xe_value...)
        end
        if len != xe_info.length
            printstyled("Model error: ", bold=true, color=:red)
            printstyled("Length of ", xe_info.x_name, " shall be changed from ",
                        xe_info.length, " to $len\n",
                        "This is currently not support in TinyModia.", bold=true, color=:red)
            return false
        end
        startIndex += xe_info.length
    end
    nx = startIndex - 1
    m.equationInfo.nx = nx

    # Temporarily remove units from x_start
    # (TODO: should first transform to the var_unit units and then remove)
    converted_x_start = convert(Vector{FloatType}, [ustrip(v) for v in x_start])  # ustrip.(x_start) does not work for MonteCarloMeasurements
    m.x_start = deepcopy(converted_x_start)

    if logStates
        # List init/start values
        x_table = DataFrames.DataFrame(state=String[], init=Any[], unit=String[], nominal=Float64[])
        for xe_info in m.equationInfo.x_info
            xe_init = get_xe(m.x_start, xe_info)
            if hasParticles(xe_init)
                xe_init = string(minimum(xe_init)) * " .. " * string(maximum(xe_init))
            end
            push!(x_table, (xe_info.x_name, xe_init, xe_info.unit, xe_info.nominal))
        end
        show(stdout, x_table; allrows=true, allcols=true, rowlabel = Symbol("#"), summary=false, eltypes=false)
        println("\n")
    end

    m.der_x = zeros(FloatType, nx)

    # Initialize model, linearEquations and compute and store all variables at the initial time
    if log
        println("      Initialization at time = ", startTime, " s")
    end
    m.nGetDerivatives = 0
    m.isInitial   = true
    m.storeResult = true
#    m.getDerivatives!(m.der_x, m.x_start, m, startTime)
    Base.invokelatest(m.getDerivatives!, m.der_x, m.x_start, m, startTime)
    m.isInitial   = false
    m.storeResult = false
    # Check vSolvedWithInitValuesAndUnit
    if length(m.vSolvedWithInitValuesAndUnit) > 0
        names = String[]
        valuesBeforeInit = Any[]
        valuesAfterInit  = Any[]

        for (name, valueBefore) in m.vSolvedWithInitValuesAndUnit
            valueAfterInit  = get_lastValue(m, name, unit=false)
            valueBeforeInit = ustrip.(valueBefore)
            if !isnothing(valueAfterInit) && abs(valueBeforeInit - valueAfterInit) >= max(abs(valueBeforeInit),abs(valueAfterInit),0.01*tolerance)*tolerance
                push!(names, name)
                push!(valuesBeforeInit, valueBeforeInit)
                push!(valuesAfterInit , valueAfterInit)
            end
        end

        if length(names) > 0
            v_table = DataFrames.DataFrame(name=names, beforeInit=valuesBeforeInit, afterInit=valuesAfterInit)
            #show(stderr, v_table; allrows=true, allcols=true, summary=false, eltypes=false)
            #print("\n\n")
            ioTemp = IOBuffer();
            show(ioTemp, v_table; allrows=true, allcols=true, rowlabel = Symbol("#"), summary=false, eltypes=false)
            str = String(take!(ioTemp))
            close(ioTemp)
            printstyled("Model error: ", bold=true, color=:red)
            printstyled("The following variables are explicitly solved for, have init-values defined\n",
                        "and after initialization the init-values are not respected\n",
                        "(remove the init-values in the model or change them to start-values):\n",
                        str, bold=true, color=:red)
            return false
        end
    end
    return true
end

"""
    terminate!(m::SimulationModel, x, time)

Terminate model.
"""
function terminate!(m::SimulationModel, x, t)::Nothing
    m.isTerminal = true
    Base.invokelatest(m.getDerivatives!, m.der_x, x, m, t)
    m.isTerminal = false
    return nothing
end

"""
    outputs!(x, t, integrator)

DifferentialEquations FunctionCallingCallback function for `SimulationModel`
that is used to store results at communication points.
"""
function outputs!(x, t, integrator)::Nothing
    m = integrator.p
    m.storeResult = true
#    m.getDerivatives!(m.der_x, x, m, t)
    Base.invokelatest(m.getDerivatives!, m.der_x, x, m, t)

    m.storeResult = false
    return nothing
end

"""
    derivatives!(derx, x, m, t)

DifferentialEquations callback function to get the derivatives.
"""
function derivatives!(der_x, x, m, t)::Nothing
#    m.getDerivatives!(der_x, x, m, t)
    Base.invokelatest(m.getDerivatives!, der_x, x, m, t)
    return nothing
end

"""
    addToResult!(simulationModel, variableValues...)

Add `variableValues...` to `simulationModel::SimulationModel`.
It is assumed that the first variable in `variableValues` is `time`.
"""
function addToResult!(m::SimulationModel, variableValues...)::Nothing
    push!(m.result, variableValues)
    return nothing
end

"""
    code = generate_getDerivatives!(AST, equationInfo, parameters, variables, functionName;
                                    hasUnits=false)

Return the code of the `getDerivatives!` function as `Expr` using the
Symbol `functionName` as function name. By `eval(code)` or
`fc = @RuntimeGeneratedFunction(code)` the function is compiled and can afterwards be called.

# Arguments

- `AST::Vector{Expr}`: Abstract Syntax Tree of the equations as vector of `Expr`.

- `equationInfo::ModiaBase.EquationInfo`: Data structure returned by `ModiaBase.getSortedAndSolvedAST
            holding information about the states.

- `parameters`: Vector of parameter names (as vector of symbols)

- `variables`: Vector of variable names (as vector of symbols). The first entry is expected to be time, so `variables[1] = :time`.

- `functionName::Function`: The name of the function that shall be generated.

# Optional Arguments

- `hasUnits::Bool`: = true, if variables have units. Note, the units of the state vector are defined in equationinfo.
"""
function generate_getDerivatives!(AST::Vector{Expr}, equationInfo::ModiaBase.EquationInfo,
                                  parameters, variables, functionName::Symbol;
                                  hasUnits=false)

    # Generate code to copy x to struct and struct to der_x
    x_info     = equationInfo.x_info
    code_x     = Expr[]
    code_der_x = Expr[]
    code_p     = Expr[]

    if length(x_info) == 1 && x_info[1].x_name == "" && x_info[1].der_x_name == ""
        # Explicitly solved pure algebraic variables. Introduce dummy equation
        push!(code_der_x, :( _der_x[1] = -_x[1] ))
    else
        i1 = 0
        i2 = 0
        for (i, xe) in enumerate(x_info)
            i1 = i2 + 1
            i2 = i1 + xe.length - 1
            indexRange = i1 == i2 ? :($i1) :  :( $i1:$i2 )
            x_name     = xe.x_name_julia
            der_x_name = xe.der_x_name_julia
            # x_name     = Meta.parse("m."*xe.x_name)
            # der_x_name = Meta.parse("m."*replace(xe.der_x_name, r"der\(.*\)" => s"var\"\g<0>\""))
            if !hasUnits || xe.unit == ""
                push!(code_x, :( $x_name = _x[$indexRange] ) )
            else
                x_unit = xe.unit
                push!(code_x, :( $x_name = _x[$indexRange]*@u_str($x_unit)) )
            end
            if hasUnits
                push!(code_der_x, :( _der_x[$indexRange] = ustrip( $der_x_name )) )
            else
                push!(code_der_x, :( _der_x[$indexRange] = $der_x_name ))
            end
        end
    end
    for (i,pi) in enumerate(parameters)
        push!(code_p, :( $pi = _m.p[$i] ) )
    end

    timeName = variables[1]
    if hasUnits
        code_time = :( $timeName = _time*u"s" )
    else
        code_time = :( $timeName = _time )
    end

    # Generate code of the function
    code = quote
                function $functionName(_der_x, _x, _m, _time)::Nothing
                    _m.time = _time
                    _m.nGetDerivatives += 1
                    simulationModel = _m
                    $code_time
                    $(code_p...)
                    $(code_x...)
                    $(AST...)
                    $(code_der_x...)

                    if _m.storeResult
                        TinyModia.addToResult!(_m, $(variables...))
                    end

                    return nothing
                end
            end
    return code
end


text:
##################################
### Particles
##################################

"""
    Particles(type::Symbol, properties::Dict{String})

Particles of a type (typically `:stars`, `:dm`, `:gas`, or `:bh`) with their properties in a `Dict`.
The `Dict` has `String` keys, which are set to all uppercase by default, with exception
of `iM` and `Zs`.

The properties are accessible via `particles.prop` (case insensitive).

# Examples
```jldoctest
julia> p = GadgetGalaxies.Particles(:stars, Dict("ID"=>[1, 2], "MASS"=>[1e11, 2e11]));

julia> p.id
2-element Vector{Int64}:
 1
 2

julia> p.pos = [1 2; 3 4; 5 6]
3×2 Matrix{Int64}:
 1  2
 3  4
 5  6

julia> p.properties["POS"] === p.pos
true
```
"""
struct Particles
    type::Symbol
    properties::Dict{String,Any}
end

function Base.getproperty(obj::Particles, sym::Symbol)
    if sym in fieldnames(Particles)
        return getfield(obj, sym)
    else
        return obj.properties[sym |> String |> uppercase |> _particle_property_exceptions]
    end
end

function Base.setproperty!(obj::Particles, sym::Symbol, val)
    if sym in fieldnames(Particles)
        setfield!(obj, sym, val)
    else
        obj.properties[sym |> String |> uppercase |> _particle_property_exceptions] = val
    end
end

Base.getindex(obj::Particles, str::String) = obj.properties[str]
Base.getindex(obj::Particles, sym::Symbol) = Base.getproperty(obj, sym)
Base.setindex!(obj::Particles, val, str::String) = Base.setindex!(obj.properties, val, str)
Base.setindex!(obj::Particles, val, sym::Symbol) = Base.setproperty!(obj, sym, val)
Base.keys(obj::Particles) = keys(obj.properties) .|> lowercase .|> Symbol
Base.haskey(obj::Particles, key) = haskey(obj.properties, key)
Base.haskey(obj::Particles, key::Symbol) = key in keys(obj)
Base.values(obj::Particles) = values(obj.properties)
Base.propertynames(obj::Particles) = [:type; keys(obj)]
Base.copy(obj::Particles) = Particles(obj.type, copy(obj.properties))

function Base.show(io::IO, ::MIME"text/plain", obj::Particles)
    printstyled(io, String(obj.type); bold=true)
    len = 0
    for val in values(obj)
        if typeof(val) <: AbstractArray && size(val, ndims(val)) > len
            len = size(val, ndims(val))
        end
    end
    println(io, ": $len Particles")
    print(io, " ")
    props = join(keys(obj.properties) .|> _particle_property_exceptions |> sort, " ")
    println(io, props)
end

function _particle_property_exceptions(key::String)
    if key == "IM"
        return "iM"
    elseif key == "ZS"
        return "Zs"
    else
        return key
    end
end

const _particle_type_id = Dict(:gas => 0, :dm => 1, :stars => 4, :bh => 5)

"""
    particle_type_id(type::Symbol)

Returns the Gadget particle type from a particle `Symbol`
(currently `:gas`, `:dm`, `:stars`, and `:bh`).
"""
particle_type_id(type::Symbol) = _particle_type_id[type]

##################################
### Galaxy
##################################

abstract type AbstractGalaxy end

"""
    Galaxy(snapshot::Snapshot,
           isub::Integer,
           subid::Union{HaloID,Nothing},
           particles::Dict{Symbol,Particles})

Galaxy of a given `snapshot` ([`Snapshot`](@ref)) with the zero-based subhalo index `isub`
with [`Particles`](@ref).
Including `subid` (containing information on the exact position in the subfind files) makes
accessing properties of the galaxy faster than when only having `isub`.

The particles are accessible via `galaxy.particletype` or `galaxy[:particletype]`
(see [`Particles`](@ref) for typical types).

# Examples
```jldoctest
julia> p = GadgetGalaxies.Particles(:stars, Dict("ID"=>[1, 2], "MASS"=>[1e11, 2e11]));

julia> g = Galaxy(Snapshot("box", 13), 1532, nothing, Dict(:stars=>p));

julia> g.stars === p
true
```
"""
struct Galaxy <: AbstractGalaxy
    snapshot::Snapshot
    isub::Integer
    subid::Union{HaloID,Nothing}
    particles::Dict{Symbol,Particles}
end

"""
    Galaxy(snapshot::Snapshot, isub[, get_id=true])

Convenience method for creating a `Galaxy` with only a `snapshot` and a subfind id.
Will try to read "MSUB" to extract the  if `get_id` is `true`.

[`Particles`](@ref) can still be added after initializing the galaxy by calling
for example `galaxy.stars = [...]`.
"""
function Galaxy(snapshot::Snapshot, isub::Integer, get_id::Bool)
    if get_id
        subbase = string(snapshot.subbase)
        subfile = GadgetIO.select_file(subbase, 0)

        n_files = read_header(subfile).num_files
        _, subid = read_halo_prop_and_id(subbase, isub, "MSUB", n_files; verbose=false)
    else
        subid = nothing
    end

    return Galaxy(snapshot, isub, subid, Dict{Symbol,Particles}())
end

Galaxy(snapshot::Snapshot, isub::Integer) = Galaxy(snapshot, isub, true)

"""
    GalaxyGroup(snapshot::Snapshot,
           igroup::Integer,
           groupid::Union{HaloID,Nothing},
           particles::Dict{Symbol,Particles})

Group of a given `snapshot` ([`Snapshot`](@ref)) with the zero-based group index `igroup`
with [`Particles`](@ref).
Including `groupid` (containing information on the exact position in the subfind files) makes
accessing properties of the galaxy faster than when only having `igroup`.

The particles are accessible via `group.particletype` or `group[:particletype]`
(see [`Particles`](@ref) for typical types).

# Examples
```jldoctest
julia> p = GadgetGalaxies.Particles(:stars, Dict("ID"=>[1, 2], "MASS"=>[1e11, 2e11]));

julia> g = GalaxyGroup(Snapshot("box", 13), 164, nothing, Dict(:stars=>p));

julia> g.stars === p
true
```
"""
struct GalaxyGroup <: AbstractGalaxy
    snapshot::Snapshot
    igroup::Integer
    groupid::Union{HaloID,Nothing}
    particles::Dict{Symbol,Particles}
end

"""
    GalaxyGroup(snapshot::Snapshot, igroup[, get_id=true])

Convenience method for creating a `GalaxyGroup` with only a `snapshot` and a subfind id.
Will try to read "MTOT" to extract the  if `get_id` is `true`.

[`Particles`](@ref) can still be added after initializing the group by calling for example
`group.stars = [...]`.
"""
function GalaxyGroup(snapshot::Snapshot, igroup::Integer, get_id::Bool)
    if get_id
        subbase = string(snapshot.subbase)
        subfile = GadgetIO.select_file(subbase, 0)

        n_files = read_header(subfile).num_files
        _, groupid = read_halo_prop_and_id(subbase, igroup, "MTOT", n_files; verbose=false)
    else
        groupid = nothing
    end

    return GalaxyGroup(snapshot, igroup, groupid, Dict{Symbol,Particles}())
end

GalaxyGroup(snapshot::Snapshot, igroup::Integer) = GalaxyGroup(snapshot, igroup, true)

function Base.getproperty(obj::T, sym::Symbol) where {T<:AbstractGalaxy}
    if sym in fieldnames(T)
        return getfield(obj, sym)
    else
        return obj.particles[sym]
    end
end

function Base.setproperty!(obj::T, sym::Symbol, val) where {T<:AbstractGalaxy}
    if sym in fieldnames(T)
        setfield!(obj, sym, val)
    else
        obj.particles[sym] = val
    end
end

Base.getindex(obj::AbstractGalaxy, sym::Symbol) = Base.getproperty(obj, sym)
Base.setindex!(obj::AbstractGalaxy, val, sym::Symbol) = Base.setproperty!(obj, sym, val)
Base.keys(obj::AbstractGalaxy) = keys(obj.particles) |> collect
Base.haskey(obj::AbstractGalaxy, key) = haskey(obj.particles, key)
Base.values(obj::AbstractGalaxy) = values(obj.particles)
Base.propertynames(obj::Galaxy) = [[:snapshot, :isub, :subid]; keys(obj)]
Base.propertynames(obj::GalaxyGroup) = [[:snapshot, :igroup, :groupid]; keys(obj)]

function Base.copy(obj::T) where {T<:AbstractGalaxy}
    particles = Dict(key => copy(p) for (key, p) in pairs(obj.particles))
    T(obj.snapshot, getind(obj), getid(obj), particles)
end

function Base.show(io::IO, ::MIME"text/plain", obj::AbstractGalaxy)
    printstyled(io, typeof(obj); bold=true)
    try
        z = round(read_redshift(obj.snapshot); digits=2)
        println(io, " at z=$z")
    catch
        println()
    end
    println(io, " $(getindname(obj)) $(getind(obj))")
    for key in sort(keys(obj))
        show(io, "text/plain", obj[key])
    end
end

getindname(g::Galaxy) = "isub"
getindname(g::GalaxyGroup) = "igroup"

getind(g::Galaxy) = g.isub
getind(g::GalaxyGroup) = g.igroup

getid(g::Galaxy) = g.subid
getid(g::GalaxyGroup) = g.groupid

get_subfind_type(g::Galaxy) = 2
get_subfind_type(g::GalaxyGroup) = 1